Inertia- and deformation-driven migration of a soft particle in confined shear and Poiseuille flow

Abstract

The cross-stream migration of dilute soft particle suspensions under simple shear flow and Poiseuille flow between two parallel plates is investigated with the lattice Boltzmann-immersed boundary method. The competition between particle elastic contraction, fluid shear forces, and fluid inertial stress drives particle migration to a particular steady state position. With a small shear rate, the migration velocity of hard and soft particles is captured by a first order analysis of the Navier–Stokes equation. With a moderate shear flow, the qualitative dependence of the migration velocity and the particle position on the shear rate for both hard and soft particles deviates from the predictions. In a moderate Reynolds number (Re) shear flow, the observed hard sphere migration velocity has a weaker dependence on the Re than predicted and also a higher order dependence on the particle distance from the channel center. For soft spheres, a migration-free zone is observed near the center at a moderate Re and Weissenberg number (Wi). In Poiseuille flow, the soft particle migrates away from the wall to an off-center position dependent on the particle deformation and inertia, in contrast to hard sphere migration where the steady state position is independent of the shear rate.

---

1 Introduction

The non-Newtonian rheological properties of soft colloidal suspensions are mainly caused by the processes of evolving particle structures under flow, such as aggregation, breakup, and particle deformation. The properties of these suspensions strongly depend on the coupling between the inter- and intra-particle structure and the fluid forces. For example, the strong shear-thinning of blood is due to red cell rouleaux break-up, red cell deformation and lateral migration away from the vessel walls.¹ Similarly, the shear-thinning and -thickening of polymer and colloidal suspensions strongly depend on the evolution of aggregation and molecular structural changes under the flow shear force.^2,3 Understanding the structural evolution of these suspensions is critical for understanding micro-circulation and other flow-related phenomena.^4–11

The complex dynamics of soft particle flow is known to depend on the particle elasticity, shape, concentration, and inter-particle repulsion/attraction. In applications such as micro-capsule drug delivery and microfluidic droplet reactors, controlling the particle shape and the dynamics is critical for ensuring delivery to the intended target via microflow.¹² In blood micro-circulation, it is known that blood viscosity decreases when the flow shear rate increases in the micro-capillaries.^1,13 This ‘shear-thinning’ effect depends strongly on the diameter of the blood vessel relative to the cell size and the flow shear rate. In addition, the dynamics of the spherical soft white blood cells (WBC) is known to be coupled to the flow of the biconcave discoid elastic red blood cells (RBC) and vice versa.^{4,6,11,14–19} To disentangle the complex interactions between the different components in blood flow, models for deformable particles (DP) and fluid flow are needed to understand how DP–fluid, DP–vessel and DP–DP hydrodynamic coupling affect the fluid and particle motion.

Cross-stream lateral migration of soft deformable particles can be attributed to the asymmetric hydrodynamic field produced by an elastic particle near a wall, and the unbalanced fluid stress at the particle surface. For a deformable particle in a suspension, the competition between the particle shape-restoring elasticity, the elastic relaxation time τ_relax, and the flow shear force, is one factor that determines the migration velocity. The deformation-driven migration force can thus be characterized by the capillary number Ca = [small gamma, Greek, dot above]ηR/G, where η is the fluid viscosity, R is the particle radius, G is the surface elastic modulus, and [small gamma, Greek, dot above] is the flow shear rate. Equivalently, the particle Weissenberg number Wi = τ_relax[small gamma, Greek, dot above] also describes the particle deformation and Wi is used in this study. For Wi > 1, the shear force overcomes the particle elastic contraction and strongly deforms the particle. The flow field around a stretched, elastically contracting soft particle near the wall results in a net migration force that pushes it away from the wall.^20–23

In addition, particle inertia in flow also leads to migration for finite particle Reynolds number Re = ρ[small gamma, Greek, dot above]R²/η ≥ 1, where ρ is the fluid density.^24–27 Re is characterized by the fluid velocity difference across the particle δU = [small gamma, Greek, dot above]R. It has been shown that the difference between the fluid velocity at the particle surface and the particle velocity results in a particle surface stress leading to cross-stream lateral migration.²⁸ In simple shear flow, the migration force pushes the particle towards the steady state position at the channel center. In Poiseuille flow, the steady state position is off-center due to the curvature of the flow velocity profile.

