Inertial focusing of a dilute suspension in pipe flow

The dynamics of rigid particle suspensions in a wall-bounded laminar flow present several non-trivial and intriguing features, including particle ordering, lateral transport, and the appearance of stable, preferential locations like the Segr\'e-Silberberg annulus. The formation of more than one annulus is a particularly puzzling phenomenon that is still not fully explained. Here, we present numerical simulation results of a dilute suspension of particles in (periodic) pipe flow based on the lattice Boltzmann and the discrete element methods (DEM). Our simulations provide access to the full radial position history of the particles while traveling downstream. This allows to accurately quantify the transient and steady states. We observe the formation of the secondary, inner annulus and show that its position invariably shifts toward the Segr\'e-Silberberg one if the channel is sufficiently long, proving that it is, in fact, a transient feature for Reynolds numbers (Re) up to 600. We quantify the variation of the channel focusing length ($L_s/2R$) with Re. Interestingly and unlike the theoretical prediction for a point-like particle, we observe that $L_s/2R$ increases with Re for both the single particle and the suspension.


I. INTRODUCTION
The idea behind inertial microfluidics stems from the pioneering work of Segré & Silberberg [1,2] on inertial focusing of neutrally buoyant particles subject to a Poiseuille flow in a straight pipe. Segré & Silberberg observed that the particles migrate across the streamlines and focus on a thin annular region. The authors attributed this phenomenon to a transverse force, which is centrifugal at the channel center and centripetal at the walls. In the absence of inertial and non-Newtonian effects, Saffman [3] pointed out that a rigid spherical particle in a creeping flow does not experience any sideways force in flows with a unidirectional velocity. Goldsmith & Mason [4] conducted similar experiments to the one of Segré & Silberberg, but under a creeping flow condition and observed that the cross-stream migration occurs only for non-rigid spherical particles. A follow-up study of the same authors [5] shows a cross-streamline inward migration of rigid spheres in non-Newtonian fluids. Bretherton showed theoretically that rigid particles with specific extreme shapes can experience a lateral drift across streamlines without inertial and non-Newtonian effects [6]. Several subsequent theoretical works, based on the method of matching asymptotic expansions, investigated systems differing from that of Segré & Silberberg in Reynolds number (Re), ratio of the particle radius to the pipe radius (a/R), bare fluid flow, or channel geometry. [7][8][9][10] Interestingly, Ho & Leal [8] found that the lateral equilibrium position of a small particle (a/R 1) in a planar two-dimensional Poiseuille flow is identical to the value measured by Segré & Silberberg along the radial direction of the cylindrical pipe. Furthermore, Asmolov [10] pointed out that the equilibrium position in a shear flow shifts further toward the walls as Re increases. In a follow-up study, Asmolov and co-workers [11] provided an expression for the lift force on a point-like particle in a confined shear flow, which is asymptotically valid independently from the lateral particle position provided that a/R 1 and Re 10.
In addition, they found a qualitative agreement with numerical simulations of a finite-size particle using the lattice Boltzmann method.
Matas, Morris & Guazzelli [12] studied the inertial migration of a dilute suspension of particles in a cylindrical channel for a wide range of Re and a/R experimentally. They observed that in addition to the outer annulus reported by Segré-Silberberg, particles accumulate into an inner annulus, which the existing theory did not predict.  [15] reported a turning point bifurcation with an unstable branch between both annuli (stable branches) that spans over a few tens of Re for a large particle (a/R = 0.15) and over a few hundreds of Re for a small particle (a/R = 0.084).
Here, we focus on the role played by the channel length by measuring the channel focusing length (also known as entry length) required to observe radial focusing of a dilute suspension of particles. We simulate a dilute suspension of finite size particles with a volume fraction φ up to 1%, a/R = 0.1 and Re up to 600. This choice of parameters is motivated by our interest in studying the existence of the unstable branch discussed in the experiment of Matas [12].

