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

## Declaration

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.

## 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 co-director. 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 lab-mates: 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 - BIA2017-83805-R - 01/01/2018 – 31/12/2020)

Finally, I want to thank my family for their extreme patience and constant love and support.

## Abstract

In this work we study the numerical simulation of particle-laden 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 particle-fluid 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 order-of-magnitude 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 history-dependent 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 backward-coupling effects according to the theory of multicomponent continua, allowing the code to deal with arbitrarily dense flow regimes.

## Nomenclature

Roman Symbols

${\displaystyle a}$ - particle's radius

${\displaystyle A_{p}}$ - particle cross-sectional area

${\displaystyle {\mathcal {B}}}$ - body (continuum mechanics)

${\displaystyle {\mathcal {B}}_{k}}$ - body belonging to component ${\displaystyle k}$ (continuum mechanics)

${\displaystyle C_{D,{\hbox{Shah}}}}$ - drag coefficient in Shah's (2007) model

${\displaystyle c_{s}}$ - local sound speed

${\displaystyle k}$ - Stokes–Einstein–Smoluchowski diffusivity

${\displaystyle \mathbb {E} ^{3}}$ - three-dimensional affine Euclidean space

${\displaystyle f}$ - friction factor (boundary layer)

${\displaystyle {\boldsymbol {F}}_{A}}$ - added mass force

${\displaystyle F_{A}}$ - characteristic scale of ${\displaystyle {\boldsymbol {F}}_{A}}$

${\displaystyle {\mathcal {F}}_{b}}$ - total body work

${\displaystyle {\boldsymbol {F}}_{\hbox{contact}}}$ - sum of contact forces

${\displaystyle {\boldsymbol {F}}_{D}}$ - steady drag force

${\displaystyle F_{D}}$ - characteristic scale of ${\displaystyle {\boldsymbol {F}}_{D}}$

${\displaystyle {\boldsymbol {F}}_{H}}$ - Boussinesq–Basset history force

${\displaystyle F_{H}}$ - characteristic scale of ${\displaystyle {\boldsymbol {F}}_{H}}$

${\displaystyle {\boldsymbol {F}}_{\hbox{hydr}}}$ - hydrodynamic force

${\displaystyle {\boldsymbol {F}}_{k}}$ - Momentum exchange between component ${\displaystyle k}$ and the rest of components

${\displaystyle {\boldsymbol {F}}_{L,{\hbox{Herr}}}}$ - Herron lift force

${\displaystyle {\boldsymbol {F}}_{L,{\hbox{Rub}}}}$ - Rubinow and Keller lift force

${\displaystyle F_{L,{\hbox{Rub}}}}$ - characteristic value of ${\displaystyle {\boldsymbol {F}}_{L,{\hbox{Rub}}}}$

${\displaystyle {\boldsymbol {F}}_{L,{\hbox{Saff}}}}$ - Saffman lift force

${\displaystyle {\mathit {Fr}}}$ - Froude number

${\displaystyle {\mathcal {F}}_{s}}$ - total surface work

${\displaystyle {\boldsymbol {F}}_{U}}$ - undisturbed-flow force

${\displaystyle F_{U}}$ - characteristic scale of ${\displaystyle {\boldsymbol {F}}_{U}}$

${\displaystyle {\boldsymbol {F}}_{W}}$ - static submerged weight

${\displaystyle {\boldsymbol {g}}}$ - gravity acceleration

${\displaystyle H^{1}(\Omega )}$ - space of functions with square integrable derivatives in ${\displaystyle \Omega }$

${\displaystyle H_{0}^{1}(\Omega )}$ - space of functions in ${\displaystyle H^{1}(\Omega )}$ vanishing on ${\displaystyle \partial \Omega }$

${\displaystyle {\boldsymbol {H}}_{D}^{1}(\Omega )}$ - space of vectors in ${\displaystyle {\boldsymbol {H}}^{1}(\Omega )}$ fulfilling the Dirichlet boundary conditions

${\displaystyle I_{p}}$ - spheric particle's moment of inertia

${\displaystyle K}$ - consistency index

${\displaystyle k}$ - Boltzmann's constant

${\displaystyle K_{B}}$ - Boussinesq–Basset kernel function

${\displaystyle {\mathit {Kn}}}$ - Knudsen number

${\displaystyle L}$ - characteristic length scale of the smallest flow structures

${\displaystyle L^{2}(\Omega )}$ - space of square integrable functions in ${\displaystyle \Omega }$

${\displaystyle L_{\hbox{filter}}}$ - characteristic filter length scale

${\displaystyle {\mathcal {M}}}$ - total momentum

${\displaystyle m_{f}}$ - displaced fluid's mass

${\displaystyle m_{p}}$ - particle mass

${\displaystyle m_{p,f}}$ - equivalent particle mass

${\displaystyle n}$ - behaviour index

${\displaystyle {\mathcal {N}}}$ - total moment of momentum

${\displaystyle n_{\hbox{dim}}}$ - number of space dimensions

${\displaystyle n_{p}}$ - number density

${\displaystyle n_{\hbox{nodes}}}$ - total number of mesh nodes

${\displaystyle N_{q}}$ - number of time steps between quadrature steps

${\displaystyle n_{\hbox{dim}}}$ - number of space dimensions

${\displaystyle N_{\hbox{Vor}}}$ - neighbourhood function

${\displaystyle {\mathcal {P}}(A)}$ - power set of set ${\displaystyle A}$

${\displaystyle p}$ - fluid pressure

${\displaystyle {\mathit {Pe}}_{g}}$ - Péclet number associated to the settling velocity

${\displaystyle {\boldsymbol {q}}}$ - heat flux

${\displaystyle {\mathcal {Q}}}$ - space of pressure scalars

${\displaystyle {\mathcal {Q}}_{h}}$ - discrete counterpart of ${\displaystyle {\mathcal {Q}}}$

${\displaystyle {\boldsymbol {r}}}$ - particle position

${\displaystyle \mathbb {R} }$ - field of the real numbers

${\displaystyle \mathbb {R} ^{+}}$ - strictly positive real numbers

${\displaystyle r}$ - internal energy source/sink

${\displaystyle {\mathit {Re}}}$ - Reynolds number

${\displaystyle {\mathit {Re}}_{\Omega }}$ - rotational Reynolds number

${\displaystyle {\mathit {Re}}_{p}}$ - particle Reynolds number

${\displaystyle {\mathit {Re}}_{p,{\hbox{Shah}}}}$ - Reynolds number for power-law fluids

${\displaystyle {\mathit {Re}}_{p,{\hbox{PL}}}}$ - particle Reynolds number for power-law fluids (Metzner and Reed, 1955)

${\displaystyle {\mathit {Re}}_{s}}$ - Shear-related particle Reynolds number

${\displaystyle \mathbb {R} ^{n}}$ - Cartesian product involving ${\displaystyle n}$ copies of ${\displaystyle \mathbb {R} }$

${\displaystyle {\boldsymbol {S}}}$ - strain rate tensor

${\displaystyle S}$ - settling coefficient

${\displaystyle S^{*}}$ - Brinkman over Oseen length scales squared

${\displaystyle {\mathit {Sc}}}$ - Schmidt number

${\displaystyle S_{L_{\hbox{filter}}}}$ - sphere centred at zero with radius equal to ${\displaystyle L_{\hbox{filter}}}$

${\displaystyle {\mathit {St}}}$ - Stokes number

${\displaystyle {\mathit {St}}_{\hbox{bubble}}}$ - Stokes number for bubble rebound

${\displaystyle {\mathit {St}}_{\hbox{coll}}}$ - collisional Stokes number

${\displaystyle {\mathit {St}}_{\Omega }}$ - rotational Stokes number

${\displaystyle {\mathit {St}}_{L}}$ - scale-dependent Stokes number

${\displaystyle {\boldsymbol {t}}_{N}}$ - surface traction on ${\displaystyle \Gamma _{N}}$

${\displaystyle T}$ - characteristic time scale of the smallest flow structures

${\displaystyle {\boldsymbol {\tau _{1}}}}$ - mass balance-related stabilization parameter

${\displaystyle {\boldsymbol {\tau _{1}}}}$ - momentum balance-related stabilization parameter

${\displaystyle {\boldsymbol {\tau _{\Delta t,e}}}}$ - matrix of stabilization parameters

${\displaystyle {\mathcal {T}}_{b}}$ - total torque due to body actions

${\displaystyle {\boldsymbol {T}}_{\hbox{contact}}}$ - sum of contact torques

${\displaystyle t_{\hbox{DEM}}}$ - discrete time variable for the particles

${\displaystyle T_{\hbox{filter}}}$ - characteristic filter time scale

${\displaystyle {\boldsymbol {T}}_{\hbox{hydr}}}$ - hydrodynamic torque

${\displaystyle {\mathcal {T}}_{s}}$ - total torque due to surface actions

${\displaystyle {\boldsymbol {u}}}$ - fluid velocity

${\displaystyle U}$ - characteristic velocity scale of the smallest flow structures

${\displaystyle u}$ - internal energy

${\displaystyle {\mathcal {U}}}$ - total energy

${\displaystyle {\boldsymbol {v}}}$ - particle velocity

${\displaystyle {\boldsymbol {\mathcal {V}}}_{D}}$ - space of velocity vectors

${\displaystyle {\boldsymbol {\mathcal {V}}}_{D,h}}$ - discrete counterpart of ${\displaystyle {\boldsymbol {\mathcal {V}}}_{D}}$

${\displaystyle V_{p}}$ - particle volume

${\displaystyle {\boldsymbol {w}}}$ - slip velocity

${\displaystyle {\mathcal {W}}_{b}}$ - power spent by the body forces

${\displaystyle {\mathcal {W}}_{s}}$ - power spent by the surface forces

${\displaystyle {\boldsymbol {\mathcal {X}}}}$ - space of solutions

${\displaystyle {\tilde {\boldsymbol {\mathcal {X}}}}}$ - space of subscales

${\displaystyle {\tilde {\boldsymbol {\mathcal {X}}}}_{0}}$ - space of subscales that vanish at the boundary

${\displaystyle {\boldsymbol {\mathcal {X}}}_{h}}$ - discrete counterpart of ${\displaystyle {\boldsymbol {\mathcal {X}}}}$

${\displaystyle y_{PL}^{+}}$ - ${\displaystyle y^{+}}$ dimensionless distance for power-law fluids

Greek Symbols

${\displaystyle \alpha _{p,\rho }}$ - particle mass fraction

${\displaystyle \beta _{\hbox{slip}}}$ - Basset slip coefficient

${\displaystyle \chi }$ - configuration (continuum mechanics)

${\displaystyle {\bar {\chi }}}$ - motion of a body (continuum mechanics)

${\displaystyle {\bar {\chi }}_{k}}$ - motion of a body of component ${\displaystyle k}$ (continuum mechanics)

${\displaystyle \chi _{t}}$ - configuration at time ${\displaystyle t}$ (continuum mechanics)

${\displaystyle \chi _{k,t}}$ - configuration of component ${\displaystyle k}$ at time ${\displaystyle t}$ (continuum mechanics)

${\displaystyle \Delta t}$ - continuous phase time step

${\displaystyle \delta t}$ - disperse phase time step

${\displaystyle \varepsilon _{ijk}}$ - Levi–Civita symbol

${\displaystyle \varepsilon _{p}}$ - amplitude of surface irregularities

${\displaystyle \Gamma }$ - radial distribution function