To understand the role of the deformation- and inertia-driven migration of soft particle suspensions, one needs to account for not only how the particle motion perturbs the fluid, but also the particle deformation due to fluid shear and hydrodynamic perturbations. Boundary integral formulations have extensively studied droplet dynamics in Stokes flow.^29–35 However, accurately modeling the interplay between particle deformation and the fluid field can require a very fine surface mesh and is very computationally expensive. Particle deformation and dynamics are coupled with the hydrodynamic field and need to be resolved simultaneously. Recent advances in multi-scale approaches such as the immersed boundary method (IB),^36–38 the multi-particle collision method (MPC),^39–41 dissipative particle dynamics (DPD)^42,43 and the lattice Boltzmann (LB) approach^44–48 have significantly improved the efficiency of complex fluid flow modeling based on individual particles. In this study, we employ a combination of LB and IB to study the lateral migration of a particle in simple shear and Poiseuille flow.

2 Method

We developed a coarse-grained model combining Lagrangian particle dynamics⁴⁹ and Eulerian LB fluid dynamics through an immersed boundary method (IB).^{36,46,47,50,51} The method has been shown to capture the physical characteristics of particle deformation in shear flow, the inter-particle hydrodynamic forces, and the near-wall hydrodynamic field around a shear-deformed particle.²⁰ In order to understand the effects of shear-induced near-wall deformation- and inertia-driven particle migration, we examine particle migration in the simple scenario of a single particle in shear flow and Poiseuille flow between two parallel plates with height H = 16R. An ensemble of one hundred particle trajectories is performed until the particle distance from the wall reaches steady state, from small to moderate Wi and Re. Particle and fluid fluctuations are not included in order to directly determine the particle migration velocities.

LB is carried out on a three-dimensional cubic lattice with 19 discrete velocities (the D3Q19 model). The lattice spacing is Δx, the kinematic viscosity is ν = η/ρ = 1/6[Δx²/Δτ], and Δτ is the LB time step. Details of the LB method can be found in recent reviews.^45,52 LB models a compressible fluid, where the fluid density can vary significantly at high velocities if the Mach number Ma = v_fluid/v_sound → 1. The error due to fluid compressibility is of the order (Ma²). In this study, the fastest fluid velocities are kept to Ma < 0.4. The systematic errors are considered by comparison with theoretical predictions.

The coarse-grained soft particle is constructed by accounting for the membrane elasticity with a two-dimensional closed network of beads and springs. The beads interact with each other via a purely repulsive Weeks–Chandler–Andersen (WCA) potential, given by

Neighboring beads are connected with a finitely extensible non-linear elastic (FENE) spring, given by

which is harmonic at a small separation but restricts the separation to l₀. The parameters of the potentials are chosen to be σ = Δx, l₀ = 3Δx, ε = k_BT. The spring constants kσ² = 0.017, 0.17, and 2.83Δε are chosen for particles with different stiffness. Δε is the simulation energy unit.

Bead positions on the spherical shell are acquired by triangle tessellation, as shown in Fig. 1. The radius of the sphere is R = 3.2Δx. A bending potential is imposed between all neighboring triangular faces

U_bend = k_bend(1 − cos θ), (3)

where θ is the angle between the surface normals of two neighboring faces, and values of k_bend = 0.056, 0.56, 2.83Δε are used for kσ² = 0.017, 0.17, 2.83, respectively. These parameters are chosen such that the ratio between the bending and the elastic moduli, characterized by the Föppl–von Kármán number Fv = GR²/k_bend, is less than 100 to ensure that the buckling modes of the coarse-grained particle are dampened.^53,54

Fig. 1 Illustration of (a) a hard sphere in simple shear flow at a position y from the wall, and its relative distance from the center is d = (1 − 2y/H). (b) A deformed particle with long axis L and short axis B in Poiseuille flow between 2 parallel plates.

The triangular bending force is exerted on the four beads with positions r₁, r₂, r₃, r₄ that constitute two adjacent triangles r₁ − r₂ − r₃ and r₂ − r₁ − r₄. The bending forces on each bead are given by

where n₁ = (r₃ − r₁) × (r₄ − r₁), n₂ = (r₄ − r₂) × (r₃ − r₂), and a = 2 × 1. The sign of θ is given by sin θ = (n̂2 × n̂1)ê, where e = r₄ − r₃, then a = sin θê. Eqn (4) has the expected directionality for a force that only affects the angle but does not affect the translation, rotation, or shape of the two triangles.

The particle volume and total surface area are conserved by applying an isotropic pressure P_V(t) and a surface dilatation tension k_A on the particle, given by

P = F_V(t)/A = [k_V(V(t) − V₀)/V₀]/A (5a)

F_A(t) = k_A(A(t) − A₀)/A₀ (5b)