II. NUMERICAL SCHEME
This section provides a summary of the numerical method and the simulation setup.
More details on the simulation method are available in the appendix.
We consider a suspension of rigid spherical particles of radius (a) subject to a parabolic flow in a pipe with a radius (R) and a length (L). The system is periodic along the z−axis, the cylindrical channel is mainly governed by the channel Reynolds number (Re = u 0 R ν ), which quantifies the importance of the inertial forces over the viscous forces (ν is the kinematic viscosity of the fluid), and the particle-to-channel radius (a/R). We set the volume fraction of the suspension to φ ≤ 1% to minimize the effect of particle-particle hydrodynamic interactions. Each simulation was carried out on 200 CPUs on the high-performance computing clusters available at the Erlangen National High Performance Computing Center (NHR@FAU) or the Jülich Supercomputing Center (JSC).
The motion of the fluid is described by the discrete lattice Boltzmann equation in velocity space that reads as where f i is the discrete probability function of finding a particle at position x and time t moving with velocity c i , i = 1, ..., 19 on a 3D Eulerian grid [16][17][18][19]. Here, is the standard Bhatnagar-Gross-Krook (BGK) collision operator [20], τ is the relaxation time, and f eq i is a truncated expansion of the Maxwell-Boltzmann distribution for the particle velocities in an ideal gas and corresponds to the local equilibrium distribution function. The lattice constant and the discrete time step are denoted by ∆x and ∆t, respectively. The LBM follows a stream and collide scheme and thus can be divided into two steps: i) the collision step wheref is the post-collision distribution function, and ii) the advection step where the discrete particle probability distributions are streamed from one lattice node to the next one according to their velocities To generate a parabolic flow, we apply a body force on the fluid along the z−axis, incorporated into Eq. 1 through the source term F i .
where f (x, t) accounts for the body force, u is the macroscopic velocity, c s = 1/3∆x/∆t is the lattice speed of sound, and ω i are the lattice weights which, for the D3Q19 LBM employed here, read as 1/3, 1/18 and 1/36 for i = 1, i = 2 . . . 7, and i = 8 . . . 19, respectively. The fluid density (ρ) and the macroscopic velocity on each lattice node are calculated from the zeroth and first moments of the discrete probability function f i such that The kinematic viscosity of the fluid is related to the relaxation time τ by with ρ 0 being the mass density. By applying a body force (f z ) mimicking a pressure gradient along the z direction and a no-slip boundary condition, we obtain a parabolic profile with a maximum velocity u 0 = fzR 2 4ρν . We restrict our investigation to a monodisperse suspension of rigid spherical particles with a density ratio between the particle (ρ p ) and the fluid fixed to unity. The particle's motion is solved numerically using the discrete element method (DEM), a widely employed numerical scheme in engineering and physics [21]. The DEM is based on integrating Newton's equations of motion for a rigid body where F tot and T tot are the net force and torque exerted on the solid particle, I and M = ρ p (4/3)πR 3 are the inertia tensor and the mass of the particle, respectively. The contributions to the total force (and torque) stem from the hydrodynamic interactions between the fluid and the particle (F H ), and the lubrication force (F lub ). Equations 7 and 8 are integrated using a leapfrog scheme [22].