${\displaystyle \gamma }$ - local strain rate

${\displaystyle \Gamma _{D}}$ - subset of ${\displaystyle \partial \Omega }$ where Dirichlet boundary conditions apply

${\displaystyle \Gamma _{N}}$ - subset of ${\displaystyle \partial \Omega }$ where Neumann boundary conditions apply

${\displaystyle \lambda }$ - mean free path of the fluid molecules

${\displaystyle \mu }$ - fluid's dynamic viscosity

${\displaystyle \mu _{0}}$ - yield viscosity

${\displaystyle \nu }$ - fluid's kinematic viscosity

${\displaystyle \nu _{m}}$ - molecular collision frequency

${\displaystyle \Omega }$ - continuum phase domain

${\displaystyle {\boldsymbol {\Omega }}_{p}}$ - particle's angular velocity

${\displaystyle \Omega _{p}}$ - characteristic value of ${\displaystyle {\boldsymbol {\Omega }}_{p}}$

${\displaystyle \Omega _{s}}$ - characteristic scale of the vorticity in pure-shear flow

${\displaystyle {\boldsymbol {\psi }}_{k}}$ - moment exchange of component ${\displaystyle k}$ with the rest of components

${\displaystyle \rho }$ - relative density of the disperse component

${\displaystyle \rho _{f}}$ - fluid's density

${\displaystyle {\boldsymbol {\sigma }}}$ - Cauchy's stress tensor

${\displaystyle \sigma _{m}}$ - tangential accommodation coefficient

${\displaystyle {\boldsymbol {\tau }}}$ - deviator stress tensor

${\displaystyle \tau _{c_{s}}}$ - sonic time scale

${\displaystyle \tau _{D}}$ - molecular diffusion time scale

${\displaystyle \tau _{\nu }}$ - viscous diffusion time scale

${\displaystyle \tau _{p,\Omega }}$ - rotational particle relaxation time

${\displaystyle \tau _{p}}$ - particle response time

${\displaystyle \tau _{w}}$ - shear stress at the wall

${\displaystyle \Theta }$ - dimensionless frequency normalized by ${\displaystyle \nu _{m}}$

${\displaystyle \theta }$ - absolute temperature

${\displaystyle \xi }$ - inverse of the particle mobility

${\displaystyle \xi _{\hbox{U}}}$ - activation coefficient for ${\displaystyle {\boldsymbol {F}}_{U}+{\boldsymbol {F}}_{A}}$ term

Other symbols

${\displaystyle \left\langle \,\cdot \,\right\rangle _{k}}$ - filtering operator over the component ${\displaystyle k}$

${\displaystyle {\overset {k}{\overline {\,\cdot \,}}}}$ - filtering operator per unit component volume over the component ${\displaystyle k}$

${\displaystyle \mathrm {D} /\mathrm {D} t}$ - material derivative of the fluid continuum

${\displaystyle \Delta }$ - Laplacian operator

${\displaystyle \partial A}$ - boundary of set ${\displaystyle A}$

${\displaystyle (\cdot {,}\cdot )}$ - inner product

${\displaystyle \langle \cdot {,}\cdot \rangle }$ - integral of the product of two functions

${\displaystyle {\overset {k}{\,.\,}}}$ - ${\displaystyle k}$-component variable over ${\displaystyle k}$-volume fraction

Acronyms / Abbreviations

ALE - arbitrary Lagrangian-Eulerian method

ALU - arithmetic logic unit

ASGS - algebraic sub-grid 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 - fluid-structure 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 - particle-in-cell

PID - particle impact drilling

PPC - particles per cell

PPR - recovery method of Zhang and Naga (2005)

PSM - pseudo-spectral method

RVE - representative elemental volume

SC - standard case

SM - streaming multiprocessors

VMS - Variational multiscale (method)

# 1 Introduction

## 1.1 Background and Motivation

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 straight-circumference arch-straight, 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 slow-moving 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 Tennis-ball 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 fluid-particle 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 action-reaction 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 fast-moving particles, in the so-called Newton-drag 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 ground2.

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, Saint-Venant (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 number3. 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 Navier-Stokes 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 non-uniformity 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 particle-laden 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 mid-eighteenth 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 squared-velocity drag law [325].

(3) This still did not resolve d'Alembert's paradox, since the fact that the Navier-Stokes 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.

### 1.1.1 Numerical simulations

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 particle-laden 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 fluid-structure interaction problems.

In fact, researchers are increasingly following this route, producing fully-resolved 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 body-fitted mesh (${\textstyle 1x10^{3}}$); in [273] the fluidization of a bed of spheres was studied numerically and experimentally (${\textstyle 1.204x10^{3}}$); in [138] the lattice-Boltzmann method was used to simulate spheres setling in a stationary fluid (${\textstyle 6.4x10^{3}}$); while in [349] an embedded-body approach was used to study particle clustering in turbulent flows (${\textstyle \approx 2x10^{4}}$). 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 non-Newontian 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 particle-laden 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 low-density 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 cut-down in the computational cost.

A quick calculation illustrates this point very convincingly: First, take a suspension of ${\textstyle N}$ 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, coarse-grained 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 ${\textstyle (10N)^{3}}$ in the first case and ${\textstyle (0.1N)^{3}}$ in the second case. That is a factor 1x106 in number of cells between the two approaches, which amounts to a significantly greater number in terms of computational resources1.

Note that such strategy corresponds to a mixed-scale 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 so-called Euler–Euler (EE) method [180], also called the Eulerian-Eulerian or two-fluid 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 (large-scale) 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 finer-scale details, and that must be closed (expressed as a function of the averaged variables) to obtain a well-posed 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 particle-laden flows for many years; especially in the chemical [332,139,192] and the nuclear [196,158] industries.

In fact, EE-type 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 particle-based 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 size-number 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 2x1012-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 particle-based 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 coarse-grained description, where a smaller-than-real number of (larger) particles is used instead of the one-to-one 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 coarse-graining, hybrid methods and, in particular, CFD-DEM methods retain an interest due to their many advantages over exclusively continuum-based methods:

• One avoids the complications associated with a continuum-based description of granular matter, with still several open theoretical questions [329,297,22]; despite recent promising advances, mainly in the area of dry granular matter [163,110].
• The non-uniform granulometry of real materials can be reproduced naturally, while the same is extremely difficult to include in a model. Often, the only practical alternative is to simplify the granulometry into a finite number of fixed sizes with different proportions and include as many phases as different sizes are present, with the associated multiplicity of laws and computational resources required. This point is actually included in the point above but deserves special mention.
• A lower degree of empiricism is expected, thanks to the lower-scale description of the solid phase.
• The dynamics of the particles are not averaged-over, so that the actual mechanisms that dictate the behaviour of the particles can be understood directly, providing a more valuable qualitative picture that does not need to be reconstructed (which would be nontrivial or impossible) as a post-process. This advantage is particularly important close to the domain boundaries, where averaged methods tend to do most poorly, due to the sharp transitions.
• The consideration of additional effects such as new types of interaction laws (e.g.electrostatic interaction) is natural and does not require the overall modification of the models in place, but usually only a straightforward additive term.
• One can take advantage of existing programs (e.g.DEM programs) and modify them in a simple way to adapt them to the interaction effects with the surrounding fluid.

So, although the number of real-world 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 continuum-only 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 EE-type equations and to understand the full-size objects of study through the analysis of simplified, small-size 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 latest-generation 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

### 1.1.2 Industrial Doctorate

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 particle-laden 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, DEM-related 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 DEM-related 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.

## 1.2 Objectives and methodology

The objectives of the research reported in this document are the following:

1. To develop an algorithm that combines the discrete element method and the finite element method to simulate particle-laden flows with the following list of requisites:
• capability of dealing with a wide range of regimes, including the possibility to have regions with dense and dispersed suspensions simultaneously
• use of the finite element method to discretize the fluid
• use of the discrete element method to model the particles
2. To study the range of applicability of the Maxey–Riley equation as a model for the motion of the individual particles submerged in a fluid, improving the current knowledge on the subject and generating, where possible, practical estimates of direct application to numerical modelling.
3. To study current alternatives for the numerical treatment of the history term in the equation of motion and compare them.
4. To improve on the method of quadrature (the history term involves a time integral, as will be explained later) proposed by van Hinsberg [353] and provide a detailed study of its efficiency and accuracy, providing convincing evidence that it is not necessary to neglect this term to have an efficient numerical method.
5. To report an account of relevant application examples of the proposed strategy with interest to the industry, as well as of the different technologies developed for their particular requirements.
6. To generalize a stabilized finite element method and use it to discretize the backward-coupled flow equations.
7. To develop a suitable inter-phase coupling strategy.

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 Multi-physics [88], or Kratos for short. Kratos is based on C++, with a Python-based external layer. It is organized according to a marked object-oriented 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 Eulerian-description, finite element method-based CFD simulator and it is in charge of solving the fluid equations of motion. DEMApplication is a general-purpose 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 inter-phase coupling tools, the hydrodynamic calculation algorithms, the modified fluid finite elements necessary to solve the modified backward-coupled 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 download2. FluidDynamicApplication has an interface usable from the Pre-/Post-processor 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 fluid-solving 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).

## 1.3 Outline of this document

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 well-established 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 state-of-the-art and compare several alternatives compatible to the finite element method. The description of the particles-fluid 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 T-shaped pipe junctions, which was only recently first studied by Vigolo et al. [356]. This test exemplifies the use of the one-way coupled strategy, with no inter-particle interactions, representative of low-Reynolds number internal flows with low-density 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 one-way coupled strategy for the fluid-particles coupling, although this time we consider the inter-particle 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 backward-coupled 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 Geldard-D 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 backward-fluid 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.

# 2 The Maxey–Riley equation

## 2.1 Introduction

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 point-mass, 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 second-order 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 warm-cloud 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 finite-sized 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 ${\textstyle {\boldsymbol {u}}:{\mathcal {D}}\times \mathbb {R} ^{+}\rightarrow \mathbb {R} ^{n}}$ define the background fluid velocity field, where ${\textstyle {\mathcal {D}}\subseteq \mathbb {R} ^{n}}$ is some spatial open domain contained in the fluid and ${\textstyle n}$ is the ambient spatial dimension (either two or three). Let ${\textstyle {\boldsymbol {r}}(t)\in {\mathcal {D}}}$ be the trajectory of a particle through this fluid, and ${\textstyle {\boldsymbol {v}}}$ its velocity. The MRE states:

 ${\displaystyle m_{p}{\frac {\mathrm {d} {\boldsymbol {v}}}{\mathrm {d} t}}=\underbrace {m_{f}{\frac {\mathrm {D} {\boldsymbol {u}}}{\mathrm {D} t}}} _{{\boldsymbol {F}}_{U}}+\underbrace {{\frac {1}{2}}m_{f}\left({\frac {\mathrm {D} }{\mathrm {D} t}}\left({\boldsymbol {u}}+{\frac {1}{10}}a^{2}\nabla ^{2}{\boldsymbol {u}}\right)-{\frac {\mathrm {d} {\boldsymbol {v}}}{\mathrm {d} t}}\right)} _{{\boldsymbol {F}}_{A}}+\underbrace {6\pi a\mu \left({\boldsymbol {u}}-{\boldsymbol {v}}+{\frac {1}{6}}a^{2}\nabla ^{2}{\boldsymbol {u}}\right)} _{{\boldsymbol {F}}_{D}}}$ ${\displaystyle +\underbrace {6a^{2}{\sqrt {\pi \rho _{f}\mu }}{\frac {\mathrm {d} }{\mathrm {d} t}}\int _{t_{0}}^{t}{{\frac {1}{\sqrt {t-\tau }}}\left({\boldsymbol {u}}-{\boldsymbol {v}}+{\frac {1}{6}}a^{2}\nabla ^{2}{\boldsymbol {u}}\right)\,{d}\tau }} _{{\boldsymbol {F}}_{H}}+\underbrace {(m_{f}-m_{p}){\boldsymbol {g}}} _{{\boldsymbol {F}}_{W}}}$