where V(t) and V₀ are the instantaneous and rest particle volumes, and A(t) and A₀ are the instantaneous and rest surface areas. The constraint forces F_V and F_A are distributed equally on all the beads that constitute the particle surface. For the softest particle examined here, it was found that with k_V = 0.56Δε/Δx and k_A = 0.56Δε/Δx, the particle volume and area can vary up to 10% at the highest shear rate examined. The total conservative forces acting on each bead is given by F^C = −∇U = F_WCA + F_FENE + F_bend + F_V + F_A. All beads have mass m_b = 36m₀, where m₀ is the simulation unit mass. The velocity Verlet method is employed to project the bead trajectories with the time step δt = 0.1Δτ.

The shape of the particle is characterized by the Taylor deformation parameter D = (L − B)/(L + B), where L and B are the longest and shortest particle axes during deformation, as shown in Fig. 1b. The particle elasticity is characterized by the surface elastic modulus G = stress/strain = δf_x/δR_x, in units of Δε/Δx², where f_x is an external stretching force in the x-direction and δR_x is the particle deformation along the direction of the applied force. The orthogonal elastic modulus δf_x/δR_z is also measured. The particle elasticity is also characterized by its deformation relaxation time τ_relax, measured from the relaxation of the particle stretch after shear deformation by fitting [L(t) − L_eq]² to exp(−2t/τ_relax). Fig. 2a shows that G ∼ k^0.7±0.05. The inset shows that Fv < 100 for the regime investigated in this study and particle buckling is not significant. G is also found to very weakly depend on k_bend, k_V and k_A. The elastic relaxation time τ_relax ∼ G^−1.1±0.1 is also observed, which gives an approximate relation between Ca and Wi, with Wi ≈ 15Ca.

Fig. 2 (a) Particle long (symbols) and short (symbols with dashed line) axis stretch elasticity as functions of the spring constant for k_bend = 0.56 (circles) and k_bend = 2.83 (squares). The solid symbols indicate the chosen parameters for the hard and soft spheres examined here. The dotted line shows the power law fit with G ∼ k^0.7. The inset shows Fv as a function of k. (b) Particle deformation as a function of Wi and Ca (inset).

The particle deformation D is known to be linearly dependent on Ca under small shear for droplets and elastic capsules.^22,29 For the elastic particles considered in this study, D is observed to be linearly proportional to Wi only for Wi < 1, and D reaches a plateau under strong shear as the particle stretches in the non-linear elasticity regime, as shown in Fig. 2b. The inset shows that D ≈ (25/4)Ca, which agrees with the theoretical prediction.²⁹ Under simple shear flow, the particle deformation is constant throughout the channel except near the wall.

Hydrodynamic interactions are captured with the exchange and propagation of frictional momentum between the fluid and the DP.^45,55 Each bead on the DP experiences a friction force

F_f = −ζ[u_b(r_b) − u_f(r_b)] (6)

where u_b(r_b) and u_f(r_b) are the bead and fluid velocity at the bead position r_b, respectively, and ζ = 3πησ = 56.4m₀/Δτ is the bead friction coefficient. The fluid velocity at position r_b is determined by a trilinear interpolation of the fluid velocity on the neighboring lattice points nn that enclose r_b, u_f(r_b) = ∑_r∈nnw_ru(r), where the weights w_r are the coefficients of the normalized linear Lagrange interpolation polynomial. The frictional momentum density Δj = −F_fδt/Δx³ is transferred to the fluid.

The beads are repelled from the walls with a repulsive potential when the beads are within 1Δx of the walls. The particle elasticity and fluid velocity are independently varied in order to investigate the contributions of particle inertia and particle deformation to particle migration. Three types of soft particles, with the two-dimensional elasticity G = 0.4, 4.4, and 44.0 Δε/Δx² are chosen to identify the hydrodynamically-induced migration fluxes due to particle inertia and deformation. With the simulation energy (Δε) matched to the thermal energy at room temperature and Δx = 1 μm, the elasticities correspond to 0.29, 2.9, and 29 μN m⁻¹. For Δx = 100 μm, the elasticities correspond to 2.9, 29, and 290 mN m⁻¹. For comparison, the elasticity of a red blood cell (diameter ≈ 6–8 μm) is ≈6 μN m⁻¹ and the surface tension of water in air is ≈72 mN m⁻¹. The elastic properties of soft particles with a size range between μm and mm can be captured with the coarse-grained particle model.

