After reviewing and organizing the literature on the problem of inertial cross-stream migration of rigid spheres in various geometries including tubes and channels, we use Dissipative Particle Dynamics (DPD) simulations to study the simplest case of migration of a single neutrally or non-neutrally buoyant sphere with diameter 20% of the gap in plane Poiseuille flow and assess the potential and limitations of DPD simulations for this and similar problems. We find that the neutrally buoyant sphere lags by up to 6% behind the surrounding fluid and is focused at a position around 50% of the distance between the channel center and the wall. With Re increasing from around 100 to 500, the sphere migrates closer to the channel center. With flow driven by gravity, a much denser non-neutrally buoyant sphere leads the surrounding fluid and is focused at a position closer to the wall, around 60% the distance from the channel center to the wall, in qualitative agreement with previous work. The lower values of the Schmidt number Sc in DPD simulations relative to real fluids, due to the relatively large diffusivity of DPD beads, are shown to not significantly affect the consistency of our DPD results with literature results although they make results noisy at low Re (i.e., 50). However, the increase in Ma and Wi with increasing Re leads to compressible flow effects and in some cases viscoelastic effects at high Re depending on the DPD parameters chosen. Even for optimally chosen parameters, we require Re500 to avoid strong compressibility effects. Thus, the relative simplicity of the DPD method for complex fluid flows is offset by the need to control the effects of unphysically high values of other parameters, such as Ma and Wi, which seriously limits the range of conditions under which DPD simulations give valid results in fluid transport problems.

Micron-sized solid particles, droplets, vesicles, and cells, when suspended in Poiseuille flow, can migrate across streamlines and even be focused at well-defined off-center positions within the channel or tube. When induced by inertial effects, this is referred to as inertial migration or inertial focusing.1–5 Inertial migration in the Poiseuille flow in straight channels and tubes results from a combination of a wall-directed shear-gradient-induced lift force due to the parabolic velocity profile and a center-directed wall-induced lift force. (Here “lift force” means simply a force acting perpendicular to the streamlines.) Particle rotation and slip-shear also can induce migration, but have a minor effect on inertial migration of solid particles that are small relative to the gap.6 The deformability of droplets, vesicles, and cells significantly changes their inertial migration relative to solid particles.7–11 Inertial focusing is significant when the particle Reynolds number (the ratio of fluid inertial to viscous force) is on the order of unity or higher, which can be produced even in micro-channels at high throughput.5,12 Thus, microfluidic devices employing inertial focusing have been used to focus, separate, and sort bio-particles, enabling applications such as high throughput cytometry, enrichment of target cells, and detection of infected or diseased cells, such as cancer cells.1,13

As first reported in 1961, Segré and Silberberg14 found that small, rigid spheres, initially uniformly dispersed at the entrance of a straight tube, focused into an annular ring at a position around 60% of the distance between the centerline and wall in a laminar flow. This became known as the “tubular pinch effect,”14–16 attributed to the inertia of the moving fluid.

This surprising phenomenon triggered many subsequent theoretical, numerical, and experimental studies, which are summarized in Tables I–IV and Figs. 1 and 2 for neutrally buoyant rigid particles flowing in tubes and channels of various cross sections: circular, square, slit, and rectangular. In these tables, α = a/Lc is the “blockage ratio,” with a being the diameter of the sphere and Lc being the diameter of a tube, or the gap of a 2D channel, or the smallest gap of a rectangular channel. The Reynolds number is defined as Re = ρfUmLc/μf for a tube or 2D channel or Re = ρfUmDh/μf for a rectangular channel, with ρf being the density of the fluid, Um being the maximum velocity of the flow, μf being the viscosity of the fluid. Dh = 2WH/(W + H) is the hydraulic diameter of a rectangular channel, with W and H being the width (of the long wall) and height (of the short wall), respectively. When W = H, i.e., for square channels, Dh = H. Notably, Re in some publications is based on the mean velocity, which we have transformed to the above definition of Re, which is based on the maximum velocity. zeq/(Lc/2) is the relative equilibrium position within the cross section of the tube or channel, where zeq is the distance from the channel center to the center of the sphere at its equilibrium position. Specifically, for tubes, zeq/(Lc/2) = zeq/R, with R being the radius of the tube; for 2D channels and square channels, zeq/(Lc/2) = zeq/(H/2). For rectangular channels, equilibrium positions can appear at the centers of the two long walls (LWCE, long-wall-center equilibrium position) and at the centers of the two short walls (SWCE, short-wall-center equilibrium position); thus, zeq/(Lc/2) = zeq/(H/2) for the LWCE and zeq/(Lc/2) = zeq/(W/2) for the SWCE. AR is the aspect ratio of the cross section of the rectangular channel. In all of these tables, the arrows indicate the direction of gradual shift in the focusing position and in other features of migration listed in the Comments column as Re increases steadily over the corresponding ranges given. Besides the brief summary in these tables, a more complete review of the literature is available in  Appendix A.

TABLE I.

Summary of the literature for inertial migration in cylindrical tubes.

Authors (year)MethodαRezeq/RComments
Segré and Silberberg 1961 and 196214–16  Experiments 0.03–0.15 4 → 1400 0.47 → 0.72 Earliest experiments showing inertially induced cross-stream migration in Poiseuille flow 
Cox and Brenner (1968)17  Theory ≪1 ≪1 … Derived formulae for the lateral migration and lift force in tubes with circular and non-circular cross sections 
Matas et al. (2004)18  Experiments 1/42-1/8 134 → 3400 Outer: 0.64 → 0.78; Inner: ≈0.5 Found not only the “outer” Segré and Silberberg focusing annulus but also an inner annulus at high Re 
Shao et al. (2008)19  Direct-forcing fictitious domain method 0.1, 0.15 100 → 1300 Outer: 0.62 → 0.84 Position shifts toward the wall with increasing Re for modest Re; an additional “inner equilibrium position” also observed for Re ≥ 640 
640 → 2200 Inner: 0.6 → 0.32 
Matas et al. (2009)20  Matched asymptotic expansion ≪1 1 → 4000 0.63 → 0.92 Lift force scaled by Um2 decreases strongly with Re above Re = 1 and is much smaller in the pipe than that in the channel 
Choi and Lee (2010)21  Experiments 1/50, 1/23, 1/12 3.2 → 154.8 0.55 → 0.7 Position shifts toward the wall with increasing Re 
Authors (year)MethodαRezeq/RComments
Segré and Silberberg 1961 and 196214–16  Experiments 0.03–0.15 4 → 1400 0.47 → 0.72 Earliest experiments showing inertially induced cross-stream migration in Poiseuille flow 
Cox and Brenner (1968)17  Theory ≪1 ≪1 … Derived formulae for the lateral migration and lift force in tubes with circular and non-circular cross sections 
Matas et al. (2004)18  Experiments 1/42-1/8 134 → 3400 Outer: 0.64 → 0.78; Inner: ≈0.5 Found not only the “outer” Segré and Silberberg focusing annulus but also an inner annulus at high Re 
Shao et al. (2008)19  Direct-forcing fictitious domain method 0.1, 0.15 100 → 1300 Outer: 0.62 → 0.84 Position shifts toward the wall with increasing Re for modest Re; an additional “inner equilibrium position” also observed for Re ≥ 640 
640 → 2200 Inner: 0.6 → 0.32 
Matas et al. (2009)20  Matched asymptotic expansion ≪1 1 → 4000 0.63 → 0.92 Lift force scaled by Um2 decreases strongly with Re above Re = 1 and is much smaller in the pipe than that in the channel 
Choi and Lee (2010)21  Experiments 1/50, 1/23, 1/12 3.2 → 154.8 0.55 → 0.7 Position shifts toward the wall with increasing Re 
TABLE II.

Summary of the literature on inertial migration in square channels.

Authors (year)MethodαRezeq/(H/2)Comments
Chun and Ladd (2006)22  Lattice-Boltzmann method 0.11 200 → 1000 → 2000 0.74 → 0.82 → 0.86 32 particles simulated. Eight equilibrium positions, four at the center of the faces and four at the corners → four face-centered positions vanish → some particles move near the center of the channel 
Di Carlo et al. (2007)12  Experiments 0.18 15 → 90 No focusing → 0.64 Four symmetric positions midway along each of the four walls, with the distance from the wall decreasing as Re increases 
Di Carlo et al. (2009)23  Finite element simulations and experiments 0.1 → 0.8 20–80 0.6 → 0.2 Four symmetric positions midway along each of the four walls. For 0.2<α<0.4,FL=CLρfUm2a3/H, at z/H < 0.5; FL=CLρfUm2a6/H4, at z/H > 0.5, where CL depends on the position 
Bhagat et al. (2009)24  Experiments 0.08 1 → 30 → 300 No focusing → 0.6 → 0.8 Focusing at sidewalls beginning at Re > 20 which progressively improves at higher Re. Method not sensitive to corner focusing 
Choi et al. (2011)25  Experiments 0.075, 0.12, 0.16 2 → 10 → 120 No focusing → 0.56 → 0.74 Focusing starts at smaller Re as α increases → “ring” around the channel perimeter → focus at four positions mid-way along the walls 
Miura et al. (2014)26  Experiments 0.108 200→2400 Centers: 0.72 → 0.44 Four face-centered positions, while four positions at corners appear at Re > 500; further increase in Re moves the corner positions toward the channel corner and the face-centered ones toward the channel center 
500 → 2400 Corners: 0.84 → 1.17 
Hood et al. (2015)27  Regular perturbation theory ≤0.3 10–80 ≈0.6 for α ≈ 0.1 Extended Ho-Leal perturbation analysis for lift force to one extra term, of order a5; use this to explain success of empirical scaling of Di Carlo et al. and show its failure for small α 
Lashgari et al. (2017)28  Immersed boundary method 0.2 100 → 400 → 800 0.57 → 0.61 → 0.56 Four face-centered positions → closer to wall → back to channel center 
Authors (year)MethodαRezeq/(H/2)Comments
Chun and Ladd (2006)22  Lattice-Boltzmann method 0.11 200 → 1000 → 2000 0.74 → 0.82 → 0.86 32 particles simulated. Eight equilibrium positions, four at the center of the faces and four at the corners → four face-centered positions vanish → some particles move near the center of the channel 
Di Carlo et al. (2007)12  Experiments 0.18 15 → 90 No focusing → 0.64 Four symmetric positions midway along each of the four walls, with the distance from the wall decreasing as Re increases 
Di Carlo et al. (2009)23  Finite element simulations and experiments 0.1 → 0.8 20–80 0.6 → 0.2 Four symmetric positions midway along each of the four walls. For 0.2<α<0.4,FL=CLρfUm2a3/H, at z/H < 0.5; FL=CLρfUm2a6/H4, at z/H > 0.5, where CL depends on the position 
Bhagat et al. (2009)24  Experiments 0.08 1 → 30 → 300 No focusing → 0.6 → 0.8 Focusing at sidewalls beginning at Re > 20 which progressively improves at higher Re. Method not sensitive to corner focusing 
Choi et al. (2011)25  Experiments 0.075, 0.12, 0.16 2 → 10 → 120 No focusing → 0.56 → 0.74 Focusing starts at smaller Re as α increases → “ring” around the channel perimeter → focus at four positions mid-way along the walls 
Miura et al. (2014)26  Experiments 0.108 200→2400 Centers: 0.72 → 0.44 Four face-centered positions, while four positions at corners appear at Re > 500; further increase in Re moves the corner positions toward the channel corner and the face-centered ones toward the channel center 
500 → 2400 Corners: 0.84 → 1.17 
Hood et al. (2015)27  Regular perturbation theory ≤0.3 10–80 ≈0.6 for α ≈ 0.1 Extended Ho-Leal perturbation analysis for lift force to one extra term, of order a5; use this to explain success of empirical scaling of Di Carlo et al. and show its failure for small α 
Lashgari et al. (2017)28  Immersed boundary method 0.2 100 → 400 → 800 0.57 → 0.61 → 0.56 Four face-centered positions → closer to wall → back to channel center 
TABLE III.