III. RESULTS
In this section, we present and discuss the role of inertia on the spatial distribution of a dilute suspension of particles for Re ≤ 600 and for a particle-to-pipe radius of a/R = 0.1. Three volume fractions are considered: φ = 0.02% (single particle), φ = 0.5% and φ = 1%. In Fig. 2   suspension of particles (red circles) exhibits a similar behavior but with a slightly higher asymptotic limit of r s /R ≈ 0.83. We present the effect of the volume fraction in Fig. 3(b).
Although we are in the dilute regime (φ ≤ 1%), we observe that r s /R is typically larger for suspensions (φ = 0.5% and φ = 1%) when compared to a single particle (φ = 0.02%) which suggests that we cannot neglect the hydrodynamic interaction between particles in the inertial regime. We also observe a similar behavior by comparing to the experiments of MMG [12] and NYYIS [23] at, respectively, φ ≈ 0.1% and φ ≈ 0.001%.
In order to exclude the effect of hydrodynamic interactions between particles, we now consider the case of an isolated particle. Looking at the radial position history in Fig. 4, we find that the particle always migrates toward the same steady-state radial position independently of its initial position. This behavior is persistent for the parameter space considered in this study and defined by Re ≤ 600 and a/R = 0.1. We now turn to the channel focusing length,  an important observable that quantifies the channel length required for the particle to reach its radial steady-state position. Matas et al. [12] estimated the channel focusing length to be inversely proportional to the channel Reynolds number in straight rectangular geometries.
In a subsequent study, Matas et al. [24] observed a similar behavior for pipe geometries, although the decrease of the focusing length with the increase of Re was found to be less pronounced when compared to rectangular geometries. The estimate of the focusing length was based under the assumption that a/R 1, and therefore Re p = Re(a/R) 2 1. Using the finite element method, Di Carlo [25] found for a finite-size particle that the scaling of the lift force depends on the position of the particle in the channel and on the particle size. To further investigate this phenomenon, we start by considering the case of an isolated particle initially located at the same radial position and measure its radial trajectory while varying Re. The results are depicted in Fig. 5 for Re = 125, 400, and 600. Interestingly, and in contrast to the analytical predictions, we observe an increase of the channel focusing length with Re. Indeed using matched asymptotic expansions, Matas [12] predicted the focusing length to decrease with Re, such that L s /2R = 6πA −1 Re −1 (R/a) 3 , where A is the magnitude of tal studies [12,14,23]. For that, we have considered three set of simulations with different initial positions. The number of particles for each simulation is 24. We have computed the probability distribution for each set of data separately and then for the combined three set of data to improve the statistics. In both cases, we have observed the same behavior that is reported in Fig. 6 for Re = 38 and Re = 600. It can be seen that the particles reach the radial steady state already after a focusing length of approximately L/2R = 300 at Re = 38.
At Re=600 on the other hand we here still observe the remnants of a bimodal distribution, corresponding to an inner and outer annulus, as initially reported in the experimental work of Matas [12] in pipes. As the particles travel further downstream, the inner peak disappears and the particles accumulate at the outer (Segré-Silberberg) annulus, which suggests that the channel focusing length for Re = 600 and a/R = 0.1 is beyond L/2R = 300 and that any measure of the radial steady-state should be done further downstream. A second inter- esting observation is that the channel focusing length (L s /2R) for Re = 600 is larger than for Re = 38, suggesting that L s /2R increases with Re, similar to the single particle case.
This unexpected results was also observed in the experimental work of Nakayama et al. [23] for Re up to 600 but not discussed further. To explain this behaviour, we report in Fig. 7 the probability of finding the particles at the outer (P out ) annulus at different downstream axial positions. The location of the outer annulus is known from the radial steady-state measurement performed in Fig.3. We observe that independently from Re, the probability of the particles accumulated around the channel center decreases as L/2R increases until the focusing length is reached, where the particles are fully migrated in the outer annulus.
Our simulations clarify that the inner annulus reported by Matas et al. [12], and in the parameter space investigated here, is a transient configuration which is in agreement with the recent experimental work of Nakayama et al. [23] using longer pipes with L/2R up to 1000.
Interestingly, the existence of the inner annulus seems to persist at even further downstream axial positions with the increase of Re, which may explain the increase of the focusing length with Re in pipe geometries. This can be seen clearly in Fig. 8, where we show snapshots of the radial particle distribution over the pipe cross-section at different axial positions and for different Re. We can observe the appearance of the inner annulus at relatively large Re at L/2R ≤ 300 and its full disappearance at a further downstream location, here chosen at L/2R = 600.