With the assumptions involved in LB and IB, systematic errors can arise due to fluid compressibility and the momentum transfer between the particle and the fluid. With the immersed boundary method, the no-slip boundary condition is imposed on the particle surface, and the friction momentum is calculated based on a first-order Lagrange interpolation of the fluid velocity at the bead position from the lattice fluid velocity. Thus, the accuracy of momentum coupling is O(Δx²). Furthermore, the fluid velocity from LB is of O(Δx²) accuracy. The systematic error may be estimated by examining the velocity field generated by a moving particle with velocity U_p in a quiescent fluid inside a periodic domain. Comparison of the LB–IB solution is made with the results using the lattice Boltzmann method with the fluid bounce-back condition (LB–BB) at the fluid–particle interface. The accuracy of LB–BB has been well characterized for colloidal hydrodynamics, and it was found to accurately capture the drag force acting on a colloid in flow up to a particle Reynolds number Re_p = ρU_pR/η ≈ 100.^56,57 The fluid velocity field along the direction orthogonal to the particle velocity is calculated using LB–IB and LB–BB. Fig. 3 shows that the results quantitatively agree for the range of particle velocities relevant to this study. This shows that the immersed boundary method captures the fluid velocity field with the same order of accuracy as the bounce-back condition. Furthermore, Fig. 3 shows that the fluid field agrees quantitatively with the prediction by Hasimoto⁵⁸ for a point particle in Stokes flow in a periodic domain for Re_p = 0.29. At higher Re_p, the reduced fluid velocity corresponds to a increase in the fluid boundary layer around the spherical particle. Fluid density variation inside and outside the particles at the highest flow rates studied were found to be less than 0.1 percent.

Fig. 3 Fluid velocity field along the velocity-normal at distance (x − R)/R from a moving particle with velocity U_p = U_pŷ. U_p is chosen with the Mach numbers Ma = 0.026 (circles), 0.26 (triangles), 0.52 (squares), corresponding to Re_p = 0.29, 2.9, 5.8. The results from LB–IB are shown by the empty symbols, and the results from LB–BB are shown by the filled symbols. The prediction for a periodic point force is shown as the dashed line.

3 Migration of hard and soft particles

We first examine the lateral migration of hard spheres. Shear-induced inertial migration of hard spheres has been well-studied experimentally and theoretically.^22,25,28,59 For small but finite Re, the unbalanced surface stress due to particle inertia leads to particle migration away from a wall boundary. In contrast, a hard sphere in Stokes flow (Re = 0) does not migrate. The inertia-driven migration effect has been exploited to separate particles of different sizes and shapes.^{22,24,25,59,60} First order theoretical analysis of the Navier–Stokes equation in the limit of Re < (R/H)² shows that the migration velocity v_lift is proportional to Re,²⁸ given bywhere v_m = v_wall, α = 0 for simple shear flow and v_m = v_max, α = −4 for Poiseuille flow. G₁(d) and G₂(d) are numerically integrated functions from the particle surface stress.²⁸ In simple shear flow, the lift velocity is expected to be linearly dependent on the particle distance from the channel center (d = 1 − 2y/H) as the particle approaches the center, with the steady state position at the center where v_lift = 0. In Poiseuille flow, the lift velocity is expected to quadratically vary with d, with zero crossings and steady state positions at d = 0.6 and −0.6.

LB–IB is employed to examine particle migration under simple shear flow and Poiseuille flow. First, the predictions for hard spheres are tested with a very stiff particle with G = 44Δε/Δx², moderate Re and negligible Wi. Under simple shear flow, Fig. 4a shows that the migration velocity is linearly dependent on d for d < 0.6 at lower Re, and that the particle steady state position is at the center. For Re < 0.1, the observed v_lift from LB–IB is in quantitative agreement with eqn (7). Quantitative agreement is also found with the results of LB–BB up to Re = 0.32. As Re increases, it is observed that the relative migration velocity does not increase linearly with Re. Furthermore, v_lift is found to become quadratically dependent on d with the steady state position at the center, which is captured by both LB–IB and LB–BB. These differences to eqn (7) may be due to higher order inertial effects since the first order analysis is strictly valid for Re < (R/H)². For comparison with experiments with 10 and 100 μm particles, the shear rates for Re = 1 are 10⁴ and 100 s⁻¹, respectively.

Fig. 4 Lift velocity of a hard sphere with G = 44Δε/Δx² in (a) simple shear for Re = 0.028 (squares), 0.07 (triangles), 0.14 (circles), 0.28 (*), 0.7 (+) and (b) pressure-driven flow for Re = 0.032 (squares), 0.16 (triangles), 0.32 (circles), 0.65 (*) from the composition of 100 initially randomly placed spheres. The corresponding solid lines are the results of LB–BB. The dashed lines show the predicted v_lift from eqn (7) ordered from bottom to top for low Re to higher Re.