Summary of the literature on inertial migration in 2D channels.

Authors (year)MethodαRezeq/(H/2)Comments
Rubinow and Keller (1961)29  Theory ≪1 ≪1 … Calculated the lift force based on rotation (slip-spin) of the particles in a uniform unbounded flow 
Saffman (1965)30  Theory ≪1 ≪1 … Proposed a formula for the lift force as a function of slip-shear in a wall bounded flow 
Ho and Leal (1974)6  Perturbation expansion ≪1 ≪1 ≈0.6 Showed that the focusing is due to a balance of lift due to the quadratic velocity profile and wall lift. At low Re, the lateral force scales as FLρfUm2a4/H2 
Schonberg and Hinch (1989)31  Matched asymptotic expansion ≪1 1 → 150 0.64 → 0.7 Position shifts toward the wall with increasing Re; inertial wall lift is limited to the distance HRe−1/2 from the wall 
Feng et al. (1994)32  2D finite element simulation 0.25 40–200 0.45 → 0.52 Position moves closer to the wall with increasing Re. The simulation is 2D, so the particle is a cylinder 
Asmolov (1999)33  Matched asymptotic expansion ≪1 1 → 3000 0.6 → 0.96 Position moves closer to the wall with increasing Re 
Lashgari et al. (2017)28  Immersed boundary method 0.2 200 0.506 Position is closer to the channel center than for the same α at the same Re in ducts 
Authors (year)MethodαRezeq/(H/2)Comments
Rubinow and Keller (1961)29  Theory ≪1 ≪1 … Calculated the lift force based on rotation (slip-spin) of the particles in a uniform unbounded flow 
Saffman (1965)30  Theory ≪1 ≪1 … Proposed a formula for the lift force as a function of slip-shear in a wall bounded flow 
Ho and Leal (1974)6  Perturbation expansion ≪1 ≪1 ≈0.6 Showed that the focusing is due to a balance of lift due to the quadratic velocity profile and wall lift. At low Re, the lateral force scales as FLρfUm2a4/H2 
Schonberg and Hinch (1989)31  Matched asymptotic expansion ≪1 1 → 150 0.64 → 0.7 Position shifts toward the wall with increasing Re; inertial wall lift is limited to the distance HRe−1/2 from the wall 
Feng et al. (1994)32  2D finite element simulation 0.25 40–200 0.45 → 0.52 Position moves closer to the wall with increasing Re. The simulation is 2D, so the particle is a cylinder 
Asmolov (1999)33  Matched asymptotic expansion ≪1 1 → 3000 0.6 → 0.96 Position moves closer to the wall with increasing Re 
Lashgari et al. (2017)28  Immersed boundary method 0.2 200 0.506 Position is closer to the channel center than for the same α at the same Re in ducts 
TABLE IV.

Summary of the literature summary on inertial migration in rectangular channels.