(2.1)

where ${\textstyle m_{p}}$ is the mass of the particle, ${\textstyle m_{f}}$ is the mass of the displaced fluid volume, ${\textstyle a}$ is the particle radius, ${\textstyle \rho _{f}}$ is the density of the fluid, ${\textstyle \mu }$ is the dynamic viscosity of the fluid (${\textstyle \nu =\mu /\rho _{f}}$ is the kinematic viscosity) and ${\textstyle {\boldsymbol {g}}}$ the acceleration due to gravity. The notation ${\textstyle \mathrm {D} /\mathrm {D} t}$ refers to the material derivative of the fluid. Eq.~2.1, together with ${\textstyle {\boldsymbol {v}}=\mathrm {d} {\boldsymbol {v}}/\mathrm {d} t}$ and the initial conditions ${\textstyle {\boldsymbol {r}}(t_{0})={\boldsymbol {r}}_{0}}$ and ${\textstyle {\boldsymbol {v}}(t_{0})={\boldsymbol {v}}_{0}}$ 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 (${\textstyle {\boldsymbol {F}}_{U}}$), the added mass or virtual mass term (${\textstyle {\boldsymbol {F}}_{A}}$), the Stokes drag term (${\textstyle {\boldsymbol {F}}_{D}}$), the Boussinesq–Basset history term (${\textstyle {\boldsymbol {F}}_{H}}$) and the term due to the weight of the particle minus its (hydrostatic) buoyancy (${\textstyle {\boldsymbol {F}}_{W}}$).

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 second-order 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).

## 2.2 Range of validity

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.

 ${\displaystyle {\frac {\mathrm {d} {\boldsymbol {v}}}{\mathrm {d} t}}={\frac {1}{\rho }}{\frac {\mathrm {D} {\boldsymbol {u}}}{\mathrm {D} t}}+{\frac {1}{2\rho }}\left({\frac {\mathrm {D} {\boldsymbol {u}}}{\mathrm {D} t}}-{\frac {\mathrm {d} {\boldsymbol {v}}}{\mathrm {d} t}}+{\frac {9}{10}}{\frac {1}{2\rho {+1}}}{\frac {1}{{\mathit {St}}{\mathit {Re}}}}{\frac {\mathrm {d} \nabla ^{2}{\boldsymbol {u}}}{\mathrm {d} t}}\right)}$ ${\displaystyle +{\frac {2\rho +1}{2}}{\mathit {St}}\left({\boldsymbol {u}}-{\boldsymbol {v}}+{\frac {3}{2}}{\frac {1}{2\rho {+1}}}{\frac {1}{{\mathit {St}}{\mathit {Re}}}}\nabla ^{2}{\boldsymbol {u}}\right)}$ ${\displaystyle +{\frac {3}{2{\sqrt {\mathrm {\pi } }}}}{\sqrt {(2\rho {+1)}{\mathit {St}}}}\int _{t_{0}}^{t}{\frac {1}{\sqrt {t-\tau }}}\left({\boldsymbol {u}}-{\boldsymbol {v}}+{\frac {3}{2}}{\frac {1}{2\rho {+1}}}{\frac {1}{{\mathit {St}}{\mathit {Re}}}}\nabla ^{2}{\boldsymbol {u}}\right)\,{d}\tau +(1-\rho ){\frac {1}{{\mathit {Fr}}^{2}}}}$
(2.2)

where

 ${\displaystyle \rho ={\frac {\rho _{p}}{\rho _{f}}}\quad {\hbox{ (relative density)}}}$ ${\displaystyle {\mathit {Re}}={\frac {UL}{\nu }}\quad {\hbox{ (Reynolds number of the ambient flow)}}}$ ${\displaystyle {\mathit {St}}={\frac {2\rho +1}{9}}{\frac {a^{2}}{\nu }}{\frac {1}{T}}\quad {\hbox{ (Stokes number)}}}$ ${\displaystyle {\mathit {Fr}}={\frac {U}{\sqrt {gL}}}\quad {\hbox{ (Froude number)}}}$
(2.3)

where ${\textstyle g:=\left\|{\boldsymbol {g}}\right\|}$ and ${\textstyle U}$, ${\textstyle L}$ and ${\textstyle T}$ are the dimensional scalars by which ${\textstyle {\boldsymbol {u}}}$ (and ${\textstyle {\boldsymbol {v}}}$), ${\textstyle {\boldsymbol {x}}}$ and ${\textstyle t}$ have been normalized. These quantities are conventionally defined such that ${\textstyle L}$ is a characteristic length scale of the flow and ${\textstyle U}$ is defined such that ${\textstyle U/L}$ is the characteristic magnitude of the gradient of the unperturbed velocity near the particle location, while ${\textstyle T}$ is defined such that ${\textstyle UT=L}$. By the term unperturbedwe refer to the flow resulting from subtracting the (Stokes) flow produced by the presence of the particle under the same far-field 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 ${\textstyle \tau _{p}:=T{\mathit {St}}}$, we obtain the relation

 ${\displaystyle {\mathit {St}}={\frac {\tau _{p}}{L}}}$
(2.4)

The quantity ${\textstyle \tau _{p}}$ 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 ${\textstyle 1-1/\mathrm {e} }$ (where ${\textstyle \mathrm {e} =2.71\dots }$) parts of the surrounding fluid velocity, under the action of a constant, uniform flow (and neglecting ${\textstyle {\boldsymbol {F}}_{H}}$). 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

 ${\displaystyle {\mathit {Re}}_{p}={\frac {aW}{\nu }}}$
(2.5)

where ${\textstyle W}$ is the characteristic magnitude of the relative (or slip) velocity ${\textstyle \left\|{\boldsymbol {v}}-{\boldsymbol {u}}\right\|}$. The quantity ${\textstyle {\mathit {Re}}_{p}}$ 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]

 ${\displaystyle {\frac {a{\frac {a}{L}}U}{\nu }},\,{\mathit {Re}}_{p},\,{\frac {a}{L}}\rightarrow 0}$
(2.6)

In practice, these asymptotic relations are interpreted as requiring the left-hand sides to be much smaller than one 1, that is

 ${\displaystyle {\frac {a{\frac {a}{L}}U}{\nu }}\ll 1}$ (2.7) ${\displaystyle {\mathit {Re}}_{p}\ll 1}$ (2.8) ${\displaystyle {\frac {a}{L}}\ll 1}$ (2.9)