Under pressure-driven flow, v_lift is found to have quadratic dependence on d for d < 0.6 for the entire range of Re examined. The LB–IB results agree quantitatively with the O(Re) analysis for Re = 0.032 and 0.16. It is observed that the particles near the central region will migrate towards the walls due to the negative v_lift, and the steady state position at d = 0.6 does not vary with Re even for Re much larger than (R/H)² as predicted. For the results of LB–BB, more velocity oscillations are found due to errors caused by the first order bounce-back condition.

In the other limiting regime where particle inertia is small and particle deformation is large, particle migration is driven by deformation-induced hydrodynamic forces near a wall. Under pressure-driven flow, DP deformation depends on the local shear force and varies linearly as the particle migrates towards the center at small shear rates, as shown in Fig. 5. At higher Wi, the particle deformation D plateaus near the wall, where the shear rate is the highest, due to the non-linear elastic spring. For Wi = 6.9, the highest shear rate examined, the particle stretches nearly 100%. Under simple shear flow, the local shear rate is constant throughout the channel and the steady state particle shape does not change except near the wall. Prior studies have measured the migration velocity.^22,23,38,61 Under simple shear flow, the migration velocity for droplets with equal inner and outer viscosity is found to be

Fig. 5 Deformation of soft particles with G = 0.4Δε/Δx² under pressure-driven flow for Re = 0.033 (squares), 0.065 (triangles), 0.165 (circles), 0.33 (*), 0.65 (+) and Wi = 0.35, 0.7, 1.7, 3.5, 6.9 ordered from bottom to top. The data set is from the composition of 50 initially randomly placed spheres. The data point spread is due to the initial particle shape change from sphere to ellipsoid under shear.

The prefactor depends on the ratio between the inner and outer fluid viscosity of the particle. To examine particle migration in the deformation-dominated regime, the migration trajectories of a softer particle with G = 0.4 are examined. The v_lift exhibits a quadratic dependence on d as shown in Fig. 6a, as expected from eqn (8). The lift velocity for Wi = 0.15 and 0.3 agrees quantitatively with eqn (8), but v_lift/v_wall does not increase proportionally with Wi as the shear rate further increases. This may be due to the non-linear dependence of particle deformation on shear force for large particle deformations, higher order Wi contributions, and the increased effect of inertia driven migration.

Fig. 6 Lift velocity of a deformable particle with G = 0.4Δε/Δx² in (a) simple shear Re = 0.014 (squares), 0.028 (triangles), 0.07 (circles), 0.14 (*), 0.28 (+) and Wi = 0.15, 0.3, 0.74, 1.5, and 3.0. The dashed lines are the predictions of eqn (8) ordered from bottom to top for low Re to higher Re. (b) In pressure-driven flow for Re = 0.033 (squares), 0.065 (triangles), 0.33 (*), 0.65 (+) and Wi = 0.35, 0.7, 1.7, 3.5, 6.9, from the composition of 50 initially randomly placed spheres. The dashed lines show eqn (10) multiplied by a factor of 2 for the same Wi, ordered from bottom to top. The dotted line shows eqn (9) for Wi = 0.35.

Under Poiseuille flow, Fig. 6b shows that the particle migration velocity has a non-linear dependence on d for Wi < 1.0, and v_lift decays to 0 as d approaches 0. First order analysis has found the lateral migration velocity to be²³

From Fig. 5, it is observed that D decays linearly as d approaches 0. This dependence combined with eqn (9) shows that v_lift should decay quadratically near the channel center.^22,23,38,61 However, Fig. 6b shows that the observed migration velocity increases faster than the quadratic dependence on d for d > 0.4. This cannot be simply explained by the non-linear elastic particle deformation near the wall. Indeed, the non-linear elastic deformation near the wall should lead to a weaker dependence on d instead of higher order dependence. The stronger measured migration velocity may be due to additional contributions to particle migration due to the particle inertia. However, simply adding the inertial migration contribution from eqn (7b) does not lead to better agreement with the measured migration velocity. Further analysis starting from the migration velocity of a droplet near a wall^22,61 and extending to a particle in Poiseuille flow between two parallel plates results in

As shown in Fig. 6b, the additional dependence on (1 + d²)/(1 − d²)² appears to account for the stronger dependence for a particle closer to the wall. The observed v_lift qualitatively agrees with eqn (10) for Wi < 1, with the observed v_lift larger than eqn (10) by a factor of 2.