IV. CONCLUSION
We have studied numerically the inertial migration of an isolated particle and a dilute suspension of particles (φ ≤ 1%) in a pipe flow for Re up to 600 and a/R = 0.1. Our simulations shed light on the particle dynamics by continuously monitoring them on their downstream propagation. We found that in the dilute regime, the particles migrate radially and form a structured ring located between the center and channel walls corresponding to the tubular pinch effect reported by Segré and Silberberg [1]. The radius of the ring increases with Re and reaches an asymptotic limit at around Re ≥ 600. Interestingly for the suspension case, the steady-state radial position r s is further away toward the wall than in the single-particle case. At the same time, we do not observe a significant difference between volume fractions of φ = 0.5% and φ = 1%. This suggests that the pair hydrodynamic  [12] at L/R = 313. Interestingly, we found that the channel focusing length increases with Re, unlike the analytical predictions for a point-like particle in rectangular [26] and pipe [24] geometries. The experiments on longer pipes [23] support our result, although the authors do not explicitly discuss it. Thus, we speculate that the disagreement with the analytical model stems from the particle being sufficiently large to disturb the flow on the one hand and the pipe geometry on the other hand.
particle such that implies that a discrete distribution function moving from a fluid node (x) toward a solid node (x s ) will bounce back, half-way from x s , on the boundary node (x b ) and travel to the opposite direction with a velocity cī = −c i . The last term in the RHS of the equation A1 accounts for a momentum correction applied when a boundary node x b is moving with a velocity u b defined as where U and Ω are the solid particle translational and angular velocities respectively, and X is the solid particle center of mass position. In the case of static solid nodes, e.g. the pipe boundaries, the BB rule is reduced to The no-slip boundary condition is satisfied on the solid boundaries (particles and walls) through the BB rule with second-order accuracy, provided that the grid resolution is high [19].
Fluid nodes without a single lattice velocity directed toward the solid nodes are typically streamed following Eq. 3.

Particle dynamics
The hydrodynamic force applied to the solid particle is obtained using the momentum exchange method (MEM). The hydrodynamic force can be calculated as the sum of the BB collisions over the boundary nodes located mid-way between the fluid and solid nodes. As a result, the hydrodynamic force and torque exerted on the solid particle can be written as When the solid particle is advected, the solid nodes on the back of the particle (opposite to the motion's direction) are cleared, and their momentum is transferred to the newly created fluid nodes. A similar procedure is applied to the fluid nodes on the front of the particle that is converted into solid nodes, but this time the momentum is transferred from the fluid to the newly created solid nodes. The covering and uncovering of fluid nodes can lead to a violation of global momentum conservation. Thus a correction term needs to be added to the hydrodynamic force and torque exerted on the particle to fulfill the global momentum conservation. We follow here the approach described by Jansen & Harting in [30,31] which the expressions of the hydrodynamic force and torque after including the corrections due to the covering and uncovering of fluid nodes during the particle's motion read as and Here ρ is the fluid density of the newly created fluid node calculated by averaging the fluid density of the first set of neighboring nodes [28]. When the surface-to-surface distance between two solid particles is below one lattice space, the flow in the gap can not be resolved anymore. To overcome this issue, we add a short-range lubrication correction force as discussed in Ref. [27,32]. The lubrication force correction exerted on particle i by the particle j is defined as where a is the radius of the particles, d ij is the mass-to-mass center distance andd ij = d ij /d ij its correspondent unit vector. The lubrication force correction is applied when the gap distance between the particles ||δ ij || = d ij − 2R is smaller than the cut off distance δ c which is fixed in this work to 2/3 as suggested by Nguyen & Ladd [27,32]. complete radial focusing of the particles at Re = 38 occurs at a shorter axial distance L/2R as compared to Re = 600, and that there is no noticeable accumulation of the particles at the inner annulus. Conversely, at Re = 600, the particles form two rings corresponding to the Segré and Silberberg annulus and the inner annulus. The latter disappears completely as the particles travel further downstream along the channel length.