Note that the LHS of Eq.~2.7 can be understood as a sort of downscaled version of the ambient ${\textstyle {\mathit {Re}}}$ to the size of the particle, since the velocity variation seen by the particle scales as ${\textstyle a/LU}$. 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 ${\textstyle F_{U}}$, 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 ${\textstyle a/L}$ (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.7-2.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

• What kind of effects are likely to appear first when moving away from its range of applicability?
• At what approximate pace are these effects expected to grow?

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 ${\textstyle \rho }$ 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 1x101 for any commonly encountered materials (for a suspension of osmium particles in gasoline one has ${\textstyle \rho \approx 30}$, 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 ${\textstyle {\mathit {St}}}$ (typically, small bubbles and low Reynolds numbers lead to spherical bubbles; see [47] for an in-depth 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

 ${\displaystyle {\mathit {St}}={\frac {2\rho +1}{9}}{\frac {a^{2}}{\nu }}{\frac {1}{T}}\sim {\frac {2\rho +1}{9}}{\frac {a^{2}}{\nu }}{\frac {U}{L}}<\max {\left\{1,{\frac {2\rho +1}{90}}{\frac {a}{\nu }}U\right\}}\lesssim }$ ${\displaystyle \lesssim \max {\left\{1,{\frac {2\rho +1}{90}}{\frac {a}{\nu }}W\right\}}=\max {\left\{1,{\frac {2\rho +1}{90}}{\mathit {Re}}_{p}\right\}}}$
(2.10)

where we have interpreted Eq.~2.9 as ${\textstyle a<0.1L}$ and used the fact that for very large ${\textstyle {\mathit {St}}}$ 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 ${\textstyle {\mathit {St}}}$ by the particle Reynolds number and in many cases (e.g.water suspended mineral particles) by an additional order of magnitude, so that if ${\textstyle {\mathit {Re}}_{p}<0.1}$, then ${\textstyle {\mathit {St}}<0.01}$, speaking in order-of-magnitude terms. In atmospheric air, a similar argument would lead to bound ${\textstyle {\mathit {St}}}$ below two to three orders of magnitude above ${\textstyle {\mathit {Re}}_{p}}$. Using the more conservative ${\textstyle a<0.01L}$ 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.7-2.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 ${\displaystyle x\ll 1}$ as meaning at least two orders of magnitude smaller than one or ${\displaystyle x<0.01}$ [197].

### 2.2.1 Inertial effects: first finite-Reₚ effects

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 no-slip boundary conditions on the particle surface and the far-field velocity conditions at infinity and expanding the stream function in powers of ${\textstyle {\mathit {Re}}_{p}}$, the hydrodynamic force on the particle can be calculated, to leading ${\textstyle {\mathit {Re}}_{p}}$ 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 ${\textstyle {\mathit {Re}}_{p}={\mathcal {O}}(1)}$, 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, re-derived Oseen's inertial, first-order correction to the steady drag force, corrected a flaw in Oseen's original reasoning and added an additional term to the expansion:

 ${\displaystyle {\boldsymbol {F}}_{D}=6\mathrm {\pi } a\mu \left(1+{\frac {3}{8}}{\mathit {Re}}_{p}+{\frac {9}{40}}{\mathit {Re}}_{p}^{2}\log {{\mathit {Re}}_{p}}+{\mathcal {O}}({\mathit {Re}}_{p}^{2})\right)\left({\boldsymbol {u}}-{\boldsymbol {v}}\right)}$
(2.11)

This way, the contribution from the outer region becomes a correction of the steady drag force, whose leading term, given by ${\textstyle 3/8{\mathit {Re}}_{p}}$, 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 ${\textstyle {\mathit {Re}}_{p}\approx 1.4x10^{-1}}$, and by this point the error in this expression is only ${\textstyle \approx 3x10^{-2}}$, taking as a reference the empirical drag coefficient by Clift et al. (1978) (which has a root mean square error of ${\textstyle \approx 1x10^{-2}}$ with respect to the empirical data gathered by Brown and Lawler [50] covering the range ${\textstyle 0<{\mathit {Re}}_{p}<4x10^{3}}$). This means that its range of validity is large enough to estimate the first-effects 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 ${\displaystyle C_{D}=2\left\|{\boldsymbol {F}}_{D}\right\|/(\mathrm {\pi \rho _{f}\mu } )}$.

The above expression of the force is only valid for long-term 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 non-zero 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 ${\textstyle {\mathit {Re}}_{p}}$. The hydrodynamic force on the particle is then computed to first order in ${\textstyle {\mathit {Re}}_{p}}$ 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 fully-resolved numerical simulations [243] that higher rate of convergence of the particle velocity to terminal conditions (i.e.it varies but is ${\textstyle ~t^{-2}}$ for smoothly accelerating motions where the slip velocity tends to increase) than that predicted by the MRE equation, which is ${\textstyle ~t^{-1/2}}$, 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 short-time-scale 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 so-called 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 ${\textstyle \tau _{\nu }=a^{2}/\nu }$, 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]

 ${\displaystyle \tau _{p}\ll {\frac {\nu }{U^{2}}}}$
(2.12)

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

 ${\displaystyle {\mathit {Re}}_{p}^{2}\ll {\frac {9}{2\rho {+1}}}}$
(2.13)

As the value of ${\textstyle {\mathit {Re}}_{p}}$ 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 ${\textstyle {\mathit {Re}}_{p}\ll 1}$, 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 non-uniformity 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 ${\textstyle w}$, 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 ${\textstyle {\mathit {Re}}_{p}}$, for the lift force that an isolated, non-rotating particle experiments under constant shear (linear non-uniformity of the flow).

 ${\displaystyle {\boldsymbol {F}}_{L,{\hbox{Saff}}}=6.45a^{2}{\sqrt {\frac {\mu _{f}\rho _{f}}{\left\|\nabla \times {\boldsymbol {u}}\right\|}}}\left(\nabla \times {\boldsymbol {u}}\right)\times {\boldsymbol {w}}}$
(2.14)

where ${\textstyle {\boldsymbol {w}}={\boldsymbol {v}}-{\boldsymbol {u}}}$ is the slip velocity. This result is valid under the restrictions ${\textstyle {\mathit {Re}}_{p}\ll 1}$, ${\textstyle {\mathit {Re}}_{s}\ll 1}$ and ${\textstyle {\mathit {Re}}_{p}\ll {\sqrt {{\mathit {Re}}_{s}}}}$; where the shear Reynolds number is defined by

 ${\displaystyle {\mathit {Re}}_{s}:={\mathit {Re}}_{p}{\frac {a\left\|\nabla \times {\boldsymbol {u}}\right\|}{\left\|{\boldsymbol {w}}\right\|}}}$
(2.15)

that is, the particle Reynolds number multiplied by a dimensionless shear rate gradient.

Unfortunately, the typical values for ${\textstyle {\mathit {Re}}_{s}}$ in turbulence are often smaller than those of ${\textstyle {\mathit {Re}}_{p}}$ [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-${\textstyle {\mathit {Re}}_{p}}$ analytical formula was provided by Herron et al. [164] and reads

 ${\displaystyle {\boldsymbol {F}}_{L,{\hbox{Herr}}}=8.6a^{2}{\sqrt {\frac {\mu _{f}\rho _{f}}{\left\|\nabla \times {\boldsymbol {u}}\right\|}}}\left(\nabla \times {\boldsymbol {u}}\right)\times {\boldsymbol {w}}}$
(2.16)

which has a coefficient about 2x101 larger compared to ${\textstyle {\boldsymbol {F}}_{L,{\hbox{Saff}}}}$. This expression is subject to the restriction ${\textstyle {\mathit {Re}}_{p}\ll 1}$. Note that ${\textstyle {\boldsymbol {F}}_{L,{\hbox{Saff}}}}$ and ${\textstyle {\boldsymbol {F}}_{L,{\hbox{Herr}}}}$ 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 solid-body 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 history-dependent contributions. Those have been studied in [60,55] for pure shear (generalizing Saffman's result) and solid-body 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.

 ${\displaystyle {\frac {F_{L,{\hbox{Saff}}}}{F_{D}}}\lesssim {\frac {6.45a^{2}\rho _{f}\nu ^{1/2}\Omega _{s}^{1/2}W}{6\mathrm {\pi } a\mu W}}\sim {\frac {1}{3}}{\mathit {Re}}_{s}}$
(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 first-effects of the force modification due to non-zero particle angular velocity were calculated by Rubinow and Keller [294], to zeroth-order accuracy in ${\textstyle {\mathit {Re}}_{p}}$. 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

 ${\displaystyle {\boldsymbol {F}}_{L,{\hbox{Rub}}}=\mathrm {\pi } a^{3}\rho _{f}{\boldsymbol {\Omega }}_{p}\times {\boldsymbol {w}}(1+O({\mathit {Re}}_{p}))}$
(2.18)

where ${\textstyle {\boldsymbol {\Omega }}_{p}}$ 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

 ${\displaystyle {\mathit {Re}}_{\Omega }:={\frac {a^{2}\Omega _{p}}{\nu }}}$
(2.19)

where ${\textstyle \Omega _{p}=\left\|{\boldsymbol {\Omega }}_{p}\right\|}$. Note that ${\textstyle {\mathit {Re}}_{\Omega }=a/L{\mathit {Re}}_{p}}$, 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 first-order rotation-induced effect of inertia is a lift force with a magnitude of order ${\textstyle {\mathit {Re}}_{\Omega }}$ [294], the effect results of order ${\textstyle a^{2}U/(L\nu )}$; more precisely, we have

 ${\displaystyle {\frac {F_{L,{\hbox{Rub}}}}{F_{D}}}\sim {\frac {\mathrm {\pi } a^{3}\rho _{f}\Omega _{p}W}{6\mathrm {\pi } a\mu W}}\sim {\frac {\mathrm {\pi } a^{3}\rho _{f}U/LW}{6\mathrm {\pi } a\mu W}}={\frac {1}{6}}{\mathit {Re}}_{\Omega }}$
(2.20)

where ${\textstyle U/L}$ give the estimation of local shear rate; ${\textstyle U}$ and ${\textstyle L}$ 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 history-dependent terms that decay much faster than the translational analogues. In the steady-state limit, the equation reads

 ${\displaystyle I_{p}{\frac {\mathrm {d} {\boldsymbol {\Omega }}_{p}}{\mathrm {d} x}}=-8\mathrm {\pi } a^{3}\mu {\boldsymbol {\Omega }}_{p}}$
(2.21)

This equation allows us to calculate the rotational relaxation time as

 ${\displaystyle \tau _{p,\Omega }={\frac {I_{p}}{8\mathrm {\pi } a^{3}\mu }}={\frac {\rho a^{2}}{15\nu }}}$
(2.22)

from which a rotational Stokes number can be defined (assuming the steady-state forces are dominant)

 ${\displaystyle {\mathit {St}}_{\Omega }={\frac {\tau _{p,\Omega }}{T}}={\frac {1}{15}}\rho {\frac {a^{2}}{\nu }}{\frac {1}{T}}\sim {\mathit {St}}}$
(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 order-of-magnitude analysis the present formulation will suffice, especially since the rate at which the steady-state formulation becomes inaccurate is quite slow, and it still yields reasonable results at ${\textstyle {\mathit {Re}}_{\Omega }}$ 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 ${\displaystyle {\frac {1}{2}}\Omega _{p}^{2}I\propto \Omega _{p}^{2}a^{5}}$ and the translational kinetic energy is ${\displaystyle {\frac {1}{2}}\left\|{\boldsymbol {v}}\right\|^{2}m_{p}\propto \left\|{\boldsymbol {v}}\right\|^{2}a^{3}}$. Thus, after a collision we have ${\displaystyle \Omega _{p}\propto \left\|{\boldsymbol {v}}\right\|/a\propto U/a}$, for light and/or small particles. Furthermore, this gives ${\displaystyle {\mathit {Re}}_{\Omega }\sim aU/L\sim {\mathit {Re}}_{p}U/W}$; and the ratio ${\displaystyle U/W}$ 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 second-order 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]

 ${\displaystyle {\boldsymbol {F}}_{D}={\frac {3\mu }{2a}}\int _{S_{p}}{\boldsymbol {u}}({\boldsymbol {x}})-{\boldsymbol {v}}({\boldsymbol {x}}_{p})\,\mathrm {d} A\approx 6\mathrm {\pi } a\mu \left({\boldsymbol {u}}({\boldsymbol {x}})-{\boldsymbol {v}}({\boldsymbol {x}}_{p})+{\frac {a^{2}}{6}}\nabla ^{2}{\boldsymbol {u}}+{\frac {a^{4}}{120}}\nabla ^{4}{\boldsymbol {u}}\right)}$
(2.24)

Where the ${\textstyle \nabla ^{2k}}$ is the ${\textstyle k}$-th successive composition of the vector Laplacian.

Now, for low ${\textstyle {\mathit {Re}}_{p}}$, the fluid is expected to present only long-wave 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 ${\textstyle \nabla ^{2}{\boldsymbol {u}}}$ across the particle diameter. Then, an adequate order-of-magnitude estimate of its value is given by (see, e.g. [245] for a discussion on the estimation of order of magnitude scales)

 ${\displaystyle \left\|{\frac {a^{4}}{120}}\nabla ^{4}{\boldsymbol {u}}\right\|\sim {\frac {a^{4}}{120}}{\frac {\left\|\nabla ^{2}{\boldsymbol {u}}\right\|}{2a^{2}}}={\frac {1}{40}}\left\|{\frac {a^{2}}{6}}\nabla ^{2}{\boldsymbol {u}}\right\|}$
(2.25)

So that the next order correction will be significantly smaller than the classical Faxén correction up to the point when ${\textstyle {\mathit {Re}}_{p}}$ is not small any more (in the limit of quasi-steady Stokes flow it holds ${\textstyle \nabla ^{4}{\boldsymbol {u}}=0}$ exactly). This means that in practice the small ${\textstyle {\mathit {Re}}}$ requirement will guarantee that the Faxén corrections will make a good enough job and thus the restriction Eq.~2.9, that is, ${\textstyle a/L\ll 1}$ by itself does not appear to be restrictive, but rather indirectly through the violation of ${\textstyle {\mathit {Re}}_{p}\ll 1}$. 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 ${\textstyle a\sim 4L}$ (where ${\textstyle L}$ is the Kolmogorov microscale). From then on, inertial (finite ${\textstyle {\mathit {Re}}_{p}}$) effects quickly kick in.

### 2.2.3 Nonsphericity effects

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 ${\textstyle {\mathit {Re}}_{p}}$; 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, ${\textstyle \varepsilon }$, 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 higher-order correction and can thus be neglected to first-order approximation.

Specifically, let the particle's surface ${\textstyle S_{p}}$ be defined in polar coordinates around its center as the set of points with coordinates ${\textstyle r}$, ${\textstyle \theta }$ and ${\textstyle \phi }$ such that

 ${\displaystyle {\frac {r}{a}}=1+\varepsilon _{p}f(\theta ,\phi )}$
(2.26)

where ${\textstyle a}$ is the radius of the equivalent-volume sphere.

Their equation for the hydrodynamic force on the translating particle reads:

 ${\displaystyle {\boldsymbol {F}}=\left({\boldsymbol {I}}-\varepsilon _{p}{\boldsymbol {H}}\right){\boldsymbol {F}}_{D}+\left({\boldsymbol {I}}-2\varepsilon _{p}{\boldsymbol {H}}\right){\boldsymbol {F}}_{H}+\left({\boldsymbol {I}}-\varepsilon _{p}{\boldsymbol {H}}\right){\boldsymbol {F}}_{A}}$
(2.27)

where

 ${\displaystyle {\boldsymbol {H}}={\frac {3}{8\mathrm {\pi } }}\int _{S_{1}}f{\boldsymbol {n}}{\boldsymbol {n}}\,\mathrm {d} S}$
(2.28)

and the integral is over the unit sphere ${\textstyle S_{1}}$, after applying a uniform scaling by its radius over the whole space. Note that ${\textstyle {\boldsymbol {n}}}$ 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 star-shaped with respect to its center of mass1. 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 ${\textstyle \left\|\,\cdot \,\right\|}$, based on the usual Euclidean norm for vectors, denoted with the same symbol. Then, for example, the corrected drag force ${\textstyle {\boldsymbol {F}}_{D,\varepsilon _{p}}=\left({\boldsymbol {I}}-\varepsilon _{p}{\boldsymbol {H}}\right){\boldsymbol {F}}_{D}}$ fulfils

 ${\displaystyle {\frac {\left\|{\boldsymbol {F}}_{D,\varepsilon _{p}}-{\boldsymbol {F}}_{D}\right\|}{\left\|{\boldsymbol {F}}_{D}\right\|}}\leq \left\|\varepsilon _{p}{\boldsymbol {H}}\right\|={\frac {3\varepsilon _{p}}{8\mathrm {\pi } }}\left\|\int _{S_{p}}f{\boldsymbol {n}}{\boldsymbol {n}}\,\mathrm {d} S\right\|\leq {\frac {3\varepsilon _{p}}{8\mathrm {\pi } }}\max _{S_{p}}{f}{\frac {4}{3}}\mathrm {\pi } }$ ${\displaystyle ={\frac {\varepsilon _{p}}{2}}\max _{S_{p}}{f}\leq \varepsilon _{p}}$
(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 steady-state motions, the relative error made by taking the Stokes expression with the average diameter is smaller than

 ${\displaystyle {\mathcal {O}}({\frac {a_{\mathrm {circ} }-a_{\mathrm {insc} }}{2a_{\mathrm {insc} }}})}$
(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 star-shaped 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 (non-star shaped) rough sphere.

### 2.2.4 Effect of neighbouring particles

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 inter-particle influences:

• Mid/long-range hydrodynamic interactions, due to the flow disturbances generated by the presence of all the other particles
• Short-range forces of various types that appear when the inter-particle gaps become very small, such as lubrication forces and inter-particle contact forces

Starting with the long-range 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 low-Reynolds 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 (${\textstyle a/{\mathit {Re}}_{p}}$) [248] from the particle's center, i.e.

 ${\displaystyle {\boldsymbol {u}}=\left\|{\boldsymbol {w}}\right\|{\begin{pmatrix}\cos \theta \left(1-{\frac {3}{2}}{\frac {1}{r}}+{\frac {1}{2}}\left({\frac {1}{r}}\right)^{3}\right)\\\sin \theta \left(1-{\frac {3}{4}}{\frac {1}{r}}-{\frac {1}{4}}\left({\frac {1}{r}}\right)^{3}\right)\\\sin \theta \left(1-{\frac {3}{4}}{\frac {1}{r}}-{\frac {1}{4}}\left({\frac {1}{r}}\right)^{3}\right)\end{pmatrix}}}$
(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 ${\textstyle r}$ the radial coordinate (normalized by the sphere's radius) and ${\textstyle \theta }$ 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 ${\textstyle 2^{k}}$-th neighbour, as a function of ${\textstyle \alpha _{p}}$, is shown for reference; for ${\textstyle k=0{..}5}$ 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 finite-size particles and is robust for all values of ${\textstyle \alpha _{p}}$ up to close packing.

 Figure 3: Modulus of the flow disturbance around a particle (black circle), calculated according to Eq.  and normalized by the modulus of the slip velocity, ${\displaystyle \left\|{\boldsymbol {w}}\right\|}$. All distances have been normalized by the radius of the particle, ${\displaystyle a}$.
 Figure 4: Average distance (normalized by the radius of the sphere) to ${\displaystyle k}$-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 multi-body in general. In the dilute limit that we are interested in (we are studying the first-effects of the presence of neighbours), it is, however, possible to give the interactions a pair-wise 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 short-range interactions, these include two types of forces: the shortest-range 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, short-range forces that can be naturally included in a contact model. Nonetheless, in some situations it may still be interesting to neglect them completely as:

1. The complicated modelling issue associated with the two types of short-range interactions is avoided.
2. The computational cost is greatly reduced, as the most expensive parts of the DEM algorithm (force calculation and search) are avoided.
3. The time step of the simulation, which is otherwise upper bounded by a fraction of the contact duration, can then be much increased, further reducing the cost.
4. It becomes easier to give a completely statistical interpretation to the particles, making it possible to alter the real concentration of particles to speed up the simulations, avoiding for instance the need to consider a realistic number of particles when that is too costly.

We wish to characterize the conditions for which these effects start to become important, causing the MRE to start loosing accuracy in modelling particle-laden flows. Such characterization has been attempted before  [116,118] in the context of general particle-laden 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 three-way interactions become important is

 ${\displaystyle \alpha _{p}\lesssim 10^{-3}}$
(2.32)

where ${\textstyle \alpha _{p}}$ is the local solid fraction. This corresponds to an average nearest-neighbour separation of about 1x101 particle radii. This is an oft-cited 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 order-of-magnitude estimates.

Before proceeding, we must make one further remark concerning the scope of the discussion. In some particle-laden 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 long-range 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 order-of-magnitude 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 first-order 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

 ${\displaystyle {\frac {V_{\hbox{settling}}-V_{0}}{V_{0}}}=S\alpha _{p}}$
(2.33)

where ${\textstyle U_{0}}$ is the terminal velocity of the target particle when isolated, ${\textstyle V_{\hbox{settling}}}$ the actual settling velocity of the particle in the suspension, ${\textstyle \alpha _{p}:=1-\alpha _{f}}$ is the solid volume fraction and ${\textstyle S}$ is the settling coefficient, which under the particular assumptions of the theorem is ${\textstyle S=-6.55}$; and ${\textstyle \alpha _{p}}$. 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 ${\textstyle S}$. For instance, ${\textstyle S\approx -2.5}$ 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 (${\textstyle S}$ 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 non-uniformity 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 first-order 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, ${\textstyle {\frac {\sqrt {3}}{2}}\phi ^{-1/2}a}$, as the simulations for higher values of ${\textstyle {\mathit {Re}}}$ were found to be prohibitively expensive, in  [167].

The derivations by Kaneda were performed for simplicity under the hypothesis ${\textstyle \alpha _{p}\sim {\mathit {Re}}_{p}^{2}}$, although the resulting expression is shown to have the correct asymptotic behaviour at the opposite limits ${\textstyle (\alpha _{p},{\mathit {Re}}_{p})\rightarrow (0,{\mathit {Re}}_{p})}$ and ${\textstyle (\alpha _{p},{\mathit {Re}}_{p})\rightarrow (\alpha _{p},0)}$. 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

 ${\displaystyle {\boldsymbol {F}}_{D}^{\hbox{Kaneda}}=C_{K}{\boldsymbol {F}}_{D}}$
(2.34)

with

 ${\displaystyle C_{K}:=\left(1+{\mathit {Re}}_{p}K^{\hbox{Kaneda}}(S^{*})\right)}$
(2.35)

where ${\textstyle S^{*}}$ 1 is a nondimensional constant such that ${\textstyle S^{*}=9/2\alpha _{p}/{\mathit {Re}}_{p}^{2}+o(\alpha _{p})}$ (when ${\textstyle \alpha _{p}\rightarrow 0}$) and where ${\textstyle {\boldsymbol {F}}_{D}}$ is the usual Stokes drag; and

 ${\displaystyle K^{\hbox{Kaneda}}(S^{*})={\frac {3}{8}}\left((2S^{*}+1){\sqrt {(4S^{*}+1)}}-4S^{*2}\log {\frac {{\sqrt {(4S^{*}+1)}}+1}{{\sqrt {(4S^{*}+1)}}-1}}\right)}$
(2.36)

Note that

 ${\displaystyle {\boldsymbol {F}}_{D}^{\hbox{Kaneda}}\sim \left(1+{\frac {3}{8}}{\mathit {Re}}_{p}+{\frac {3}{2}}{\mathit {Re}}_{p}S^{*}\right){\boldsymbol {F}}_{D}}$ ${\displaystyle \sim \left(1+{\frac {3}{8}}{\mathit {Re}}_{p}+{\frac {27}{4}}\alpha _{p}{\mathit {Re}}_{p}^{-1}\right){\boldsymbol {F}}_{D}\quad {\hbox{ for }}S^{*}\rightarrow 0}$ ${\displaystyle {\boldsymbol {F}}_{D}^{\hbox{Kaneda}}\sim \left(1+{\mathit {Re}}_{p}{\sqrt {S^{*}}}+{\frac {3}{40}}{\frac {{\mathit {Re}}_{p}}{\sqrt {S^{*}}}}\right){\boldsymbol {F}}_{D}}$ ${\displaystyle \sim \left(1+{\frac {3}{\sqrt {2}}}\alpha _{p}^{1/2}+{\frac {\sqrt {2}}{40}}{\mathit {Re}}_{p}^{2}\alpha _{p}^{-1/2}\right){\boldsymbol {F}}_{D}\quad {\hbox{ for }}S^{*}\rightarrow +\infty }$
(2.37)

Which, as stated, recover the first order correction due to Oseen for the drag force at ${\textstyle S^{*}\rightarrow 0}$ (and thus ${\textstyle \alpha _{p}\rightarrow 0}$), see Eq.~2.11, and also the classic expression by Brinkman [49], at ${\textstyle {\mathit {Re}}_{p}\rightarrow 0}$.

Figure 5 shows the magnitude of the correction coefficient ${\textstyle C_{K}}$ for a wide range of solid volume fractions, for different values of ${\textstyle {\mathit {Re}}_{p}}$. 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.
 Figure 5: Correction factor of the steady drag force as a function of the average solid fraction for a randomly distributed array of spheres. The different curves are labelled according to their corresponding ${\displaystyle {\mathit {Re}}_{p}}$.

Although Kaneda's result applies to arrays of fixed spheres under a uniform background flow, it clearly reflects the inextricable relationship between the low-Reynolds number and low-concentration hypotheses through the dimensionless number ${\textstyle S^{*}}$. This number can be understood as the quotient between two length scales:

 ${\displaystyle S^{*}={\frac {9}{2}}\left({\frac {L_{\hbox{Oss}}}{L_{\hbox{Brink}}}}\right)^{2}}$
(2.38)

where ${\textstyle L_{\hbox{Oss}}=a/{\mathit {Re}}_{p}}$ is the Oseen length (distance from the particle at which inertia becomes as important as viscosity) and ${\textstyle L_{\hbox{Brink}}=a/{\sqrt {\alpha _{p}}}}$ (distance from the particle at which the Stokesian disturbance caused by it switches from a linear to a cubic-order decay due to the presence of neighbours). The physical interpretation of the two regimes in Eq.~2.36 goes as follows: when ${\textstyle S^{*}}$ is very small (${\textstyle L_{\hbox{Brink}}\gg L_{\hbox{Oss}}}$), inertial effects enter much before the Stokesian disturbance has been screened by other particles and thus their influence is dictated by finite-${\textstyle {\mathit {Re}}_{p}}$ effects. On the contrary, if ${\textstyle L_{\hbox{Brink}}\ll L_{\hbox{Oss}}}$, 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 ${\textstyle S^{*}\ll 1}$, 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 ${\textstyle L_{\hbox{cloud}}/L_{\hbox{Oss}}<0.1}$, ${\textstyle {\mathit {Re}}_{\hbox{cloud}}<0.1}$. 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 filament-like 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 well-known of these mechanisms consists in the progressive expulsion of heavy particles from high-vorticity regions into straining-flow 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 high-strain 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 ${\textstyle {\mathit {St}}}$ grows. For instance, the sweep-stick mechanism  [149] consists in the tendency of particles to move away from regions of high accelerations and stick to the low-acceleration zones. The sweep-stick mechanism becomes important at around ${\textstyle {\mathit {St}}=1}$ and higher (inertial range).

Another effect that becomes important at finite Stokes numbers is the so-called 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 sub-Kolmogorov 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 length-scale, 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 ${\textstyle {\mathit {St}}}$ 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 ${\textstyle {\mathit {St}}\sim 1}$  [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 single-particle 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 three-way 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 single-particle 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

 ${\displaystyle S=K_{T}\alpha _{p}L_{\hbox{cluster}}^{2}}$
(2.39)

where ${\textstyle K_{T}}$ is the volumetric shape factor [285], an order-one constant that depends on the shape of the cluster and is equal to one for spherical clusters and ${\textstyle L_{\hbox{cluster}}}$ is the characteristic size of the cluster. This formula was derived under the hypothesis of quasi-steady Stokesian flow around the cluster. That is, this work is concerned with the low-Reynolds 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

1. Rigorously proving the following relation
 ${\displaystyle V_{\hbox{settling}}\sim {\frac {(\rho -1)}{\nu }}CN_{\hbox{cluster}}^{2/3}{\frac {a^{3}}{\Delta d_{\hbox{cluster}}}}g}$
(2.40)
2. where ${\textstyle C}$ is an unknown constant and where ${\textstyle \Delta d_{\hbox{cluster}}}$ is the average separation between the particles in the cluster. This expression is equivalent to Eq.~2.39 in the regime (${\textstyle V_{0}\ll V_{\hbox{settling}}}$), as it can be immediately shown using ${\textstyle \alpha _{p}\sim a^{3}/\Delta d_{\hbox{cluster}}^{3}}$ and ${\textstyle L_{\hbox{cluster}}^{2}\sim N_{\hbox{cluster}}^{2/3}\Delta d_{\hbox{cluster}}^{2}}$.

3. Establishing the existence of a critical, minimum number of particles (${\textstyle N_{\hbox{crit}}=C\left({\frac {d}{R}}\right)^{3/2}}$) that must be contained in the cluster, so that it behaves as a macro-particle and where Eq.~2.39 applies. The macro-particle argument was precisely the one used in [4].
4. Establishing a lower bound to the time scale in which the particles will not come close enough to invalidate the model, which is the time scale of the time taken by the cluster to settle a distance equivalent to ${\textstyle L_{\hbox{cluster}}}$. This point is fundamental in proving that the macro-particle model is stable for a sufficiently long time to apply the theory effectively.

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 non-uniform spatial distribution of particles, to provide an expression for settling coefficient ${\textstyle S}$ in the much more general setting of isotropic turbulence. The general expression reads:

 ${\displaystyle S=S_{1}+S_{2}-5}$ ${\displaystyle S_{1}=\int _{2}^{\infty }{\left(A_{11}+A_{12}+2B_{11}+2B_{12}-3-{\frac {3}{r_{a}}}\right)\Gamma r_{a}^{2}\,\mathrm {d} r_{a}}}$ ${\displaystyle S_{2}=3\int _{2}^{\infty }{(\Gamma -1)r_{a}\,\mathrm {d} r_{a}}}$
(2.41)

where ${\textstyle r_{a}}$ is the radial distance normalized by the particle radius ${\textstyle a}$; ${\textstyle A_{11}}$, ${\textstyle A_{12}}$, ${\textstyle B_{11}}$ and ${\textstyle B_{12}}$ are the mobility functions of the particle pair [24] and ${\textstyle \Gamma }$ 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

 ${\displaystyle N(r_{a})=\underbrace {\frac {\alpha _{p}}{V_{p}}} _{n_{P}}\underbrace {4\mathrm {\pi } r_{a}^{2}\mathrm {d} r_{a}a^{3}} _{\left|B(r_{a}+\mathrm {d} r_{a})-B(r_{a})\right|}}$
(2.42)

where the number density ${\textstyle n_{p}}$ is the number of particles per unit volume and ${\textstyle B(r_{a})}$ is the ball of (nondimensional) radius ${\textstyle r_{a}}$ 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 particle-pair interactions (${\textstyle S_{1}}$), the effect of preferential concentration (${\textstyle S_{2}}$) 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 time-averaged sense both notions are one-to-one related) of the presence of neighbours in remarkably general conditions, and in a time-averaged 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 order-of-magnitude 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 ${\textstyle \Gamma }$ to calculate the settling coefficient. ZA distinguish between three different regimes, valid for heavy particles (${\textstyle \rho \gg 1}$):

Vanishing inertia (${\displaystyle {\mathit {St}}<0.01}$) with no turbulence In this case there is no clustering and thus ${\textstyle \Gamma \equiv 1}$, 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 (${\displaystyle {\mathit {St}}<0.01}$) In general, in a turbulent flow field the particles interactions produce a drift velocity of the non-inertial particles towards one another [382,73]. Taking this drift into account modifies the RDF (see [381]), leading to ${\textstyle S=-4}$. 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 4x101 less than that corresponding to a laminar settling situation. In terms of order-of-magnitude arguments this is in any case a negligible correction.

Finite but small inertia (${\displaystyle 0.01<{\mathit {St}}<0.1}$) In this case, the RDF is derived from the equations presented in [381] for low inertia particles, based on a power-law model for the RDF:

 ${\displaystyle \Gamma (r_{L})={\begin{cases}&c_{0}r_{L}^{-c_{1}}\quad {\hbox{ if }}r_{L}
(2.43)

where ${\textstyle L_{\hbox{cluster}}}$ is defined so as to guarantee continuity and ${\textstyle c_{0}>0}$ and ${\textstyle c_{1}>0}$ depend only on the Stokes number. ${\textstyle L_{\hbox{cluster}}}$ 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 micro-length scale ${\textstyle L}$, 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 two-way 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 ${\textstyle r_{L}\in [0.75,2.75]}$, we will consider it reasonable for order-of-magnitude arguments.

Moreover, ZA point out that for heavy particles, with ${\textstyle \rho \geq 1x10^{3}}$ (although, judging from their [381], this argument could be extended to ${\textstyle \rho \geq 1x10^{2}}$), one has ${\textstyle S_{2}\gg S_{1}}$ and so using Eq.~2.43 we can write

 ${\displaystyle S\approx S_{2}={\frac {3}{(a/L)^{2}}}{\frac {c_{1}c_{0}^{1/c_{1}}}{2-c_{1}}}}$
(2.44)

which can be rewritten as

 ${\displaystyle S\approx b\left({\frac {L_{\hbox{cluster}}}{a}}\right)^{2}}$
(2.45)

where (${\textstyle c_{1}}$ corresponds to ${\textstyle -\chi }$ in [381])

 ${\displaystyle L_{\hbox{cluster}}=c_{0}^{1/c_{1}};\quad b={\frac {3c_{1}}{4-2c_{1}}}}$
(2.46)

The values of ${\textstyle c_{0}}$ and ${\textstyle c_{1}}$ are reasonably approximated in ${\textstyle {\mathit {St}}\in [0.01,0.1]}$ by [381]

 ${\displaystyle c_{0}=1+12St^{2}}$ ${\displaystyle c_{1}=6St^{2}-9.8St^{3}+6.2St^{4}}$
(2.47)

Large inertia (${\displaystyle {\mathit {St}}>0.1}$) In this case the ${\textstyle S_{2}}$ outweighs ${\textstyle S_{1}}$ 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

 ${\displaystyle S\approx S_{2}={\frac {3}{(a/L)^{2}}}\int _{2a/L}^{\infty }{\left(\Gamma -1\right)r_{L}\,\mathrm {d} r_{L}}}$
(2.48)

where now the radial coordinate (${\textstyle r_{L}}$) is normalized by the Kolmogorov micro-scale. ZA provide a set of stochastic, probability conservation equations that can be numerically solved to determine ${\textstyle \Gamma }$.

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 non-local diffusion effect, as had already been applied by [73]. This mechanism is due to the over-sampling 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 ${\textstyle {\mathit {St}}<1x10^{-1}}$. For this range of values of ${\textstyle {\mathit {St}}}$ the modification is reduced to modifying the value of ${\textstyle c_{1}}$ in the model of Chun et al. by dividing it by the value ${\textstyle B_{\hbox{NL}}\approx 4.2x10^{-1}}$. And since in this regime the values for ${\textstyle c_{1}}$ 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 ${\textstyle c_{1}}$ at small values of ${\textstyle {\mathit {St}}}$ (see Fig. 6 to appreciate the corrected level of matching).

Still, for larger values of ${\textstyle {\mathit {St}}}$ 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 ${\textstyle {\mathit {St}}\approx 0.1}$, even when the non-local diffusion correction is applied due, they argue, to its poor prediction of the relative velocity statistics at ${\textstyle {\mathit {St}}\gtrsim 0.1}$.

We have thus taken a different path at ${\textstyle {\mathit {St}}>0.1}$. We use the fitted expression for ${\textstyle \Gamma }$ given by the former authors in this regime. It is again based on the power-law model of Eq.~2.43, where the ${\textstyle c_{0}}$ and ${\textstyle c_{1}}$ are empirically derived. Evidence supporting the power-law scaling of the RDF has been presented for ${\textstyle {\mathit {St}}\lesssim 1}$ [288,178].

Once again, the coefficients ${\textstyle c_{0}}$ and ${\textstyle c_{1}}$ are functions of the ${\textstyle {\mathit {St}}}$ only. This hypothesis can be seen to be reasonable up to ${\textstyle {\mathit {St}}\sim 1}$ 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 ${\textstyle c_{0}}$ and ${\textstyle c_{1}}$, as a function of ${\textstyle {\mathit {St}}}$ 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

 ${\displaystyle f(a,b,c;x)=\quad {\frac {a}{x}}\mathrm {e} ^{-{\frac {\log {x}-b}{2c^{2}}}}}$
(2.49)

which is proportional to a lognormal distribution probability density function. The optimal parameters ${\textstyle c_{0}}$ and ${\textstyle c_{1}}$ were found to be

 ${\displaystyle c_{0}(St)=1+f(1.0713,0.985322,1.09278;St)}$ ${\displaystyle c_{1}(St)=f(11.1892,1.19499,0.93432;St)}$
(2.50)
 (a) (b) Figure 6: Fits for the ${\displaystyle c_{0}}$ and ${\displaystyle c_{1}}$ as a function of ${\displaystyle {\mathit {St}}}$. 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 ${\displaystyle \sigma }$ 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 ${\textstyle c_{0}}$ and ${\textstyle c_{1}}$ we are able to produce an explicit analytic formula for the settling coefficient ${\textstyle S}$ that can be used to establish the first-effects of the influence of the neighbouring particles in a situation were non-uniform 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, ${\textstyle \alpha _{p}\ll S^{-1}}$ or, assuming a maximum value for the acceptable influence at 1 change in the settling velocity, the maximum solid fraction can be expressed as

 ${\displaystyle \alpha _{p}<{\frac {0.01}{S({\mathit {St}})}}}$
(2.51)

where ${\textstyle S}$ has the expression of Eq.~2.45 and where the expression for ${\textstyle c_{0}}$ and ${\textstyle c_{1}}$ can be taken from Eq.~2.47 (dividing ${\textstyle c_{1}}$ by 4.2x10-1) for small values of ${\textstyle {\mathit {St}}}$ and using Eqs.~2.49 and 2.50 for larger values. The behaviour of ${\textstyle S}$ as a function of ${\textstyle {\mathit {St}}}$ 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 ${\displaystyle {\mathit {St}}<0.1}$ (with the non-local 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 ${\textstyle {\mathit {St}}=1}$ toward a positive one. This is unexpected, as it has been empirically observed that ${\textstyle S}$ 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 ${\textstyle {\mathit {St}}=1}$, as it is at around this point that the first signs of Reynolds number dependence in the ${\textstyle c_{0}}$ and ${\textstyle c_{1}}$ 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 power-law model for the RDF is supposed to be valid all the way down, below sub-Kolmogorov scales, in accordance with DNS results for monodisperse suspensions. However, such model might suffer from the effects of a certain level of over-idealization. Chun et al.  [73] studied this issue theoretically for small-inertia particles and concluded that arbitrarily small differences between the diameters of the particle pairs lead to the appearance of a cut-off 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 real-world 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 CFD-DEM 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 inter-particle interactions:

• Removing interactions eliminates the smallest scales in the system (as long as the particle-wall interactions are also eliminated), which leads to the possibility of using much larger time steps
• It can be difficult to accurately model inter-particle interactions, due to uncertainties in the rheology involved
• The number of particles is sometimes made smaller in the simulation compared to the real system to reduce the overall cost, which alters the collision frequency artificially, rendering the effort pointless
• The overall computations become simpler and, thus, cheaper

Therefore, it is interesting to study under what conditions inter-particle 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 history-dependent 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 inter-particle collisions have a negligible effect to the overall trajectories of the particles if the average inter-collision time is substantially larger than the particle relaxation time. This criterion naturally gives rise to the definition of a collisional Stokes number as

 ${\displaystyle {\mathit {St}}_{\hbox{coll}}:={\frac {\tau _{p}}{\tau _{\hbox{coll}}}}\ll 1}$
(2.52)

where ${\textstyle \tau _{\hbox{coll}}}$ is the reciprocal of the average collision rate, ${\textstyle f_{\hbox{coll}}}$. 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 ${\textstyle i,j}$, with their expected relative velocity:

 ${\displaystyle f_{\hbox{coll}}=A_{\hbox{coll}}{\boldsymbol {v}}_{ij}n_{j}}$
(2.53)

where

 ${\displaystyle A_{\hbox{coll}}=\mathrm {\pi } (a_{i}+a_{j})^{2}}$
(2.54)

where ${\textstyle A_{\hbox{coll}}}$ is the swept area formed by the projected silhouette of the two particles onto the plane orthogonal to ${\textstyle {\boldsymbol {v}}_{ij}}$, ${\textstyle a_{i}}$ and ${\textstyle a_{j}}$ are the respective radii and ${\textstyle n_{j}}$ is the number density of the species of the particle ${\textstyle j}$. 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 ${\textstyle A_{\hbox{coll}}}$ 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 gravity-driven 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 ${\textstyle \delta N}$ in a given time interval ${\textstyle \delta t}$ can be calculated as ${\textstyle \delta N=4n\mathrm {\pi } a^{2}v\delta t}$, corresponding to a cylindrical volume of base ${\textstyle 4\mathrm {\pi } a^{2}}$ and height ${\textstyle vt}$, where ${\textstyle v}$ 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

 ${\displaystyle {\mathit {St}}_{\hbox{coll,sett}}={\frac {2}{27}}{\frac {g}{\nu }}\alpha _{p}\left(2\rho +1\right)\left(\rho -1\right)\varepsilon a^{3}}$
(2.55)

and the equation is valid to first order in ${\textstyle \varepsilon }$, where ${\textstyle \varepsilon }$ is the range in diameters normalized by the average diameter.

The turbulence-driven 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 ${\textstyle {\mathit {St}}}$ Staffman and Turner [301] derived the following expression, which has been validated in numerical simulations in the limit of very small Stokes numbers [331].

 ${\displaystyle {\mathit {St}}_{\hbox{coll}}\approx 1.24\alpha _{p}{\mathit {St}}}$
(2.56)

where ${\textstyle {\mathit {St}}}$ represents the usual Stokes number based on the Kolmogorov microscales.

At the opposite extreme, in the so-called uncorrelated regime, the following expression can be derived [220]

 ${\displaystyle {\mathit {St}}_{\hbox{coll, turb}}={\frac {13}{2}}\alpha _{p}{\frac {{\mathit {St}}_{\Lambda {\mathcal {L}}}}{\sqrt {1+{\mathit {St}}_{\Lambda {\mathcal {L}}}}}}{\frac {\Lambda _{\perp }}{a}}}$
(2.57)

where ${\textstyle {\mathit {St}}_{\Lambda {\mathcal {L}}}}$ is the Stokes number based on the eddy turnover time and where ${\textstyle \Lambda _{\perp }}$ 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 particle-laden simulations. The result is

 ${\displaystyle {\mathit {St}}_{\hbox{coll, turb}}={\frac {\alpha _{p}}{1+{\mathit {St}}_{\Lambda {\mathcal {L}}}}}\left(1.2{\mathit {St}}+{\frac {13}{2}}{\frac {{\mathit {St}}_{\Lambda {\mathcal {L}}}^{2}}{\sqrt {1+{\mathit {St}}_{\Lambda {\mathcal {L}}}}}}{\frac {\Lambda _{\perp }}{a}}\right)}$
(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

 ${\displaystyle {\mathit {St}}_{\hbox{coll}}=\min({\mathit {St}}_{\hbox{coll,sett}},{\mathit {St}}_{\hbox{coll, turb}})\ll 1}$
(2.59)

(1) The asterisk has been added to avoid a symbolic clash with the settling coefficient, ${\displaystyle S}$.

(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 ${\displaystyle 2a}$ because no particle centres can be found closer than ${\displaystyle 2a}$.

### 2.2.5 Small-size effects

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.7-2.9 can all be fulfilled by making ${\textstyle a}$ 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 quasi-equilibrium) that lead to the linear stress-strain rate relation [135]. The progressive deterioration of the validity of this assumption as the length scales shrink manifests itself as the so-called 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, ${\textstyle {\mathit {Kn}}}$, specifically quantifies this phenomenon and is defined as

 ${\displaystyle {\mathit {Kn}}={\frac {\lambda }{a}}}$
(2.60)

with ${\textstyle \lambda }$ the mean free path (MFP) and ${\textstyle a}$ 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 no-slip 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 molecular-scale 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 out-of-equilibrium region is known as the Knudsen layer [193]. We can conclude that the assumption of no-slip boundary condition at the solid wall must introduce an error of the order of ${\textstyle \lambda }$, that is, proportional to ${\textstyle {\mathit {Kn}}}$.

While such an error is the manifestation of the onset of non-continuum 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 ${\textstyle {\mathit {Kn}}}$ is small enough. The set of low ${\textstyle {\mathit {Kn}}}$ 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 no-slip 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 low-Knudsen number limit (weak non-continuum effects), neglecting for now the other forces.

As previously mentioned, the Navier-Stokes 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:

 ${\displaystyle {\frac {a}{\mu \beta _{\hbox{slip}}}}\tau }$
(2.61)

where ${\textstyle \tau }$ is the local tangential stress and ${\textstyle \beta _{\hbox{slip}}}$ is known as the Basset slip coefficient, in general ranging from zero (perfect slip) to infinity (no-slip conditions).

The corresponding drag force on a sphere moving in uniform motion through an otherwise quiescent infinite fluid is (see, e.g. [283])

 ${\displaystyle {\frac {\left|{\boldsymbol {F_{D,{\hbox{1st order}}}}}\right|}{\left|{\boldsymbol {F_{D,{\hbox{Stokes}}}}}\right|}}={\frac {1+2\beta _{\hbox{slip}}^{-1}}{1+3\beta _{\hbox{slip}}^{-1}}}\Rightarrow {\frac {\left|{\boldsymbol {F_{D,{\hbox{Stokes}}}-F_{D,{\hbox{1st order}}}}}\right|}{\left|{\boldsymbol {F_{D,{\hbox{Stokes}}}}}\right|}}={\frac {2\beta _{\hbox{slip}}^{-1}}{1+3\beta _{\hbox{slip}}^{-1}}}}$
(2.62)

since ${\textstyle \beta _{\hbox{slip}}}$ 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 ${\textstyle \beta _{\hbox{slip}}}$ with the mean free path from the kinetic theory of monatomic gases in isothermal conditions.

 ${\displaystyle \beta _{\hbox{slip}}^{-1}={\frac {2-\sigma _{m}}{\sigma _{m}}}{\mathit {Kn}}}$
(2.63)

where ${\textstyle \sigma _{m}}$ 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 (no-slip). 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 ${\textstyle \sigma _{m}}$ depends on the gas and the surface material and finish, but experimental evidence shows ${\textstyle \sigma _{m}}$ to be between 2x10-1 and 8x10-1 (or ${\textstyle \beta _{\hbox{slip}}^{-1}\in (0.11{\mathit {Kn}},0.67{\mathit {Kn}})}$), with the former corresponding to specially treated, ultra-smooth surfaces and the latter to most practical surfaces [133]. Maxwell model is derived, by taking the first order (in ${\textstyle {\mathit {Kn}}}$) effects into account, and thus its validity is bounded below by ${\textstyle {\mathit {Kn}}<0.1}$; the so-called 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 ${\textstyle {\mathit {Kn}}\ll 1}$, which imposing our accustomed one percent error bound reads

 ${\displaystyle {\mathit {Kn}}<1x10^{-2}}$
(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 ${\textstyle {\mathit {Kn}}\lesssim 1x10^{-3}}$ for ${\textstyle \sigma _{m}=0.2}$, putting the minimum radius at a ten times larger value.

While small-frequency conditions are the most relevant for a majority of applications with low ${\textstyle {\mathit {Kn}}}$, 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 ${\textstyle {\mathit {Kn}}}$, 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 ${\textstyle {\mathit {Kn}}}$. We consider this case to be representative for unsteady motion in low ${\textstyle {\mathit {Kn}}}$ conditions, at least for order-of-magnitude arguments. In the referred paper, expressions are provided for the asymptotic limits in which we are interested, that is, the limit of ${\textstyle {\mathit {Kn}}\ll 1,\Theta \ll 1}$, where ${\textstyle \Theta }$ is a nondimensional frequency defined as (see [376])

 ${\displaystyle \Theta ={\frac {\omega }{\nu _{m}}}={\mathit {Kn}}^{2}\beta }$
(2.65)

where ${\textstyle \nu _{m}}$ is the molecular collision frequency, ${\textstyle \omega }$ is the characteristic frequency of the relative motion and ${\textstyle \beta }$ is a Stokes number 3, defined as

 ${\displaystyle \beta ={\frac {a^{2}\omega }{2\nu }}}$
(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 no-slip limit (which can be found, for example, as a solved exercise in [202]). The ratio of the amplitudes of both forces (with the no-slip case in the denominator), which we refer to as ${\textstyle \alpha _{\mathit {Kn}}}$, provides an adequate measure of the distance from the MRE range of applicability, to first order in the ${\textstyle {\mathit {Kn}}}$. 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 ${\textstyle \Theta }$ and ${\textstyle {\mathit {Kn}}}$ grow, as the momentum transfer at the particle's surface weakens.
 Figure 8: Ratio of force calculated as in [376], valid for small ${\displaystyle {\mathit {Kn}}}$ and ${\displaystyle \beta }$, to the magnitude of the force calculated with the MRE, valid under the continuum, no-slip conditions on the surface of the sphere.

Note that the intersection of the graph with the plane ${\textstyle \beta =0}$ recovers the steady-state asymptotic behaviour discussed above. For the case shown, ${\textstyle \beta _{\hbox{slip}}^{-1}=1.016191}$, consistent with [376]. Such value comes as a result of refined gas-kinetic simulations, which correct expression Eq.~2.63 and taking ${\textstyle \sigma _{m}=1}$, 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 ${\textstyle {\mathit {Kn}}}$ and ${\textstyle \beta }$ so that ${\textstyle \Theta }$ is bounded above by 1x10-1 since, as mentioned, the validity of this formulation is confined to small values of ${\textstyle {\mathit {Kn}}}$ and ${\textstyle \Theta ={\mathit {Kn}}^{2}\beta }$. Nonetheless, within these restrictions, we have chosen the range of values for ${\textstyle \beta }$ 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 ${\textstyle \beta }$ can be related to Eq.~2.4 by using the estimate ${\textstyle \omega \sim 2\mathrm {\pi } /T}$, with ${\textstyle T}$ the characteristic time of the relative motion (which in turbulence-driven flows will be governed by the fluid's small scales motion). This gives

 ${\displaystyle {\mathit {St}}\sim {\frac {2\rho +1}{9\mathrm {\pi } }}\beta }$
(2.67)

The order of magnitude of the coefficient of ${\textstyle \beta }$ can reach values of over 1x103 for some materials (e.g.heavy metals suspended in air). Thus, roughly speaking, the graph covers ${\textstyle {\mathit {St}}\in (0,10^{4})}$, well beyond our range of interest (see Table 1).

In summary, at the small-size 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 ${\textstyle {\mathit {Kn}}}$ 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 ${\textstyle {\mathit {Kn}}}$, 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 no-slip boundary condition is not available yet [134]. Another consequence of the higher density of gases is the change in the order of appearance of small-size effects (see [249]; compare Figures 1 and 4). Indeed, here non-continuum effects come before thermodynamic non-equilibrium effects. That is, statistical fluctuations become significant before the no-slip 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 no-slip boundary condition was employed by Einstein in his famous analysis of Brownian motion. Thus, the first small-size 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

 ${\displaystyle {2}\tau _{c_{s}}={\frac {a}{c_{s}}}\quad {\hbox{ (sonic time)}}}$ ${\displaystyle \tau _{p}={\frac {m_{p,f}}{\xi }}\quad {\hbox{ (particle's response time)}}}$ ${\displaystyle \tau _{\nu }={\frac {a^{2}}{\nu }}\quad {\hbox{ (viscous diffusion time)}}}$ ${\displaystyle \tau _{D}={\frac {a^{2}}{D}}\quad {\hbox{ (molecular diffusion time)}}}$
(2.68)

where ${\textstyle m_{p,f}=m_{p}+{\frac {1}{2}}m_{f}}$ is the equivalent mass (particle mass plus added mass), ${\textstyle c_{s}}$ is the local sound speed in the fluid, ${\textstyle \xi }$ is the inverse of the mobility (that is, the constant quantity in the expression of ${\textstyle {\boldsymbol {F}}_{D}}$ drag coefficient 5); and ${\textstyle D}$ is the Stokes–Einstein–Smoluchowski diffusivity, given by the formula

 ${\displaystyle D={\frac {k\theta }{\xi }}}$
(2.69)

where ${\textstyle k}$ is Boltzmann's constant and ${\textstyle \theta }$ is the absolute temperature.

The meaning of the different time scales is as follows: ${\textstyle \tau _{c_{s}}}$ is the time it takes a sound wave to travel through the fluid a distance equivalent to the particle's radius, ${\textstyle a}$; ${\textstyle \tau _{p}}$ is the time it takes for a particle to lose a ${\textstyle 1-1/\mathrm {e} }$ fraction of its initial slip velocity, due to the action of the steady drag force and added mass forces alone; ${\textstyle \tau _{\nu }}$ is the time it takes vorticity to diffuse over a length ${\textstyle a}$ from the particle's surface, while the time scale ${\textstyle \tau _{D}}$ is the long-term average time it takes for a Brownian particle to diffuse, again, to a distance ${\textstyle a}$ (see the Langevin equation Eq.~2.79). Note that the definitions of ${\textstyle \tau _{p}}$ (see beginning of this Section 2.2) and of ${\textstyle \tau _{\nu }}$ (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, ${\textstyle \tau _{cs}}$ is much smaller than the rest, ${\textstyle \tau _{\nu }}$ is typically smaller than ${\textstyle \tau _{p}}$ for gases, while both are comparable for liquids and ${\textstyle \tau _{D}}$ is much larger than the rest for all but the tiniest length-scales, where the continuum hydrodynamic theory has started to break down anyway. Manifestly, such ordering is robust, at least for particles larger than 1x102 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 time-scale, 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 ${\textstyle {\mathit {Re}}\sim {\frac {L^{2}/\nu }{L/U}}}$. At the Kolmogorov scales the Reynolds number is by definition of order one, so

 ${\displaystyle Re\sim {\frac {\eta ^{2}/\nu }{\tau _{\eta }}}\sim 1\Rightarrow a^{2}/\nu =\tau _{\nu }\ll \tau _{\eta }}$
(2.70)

since we have by assumption that ${\textstyle a\ll \eta }$. In fact, by the same argument, we have the even stronger

 ${\displaystyle {\sqrt {\tau _{\nu }}}\ll {\sqrt {\tau _{\eta }}}}$
(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 ${\textstyle \tau _{D}}$, due to its pseudo-fractal 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 ${\textstyle \tau _{p}}$ [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:

 ${\displaystyle \delta (t):={\sqrt {\hbox{MSD(t)}}}}$
(2.72)

where

 ${\displaystyle {\hbox{MSD}}(t)=\left\langle \left\|{\boldsymbol {r}}(t)-{\boldsymbol {r}}(0)\right\|^{2}\right\rangle }$
(2.73)

That is, the ensemble average of the Euclidean norm of the displacement after a given time lapse ${\textstyle t}$. In this way, ${\textstyle \delta (t)}$ can be interpreted as the amplitude of the Brownian motion at the time scale ${\textstyle t}$. So that, instead of comparing forces, we can compare the characteristic displacements ${\textstyle \delta (\tau _{\eta })}$ and ${\textstyle \delta _{f}(\tau _{\eta })}$ 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 inertia7. 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 ${\textstyle \tau _{\eta }}$, 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:

 ${\displaystyle W\sim (1-\beta _{\rho }){\frac {\tau _{p}}{\tau _{\eta }}}U=(1-\beta _{\rho }){\mathit {St}}U}$
(2.74)

where

 ${\displaystyle \beta _{\rho }={\frac {3}{2\rho +1}}}$
(2.75)

and ${\textstyle \tau _{\eta }}$ is the Kolmogorov time scale, while ${\textstyle U}$ 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 ${\textstyle {\mathit {St}}<1}$, according to Balachandar, as we will discuss in the following section8).

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 (${\textstyle \equiv \epsilon ^{3/4}/\nu ^{1/4}}$, with ${\textstyle \epsilon }$ the turbulent energy dissipation rate) is smaller than ${\textstyle \left\|{\boldsymbol {g}}\right\|}$. In order to derive our order-of-magnitude relation, it is sufficient to separately consider the steady settling velocity, ${\textstyle W_{g}}$, which is obtained by balancing the drag force with the submerged weight, giving

 ${\displaystyle W_{g}=(1-\beta _{\rho }){\mathit {St}}\tau _{\eta }g}$
(2.76)

Once the order of magnitude of the slip velocity has been identified, we can construct the following time-dependent distance scale

 ${\displaystyle \delta _{f}(t)=Wt}$
(2.77)

which is a measure of the inertial deviation from the perfect tracer particle trajectory of a given particle. Similarly, we define ${\textstyle \delta _{g}(t)=W_{g}t}$.

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 ${\textstyle \delta (\tau _{\eta })/\delta _{f}(\tau _{\eta })}$ (or ${\textstyle \delta (\tau _{\eta })/\delta _{g}(\tau _{\eta })}$); 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

 ${\displaystyle {\frac {\delta (\tau _{\eta })}{\delta _{f}(\tau _{\eta })}}\,,{\frac {\delta (\tau _{\eta })}{\delta _{g}(\tau _{\eta })}}\ll 1}$
(2.78)

The Langevin theory [345] provides the right tool to calculate ${\textstyle \delta (t)}$. 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 often-called 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 ${\textstyle \tau ^{*}=t/\tau _{\nu }}$:

 ${\displaystyle {\hbox{MSD}}(\tau ^{*})=2D\tau _{\nu }\tau ^{*}\left(1-{\sqrt {\frac {4}{\mathrm {\pi } }}}\tau ^{*-{\frac {1}{2}}}+\left(1-{\frac {\tau _{p}}{\tau _{\nu }}}+\varepsilon (\tau ^{*})\right)\tau ^{*-1}\right)}$
(2.79)

where

 ${\displaystyle \varepsilon (\tau ^{*})={\frac {1}{\sqrt {1-4{\frac {\tau _{p}}{\tau _{\nu }}}}}}\left({\frac {1}{\alpha _{+}^{3}}}\mathrm {e^{\alpha _{+}^{2}\tau ^{*}}} {\hbox{erfc}}(\alpha _{+}{\sqrt {\tau ^{*}}})-{\frac {1}{\alpha _{-}^{3}}}\mathrm {e^{\alpha _{-}^{2}\tau ^{*}}} {\hbox{erfc}}(\alpha _{-}{\sqrt {\tau ^{*}}})\right)}$
(2.80)

with

 ${\displaystyle \alpha _{\pm }={\frac {1\pm {\sqrt {1-4{\frac {\tau _{p}}{\tau _{\nu }}}}}}{2{\frac {\tau _{p}}{\tau _{\nu }}}}}}$
(2.81)

The parameter ${\textstyle \varepsilon }$ is a small correction, relevant only for ${\textstyle \tau ^{*}\ll 1}$ [277]. The function defined in Eq.~2.79 has two extreme regimes; those corresponding to ${\textstyle \tau ^{*}\ll 1}$ (ballistic regime), and ${\textstyle \tau ^{*}\gg 1}$ (diffusive regime). It is a direct calculation to check that their corresponding asymptotic relations are given by

 ${\displaystyle {2}{\hbox{MSD}}(\tau ^{*})\sim {\frac {k\theta }{m_{p,f}}}\tau _{\nu }^{2}\tau ^{*2}\quad {\hbox{ as }}\tau ^{*}\rightarrow 0}$ ${\displaystyle {\hbox{ MSD}}(\tau ^{*})\sim 2D\tau _{\nu }\tau ^{*}\quad {\hbox{ as }}\tau ^{*}\rightarrow \infty }$
(2.82)

Eq.~2.79 can be used to compute ${\textstyle \delta (\tau _{\eta })}$. 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 ${\textstyle \delta }$ for ${\textstyle \tau ^{*}<10}$ and ${\textstyle \tau ^{*}\geq 10}$ as the ballistic and diffusive regimes. By gluing together the two solutions at ${\textstyle \tau ^{*}=10}$ we generate the function ${\textstyle \delta _{b,d}}$