At high Wi where the particle is strongly deformed near the wall, it is found that closer to the wall, the dependence of v_lift on d decreases. This is consistent with the smaller changes in the particle deformation at higher shear rates near the wall, and leading to a weaker dependence of v_lift on d. In addition, v_lift becomes negative as the particle approaches the channel center, showing a migration flux away from the center for particles near the center. This may be attributed to inertia-driven particle migration away from the center at the higher shear rates, where a negative v_lift is found near the channel center as shown for hard spheres. Furthermore, as the particle moves cross-stream towards the center, the projected cross-stream surface area (A_xz = ∑n_i(x̂ + ẑ)) decreases as the particle shape becomes less deformed. This leads to a gradient in the migration friction and a migration flux away from the center. The position at which v_lift = 0 also shifts towards the channel center as the shear rate increases, which is due to the change in the particle deformation for different Wi. This is in contrast to hard sphere migration for which the steady state position is independent of Re in Poiseuille flow.

We further examine particle migration under simple shear and Poiseuille flow chosen such that Re ∼ Wi < 1 for soft particles with G = 4.4Δε/Δx². In simple shear flow, the linear combination of eqn (7) and (8) is in quantitative agreement with the observed v_lift at low shear rates as shown in Fig. 7a. Quadratic dependence of v_lift on d is found for Wi > 0.3, where the particle is stretched by ≈10%. For Re = 0.7 and Wi = 0.74, v_lift is found to quadratically decrease to zero at y/H ≈ 0.3 and v_lift = 0 in the central region for d < 0.35. In this region, there is no observed net lift force about the particle and the particle migrates neither towards nor away from the center. Indeed, for Wi > 0.3, the migration-free region is found to expand from the channel center as the shear rate increases, as shown in Fig. 8. This is in stark contrast to the inertia-dominated or deformation-dominated lateral migration in simple shear flow, where the particle migrates to the channel center. It also differs from purely inertia- or deformation-driven migration results found in Poiseuille flow, where the particles migrate to a specific position in the channel. One possible explanation for the “migration-free” zone may be the rotating motion of the particle surface. At all Wi examined here, the soft particle appears to undergo a tank-treading motion in simple shear flow. It has been previously suggested that the tank-treading motion of an ellipsoidal particle also induces a migration flux.^62,63 For the higher Wi where the particle does not migrate, they appear to simply tank-tread at an inclined angle, without lateral migration. The tank-treading frequency is proportional to the shear rate, reaching as high as 10⁻²[1/Δτ]. The fluid perturbations due to the tank-treading motion at high frequency may overcome the inertia- and deformation-driven fluid stress gradient.

Fig. 7 Lift velocity of a deformable particle with G = 4.4Δε/Δx² in (a) simple shear for Re = 0.028 (squares), 0.07 (triangles), 0.14 (circles), 0.28 (*), 0.7 (+) and Wi = 0.03, 0.074, 0.15, 0.3, and 0.74. The dashed lines show the linear combinations of eqn (7) and (8), ordered from bottom to top. (b) Pressure-driven flow for Re = 0.033 (squares), 0.065 (triangles), 0.33 (circles), 0.65 (*), 0.98 (+) and Wi = 0.035, 0.17, 0.35, 0.70, 1.05. The figures are an overlaid composition of 100 trajectories of initially randomly placed spheres.

Fig. 8 Lift velocity of a deformable particle with G = 4.4Δε/Δx² in simple shear for Re = 0.28 (*), 0.42 (solid triangles), 0.56 (solid squares), and 0.74 (+) and Wi = 0.3, 0.44, 0.59, 0.74.

Under Poiseuille flow, soft sphere migration velocity is quadratically dependent on d, similar to hard sphere migration at small Wi and Re, as shown in Fig. 7b. The position at which v_lift crosses zero shifts towards the center as the shear rate increases in contrast to the inertia-driven migration of hard spheres. For Wi = 1.05, v_lift appears to decrease linearly near the wall for d > 0.5, and quadratically as the particle migrates closer to the channel center. These differences may be attributed to particle deformation due to the shear force, and also the change of the particle shape as the local fluid shear force decreases as the particle moves from the wall towards the center. Near the wall, the shear force is the strongest and the particle is deformed to an ellipsoid. As the particle approaches the center, it becomes more spherical. The change of particle shape introduces a deformation-dependent particle drift away from the channel center, which is more pronounced for softer particles.

4 Conclusion

For droplets, vesicles, and capsules, inertial- and deformation-driven migration away from the channel walls in microfluidic flow may be utilized to separate particles based on size and deformability. It is interesting to observe that particle migration velocity in Poiseuille flow is 3 to 4 times faster than in simple shear flow for hard and soft particles, which is due to the larger fluid velocity change across the particle. This study also found that for small Reynolds and Weissenberg numbers flow, the particle migration velocities are quantitatively captured by a first-order analysis of the Navier–Stokes equation in shear flow. For hard spheres, theoretical analysis also captures the migration velocity in Poiseuille flow. For very soft particles, the predictions of the deformation-driven particle migration velocity in pressure-driven flow show a weaker dependence on d than observed, but analysis based on the migration velocity near a surface appears to be qualitatively consistent with the observed migration velocity. For soft particles in shear flow, the linear combination of the O(Re) and O(Wi) predictions capture the particle migration velocities under conditions where the particle inertia and deformation are small but not negligible.