Authors (year)MethodαRezeq/(H/2) for LWCE; zeq/(W/2) for SWCEComments
Bhagat et al. (200834 and 200924Experiments AR = 0.5 ≈0.08 1 → 20 → 200 SWCE: No focusing → 0.76 → 0.86 Focusing at sidewalls beginning at Re > 10 which progressively improves at higher Re. Method not sensitive to corner focusing 
   
AR = 2.5 ≈0.1 1 → 20 → 80 LWCE: No focusing → 0.56 → 0.7 
Masaeli et al. (2012)35  Experiments (AR ≈ 1.4–1.8) 0.09–0.24 10–70 LWCE: ≈0.5 for AR = 1.6 and α = 0.2 Position is influenced by the size of the particle (particle rotational diameter) 
Ciftlik et al. (2013)36  Experiments (AR = 1.6) 0.2 1501500216030004000 LWCE: 0.60.360.450.4no focusing Focusing at the centers along the two long walls; migration toward the channel center at Re ≈ 600; for 300 < Re < 600, focusing also at the centers of the two short walls; focusing along one long wall lost at Re = 3000; focusing lost altogether at Re ≈ 4000 
Zhou and Papautsky (2013)37  Theory and experiments (AR = 2) 0.11–0.3 30–240 LWCE: 0.6 for α ≈ 0.15; 0.52 for α ≈ 0.29 Proposed and validated a two-stage migration model in which shear-gradient lift pushes particles to the wall and then slip-spin lift pushes them to the center of the wall 
Liu et al. (2015)38  Direct numerical simulations (AR = 1, 2, 4, 6) 0.1–0.3 30–468 For AR = 2, SWCE: 0.75–0.85; LWCE: 0.5–0.6 At high Re, multiple focusing positions along long walls and one at center of short wall 
Liu et al. (2016)39  Finite element simulations and experiments (AR = 1, 2, 4, 6) 0.1–0.3 5–300 LWCE: ≈0.6 Proposed a correction to the Ho-Leal lift formula 
Lashgari et al. (2017)28  Immersed boundary method (AR = 2) 0.2 200 LWCE: 0.55 SWCE: 0.82 Focusing at the centers along the two long walls and the tow short walls at Re = 200 
Authors (year)MethodαRezeq/(H/2) for LWCE; zeq/(W/2) for SWCEComments
Bhagat et al. (200834 and 200924Experiments AR = 0.5 ≈0.08 1 → 20 → 200 SWCE: No focusing → 0.76 → 0.86 Focusing at sidewalls beginning at Re > 10 which progressively improves at higher Re. Method not sensitive to corner focusing 
   
AR = 2.5 ≈0.1 1 → 20 → 80 LWCE: No focusing → 0.56 → 0.7 
Masaeli et al. (2012)35  Experiments (AR ≈ 1.4–1.8) 0.09–0.24 10–70 LWCE: ≈0.5 for AR = 1.6 and α = 0.2 Position is influenced by the size of the particle (particle rotational diameter) 
Ciftlik et al. (2013)36  Experiments (AR = 1.6) 0.2 1501500216030004000 LWCE: 0.60.360.450.4no focusing Focusing at the centers along the two long walls; migration toward the channel center at Re ≈ 600; for 300 < Re < 600, focusing also at the centers of the two short walls; focusing along one long wall lost at Re = 3000; focusing lost altogether at Re ≈ 4000 
Zhou and Papautsky (2013)37  Theory and experiments (AR = 2) 0.11–0.3 30–240 LWCE: 0.6 for α ≈ 0.15; 0.52 for α ≈ 0.29 Proposed and validated a two-stage migration model in which shear-gradient lift pushes particles to the wall and then slip-spin lift pushes them to the center of the wall 
Liu et al. (2015)38  Direct numerical simulations (AR = 1, 2, 4, 6) 0.1–0.3 30–468 For AR = 2, SWCE: 0.75–0.85; LWCE: 0.5–0.6 At high Re, multiple focusing positions along long walls and one at center of short wall 
Liu et al. (2016)39  Finite element simulations and experiments (AR = 1, 2, 4, 6) 0.1–0.3 5–300 LWCE: ≈0.6 Proposed a correction to the Ho-Leal lift formula 
Lashgari et al. (2017)28  Immersed boundary method (AR = 2) 0.2 200 LWCE: 0.55 SWCE: 0.82 Focusing at the centers along the two long walls and the tow short walls at Re = 200 
FIG. 1.

Equilibrium positions of neutrally buoyant spherical particles flowing in tubes (a) and square channels (b) from experiments and predictions in the literature. The numbers beside the lines are blockage ratios. In (a), solid symbols are for outer equilibrium positions, while open symbols are for the “inner equilibrium positions.” The inner equilibrium positions for Matas et al.18 and Shao et al.19 fall on nearly the same line. In (b), solid symbols are for the equilibrium position at the centers of the walls, while open symbols are for corner positions. Lines are guides to the eye.

FIG. 1.

Equilibrium positions of neutrally buoyant spherical particles flowing in tubes (a) and square channels (b) from experiments and predictions in the literature. The numbers beside the lines are blockage ratios. In (a), solid symbols are for outer equilibrium positions, while open symbols are for the “inner equilibrium positions.” The inner equilibrium positions for Matas et al.18 and Shao et al.19 fall on nearly the same line. In (b), solid symbols are for the equilibrium position at the centers of the walls, while open symbols are for corner positions. Lines are guides to the eye.

Close modal
FIG. 2.

Equilibrium positions of neutrally buoyant spherical particles flowing in 2D channels (a) and rectangular channels (b) from experiments and theory in the literature. The numbers beside the lines are blockage ratios. In (a), Feng et al. simulated inertial migration of circular particles (equivalent to cylinders in 3D) in a 2D plane Poiseuille flow using a finite element simulation. In (b), solid symbols are for the long-wall-center equilibrium position, LWCE (zeq/(H/2)), while open symbols are for the short wall, SWCE (zeq/(W/2)). Lines are guides to the eye.

FIG. 2.

Equilibrium positions of neutrally buoyant spherical particles flowing in 2D channels (a) and rectangular channels (b) from experiments and theory in the literature. The numbers beside the lines are blockage ratios. In (a), Feng et al. simulated inertial migration of circular particles (equivalent to cylinders in 3D) in a 2D plane Poiseuille flow using a finite element simulation. In (b), solid symbols are for the long-wall-center equilibrium position, LWCE (zeq/(H/2)), while open symbols are for the short wall, SWCE (zeq/(W/2)). Lines are guides to the eye.

Close modal

The results obtained so far from the above studies on the inertial migration of neutrally buoyant rigid spherical particles in Poiseuille flow in straight tubes or channels are summarized in the following few paragraphs, with a more complete summary given in  Appendix A.

In pipes, rigid spherical particles are focused in a single “outer” annulus at a small blockage ratio α and at an additional “inner annulus” at α ≥ 0.06 and Re ≥ 600.14–16,18,21 In square channels, four face-centered equilibrium positions appear at a modest Re and small blockage ratio12,22–26,28 and four corner equilibrium positions sometimes appear upon increasing Re to 200 or higher, for α ≥ 0.1.22,26 Moreover, particles focused at the face-centered equilibrium positions migrate to the corners with increasing Re up to approximately 400 or higher, depending on α.22,26 In rectangular channels, there are always equilibrium positions at the centers of each of the two long channel faces,24,28,34–39 and at moderate Re (20–600), for a modest aspect ratio of the channel cross section (AR = 1.6, 2.0), there can also be two additional equilibrium positions, one at the center of each short channel face.24,28,35,39 The high variability of results in Fig. 2(b) is partly due to the various blocking ratios α and various aspect ratios of the channel cross section.

At low Re, the equilibrium positions of small particles, α ≪ 1, are about 60% of the distance between the centerline and wall in both plane and circular Poiseuille flow and are closer to the center as the blockage ratio α increases. The equilibrium position moves toward the wall with increasing Re for Re O(102).6,12,14–16,18–25,28,31–34 However, when Re exceeds a critical value on the order of hundreds, with smaller values for higher blockage ratios α, experimental18,26,36 and numerical19,28 studies show that the particles begin to shift back toward the centerline, which is the so-called “inner equilibrium position” which can coexist with the outer equilibrium position discussed above. Matas et al.18 first observed the “inner equilibrium position” experimentally in pipe flow; at high Re > 1000, particles focused at both the outer and inner equilibrium positions. However, they found that the inner annulus is much broader and attracts more particles at high Re, while the outer annulus gets sparser and thin. No such reversal, however, was predicted in the perturbation studies of Asmolov33 and Matas et al.,20 despite the fact that these predictions reached values of Re > 1000 where the reversal in the preferred particle position from the outer to inner position was observed experimentally. However, their theoretical studies considered only asymptotically small particle sizes, relative to the gap, and the reverse migration toward the center probably depends on higher order terms in α.

Scaling laws for the lateral force on particles have been extensively studied. After the seminal work of Rubinow and Keller29 and Saffman,30 who proposed formulae based on the rotation and slip-shear, respectively, asymptotic theory became the dominant method for subsequent studies. The work of Ho and Leal,6 Schonberg and Hinch,31 and Asmolov33 and Matas et al.20 all arrived at the well-established expression for the magnitude of the lift force: FL=CLρfUm2a4/H2 for small blockage ratios. Starting from this, Hood et al.27 derived a higher order regular perturbation expansion giving the next order term for the lift force, proportional to ρfUm2a5/H3, valid for larger ranges of particle sizes and Re, while Di Carlo et al.23 and Liu et al.39 supplemented their models with numerical and experimental results and developed empirical scaling laws based on these. The scaling law of Di Carlo et al.23 could be explained using the higher-order perturbation theory of Hood et al.27 

In addition to neutrally buoyant spheres, non-neutrally buoyant spheres in Poiseuille flow have also drawn some attention in the past decades. Repetti and Leonard40 experimentally studied the inertial migration behavior of non-neutrally spheres (α = 0.05, 0.1) suspended in laminar flow through a vertical rectangular channel at moderate Re (10–200), while Jeffery and Pearson41 did a similar study in a vertical pipe. Using matched asymptotic theory, theoretic predictions of the inertial migration of small (α ≪ 1) non-neutrally buoyant spheres in 2D plane Poiseuille flow were made by Vasseur and Cox42 in the near-Stokes regime (Re ≪ 1) and by Hogg43 for Re = 1–100. All these studies concluded that, relative to the neutrally buoyant case, the sphere migrates toward the wall when the sphere leads the undisturbed flow and toward the channel center if it lags behind the fluid. Later 2D finite element simulations of a circular particle (α = 0.25) in plane Poiseuille flow performed by Feng et al.32 showed that the above conclusion is valid when the difference between the particle and the fluid velocity is not large, meaning that the slip Reynolds number is less than 4 or so; otherwise, the particle will migrate toward the channel center whether the particle leads or lags behind the fluid. In addition, Asmolov33 characterized the non-dimensional lift on a non-neutrally buoyant sphere using three dimensionless groups: the channel Reynolds number, the fractional distance from the wall, and the slip Reynolds number, using an orthonormalization method to integrate the equations of the matched asymptotic expansion.

Although the physics of the motion of rigid spheres in Poiseuille flow has become clearer in recent years, there are still observations with incomplete explanations, especially the reverse migration toward the channel center at high Re. In addition, the effects of the rotation (slip-spin) and slip-shear on the migration of rigid spheres are less discussed. Given the difficulty of both well-designed experiments, with neutrally buoyant particles, and of traditional methods of simulation of the continuum Navier-Stokes equations in three dimensions with moving boundaries, efforts have recently been made to employ mesoscopic particle-based methods to solve suspension problems such as these44–48 since such methods are readily adaptable to different geometries and to soft particles, droplets, capsules, etc. Numerical packages for such methods are now available in open-source forms, such as in LAMMPS or HOOMD-Blue, making their use even more attractive. Here, we explore the use of a common particle-based simulation method, dissipative particle dynamics (DPD), to assess the inertial migration behavior of a single rigid sphere in plane Poiseuille flow over a range of Re. The effect of rotation should be captured in this method when Re becomes large. We are particularly interested in assessing limitations of the DPD method, due to the use of fluid beads to transport momentum, including the effect on particle diffusivity, fluid compressibility, and fluid viscoelasticity.

The validity and accuracy of the DPD method for fluid mechanical problems has been assessed by a number of authors. Particularly noteworthy is the thorough study of Kim and Phillips,49 who solved a set of test problems involving flow around spheres and cylinders, addressing in particular the effect of the Reynolds number. In addition to the finding that DPD simulations can be accurate, to be within 10% or better for the drag coefficient, for Reynolds numbers up to around 100, they also found the range of DPD parameters needed to achieve this level of accuracy and helped define what restricts the method to this range of Reynolds numbers. Of particular significance in this regard is the Mach number, the ratio of the characteristic velocity to the speed of sound of the fluid, which is rather low because of compressibility of DPD particles. Zhao and Larson50 also found significant compressibility effects in their test of the DPD method, which, however, did not lead to errors in the velocity greater than about 10% until Re reached 200. The speed of sound can be increased by increasing the repulsive interactions between DPD beads, but Kim and Phillips49 found that this becomes limited by a transition of the fluid rheology to shear thinning. Evidently, strong repulsion turns the DPD fluid glassy, giving it an appreciable relaxation time and hence viscoelasticity. As described below, we find very similar limitations in applying the method to the problem of inertial lift. In addition to problems with accuracy at high Re due to the limited speed of sound, for the problem of inertial migration, we find severe problems with noise due to Brownian motion at low Re, as detailed below.

A description of the DPD methodology will be presented in Sec. II, following by our results and discussion, and ending with our concluding remarks and perspectives regarding both the migration phenomenon and the suitability of the DPD method for tackling this problem.

DPD is a particle-based algorithm first proposed by Hoogerbrugge and Koelman,51 which is applicable to the simulation of simple and complex fluids at the mesoscopic scale. In a DPD simulation, a single DPD bead could represent a very large number of molecules, depending on the degree of coarse-graining chosen. The motion of each bead follows Newton’s second law,

dridt=vi,midvidt=Fij+Fext,
(1)

where bead i located at ri with mass mi and velocity vi is subjected to the pairwise additive inter-particle force Fij from bead j and external force Fext. Fij is composed of three center-to-center directed forces, namely, a conservative force FijC, a dissipative force FijD, and a random force FijR,

Fij=FijC+FijD+FijR.
(2)

The conservative force is a soft repulsion which allows a larger time step relative to that permitted when using atomistic molecular potentials employed in standard molecular dynamics. Since it is purely repulsive, effective attractive interactions between species can only be modeled by making those interactions less repulsive than others, relying on the system pressure to force less repulsive species into preferred contact. Thus, the method, in its standard form, is appropriate for fully condensed states of matter and not for problems involving vapor-liquid equilibrium. The dissipative force slows down the motion of the beads, allowing for frictional forces lost by the coarse-graining to be recovered, while the random force speeds up particles to compensate for the loss of kinetic energy introduced by the dissipative force. The three forces are given as

FijC=aij(1rij/rcC)r^ij,rij<rcC0,rijrcC,FijD=γwDrijr^ijvijr^ij,FijR=σwRrijθijΔt12r^ij,
(3)

where aij is the bead-bead repulsion coefficient, rij = rirj, r^ij=rij/rij, vij = vivj, rcC is the cutoff radius for the conservative inter-particle interactions, γ is the friction coefficient, σ is the amplitude of noise, Δt is the time step, wD(rij) and wR(rij) are dimensionless weighting functions that depend on rij, and θij is the Gaussian white noise satisfying θij(t)=0 and θij(t)θkl(t)=(δikδjl+δilδjk)δ(tt). According to the fluctuation-dissipation theorem,52 we must have wDr=wR(r)2 and σ2 = 2γkBT, which link the dissipative and random forces, where kBT is the Boltzmann temperature of the system. Typically, these functions are given by

wDr=wR(r)2=(1r/rc)s,r<rc0,rrc,
(4)

where rc is the cutoff radius for dissipative and random inter-particle interactions, s is the exponent of the dissipative weight function which is most often taken as s = 2, but we will sometimes take it to be s = 1/2.44 Note that rc without superscript is the cut-off radius for the dissipative and random forces and can be made distinct from the cut-off radius rcC for the conservative force. A modified version of the velocity-Verlet algorithm is used to advance the set of positions and velocities over the short time step Δt.53 

In DPD simulations, the unit of length can be taken to be rcC and that of time is defined as thermal time τ=rcCm/kBT. Besides the exponent s above, the most critical DPD parameters are aff, m, n, rc, rcC,γ, and kBT, where n = ρ/m is the number density of DPD beads and ρ is the mass density. Of these seven parameters, dimensional analysis implies that only four are independent; three of them, typically m, rcC, and kBT, can therefore be set to unity without loss of generality,52 which is done here. [Or, we could define four dimensionless DPD parameters: ãffaffrcC/kBT, ñn(rcC)3, r̃crc/rcC, and γ̃γrcC/mkBT.] The four remaining quantities, and exponent s, determine the physically relevant dimensionless groups, which are the Schmidt number Sc, the ratio of the Mach to Reynolds number Ma/Re, and the ratio of the Weissenberg number to the Reynolds number Wi/Re, which are dependent only on material properties and the gap (which is fixed) and not on the velocity. We discuss these groups later.

An array of aij parameters, which has previously been employed to simulate Poiseuille flow in a plane channel,45 is adapted to our system. The array is given as

aij=wfpw310f3253p103,
(5)

where w, f, and p represent DPD beads that compose the wall, the fluid, and the suspended particle, respectively. While most of the work here uses these repulsive parameters, later we will consider increased values to expand the range of DPD parameter explored in this work.

For pure fluids in the channel without the particle, the velocity profiles from our DPD simulations (examples are shown in Fig. 3) are compared to the Hagen-Poiseuille equation,

vx(z)=12μfdpdxz2H22=ρffxf2μfz2H22,
(6)

where vx(z) is the velocity in the x direction (flow direction) as a function of the z coordinate, −dp/dx is the pressure gradient in the x direction, and fxf is the external driving force in the x direction. The external driving forces in y and z directions are set to zero, i.e., Ffext=fxf,fyf,fzf=[fxf,0,0]. The good match shown in Fig. 3 of the simulated velocity profile to this form demonstrates that our parameter choice yields the Hagen-Poiseuille flow with a no-slip boundary condition and allows us to calculate the viscosity of the fluid.

FIG. 3.

Example DPD velocity profiles of pure fluid with different γ = γff (with the lowest γ giving the highest velocity) and driven by a force fxf=0.005kBT/rcC, fitted by the Hagen-Poiseuille equation (the solid black lines).

FIG. 3.

Example DPD velocity profiles of pure fluid with different γ = γff (with the lowest γ giving the highest velocity) and driven by a force fxf=0.005kBT/rcC, fitted by the Hagen-Poiseuille equation (the solid black lines).

Close modal

Simulations are performed in a box with size Lx×Ly×Lz=60rcC×40rcC×40rcC. A single wall, constructed using three layers of beads in an HCP lattice, was placed at z = −10rc parallel to the x-y plane (see Fig. 4). Due to the periodic boundary condition, this single wall models a pair of walls perpendicular to the z direction. Wall beads are frozen in position and packed at a high density ρw = 61.35, to prevent the fluid beads from penetrating across the wall. Momentum is thereby lost from the fluid beads as it is transferred to the wall beads and is then replaced by the momentum from the external force Ffext applied to the fluid beads, similar to the momentum transfer occurring in conventional application of the no-slip condition of macroscopic fluid mechanics. The thickness of the wall is negligible; thus, the height of the channel is H=Lz=40rcC. The box length was verified to be long enough for well-developed plane Poiseuille flow upstream and downstream of the particle, and the width (perpendicular to both gap and flow direction) was large enough to minimize interactions of periodic images of the sphere.11 

FIG. 4.

Model for a rigid sphere (yellow beads and red bonds) initially in the center of the plane channel with a periodic rigid wall (in pink), where solvent beads are omitted for clarity. The white dashed lines show the box periodic boundary.

FIG. 4.

Model for a rigid sphere (yellow beads and red bonds) initially in the center of the plane channel with a periodic rigid wall (in pink), where solvent beads are omitted for clarity. The white dashed lines show the box periodic boundary.

Close modal

Rigid spheres are constructed by assembling beads into a dense FCC lattice and connecting them by harmonic bonds. The density of beads in the spheroid is ρs = 32 which is less than that in the wall although it still effectively prohibits the penetration of fluid beads into it due to its large size (compared to the effective size of a fluid bead), while keeping the computational cost under control. To render the sphere neutrally buoyant, we used eight-fold lower bead mass for particles in the sphere (mp = 1/8) and applied an eight-fold lower force to the DPD particles in the sphere relative to the fluid (i.e., fxp = fxf/8) so that the sphere has the same inertia and density as an equivalent volume of fluid. In addition, to study the inertial migration of non-neutrally buoyant spheres, we also performed simulations with a sphere of bead mass mp = 1 and driven by the same force on each bead as that applied to the fluid bead (fxp = fxf). This produces an 8-fold higher density of the sphere than of the fluid, allowing the effect of strong sedimentation to be examined, with both sphere and fluid flow driven in the same direction by gravity, so that the sphere tends to move faster than the surrounding fluid. Each bead in the sphere is connected with its nearest, next-nearest, and next-next-nearest beads through a bond with a potential energy

Vr=12krr02,
(7)

where the force constant k is set to 125 for mp = 1/8 and 1000 for mp = 1, a value high enough to restrain stretching or retraction, making the particles non-deformable. r is the distance between two beads, and the bond rest lengths are r0=2b/2,b,3b for bonding to the nearest, next-nearest, and next-next-nearest beads, respectively, with b=0.5rcC being the size of the unit lattice. The spheres are initially placed at the position (0,0,10rcC), which is at the center of the channel.

The time step is Δt/τ ≤ 0.01, and it must also satisfy Δt = 0.1τ/γ to keep the equilibrium temperature to within 2%.54,55 All simulations with γ = 4.5 are first run for one million steps to achieve a confirmed steady state position of the particle and then for nine million more steps to collect fully developed data. For larger γ, where the time step was smaller, we took a time of 104τ to equilibrate the position and 9 × 104τ to collect data. We found that smaller time steps than this or longer runs did not lead to changes in the predictions shown here. Since the velocities are of order unity in DPD units (rcC/τ), over a run time of 104τ, particles travel a distance of around 104rcC over the period used to equilibrate the particle position, which is around 200 channel heights, before averaging starts. Since it has been reported that for low Re = O(10), a distance of around 100 channel heights is necessary for the particle to reach its equilibrium position,24 and less than this for higher Re, our run times are sufficient to equilibrate the position and obtain good averages. Note also that the strain experienced by the fluid is no higher than about 10% per time step, implying that the particle will only travel a fraction of its diameter per time step and, again, we observed no changes on reducing the time step. Simulations are performed with DPD using the Highly Optimized Object-oriented Molecular Dynamics (HOOMD)-Blue software package,56,57 a GPU-based simulation toolkit.

As described in numerous publications,44,53,58–62 the DPD method can exhibit unrealistic behavior owing to unphysical values of transport coefficients, compressibility, and even solvent viscoelasticity in certain cases due to the use of soft beads to represent fluid packets, leading to distortion of their spatial structure under fast flow. For the inertial migration of a non-Brownian sphere in a Newtonian fluid, for a fixed ratio of sphere diameter to gap (i.e., blockage ratio α), one expects the Reynolds number Re to be the only parameter controlling the final position of the particle. However, in our DPD simulations, besides Re, other dimensionless groups that reflect the sphere diffusivity, the fluid compressibility, and fluid viscoelasticity may have an influence on the particle migration.

Among the transport coefficients in DPD simulations, the fluid bead diffusivity (D) is of key interest, where the Schmidt (Sc) number measures the ratio of the fluid momentum diffusivity to bead diffusivity, Sc = μf/ρfD. Groot58 pointed out that for liquids, the diffusivity of a DPD bead is often too fast relative to the momentum diffusivity because of the soft repulsion between DPD beads, allowing beads to pass through each other. The large diffusivity coupled with the slow speed of momentum transfer,61 both of which can be of the same order, leads to the result Sc = Ο(1), while for a real liquid such as water, Sc is Ο(103) or higher. According to the equations (shown in Table V) for D, μf, and Sc of a conventional DPD fluid51 derived by Marsh et al.,63 decreasing kBT and increasing γ, n, and/or rc can result in a smaller diffusivity and a larger Sc. We note that the estimate of the viscosity in Table V ignores any influence of the conservative repulsive interactions between DPD particles, and we will therefore obtain the viscosity directly from the wall shear stress in Poiseuille flow, as discussed below. Decreasing kBT at a fixed DPD friction coefficient γ will decrease the Brownian force coefficient σ according to the fluctuation-dissipation theorem, which thus gives rise to shear-thinning behavior, as the random motion then becomes too weak relative to the conservative repulsion.49 Increasing either n or rc, on the other hand, increases the number of degrees of freedom and the number of pair interactions among particles and sacrifices some of the speed advantages of DPD. We therefore initially choose to vary γ but with an upper limit of 80 to study the effect of the Schmidt number on our results. Our upper limit somewhat exceeds that of Groot and Warren,53 who noted that γ > 32 leads to an instability in the temperature control. Our larger maximum value of γ is due to the much smaller time step Δt = 0.1/γ than used by Groot and Warren53 did, although at the cost of longer run times.

TABLE V.

Transport coefficients62,63 of conventional DPD and modified DPD by Fan et al.44 

CoefficientsConventional (s = 2)Modified (s = 1/2)
Diffusivity, D 45kBT2πγnrc3 315kBT64πγnrc3 
Viscosity, μ 45mkBT4πγrc3+2πγn2rc51575 315mkBT128πγrc3+512πγn2rc551975 
Schmidt number, Sc 12+2πγnrc4270875mkBT 12+32768πγnrc4216372125mkBT 
CoefficientsConventional (s = 2)Modified (s = 1/2)
Diffusivity, D 45kBT2πγnrc3 315kBT64πγnrc3 
Viscosity, μ 45mkBT4πγrc3+2πγn2rc51575 315mkBT128πγrc3+512πγn2rc551975 
Schmidt number, Sc 12+2πγnrc4270875mkBT 12+32768πγnrc4216372125mkBT 

Fan et al.44 proposed to increase Sc to a physically realistic level Ο(103) by altering the weight function exponent of the dissipative force to s = 1/2 and increasing the cutoff radii of both the dissipative force and random force (rc), while keeping the cutoff radius of the conservative force (rcC) at unity. In this way, momentum transfer is increased but without changing the fluid density or the rate of mass transfer. Here we also employed this modified dissipative and random cut-off radius within DPD to study the effect of Sc on inertial migration of rigid spheres.

We note here that the DPD bead diffusivity is indirectly relevant to the migration of a spherical particle through the value of the particle Peclet number (Pep). The particle Peclet number Pep measures the ratio of the rates of advective transport to diffusive transport by Brownian motion and is defined as Pep=UmHDp, where Dp=kBT3πμfa is the diffusivity of the rigid sphere according to the Stokes-Einstein formula. The particle Peclet number can be related to the Schmidt number through the relationship Pep=DpDPe=DpDReSc where again D is the DPD bead diffusivity and Pe is the DPD bead Peclet number. Thus, for a given Reynolds number, a low value of the Schmidt number, around unity, implies that the particle Peclet number will be abnormally low, relative to a real fluid at the same Reynolds number. This, in turn, implies that the sphere diffusivity will be more important than it should be in the migration phenomenon, relative to convection, in DPD simulations. Since most experimental inertial migration studies are carried out in the non-Brownian limit of very high particle Peclet number, Pep > 105, a smaller particle Peclet number, Pep < 1000, for instance, implies that Brownian motion will be significant enough that long simulation runs will be needed to average out diffusive drift of the particle. In fact, for Re < 50, we find that it can be impossible to obtain particle focusing in DPD simulations.

Another issue relevant for a fluid composed of DPD beads, arising from their soft repulsion, is fluid compressibility. The soft repulsion between the beads results in a low sound velocity, which means that the DPD fluid flow can become supersonic when Re is high.58 The relevant dimensionless group is the Mach number Ma = Um/cs, where cs is the speed of sound, which is given by cs2=p/ρf,64 where p is the pressure. For a DPD fluid, p=nkBT+πaijn2rcC430 for n > 2.53,65 Therefore,

cs2=kBTm+πaijnrcC415m.
(8)

We note that the classical values of DPD parameters, including the value aff = 18.75, were chosen so that the compressibility of the DPD fluid, relative to thermal energy, matches that of typical condensed-phase fluids. However, we shall see that this choice, while being suitable for static, thermodynamic calculations, does not mean that the DPD fluid is effectively incompressible in fluid flow.

The soft repulsion of DPD beads can also give rise to viscoelasticity of the fluid when Re becomes high. This can be quantified by the Weissenberg number (Wi), which can be defined as Wi=γ̇λ, where γ̇ is the shear rate and λ is the relaxation time, estimated by Mai-Duy and co-workers65 to be λγrcC/aff. Below, we provide a more relevant estimate of λ from the shear rate at which shear thinning becomes pronounced. In plane Poiseuille flow, γ̇ is roughly the average shear rate γ̇=Um/(H/2).

The hard suspension particles here are chosen to have a diameter a = 8rc, and thus the blockage ratio is α = a/H = 0.2. Smaller spheres suffer too much diffusion to reach to equilibrium positions in our simulations since the asymptotic theory6,18,20,31,33 predicts that lift force becomes very small for small particles, scaling as FLa4/H2. The behavior of a larger blockage ratio is less characteristic of experimental systems and requires too many DPD beads to be readily simulated with our methods.

Before presenting results, we note that at Re above around 500, we cannot keep the Mach number Ma below around unity, consistent with the finding of Kim and Phillips.49 Figure 5(a) shows that the Mach number Ma, derived from the formula for the speed of sound in Eq. (8), exceeds unity as the Reynolds number becomes large with a higher ratio of Ma/Re as the Schmidt number Sc increases. The Weissenberg number Wi=λγ̇, defined using the damping time λγrcC/aff by Mai-Duy et al.,65 exceeds 0.5 at large Re for fluid with large friction parameter [γ = 40, given by the blue triangles, and γ = 80, given by the green squares in Fig. 5(b)]. We expect that at Wi > 1, there should be significant shear thinning, which is not evident in Fig. 5(c) despite the large values of Wi reported in Fig. 5(b). Thus, the definition of the relaxation time suggested by Mai-Duy et al.65 seems not to be relevant if we are interested in the conditions under which the DPD becomes shear thinning. We do find significant shear-thinning in our simulations when fluid beads have a highly repulsion parameter (aff = 100); see Fig. 5(d). (Equivalently, shear thinning is obtained if kBT is made small so that the dimensionless group ãffaffrcC/kBT becomes large.) Taking γ = 40, Fig. 5(d) shows that if aff is increased to a high value of around 100, while increasing the other repulsion parameters in Eq. (5) in the same ratio (to minimize fluid penetration into the sphere or the wall), shear-thinning occurs. We find that a high value of aff does increase the speed of sound in a DPD fluid, and it also increases the viscosity, which decreases Re for a fixed fluid velocity. Thus, to obtain the same Re, the velocity must be increased, which again increases Ma, counteracting the effect of high aff on the speed of sound. Equations for the viscosity of a DPD fluid in Table V contain no dependence on aff since they were derived under the assumption of negligible aff. However, the effect of aff on DPD fluid viscosity has been pointed out by Azarnykh et al.,66 and Zhao et al.67 found that for γ > 10, the shear viscosity increases with increasing aff. Since increased aff does not help us reach a higher Re and introduces shear thinning, we pursue that option no further. The onset of shear thinning with increased repulsive parameter aff is consistent with the findings of Kim and Phillips.49 

FIG. 5.

Dependence on Reynolds number Re of (a) the Mach number Ma and (b) the Weissenberg number Wi for DPD simulations. (c) Dependence of viscosity μf of the fluid on the shear rate at the wall for the DPD simulations, where μf is calculated by fitting the Hagen-Poiseuille equation [Eq. (6)] to our simulation results; additional simulation details can be found in the caption of Fig. 6. (d) Dependence of viscosity μf of the DPD fluid on the repulsion parameter aff, with γ = 40, and with the dashed line representing the value estimated by the equation in Table V derived by Marsh et al.63 

FIG. 5.

Dependence on Reynolds number Re of (a) the Mach number Ma and (b) the Weissenberg number Wi for DPD simulations. (c) Dependence of viscosity μf of the fluid on the shear rate at the wall for the DPD simulations, where μf is calculated by fitting the Hagen-Poiseuille equation [Eq. (6)] to our simulation results; additional simulation details can be found in the caption of Fig. 6. (d) Dependence of viscosity μf of the DPD fluid on the repulsion parameter aff, with γ = 40, and with the dashed line representing the value estimated by the equation in Table V derived by Marsh et al.63 

Close modal

In Fig. 6(a), by limiting data to those for which Ma < 0.8, i.e., keeping the DPD fluid subsonic, we obtain equilibrium positions (colored symbols) of the neutrally buoyant sphere with α = 0.2 that are around 50% of the distance between the channel center and wall over the range Re ≈ 50–500, which is near the position for spheres with the same or slightly larger blockage ratio in 2D28,32 (α = 0.25 in 2D simulations by Feng et al.32) and in rectangular28,35 channels. However, the equilibrium position of a sphere with α = 0.2 in a square channel observed numerically by Lashgari et al.28 and that in a rectangular channel (AR = 2) observed experimentally by Ciftlik et al.36 are around 60% the distance from the channel center to the wall. Thus, it appears that the equilibrium position of a sphere with α = 0.2 in a 2D channel (AR → ∞) is somewhat closer to the channel center than that in a square or rectangular channel with AR ≤ 2. In addition, according to our simulations, the equilibrium position moves toward the channel center with increasing Re above 100. From our literature review above, the equilibrium positions of neutrally buoyant spheres move toward the wall with increasing Re for ReO(102);6,12,14–16,18–25,28,31–34 however, when Re exceeds a critical value on the order of 200, with smaller values for higher blockage ratios α, experimental18,26,36 and numerical19,28 studies show that the particles begin to shift back toward the centerline, in accord with our simulations. The velocity of the sphere lags that of the fluid, as shown in Fig. 6(b), which is consistent with findings for neutrally buoyant spheres by previous studies.6,32

FIG. 6.

Equilibrium positions (colored and gray symbols) of neutrally buoyant (mp = 1/8, mf = 1, fxp = fxf/8) (a) and non-neutrally buoyant (mp = mf = 1, fxp = fxf) (c) rigid spheres in Poiseuille flow from DPD simulations with Ma < 0.8, compared to those from the literature [black symbols in (a)]. DPD results for s = 2.0, rcC=rc=1.0, are given for γ = 4.5 (red circles), γ = 40 (blue triangles), and γ = 80 (green squares), while for s = 0.5, rcC=1.0, γ = 4.5, results are given for rc = 1.0 [gray diamonds in (c)], rc = 1.5 [orange pentagons in (c)], rc = 1.8 [purple hexagons in (c)], and rc = 2.0 [cyan asterisks in (c)]. In all cases (including results from the literature), the blockage ratio α is 0.2, except for that of Feng et al.25 for which α = 0.25, where these simulations are for a disk in a 2D flow, or, equivalently, for a cylinder in plane Poiseuille flow. The viscosity μf is calculated by fitting the Hagen-Poiseuille equation [Eq. (6)] to our simulation results, while Sc is calculated from the equations in Table V. Note that the purple hexagons and cyan asterisks in (c) are confined to the lower left corner of (c). Velocities of neutrally (b) and non-neutrally buoyant (d) spheres at the equilibrium position are given for μf = 1.13 and Sc = 0.7.

FIG. 6.

Equilibrium positions (colored and gray symbols) of neutrally buoyant (mp = 1/8, mf = 1, fxp = fxf/8) (a) and non-neutrally buoyant (mp = mf = 1, fxp = fxf) (c) rigid spheres in Poiseuille flow from DPD simulations with Ma < 0.8, compared to those from the literature [black symbols in (a)]. DPD results for s = 2.0, rcC=rc=1.0, are given for γ = 4.5 (red circles), γ = 40 (blue triangles), and γ = 80 (green squares), while for s = 0.5, rcC=1.0, γ = 4.5, results are given for rc = 1.0 [gray diamonds in (c)], rc = 1.5 [orange pentagons in (c)], rc = 1.8 [purple hexagons in (c)], and rc = 2.0 [cyan asterisks in (c)]. In all cases (including results from the literature), the blockage ratio α is 0.2, except for that of Feng et al.25 for which α = 0.25, where these simulations are for a disk in a 2D flow, or, equivalently, for a cylinder in plane Poiseuille flow. The viscosity μf is calculated by fitting the Hagen-Poiseuille equation [Eq. (6)] to our simulation results, while Sc is calculated from the equations in Table V. Note that the purple hexagons and cyan asterisks in (c) are confined to the lower left corner of (c). Velocities of neutrally (b) and non-neutrally buoyant (d) spheres at the equilibrium position are given for μf = 1.13 and Sc = 0.7.

Close modal

We also show in Fig. 6(c) the equilibrium positions of non-neutrally buoyant spheres, with the sphere 8 times as dense as the surrounding fluid, as discussed in Sec. II. In this case, when Sc is of order Ο(1)–Ο(10), the sphere’s positions are around 60% of the distance from the channel center to the wall, closer to the wall than that of neutrally buoyant spheres at the same Re. A drift of the equilibrium position toward the wall with increasing Re is observed when Re ≈ 50–200, above which the direction of drift reverses and the particle migrates to the channel center with further increases in Re. Figure 6(d) shows that the sphere leads the fluid, the more so with increasing Re. The dependence of the equilibrium position of the neutrally and non-neutrally buoyant spheres on the difference in slip velocity relative to the surrounding fluid in our simulations agrees with previous experimental40,41 and theoretical42,43 studies, in which the sphere migrates closer to the wall when it leads the undisturbed flow than if it lags behind the fluid. Moreover, Feng et al.32 observed that when the slip between the denser non-neutrally buoyant sphere and the fluid is relatively low, the sphere migrates toward the wall; however, when the slip is high, the sphere reversely migrates toward the channel center. This may be the reason for the reversal of the equilibrium position of non-neutrally buoyant spheres in our simulations.

Note in Fig. 6(c) that the data for s = 0.5 are closer to the centerline, or even on the centerline, even for small Re, if the Schmidt number is very high, even though the Mach number, as defined using the speed of sound given above, is less than 0.8. We suspect that the simple expression for the speed of sound is no longer relevant when Sc is large. Azarnykh et al.66 found that the viscosity and speed of sound depend on wave number and are controlled by a dimensionless parameter Ω0 that can be written (in our notation) as Ω0=2π/45nrc4γ2/mkBT, for s = 2, and a value more than 20 times higher than this, Ω0=π/5nrc4γ2/mkBT for s = 1/2. At small distance scales, the fluid properties cease to be those of the bulk, and this small-scale breakdown occurs at larger length scales as Ω0 increases. This results from the finite range of interaction of the DPD particles, which leads to a breakdown of the continuum Navier-Stokes fluid dynamics at small length scales. The DPD fluid particles form a fluid that is “lumpy” like porridge or oatmeal, and its “lumpiness” increases with increased γ, rc, or decreased s. It acts as a continuum fluid at large length scales, but at smaller scales, the “lumpiness” likely affects the flow. This is relevant to flow around a sphere and to sphere migration since the viscous boundary layer around the sphere, which influences migration, scales as Rep1/2 units of sphere diameters, where the particle Reynolds number is Rep = Reα2, with a blockage ratio α = 0.2 for our simulations. Our spheres have a diameter of a=8rcC so that at Re = 500, the boundary layer around a sphere becomes comparable to the size of a single DPD fluid particle even when rc=rcC, and this occurs for much smaller values of Re when rc>rcC and s = 0.5, as is the case for many of the high-Sc simulations with s = 0.5 in Figs. 5(a)–5(c). Note the very similar expressions for the Schmidt number in Table V and the expression for Ω0 above. Thus, a high Schmidt number practically guarantees that DPD simulations will not resolve thin boundary layers around the spheres that occur at high Rep. Hence, high Sc may help suppress Brownian motion and diffusion at low Rep but worsen resolution of fine-scale fluid structure at high Rep. This is not surprising; although DPD simulations have no grid, they are still affected by discretization error as just as other numerical methods are, and one must remain mindful of this limitation. When boundary layers become thin (as they do at high Rep), finer resolution is needed, which requires appropriate choice of DPD parameters.

Having established in Figs. 5 and 6 the range of values of Re for which the DPD method gives results in semi-quantitative agreement with experiments and other simulation methods for flow in a slit, the migration of spheres and other particles or droplets in other geometries, such as rectangular or square ducts, curved channels, and other geometries, could be carried out for single or multiple spheres, and the results compared with experimental findings discussed in Sec. II above.

We show in  Appendix B the dependencies of the ratio Ma/Re on DPD parameters based on the formulae in Table V and find that there is only limited scope for increasing Re beyond 500 when using “classical” values of DPD parameters, without increasing Ma to unacceptably high values. The exception to this seems to be in the fluid-fluid bead repulsive parameter aff; increasing this value, using the formulae in Table V, in principle monotonically reduces both Wi/Re and Ma/Re. However, we find that the simple formulae for both viscosity and relaxation time are highly inaccurate for large aff and the increased viscosity and increased shear thinning for large aff make both Ma and Wi high for large Re.

The unusually high ratio of Wi/Re is ultimately a consequence of using a limited number of DPD beads to represent a vastly larger number of real molecules. At the same time, to keep the noise in the velocity field to within a reasonable level, the velocity needs to be reasonably high, relative to the thermal motion, and thus the Mach number tends to be high. In short, it is difficult to make all dimensionless groups: Wi, Sc, and Ma attain physically realistic values, especially at high Re. We are also limited at the lower end of the range of Re since Re below 50 produces large noise as the Peclet number Pes decreases with Re, and this leads to large statistical noise due to diffusion of the sphere, especially when Sc (and therefore Pes) is reduced to keep Wi and Ma from becoming too high. Thus, the best one can do is to work in regimes in which the relevant groups that are vanishingly small in experiments (Wi, Ma, and Pes−1) are small enough in DPD simulations that their values have small enough impact on the results. DPD is a mesoscopic method, and all dimensionless groups in DPD have “mesoscopic” values, in contrast to continuum methods where some parameters, such as Ma, Wi, and Re, can be set rigorously to zero, or, in the case of Pes, to infinity. This limits the range of conditions under which DPD gives realistic results, in our case to only a decade or so in Re, and forces one to consider carefully the magnitudes and influences of dimensionless groups that are not of experimental relevance.

We are also forced to take the sphere to be of appreciable size compared to the gap since we need the sphere to be much larger than a DPD bead. Keeping it significantly smaller than the gap limits us to a very narrow range of blockage ratios. This again arises from the mesoscopic character of the DPD method and the limit that computational resources place on the number of DPD beads that can be used. In addition, even when we stay within the limits of the DPD method, the results shown in Fig. 6 are only accurate to within 10% or so, in part because the Mach number, while limited to lie below 0.8, is certainly not vanishingly small. Thus, while the DPD method is flexible and relatively easy to use, for problems involving transport, care is required to establish the range of its validity, which can be quite limited, and its accuracy within this limited range is semi-quantitative. We note that in our previous paper,11 in which DPD simulations were used to determine the focusing position of droplets, the Reynolds numbers were low enough that the Mach number rarely exceeded unity, and so the difficulties described here were much less of a concern. We also note that for fluid flow problems without suspended particles, we have found that DPD simulations can give complex velocity fields, including secondary flows that are within 10% of the exact solutions over a range of Re from 5 to 200.50 

DPD simulations have been performed to study the inertial focusing position of a single neutrally or non-neutrally buoyant rigid sphere with diameter equal to 20% of the gap in plane Poiseuille flow. The exaggerated effects of sphere diffusivity, fluid compressibility, and viscoelasticity in DPD simulations relative to experiments are assessed for various values of DPD parameters. Depending on the choice of the exponent s in the radial dependence of the dissipation function wDrij1rs, both the Schmidt number Sc of the fluid and the particle Peclet number Pep of the sphere can be less than for typical fluids at the same Reynolds number Re by a factor of up to 1000. Nevertheless, these relatively low values do not much affect our results or their agreement with literature results. In particular, a low Schmidt number (Sc), which is commonly deemed a deficiency of DPD fluids, gives physically realistic results for particle migration for Reynolds numbers up to 500. The particle Peclet number Pep remains asymptotically high for reasonably sized spheres and Re > 50 even for low Sc. Attempts to raise Sc of DPD fluids to values closer to realistic fluids by increasing the range of the viscous interactions of DPD particles significantly erodes the ability of the method to resolve viscous boundary layers in flows at an elevated Reynolds number. In particular, when the exponent s is chosen to be 1/2 so that Sc is also high, the Mach number (Ma) is found to reach order unity or higher for moderate or high Re (=O(10)–O(102)). Because of these high values of Ma, we find a dependence of the sphere’s focusing position on parameters other than Re at a high Reynolds number (Re) and a migration of the sphere all the way to the centerline even at modest Re, evidently because of the inability of DPD fluid to resolve thin boundary layers when Sc is high.

If the Mach number is kept within the range Ma < 0.8, and Sc < 100, the neutrally buoyant sphere is focused at a position around midway between the channel center and the wall, and with increasing Re from around 100 to 500, it migrates toward the channel center. The equilibrium position of a non-neutrally buoyant sphere is around 60% of the distance from the channel center to the wall, closer to the wall than that of neutrally buoyant sphere. The neutrally buoyant sphere lags behind the fluid, pushing the sphere closer to the center, while the non-neutrally buoyant sphere leads the fluid, pushing it toward the wall. Our DPD results give spurious predictions for the particle position when Ma exceeds unity or so.

Overall, our results carry implications for the use of DPD simulations for fluid mechanics and transport. In particular, because all dimensionless groups are potentially in play in DPD, and none are asymptotically small or large, one must check carefully that normally irrelevant parameters, such as the Mach number or the Weissenberg number, are not too high to be considered “small enough” for reasonable predictions to be made. This can be treacherous since the fluid parameters such as the speed of sound and the viscosity are not controllable directly but only through artificial DPD parameters such as friction coefficient and repulsive interaction strength, each of which affects more than one physical property. Because changing one DPD parameter affects multiple physical constants, and all dimensionless groups are in play, a reduction in the magnitude of one dimension group typically means increasing that of some other group. Hence, quantification of error has proven elusive so far for DPD.

One may then ask: given these serious shortcomings, why DPD should be used at all for transport problems? In general, the motivation for using DPD is its versatility. It is quite easy to apply DPD to complex fluids which would be highly challenging to study by more conventional numerical methods, fluids for which the standards for precision, in measurement and prediction, are much less demanding than for conventional fluid mechanics. For example, DPD could readily be applied to a solution containing mixtures of polymers and particles or droplets, to asymmetric, deformable colloidal particles in flow, and to non-ideal geometries, such as channels with wavy walls, and to combinations of these. Of course, the uncertain accuracy of DPD and its ability to be applied to almost arbitrarily complex fluid mixtures mean that the predictions provided by DPD for complex fluids are not always easy to validate by other, more reliable methods. The example considered here of migration of a single rigid sphere is therefore valuable for calibrating DPD and telling us the limits of the method, both in the range of Reynolds numbers (around Re = 50–500) over which the DPD works at all and the limits on its accuracy over this range. The error appears to be about a 10% deviation of the focusing position of the particle from the position predicted by a traditional numerical method, as shown in Fig. 6. (We note that the only numerical methods that have gone to higher Re than we have attained with DPD are the direct-forcing fictitious domain method and lattice Boltzmann method, as listed in Tables I–IV, and the errors in these methods were also not clearly assessed.)

Another example67 for which the accuracy of DPD predictions of polymer rheology was assessed recently showed roughly a 20% error in the polymer relaxation time in the presence of hydrodynamic interactions. This may be acceptable for such fluids, whose experimental characterization rarely exceeds accuracy of this order. Similar calibration tests need to be performed for other classes of problems using DPD.

Y.H. was supported by a grant from the China Scholarship Council. This work used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by the National Science Foundation (NSF), Grant No. TG-CHE140009. R.G.L.’s contribution is supported by the NSF under Grant No. CBET1602183. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the National Science Foundation. This research was also supported through computational resources and services provided by Advanced Research Computing at the University of Michigan, Ann Arbor.

1. Spheres in 2D channels and circular tubes

a. Theory and simulations

Soon after the “tubular pinch effect” found by Segré and Silberberg,14 this surprising phenomenon quickly triggered many subsequent studies, including attempts to provide a theoretical explanation. Rubinow and Keller29 calculated the lift force based on the rotation (slip-spin) of the particles in a uniform unbounded flow, while Saffman30 proposed a formula for the lift force as a function of slip-shear in a wall bounded flow, where this force was taken to be independent of whether or not the particles were allowed to rotate. Both of these theories are valid only at asymptotically small channel Reynolds numbers. “Asymptotically small” here means that Re ≪ α ≪ 1. Cox and Brenner17 derived implicit analytic formulae for the lateral migration (velocity and lift force) of spheres in laminar flows through tubes with circular and non-circular cross sections, taking account of the effects of the base flow, solid boundaries, inertia, and particle buoyancy, using a regular perturbation solution of the Navier-Stokes equations. However, the formulae they derived are integrals that do not explicitly give the direction and magnitude of the lift force. Subsequently, Ho and Leal6 used an asymptotic theory based on the approach of Cox and Brenner17 for the inertial migration of rigid spheres in plane Couette and plane Poiseuille flows, valid when Re ≪ 1, α2 ≪ 1. They concluded that under these conditions, the lateral force originates from the shear flow, rather than from wall-induced lag (slip-shear) and wall-induced particle rotation (slip-spin), which had been proposed earlier to explain the migration observed by Segré and Silberberg.16 Moreover, Ho and Leal6 found an explicit expression for the lift force that scales as FLρfUm2a4/H2, by retaining only leading-order terms. Their predicted equilibrium position for a particle in plane Poiseuille flow precisely agreed with the position found by Segré and Silberberg16 at their lowest Re in laminar flows in pipes. Although Segré and Silberberg16 acquired their data for Re ≈ 4–1400, Ho and Leal6 assumed Re ≪ 1 so that inertia could be taken to be small throughout the channel. To solve the problem at Re > 1, Schonberg and Hinch31 solved the resulting singular perturbation problem using the method of matched asymptotic expansions with Oseen’s linearized Navier-Stokes equation and were able to extend the predictions to Re = 1–150, for a sphere small enough that α ≪ 1 and Rep ≪ 1 in plane Poiseuille flow, where Rep = Reα2 ≪ 1. The equilibrium position was found to shift toward the wall with increasing Re, which is in agreement with the experiments of Segré and Silberberg.16 Furthermore, Schonberg and Hinch31 showed that the inertial wall lift resulting from the wall-induced asymmetry of the flow around the particle, which pushes the particle away from the wall, was limited to the distance HRe−1/2 from the wall. This scaling explains why the particle is driven closer to the wall at higher Re by the shear-gradient lift. Later, Asmolov33 extended the study of Schonberg and Hinch31 to the range Re = 1–3000 using the same theoretical approach and proposed an empirical equation for the lateral lift force FL=CLρfUm2a4/H2, where CL is the lift coefficient, which is dependent on both the Reynolds number and lateral position. The above scaling of the force with velocity, particle radius, and gap is consistent with the scaling law of Ho and Leal.6 Matas et al.20 predicted the lift force exerted on a rigid sphere in a pipe under Poiseuille flow for Re = 1–4000 using the same method and found that its magnitude decreases strongly with Re and is much smaller in the pipe than that in the channel. Recently, Hood et al.27 re-assessed and improved on the regular perturbation expansion of Ho and Leal6 and showed how it remains valid for Re up to 80, for small particles despite its formal limitation to Re < 1. They also derived a higher order asymptotic model with FL=CL,4ρfUm2a4/H2+CL,5ρfUm2a5/H3 to predict the lateral force and equilibrium particle position based on the particle size, valid for a larger range of particle sizes and Reynolds numbers than for the theory of Ho and Leal.6 

With the advance of computing capabilities, numerical studies have become powerful enough to investigate the inertial migration of particles in 2D or 3D Poiseuille flow. Feng et al.32 simulated the inertial migration of circular particles (equivalent to cylinders in 3D) with α = 0.25 in 2D plane Couette and Poiseuille flow using a finite element simulation at Re = 40–200. They observed that particles migrate to an equilibrium position around 0.5 times the half height of the channel from the centerline in plane Poiseuille flow, and this position moves closer to the wall with increasing Re, which is in qualitative agreement to the previous theory. Finally, Shao et al.19 used a direct-forcing fictitious domain method to simulate the inertial migration of spherical particles with blockage ratios of 0.1 and 0.15 in circular Poiseuille flow at Re = 100–2200. In addition to finding a shift of the equilibrium position of the sphere toward the wall with increasing Re for modest Re, an additional “inner equilibrium position” was also observed for Re ≥ 640 for these spheres of a finite blockage ratio α in a relatively long tube with a periodic boundary condition.

b. Experiments

Matas et al.18 experimentally studied the inertial migration of diluted spherical particles (α = 1/42–1/8) in circular Poiseuille flows at Re = 134–3400 and confirmed the occurrence of the Segré and Silberberg annulus. Their results for large particles (α = 1/17–1/8) agree with those of experiments with similar α performed by Segré and Silberberg,16 while for smaller particles, α = 1/42, their results agree with the predictions of asymptotic theory6,31,33 up to Re ≈ 2400, in which the annulus moves toward the wall with increasing Re. However, in addition to this “outer annulus,” another “inner annulus” with a radial position 50% of the distance from the center of the pipe radius to the wall was observed when Re > 1200 for “large” particles with blockage ratios α greater than around 1/20 and at Re = 2400 for the smallest particles studied, with α = 1/42. This inner annulus coexisted with the outer annulus from Re about 800 to around 2400 although the percentage of particles in the outer annulus diminished toward zero with increasing Re. They tried to explain this result using asymptotic (α ≪ 1 and Rep ≪ 1) theory but failed, finding that although there is a local minimum in lateral force near the position of the inner annulus seen in their experiments, there is no new zero of the lateral force, i.e., no new equilibrium positions in the lateral direction, in the pipe geometry.18,20 They therefore attributed the inner equilibrium position in their experiments to finite-particle-size effects, i.e., higher order terms in the blockage ratio α. Choi and Lee21 measured the inertial migration of neutrally buoyant spherical particles (α = 1/50, 1/23, 1/12) in a micron-scale-diameter pipe at Re = 3.2–154.8 using a digital holography technique. The Segré and Silberberg annulus was observed at Re = 154.8 for α = 1/23, and at a smaller Re = 19.4 for α = 1/12, qualitatively consistent with the predicted decrease in lift force with decreasing blockage ratio α. Given the drift of the outer annulus toward the wall and the appearance of an inner annulus for blockage ratios that are not asymptotically small, one might infer that for still larger blockage ratios (α > 1/12), the outer annulus might be lost at Re large enough for an inner annulus to appear, and so there might in that case be a single annulus that drifts toward the wall up to a critical Re around 1000, followed by a drift back toward the center for still higher Re. Such a non-monotonic trend is in fact observed in studies discussed below.

2. Spheres in channels with square or rectangular cross section

a. Experiments

With the advent of microfluidics, more experiments on lateral drift have been performed in square or rectangular channels than in circular tubes. Di Carlo et al.12 seem to have been the first to study experimentally the inertial migration of spherical particles in square channels and observed that dilute spherical particles with a blockage ratio α = 0.18 migrated to four symmetric equilibrium positions midway along each of the four walls, with the distance from the wall decreasing as Re increased from 15 to 90. In the same geometry, Choi et al.,25 using particles with smaller blockage ratios from 0.07 to 0.16, observed using digital holography that the particles do not focus over the length of the channel (around 640 channel gaps) at the lowest Re, less than Re = 20 or so, with focusing starting at lower Re as α increases. Above this Re, the particles drift to a distance from the wall similar to that predicted for slits by Ho and Leal6 (i.e., 60% of the distance from the center of the channel). At such low Re, there was no tendency of particles to center themselves along the wall, but they spread into a square “ring” following the channel perimeter. With increasing Re (>20), the particles began gradually to focus toward the four equilibrium positions mid-way along the walls, to a gradually increasing extent with Re increasing up to 240, and with increasing α. Bhagat et al.24 had found earlier behavior similar to that of Choi et al.25 over a similar range of Re in square channels with α > 0.07, except that rather than using a three dimensional imaging technique as Choi et al., Bhagat et al. imaged fluorescent particles through a horizontal side wall, making it impossible to distinguish particles accumulating at the corners from those spread out along the vertical side walls. Extending the study of migration in square channels to the range Re = 200–2400, Miura et al.26 again observed the four equilibrium positions at the centers of the four walls at Re < 200 for particles with α = 0.108. However, four additional equilibrium positions appeared near the corners for higher Re, reaching the corners at Re ≈ 500. The particles at the wall-midpoint equilibrium positions drifted away from the walls toward the channel center, and the four corner equilibrium positions were driven slightly closer to the corners along the diagonals with increasing Re.

In channels with a rectangular cross section, Bhagat et al.24,34 found experimentally that particles in high-aspect-ratio (AR = 2.5) rectangular channels with a blockage ratio (using the thinnest channel dimension) α ≈ 0.1 were spread in a band along the two long walls at Re < 40 similar to the behavior in a slit, with additional focusing of these bands toward the midpoints of the long walls for higher Re up to 200 or so, except that particles also appeared in the corners at these higher Re, similar to the behavior they observed for square cross sections. Masaeli et al.35 investigated the inertial migration of spherical (α ≈ 0.09–0.24) and prolate particles in rectangular channels with AR≈1.4–1.8 at Re = 10–70. The spheres focus only at the midpoints of the two long walls, and the equilibrium positions are about the midway between the channel center and the long wall for spheres with α ≈ 0.2 within the Re range they studied. Ciftlik et al.36 also investigated the inertial migration of spherical particles in a rectangular channel with a similar aspect ratio (AR = 1.6). For a blockage ratio of 0.2, they again found focusing at the midpoints along the two long walls when 150 < Re < 2160 with no focusing at the corners although they did find particles at the corners for smaller blockage ratios. At Re = 3000, focusing along one wall was lost apparently due to some asymmetry in the inlet particle distribution, and finally focusing was lost altogether when Re ≈ 4000 because the flow became turbulent. Notably, again for α = 0.2, the two equilibrium positions started migrating toward the center of the channel to the inner equilibrium positions around Re ≈ 600, which is consistent with the behavior in circular pipes discussed above. In addition, for α = 0.2, for 300 < Re < 600, particles also focused at the midpoints of the short walls, creating four equilibrium positions, similar to that in a square cross section. Bhagat et al.24 summarized their results in a “phase plot” with the blockage ratio and particle Reynolds number as axes, showing that focusing could be attained at blockage ratios above about 0.07 for particle Reynolds numbers Rep above about 0.1. Note, however, that this “phase plot” depends on the length of the channel; since Brownian motion is generally negligible for particles larger than 5 µm or so, even miniscule drift velocities will eventually focus very dilute particles if the channel is long enough. Ciftlik et al.36 worked out the practical limitations on inertial-focusing separation strategies posed by considerations of degree of focusing, pressure drop limitations, onset of turbulence, and other factors.

b. Theory and simulation

Zhou and Papautsky37 analyzed theoretically the physics of inertial particle migration in rectangular channels and validated through experiments a two-stage migration model: the shear-gradient-induced lift drives the particles to the channel wall until stopped by the wall lift (stage I); once these dominant forces come into balance and the particle establishes an equilibrium distance from the wall, a weaker rotation-induced (slip-spin) lift force can exert its influence and drives the particles along the wall toward its midpoint (stage II). Experiments by Zhou and Papautsky also allowed measurement of the lift coefficient during the two stages of migration. They also inferred that the rotation-induced lift scales more strongly with particle size than do the dominant shear and wall lift forces, leading to particle-size-dependent differences in focusing. Liu et al.38 used direct numerical simulations (DNSs) to determine the equilibrium positions of spherical particles in rectangular channels with AR = 2 and found that in addition to the two equilibrium positions at the centers of the long channel faces at low Re that: (a) increasing Re beyond a critical value of around 100 can stabilize positions near the center of the short channel faces, with the critical Re for stabilization dependent on the blockage ratio, and (b) a further increase in Re introduces new equilibrium positions along the long faces due to a Bernoulli-type effect, explained in more detail below.

Chun and Ladd22 investigated the inertial migration of a collection of 32 spherical particles with a blockage ratio α ≈ 0.11 at a volume fraction of 1% in a square channel using a lattice-Boltzmann method at Re = 200–2000. Eight equilibrium positions, four at the center of the faces and four at the corners, were observed in the cross section at Re = 200, similar to observations in experiments. However, at Re = 1000, the four center equilibrium positions vanished, and some particles moved near the center of the channel at Re = 2000. The focusing behavior of particles in these high-Re simulations differs significantly from predictions for isolated or very dilute particles, and the differences can probably be attributed to multi-body hydrodynamic interactions, especially since the 32 particles organized into regularly spaced trains (also discussed by Di Carlo1). A numerical study using the immersed boundary method performed by Lashgari et al.28 showed that spherical particles with α = 0.2 in a channel with a square cross section focus at face-centered equilibrium positions that move closer to the wall with increasing Re until Re ≥ 600, above which they reverse, moving closer to the channel center. In addition, they showed that in plane Poiseuille flow, the particle reaches an equilibrium position (zeq/(H/2) ≈ 0.51), which is closer to the channel center than that in the rectangular channel with AR = 2 (zeq/(H/2) = 0.55) and in the square channel (zeq/(H/2) = 0.6) at Re = 200.

Di Carlo et al.23 performed finite element simulations as well as experiments on inertial migration of spherical particles with blockage ratios α > 0.2 in square channels. Their results for the lateral force at Re = 80 were better matched by empirical scaling formulae than by formula from the matched asymptotic expansion, valid for small α.18,20,31,33 For blockage ratios between 0.2 and 0.4, Di Carlo et al.23 found that FL=CLρfUm2a3/H for particles near the channel centerline (z/H < 0.5, where z is the distance from the center line to the center of the particle) and that FL=CLρfUm2a6/H4 for particles in the near-wall region (z/H > 0.5), where CL depends on the position. The success of these empirical scaling laws over a limited range of α was explained by Hood et al.,27 who carried out the regular perturbation analysis of Ho and Leal6 to one extra order in particle radius, yielding a term of order a5, which is one order higher than in the Ho-Leal analysis,6 who showed that the lowest order term is order a4. The analysis of Hood et al.,27 which was applied to both slits and square channels, shows that the scaling of Di Carlo et al.23 will fail as α drops below 0.1, and the Ho-Leal scaling is recovered. Also based on finite element simulations, corroborated by experiments for a channel with a square cross section, Liu et al.39 proposed a correction to the Ho-Leal lift formula consisting of a sum of four lift terms that depend on the local characteristics of the flow as well as the slip velocity: the shear-gradient-induced lift, the wall-induced lift, the slip-shear lift, and a correction of the shear-gradient-induced lift due to the finite particle size relative to curvature of the Poiseuille flow profile. Values of the coefficients of these terms were obtained by fits to the finite element simulations for a blockage ratio of α = 0.3. The findings of Di Carlo et al.23 imply that these four terms should be sensitive to α. Liu et al.39 then applied these formulae to predict successfully the trajectory of a particle in a serpentine channel. Such an approach may one day allow the computational design of complex channel geometries for separations of particles, droplets, and cells.

We here show the dependencies of the Schmidt number Sc and ratios of the Mach number to the Reynolds number Ma/Re on DPD parameters aff, n, rc, and γ, where we use the formulae in Table V, as well as Eq. (8) for the speed of sound and set m=rcC=kBT=1 as discussed in the main text. The four remaining quantities, and exponent s, determine the physically relevant dimensionless groups: the Schmidt number Sc, the ratio of the Mach to Reynolds number Ma/Re, and the ratio of the Weissenberg number to Reynolds number Wi/Re, which are dependent only on material properties and the gap (which is fixed) and not on the velocity. The formula for the viscosity is shown in the text, which we use to compute Sc, and Ma/Re, to be suited only for small enough aff. In each plot, all variables other than the one plotted on the x axis are held fixed at their reference values, which is unity, except for the reference value of aff = 25 and γ = 4.5 and n = 4. Note that according to Figs. 7 and 8 increased aff and n can decrease Ma/Re, but increased n leads to an increase in the Schmidt number and increased aff produces shear thinning, both of which lead to unphysical results.

FIG. 7.

Dependence of the Schmidt number (Sc) on (a) number density n, (b) cutoff radius for the dissipative and random forces rc, and (c) friction coefficient γ for the DPD simulations with dissipation exponent s = 2.0 and 0.5. The left y label is for s = 2.0 (the red solid line), while the right y label is for s = 0.5 (the blue dashed line).

FIG. 7.

Dependence of the Schmidt number (Sc) on (a) number density n, (b) cutoff radius for the dissipative and random forces rc, and (c) friction coefficient γ for the DPD simulations with dissipation exponent s = 2.0 and 0.5. The left y label is for s = 2.0 (the red solid line), while the right y label is for s = 0.5 (the blue dashed line).

Close modal
FIG. 8.

Dependence of the ratio of the Mach number to the Reynolds number Ma/Re on (a) repulsion coefficient aff, (b) number density n, (c) cutoff radius for the dissipative and random forces rc, and (d) friction coefficient γ for the DPD simulations. The left y label is for s = 2.0 (the red solid line), while the right y label is for s = 0.5 (the blue dashed line).

FIG. 8.

Dependence of the ratio of the Mach number to the Reynolds number Ma/Re on (a) repulsion coefficient aff, (b) number density n, (c) cutoff radius for the dissipative and random forces rc, and (d) friction coefficient γ for the DPD simulations. The left y label is for s = 2.0 (the red solid line), while the right y label is for s = 0.5 (the blue dashed line).

Close modal
2.
P.
Sajeesh
and
A. K.
Sen
,
Microfluid. Nanofluid.
17
,
1
(
2014
).
3.
H.
Amini
,
W.
Lee
, and
D.
Di Carlo
,
Lab Chip
14
,
2739
(
2014
).
4.
J. M.
Martel
and
M.
Toner
,
Annu. Rev. Biomed. Eng.
16
,
371
(
2014
).
5.
J.
Zhang
,
S.
Yan
,
D.
Yuan
,
G.
Alici
,
N. T.
Nguyen
,
M. E.
Warkiani
, and
W.
Li
,
Lab Chip
16
,
10
(
2016
).
6.
B. P.
Ho
and
L. G.
Leal
,
J. Fluid Mech.
65
,
365
(
1974
).
7.
S. C.
Hur
,
N. K.
Henderson-MacLennan
,
E. R.
McCabe
, and
D.
Di Carlo
,
Lab Chip
11
,
912
(
2011
).
8.
C. A.
Stan
,
L.
Guglielmini
,
A. K.
Ellerbee
,
D.
Caviezel
,
H. A.
Stone
, and
G. M.
Whitesides
,
Phys. Rev. E
84
,
036302
(
2011
).
9.
X.
Chen
,
C.
Xue
,
L.
Zhang
,
G.
Hu
,
X.
Jiang
, and
J.
Sun
,
Phys. Fluids
26
,
112003
(
2014
).
10.
S. K.
Doddi
and
P.
Bagchi
,
Int. J. Multiphase Flow
34
,
966
(
2008
).
11.
R. L.
Marson
,
Y.
Huang
,
M.
Huang
,
T.
Fu
, and
R. G.
Larson
,
Soft Matter
14
,
2267
(
2018
).
12.
D.
Di Carlo
,
D.
Irimia
,
R. G.
Tompkins
, and
M.
Toner
,
Proc. Natl. Acad. Sci. U. S. A.
104
,
18892
(
2007
).
13.
A. J.
Mach
,
O. B.
Adeyiga
, and
D.
Di Carlo
,
Lab Chip
13
,
1011
(
2013
).
14.
G.
Segré
and
A.
Silberberg
,
Nature
189
,
209
(
1961
).
15.
G.
Segré
and
A.
Silberberg
,
J. Fluid Mech.
14
,
115
(
1962
).
16.
G.
Segré
and
A.
Silberberg
,
J. Fluid Mech.
14
,
136
(
1962
).
17.
R. G.
Cox
and
H.
Brenner
,
Chem. Eng. Sci.
23
,
147
(
1968
).
18.
J. P.
Matas
,
J. F.
Morris
, and
É.
Guazzelli
,
J. Fluid Mech.
515
,
171
(
2004
).
19.
X.
Shao
,
Z.
Yu
, and
B.
Sun
,
Phys. Fluids
20
,
103307
(
2008
).
20.
J.
Matas
,
J.
Morris
, and
É.
Guazzelli
,
J. Fluid Mech.
621
,
59
(
2009
).
21.
Y. S.
Choi
and
S. J.
Lee
,
Microfluid. Nanofluid.
9
,
819
(
2010
).
22.
B.
Chun
and
A. J.
Ladd
,
Phys. Fluids
18
,
031704
(
2006
).
23.
D.
Di Carlo
,
J. F.
Edd
,
K. J.
Humphry
,
H. A.
Stone
, and
M.
Toner
,
Phys. Rev. Lett.
102
,
094503
(
2009
).
24.
A. A. S.
Bhagat
,
S. S.
Kuntaegowdanahalli
, and
I.
Papautsky
,
Microfluid. Nanofluid.
7
,
217
(
2009
).
25.
Y. S.
Choi
,
K. W.
Seo
, and
S. J.
Lee
,
Lab Chip
11
,
460
(
2011
).
26.
K.
Miura
,
T.
Itano
, and
M.
Sugihara-Seki
,
J. Fluid Mech.
749
,
320
(
2014
).
27.
K.
Hood
,
S.
Lee
, and
M.
Roper
,
J. Fluid Mech.
765
,
452
(
2015
).
28.
I.
Lashgari
,
M. N.
Ardekani
,
I.
Banerjee
,
A.
Russom
, and
L.
Brandt
,
J. Fluid Mech.
819
,
540
(
2017
).
29.
S. I.
Rubinow
and
J. B.
Keller
,
J. Fluid Mech.
11
,
447
(
1961
).
30.
P. G.
Saffman
,
J. Fluid Mech.
22
,
385
(
1965
).
31.
J. A.
Schonberg
and
E. J.
Hinch
,
J. Fluid Mech.
203
,
517
(
1989
).
32.
J.
Feng
,
H. H.
Hu
, and
D. D.
Joseph
,
J. Fluid Mech.
277
,
271
(
1994
).
33.
E. S.
Asmolov
,
J. Fluid Mech.
381
,
63
(
1999
).
34.
A. A. S.
Bhagat
,
S. S.
Kuntaegowdanahalli
, and
I.
Papautsky
,
Phys. Fluids
20
,
101702
(
2008
).
35.
M.
Masaeli
,
E.
Sollier
,
H.
Amini
,
W.
Mao
,
K.
Camacho
,
N.
Doshi
,
S.
Mitragotri
,
A.
Alexeev
, and
D.
Di Carlo
,
Phys. Rev. X
2
,
031017
(
2012
).
36.
A. T.
Ciftlik
,
M.
Ettori
, and
M. A.
Gijs
,
Small
9
,
2764
(
2013
).
37.
J.
Zhou
and
I.
Papautsky
,
Lab Chip
13
,
1121
(
2013
).
38.
C.
Liu
,
G.
Hu
,
X.
Jiang
, and
J.
Sun
,
Lab Chip
15
,
1168
(
2015
).
39.
C.
Liu
,
C.
Xue
,
J.
Sun
, and
G.
Hu
,
Lab Chip
16
,
884
(
2016
).
40.
R. V.
Repetti
and
E. F.
Leonard
,
Nature
203
,
1346
(
1964
).
41.
R. C.
Jeffrey
and
J. R. A.
Pearson
,
J. Fluid Mech.
22
,
721
(
1965
).
42.
P.
Vasseur
and
R. G.
Cox
,
J. Fluid Mech.
78
,
385
(
1976
).
43.
A. J.
Hogg
,
J. Fluid Mech.
272
,
285
(
1994
).
44.
X.
Fan
,
N.
Phan-Thien
,
S.
Chen
,
X.
Wu
, and
T.
Yong Ng
,
Phys. Fluids
18
,
063102
(
2006
).
45.
J. A.
Millan
,
W.
Jiang
,
M.
Laradji
, and
Y.
Wang
,
J. Chem. Phys.
126
,
124905
(
2007
).
46.
N.
Phan-Thien
,
N.
Mai-Duy
, and
B. C.
Khoo
,
J. Rheol.
58
,
839
(
2014
).
47.
D.
Pan
,
Y.
Lin
,
L.
Zhang
, and
X.
Shao
,
J. Hydrodyn., Ser. B
28
,
702
(
2016
).
48.
A. L.
Blumers
,
Y. H.
Tang
,
Z.
Li
,
X.
Li
, and
G. E.
Karniadakis
,
Comput. Phys. Commun.
217
,
171
(
2017
).
49.
J. M.
Kim
and
R. J.
Phillips
,
Chem. Eng. Sci.
59
,
4155
(
2004
).
50.
T.
Zhao
,
X.
Wang
,
L.
Jiang
, and
R. G.
Larson
,
J. Chem. Phys.
139
,
084109
(
2013
).
51.
P. J.
Hoogerbrugge
and
J. M.
Koelman
,
Europhys. Lett.
19
,
155
(
1992
).
52.
P.
Espanol
and
P.
Warren
,
Europhys. Lett.
30
,
191
(
1995
).
53.
R. D.
Groot
and
P. B.
Warren
,
J. Chem. Phys.
107
,
4423
(
1997
).
54.
D. C.
Visser
,
H. C.
Hoefsloot
, and
P. D.
Iedema
,
J. Comput. Phys.
214
,
491
(
2006
).
55.
C. A.
Marsh
and
J. M.
Yeomans
,
Europhys. Lett.
37
,
511
(
1997
).
56.
J. A.
Anderson
,
C. D.
Lorenz
, and
A.
Travesset
,
J. Comput. Phys.
227
,
5342
(
2008
).
57.
J.
Glaser
,
T. D.
Nguyen
,
J. A.
Anderson
,
P.
Lui
,
F.
Spiga
,
J. A.
Millan
,
D. C.
Morse
, and
S. C.
Glotzer
,
Comput. Phys. Commun.
192
,
97
(
2015
).
58.
R. D.
Groot
,
Novel Methods in Soft Matter Simulations
(
Springer
,
Berlin, Heidelberg
,
2004
), p.
5
.
59.
E. A.
Peters
,
Europhys. Lett.
66
,
311
(
2004
).
60.
V.
Symeonidis
,
G. E.
Karniadakis
, and
B.
Caswell
,
J. Chem. Phys.
125
,
184902
(
2006
).
61.
N.
Mai-Duy
,
D.
Pan
,
N.
Phan-Thien
, and
B. C.
Khoo
,
J. Rheol.
57
,
585
(
2013
).
62.
N.
Phan-Thien
,
Understanding Viscoelasticity: An Introduction to Rheology
, Graduate Texts in Physics (
Springer
,
2012
).
63.
C. A.
Marsh
,
G.
Backx
, and
M. H.
Ernst
,
Europhys. Lett.
38
,
411
(
1997
).
64.
D.
Pan
,
N.
Phan-Thien
,
N.
Mai-Duy
, and
B. C.
Khoo
,
J. Comput. Phys.
242
,
196
(
2013
).
65.
N.
Mai-Duy
,
N.
Phan-Thien
, and
T.
Tran-Cong
,
Comput. Phys. Commun.
221
,
290
(
2017
).
66.
D.
Azarnykh
,
S.
Litvinov
,
X.
Bian
, and
N. A.
Adams
,
Phys. Rev. E
93
,
013302
(
2016
).
67.
T.
Zhao
,
X.
Wang
,
L.
Jiang
, and
R. G.
Larson
,
J. Rheol.
58
,
1039
(
2014
).