At moderate Re and Wi, it is found that the particle migration velocities for hard and soft particles are smaller than expected and do not increase linearly with Re and Wi. The qualitative dependence of the migration velocity on the particle distance from the channel center also changes as Re approaches 1 and as Wi approaches 1. Under simple shear flow, the migration velocity of a hard sphere has non-linear dependence on d. For soft particles, the migration velocity becomes zero in a central migration-free region. The steady state position of the particles remain at the channel center or in the central region. It is suspected that this may be due to a combination of (1) contributions from higher order inertia- and deformation-driven flux, (2) the fluid boundary layer around the particle increasing at higher Re, (3) fluid field perturbation due to the particle tank-treading motion and (4) systematic errors due to high fluid velocity. The cause will be investigated in future studies.

Under pressure-driven flow, the steady state position of a hard sphere is found to be weakly dependent on or independent of Re. For soft particles, the steady state position moves towards the center as the shear rate and particle deformation increases. This indicates that particles with different elasticities can be separated by their distance from the center-line. For hard spheres, different species with the same size will not be separated in the stream, but particles of different sizes and shapes may be separated. For soft spheres of the same size but different elasticities, increased particle deformation results in a steady state position closer to the center. With the parabolic flow profile, this means that the soft particles would travel at different velocities, which is being exploited for cell and particle flow-based fractionation in microfluidic devices.

Acknowledgements

Funding for this study comes from the Academia Sinica Career Development Award 100-CDA-M01 and the Republic of China National Science Council project 101-2112-M01-03-MY3. Support also came from NCTS Taiwan. I would also like to thank Dr D. Leighton, Dr H. C. Chang, Dr C. S. Hur, Dr M. Dao, and Dr Z. Peng for their help and support through many discussions.

References

1. R. Fåhraeus and T. Lindqvist, Am. J. Physiol., 1931, 96, 562.
2. J. Vermant and M. Solomon, J. Phys.: Condens. Matter, 2005, 17, R187.
3. E. Nazockdast and J. F. Morris, Soft Matter, 2012, 8, 4223.
4. S. Chien, S. Usami and R. Skalak, in Handbook of Physiology - The Cardiovascular System, American Physiological Society, 1984, vol. IV, p. 217.
5. A. S. Popel and P. C. Johnson, Annu. Rev. Fluid Mech., 2005, 37, 43–69.
6. H. H. Lipowsky, Microcirculation, 2005, 12, 5.
7. C. Pozrikidis, Modeling and Simulation of Capsules and Biological Cells, Chapman & Hall/CRC, London, 1st edn, 2003, vol. 2.
8. A. Kumar and M. D. Graham, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 2011, 84, 066316.
9. A. Kumar and M. D. Graham, Phys. Rev. Lett., 2012, 109, 108102.
10. H. Zhao, A. P. Spann and E. S. G. Shaqfeh, Phys. Fluids, 2011, 23, 121901.
11. H. Zhao, E. S. G. Shaqfeh and V. Narsimhan, Phys. Fluids, 2012, 24, 011902.
12. P. Garstecki, I. Gitlin, W. DiLuzio, G. M. Whitesides, E. Kumacheva and H. A. Stone, Appl. Phys. Lett., 2004, 85, 2649–2651.
13. H. L. Goldsmith, G. R. Cokelet and P. Gaehtgens, Am. J. Physiol., 1989, 257, H1005–H1015.
14. H. Goldsmith and R. Skalak, Annu. Rev. Fluid Mech., 1975, 7, 213.
15. J. J. Bishop, P. R. Nance, A. S. Popel, M. Intaglietta and P. C. Johnson, Am. J. Physiol., 2001, 280, H222–H236.
16. B. P. Helmke, M. Sugihara-Seki, R. Skalak and G. W. Schmid-Schoenbein, Biorheology, 1998, 35, 437.
17. L. L. Munn, R. J. Melder and R. K. Jain, Biophys. J., 1996, 71, 466.
18. A. Kumar and M. D. Graham, Soft Matter, 2012, 8, 10536.
19. J. B. Freund and M. Orescanin, J. Fluid Mech., 2011, 671, 466–490.
20. C. Hsu and Y.-L. Chen, J. Chem. Phys., 2010, 133, 034906.
21. J. R. Smart and D. T. Leighton, Phys. Fluids, 1991, 3, 21.
22. A. Karnis and S. G. Mason, J. Colloid Interface Sci., 1967, 24, 164.
23. P. C. H. Chan and L. G. Leal, J. Fluid Mech., 1979, 92, 131–170.
24. S. C. Hur, H. T. Tse and D. Di Carlo, Lab Chip, 2010, 10, 274.
25. G. Segre and A. Silberberg, Nature, 1961, 189, 209.
26. A. Karnis, S. G. Mason and H. L. Goldsmith, Nature, 1963, 200, 159.
27. E. J. Lim, T. J. Ober, J. F. Edd, G. H. McKinley and M. Toner, Lab Chip, 2012, 12, 2199.
28. B. P. Ho and L. G. Leal, J. Fluid Mech., 1974, 65, 365.
29. D. Barthès-Biesel and J. M. Rallison, J. Fluid Mech., 1981, 113, 251.
30. E. Lac, A. Morel and D. Barthès-Biesel, J. Fluid Mech., 2007, 573, 149.
31. C. Pozrikidis, Ann. Biomed. Eng., 2003, 31, 1194.
32. C. Pozrikidis, Phys. Fluids, 2005, 17, 031503.
33. C. Pozrikidis, Ann. Biomed. Eng., 2005, 33, 165–178.
34. A. Kumar and M. D. Graham, J. Comp. Physiol., 2012, 231, 6682.
35. C. Pozrikidis, Computational Hydrodynamics of Capsules and Biological Cells, Chapman & Hall/CRC, London, 1st edn, 2010.
36. C. Peskin, Acta Numerica, 2002, 1.
37. S. K. Doddi and P. Bagchi, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 2009, 79, 046318.
38. S. K. Doddi and P. Bagchi, Int. J. Multiphase Flow, 2008, 34, 966–986.
39. G. Gompper, T. Ihle, D. Kroll and R. Wrinkler, Adv. Polym. Sci., 2009, 221, 1.
40. H. Noguchi and G. Gompper, Phys. Rev. Lett., 2007, 98, 128103.
41. H. Noguchi and G. Gompper, Proc. Natl. Acad. Sci. U. S. A., 2005, 102, 14159–14164.
42. I. V. Pivkin and G. E. Karniadakis, Phys. Rev. Lett., 2008, 101, 118105.
43. D. A. Fedosov, W. Pan, B. Caswell, G. Gompper and G. E. Karniadakis, Proc. Natl. Acad. Sci. U. S. A., 2011, 108, 11772–11777.
44. J. F. Zhang, P. C. Johnson and A. S. Popel, Phys. Biol., 2007, 4, 285.
45. B. Dünweg and A. Ladd, Adv. Polym. Sci., 2008, 221, 89.
46. S. Succi, The Lattice Boltzmann Equation for Fluid Dynamics and Beyond, Oxford University Press, Oxford, 1st edn, 2001.
47. D. Wolf-Gladrow, Lattice-Gas Cellular Automata and Lattice Boltzmann Models: An Introduction, Springer, 1st edn, 2000.
48. J. R. Clausen and C. K. Aidun, Phys. Fluids, 2010, 22, 123302.
49. M. Allen and D. Tildesley, Computer Simulation of Liquids, Clarendon Press, 1987.
50. P. Ahlrichs and B. Dünweg, Int. J. Mod. Phys. C, 1998, 9, 1429.
51. C. K. Aidun and J. R. Clausen, Annu. Rev. Fluid Mech., 2010, 42, 439.
52. R. Adhikari, K. Stratford, M. Cates and A. Wagner, Europhys. Lett., 2005, 71, 473.
53. M. Widom, J. Lidmar and D. R. Nelson, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 2007, 76, 031911.
54. D. R. Nelson and L. Radzihovsky, Phys. Rev. A, 1992, 46, 7474–7479.
55. P. Ahlrichs and B. Dünweg, J. Chem. Phys., 1999, 111, 8225.
56. A. Ladd, J. Fluid Mech., 1994, 271, 285.
57. A. Ladd, J. Fluid Mech., 1994, 271, 311.
58. H. Hasimoto, J. Fluid Mech., 1959, 5, 317.
59. A. Karnis, H. L. Goldsmith and S. G. Mason, Can. J. Chem. Eng., 1966, 44, 181.
60. S. C. Hur, A. J. Mach and D. Di Carlo, Biomicrofluidics, 2011, 5, 22206.
61. C. E. Chaffey, H. Brenner and S. G. Mason, Rheol. Acta, 1965, 4, 56.
62. P. Olla, J. Phys. II, 1997, 7, 1533.
63. M. Abkarian and A. Viallat, Biophys. J., 2005, 89, 1055.

---