Inertial migration of oblate spheroids in a plane channel

We study the inertial migration of neutrally buoyant oblate spheroids in a plane channel at moderate Reynolds numbers using lattice Boltzmann simulations. Spheroids reorient to perform a log-rolling motion with their minor axis in the vorticity direction. We demonstrate that, for moderate aspect ratios, the equilibrium positions relative to the channel walls for such a stable motion depend only on the equatorial radius $a$ of the spheroid, and the inertial lift force on the spheroid is proportional to $a^3b$, where $b$ is the polar radius. Therefore, the lift force on the spheroid can be expressed in terms of that for the sphere of the same $a$.


I. INTRODUCTION
The inertial migration of particles in microchannels has an important impact on the efficiency of microfluidic focusing and separation devices 1,2 . At finite Reynolds numbers an inertial lift force causes suspended particles to migrate across flow streamlines to equilibrium positions at some distance from the channel walls where the force is effectively zero. In microfluidic applications this force is often balanced by some other force, which allows to control particle positions. So far, studies mainly focused on spherical particles. In their pioneering experiments, Segrè and Silberberg found that small spheres focus to a narrow annulus at a radial position of about 0.6 of a pipe radius 3 . Later, theoretical [4][5][6][7] and numerical 8 papers proposed several scaling formulas and approximations for the lift force in a channel flow which are widely applied in practical applications. While microfluidic devices often operate with non-spherical objects (biological cells etc.), the effect of particle shape on inertial migration remains largely unexplored, and a shape-based separation of neutrally buoyant particles is still a challenging problem in microfluidics 9 .
First experiments on inertial focusing of non-spherical particles 10,11 demonstrated the possibility of inertial separation of particles of equal volume (spheres and rods of different aspect ratios) in a straight channel at moderate channel Reynolds numbers Re ≤ 100. More recently, the separation of spheres, ellipsoids and peanut-shaped particles was attained in a spiral microfluidic device 12 , where the inertial lift force is balanced by a drag force with respect to the so-called "Dean flow" resulting from the channel curvature 1 . It was found in both configurations that the rotational diameter of particles is a key parameter which controls their equilibrium positions. However, variation of the lift force across the channel and its dependence on particle shape is unaccessible in experiment and can be obtained only be means of numerical simulation.
Numerical simulation of particle motion in microfluidic devices is a large area of research 13 . One of the most popua) Corresponding author: oivinograd@yahoo.com lar methods for such problems is Lattice Boltzmann method (LBM), which is well-suited for parallel processing and allows efficient tracking of particle-fluid interface 14 . The main difficulty in theoretical description of non-spherical particles is that their migration is strongly coupled with their orientation and rotational regime. The focus of most simulation studies of spheroids is on their rotational behavior in shear flows [15][16][17][18][19] . Typically, at moderate Reynolds numbers oblate particles reorient to perform a log-rolling motion with their minor axis in the vorticity direction, while prolate particles tumble in the flow plane. However, transitions between these rotational regimes may occur at high Reynolds numbers Re ≥ 100, in highly confined flows or under external forces. The regime changes were successfully used for e.g. the separation of elongated biological objects 20 . Simulations of inertial migration have received less attention, however. Variations of stable equilibrium positions and orientations of oblate spheroids in rectangular ducts with varying Reynolds numbers and particle aspect ratios were evaluated in 21 . Su et al. 19 evaluated the lift force on cylindrical particles in rectangular ducts. Particles execute a periodic tumbling motion, so that the lift force is unsteady. However, its average dependence on the particle position is similar to that for a sphere. Dissipative particle dynamics (DPD) simulations were used in Ref. 22 to find the equilibrium positions for prolate and oblate spheroids in a plane Poiseuille flow. Despite these numerous contributions, there is still a lack of understanding of how the lift force depends on the shape of a nonspherical particle.
In this paper we study the inertial migration of oblate spheroids by using the LBM in a plane Poiseuille flow at moderate Reynolds numbers, 11 ≤ Re ≤ 44, relevant for practical applications. We consider the rotational behavior and its effect on particle inertial migration. The lift force is measured for the stable log-rolling regime. We propose a scaling formula which expresses the force on the spheroids in terms of the lift coefficient for spherical particles of the same rotational radii. This formula allows us to predict equilibrium positions of spheroids, either neutrally-buoyant or subject to an external force.
The remainder of the paper is arranged as follows. In sec- tion II we formulate the physical problem and in Section III we deduce the scaling formula for the lift force based on the results of our simulations. Finally, we conclude in Section IV.

II. PROBLEM SETUP
We consider an oblate spheroid with an equatorial radius a and a polar radius b < a in a pressure-driven flow between two parallel walls, separated by a distance H (see Fig.1(a)). The position of the spheroid in the channel is defined by coordinates of its center x = (x, y, z) and by a unit vector directed along the spheroid's axis of symmetry (orientation vector), n = (n x , n y , n z ), |n| = 1. No-slip boundary conditions are applied at the channel walls and the particle surface.
The velocity profile in the channel in the absence of the particle is parabolic, where U m = |∇p|H 2 /(8µ) is the fluid velocity at the channel center, ∇p is a pressure gradient and µ is the dynamic viscosity. The channel Reynolds number Re = ρU m H/µ, where ρ is the fluid density, varies in our simulations in the range Re ≃ 11 − 44 which is typical for inertial microfluidics. The inertial lift force in the channel flow at finite Reynolds numbers drives particles across the flow streamlines. For spherical particles it can be written as 5,7 where G m = 4U m /H is the shear rate at the wall and c l is the lift coefficient, given by Here, the coefficients c li , i = 0, 1, 2 depend on the dimensionless particle position z/H and size a/H as well as on the Reynolds number Re; V s = (V x p − U(z))/U m is a dimensionless slip velocity, where V x p is the x-component of the particle velocity and U is the fluid velocity at the particle center. The slip velocity can arise under an external force parallel to the streamwise direction, e.g. gravity on non-neutrally buoyant particles in vertical channels. In this case it is finite over the entire channel. For neutrally buoyant particles or when the external force is perpendicular to the flow direction, the slip velocity is finite close to the wall and it is very small in the central region 7 .
Equation (2) is widely used to estimate the velocity of migration 1 of neutrally buoyant particles. It also enables us to predict the equilibrium positions, when the lift force is balanced by some external force (e.g. gravity, electric 23 or magnetic 2 fields), F l (z eq ) + F ex = 0, or by the Dean drag in curved channels 1 . To apply a similar approach to shape-based separation one needs to understand how the lift force scales with particle size and aspect ratio.
Obtaining a scaling formula like Eq. (2) for non-spherical particles is not straightforward. The main difference of such particles as compared to spherical ones is a change of orientation due to rotation in a shear flow. Non-inertial spheroids perform a periodic kayaking motion in a shear flow at small particle Reynolds numbers Re p = ρG m a 2 /µ along one of the Jeffrey orbits depending on initial conditions 24 . Such rotation induces unsteady flow disturbances, and the lift becomes time-dependent 19 . The orientation of oblate spheroids at finite Re p eventually tends to a stable state due to the inertia of the fluid and the particle 15,25 . For neutrally-buoyant particles, the log-rolling motion is stable with the symmetry axis parallel to the vorticity direction, n eq = (0, 1, 0) (see Fig. 1(c)). Therefore, a steady lift force for this stable position should be found to predict the long-term migration and equilibrium position of the particle. We extend the scaling formula Eq. (2) to oblate spheroids with a stable orientation in the following section.

III. SIMULATION SETUP
To simulate fluid flow in the channel we use a 3D, 19 velocity, single relaxation time implementation of the lattice Boltzmann method (LBM) with a Batnagar Gross Krook (BGK) collision operator 26,27 . Spheroids are are discretized on the fluid lattice and implemented as moving no-slip boundaries following the pioneering work of Ladd 14 . Details of our implementation can be found in our previous publications 7, [27][28][29][30][31] .
The size of the simulation domain is (N x , N y , N z ) = (128, 128, 81), with corresponding channel height H = 80 (all units are in simulation units). No-slip boundaries are implemented at the top and bottom channel walls using mid-grid bounce-back boundary conditions and all remaining boundaries are periodic. The kinematic viscosity is ν = 1/6 and the fluid is initialized with a density ρ f = 1. A body force directed along x with volumetric density g = 0.5 . . . 2 × 10 −5 is applied both to the fluid and the particle, resulting in a Poiseuille flow with Re ≈ 11...44.
We use neutrally buoyant spheroids of equatorial radii a = 6, 8 and 12 and various aspect ratios, 0.33 ≤ b/a ≤ 1. The particles are initialized close to the expected equilibrium with zero initial velocity and in log-rolling orientation. We assume that the equilibrium is reached when the change of the running average of particle z-coordinate over 10 steps does not exceed 10 −5 . To measure the lift force acting on a particle at some distance z * from the wall, we constrain its z-coordinate to z = z * and launch the particle with zero initial velocity and n = (0, 1, 0) which corresponds to the stable log-rolling state. Once a stationary velocity is obtained, the vertical component of the force is measured and averaged over 10 4 simulation steps.
To check if the results depend on the box size due to periodic boundary conditions in x− and y−directions, we also simulate the migration of the largest sphere (a = b = 12) in a larger simulation box with (N x , N y , N z ) = (256, 256, 81). The difference in equilibrium positions for the two box sizes is 100 times smaller than the typical separation of equilibrium positions of different particles.
To test the resolution of the method in the near-wall zone, we measure the slip velocity of a neutrally-buoyant sphere with a = 8, fixed at various distances from the wall, and compare it to an analytical solution for wall-bounded shear flow 32 . Sufficient accuracy is attained for separations as small as two lattice nodes, similar to previous results for the sphere approaching a rough wall 27 .

IV. NUMERICAL RESULTS AND DISCUSSION
We first study trajectories and orientations of freely moving neutrally buoyant spheroids of different sizes and aspect ratios in a flow with Re = 22. Regardless of the initial position and orientation, oblate spheroids eventually reorient to the stable log-rolling motion with the axis of symmetry and the angu- lar velocity parallel to the y axis, n = (0, 1, 0), ω = (0, ω y , 0). They also focus at some distance z eq from the wall due to inertial migration. The rates of reorientation and migration depend on the particle size and the aspect ratio. We compare in Fig. 2 the rotational behavior and trajectories of spheroids with various aspect ratios, b/a = 1 (sphere), 0.8 and 0.5 and the same equatorial radius a/H = 0.15, initial position z 0 /H = 0.2 and orientation n 0 = (0.66, 0.75, 0). The x-component of the orientation vector n x experiences decaying oscillations (see Fig.2(a)), while n y converges to unity. This means that particles start with ia kayaking motion as shown in Fig. 1(b) and slowly converge to a log-rolling motion (see Fig. 1(c)). The kayaking motion is responsible for oscillations of the trajectories as can be seen in Fig. 2(b). Oscillations in the orientation of a less oblate spheroid with b/a = 0.8 decay much slower than those for b/a = 0.5 (more oblate), so that the migration to the equilibrium position in this case is faster than the reorientation. However, the particle trajectory with b/a = 0.8 in Fig. 2(b) is much less affected by the kayaking motion.
The equilibrium positions for all spheroids with given a/H are very close (see Fig. 2(b)). This means that they are controlled by the equatorial radius a, in full agreement with the experimental observations 10,11 . To validate further this conclusion we evaluate z eq for spheroids of different sizes and aspect ratios.
The results are plotted in Fig. 3 as functions of b/a. We only show a lower equilibrium position since an upper one is symmetric with respect to the channel axis z = H/2 for neutrally buoyant particles. We can see that z eq indeed depends on the equatorial radius a, while the dependence on the aspect ratio is weak. We also find that, while the magnitude of the lift force and the shape of the lift curve c l (z) depend on the Reynolds number, equilibrium positions are almost independent of it for Re = 11 − 44, as has been previously shown for spheres 7 .
Based on these observations, we suppose that the lift coefficient at any position z, not only at z eq , is controlled by the equatorial radius. Thus, we can present the lift force in a way similar to that for a sphere (Eq. (2)), Here, the force depends on the particle position via the lift coefficient c l only (the same as for a neutrally-buoyant sphere), while the function f only includes the dependence on the aspect ratio. To verify Eq. (4) we measure the lift force on spheroids with fixed z position in a stable log-rolling orientation but free to rotate and translate in other directions. We evaluate the ratios of the forces on spheroids and on spheres, F l / ρa 4 G 2 m c l which should be equal to f (b/a). The results shown in Fig. 4 as functions of particle positions confirm our formula (4) since the ratios are nearly constant for given b/a. Moreover, we can deduce from Fig. 4 that f = b/a. This conclusion is verified for three Reynolds numbers, Re = 11, 22, 44. Therefore, we can rewrite Eq. (4) as We evaluate the ratio F/(ρa 3 bG 2 m ) to obtain c l , which follows from (5) and is expected to depend only on a/H. Figure  5 shows c l in the central region, 0.2 ≤ z/H ≤ 0.5. We observe that our results for a fixed a/H and for various b/a do collapse onto the same curves corresponding to the lift coefficient for spheres.
Our scaling formula (5) predicts correctly the lift force for the central region of the channel. However, it is less accurate in a small near-wall region where the gap between the spheroid and the wall, z − a, is much less than a. We plot the lift coefficient for this region in Fig. 6(a) by using Eq. (5). At finite gaps, the results for different spheroids collapse again onto the same curves, but they diverge at small gaps, z − a ≤ 0.2a. The possible explanation for this discrepancy is that the hydrodynamic interaction of spheroids in the near-wall region depends significantly on both particle radii 33 . This effect is illustrated in Fig. 6(b). The slip velocities are finite at small gaps and vary with both a and b/a. Thus their contribution to the lift force is also finite in accordance with Eq. (3). For this reason the behavior of V s and c l is consistent. More oblate particles have a smaller slip velocity and therefore experience a smaller lift force. Figure 6 suggests that the scaling Eq. (5) is valid for practical estimates when the gap between the particle and the wall is greater than 0.2a. Equation (5) can be used to predict equilibrium positions of spheroids with different aspect ratios. Experiments often study the separation of particles of equal volume, V = 4 3 πa 2 b = const, based on a variation of equilibrium positions for such particles 10,11 . Figure 7 shows z eq (b/a) (black symbols) for such oblate spheroids. More oblate particles focus at larger distances from the wall, however, the variation of equilibrium positions is rather small. It follows from Eq. (5) that the lift coefficient and z eq (the zero of c l ) are both independent of b, and that they differ since a changes at constant volume. For comparison we also plot the data by Lashgari et al. 21  at moderate values. Our data for the lift coefficient of spheres with a = 0.075H, 0.1H, 0.15H shown in Fig. 5 can be used to evaluate c l (z/H, a/H) and z eq (b/a) by data interpolation. The equilibrium positions calculated in such manner are plotted in Fig. 7 (solid curve) and agree well with the simulation data. Using Eq. (5) we are also able to predict the equilibrium positions for more complex situations when spheroidal particles are moving under the influence of an external force F ex directed normal to the walls. We naturally assume this force to be proportional to the particle volume, so that F ex = V f ex , where f ex is a volume force. For example, for nonneutrally buoyant particles in a horizontal channel we have f ex = − (ρ p − ρ)g, where ρ p and ρ are the particle and fluid densities and g is the gravitational acceleration. Balancing the two forces, F l (z eq ) + F ex = 0, and using Eq. (5), we obtain the equilibrium position z eq from the following equation: Here, the dimensionless parameter c ex characterizes the relative value of the external force. Equation (6) does not involve b, and this means that the equilibrium positions for particles with the same a under a constant volume force coincide again, similar to neutrally buoyant particles. Using interpolated curves c l (z/H, a/H) we can solve Eq. (6) and predict equilibrium positions of spheroids under an external force for a range of aspect ratios. Figure 7 also shows the equilibrium positions for nonneutrally buoyant spheroids of fixed volume equal to that of a sphere of radius a = 0.1H and for a dimensionless external force c ex = −0.045 (open symbols). This value of c ex corresponds, for example, to spheroids slightly heavier than the fluid, ρ p = 1.1ρ, in horizontal water flow with Re = 22 and H = 200µm. We stress that Eq. (6) has only one zero, corresponding to one equilibrium position, in this case. Indeed, the dependence c l (z/H) is antisymmetric with respect to the channel axis z = H/2, and its maximum, c max l ≃ 0.18 (see also Fig. 5 in 7 ), is smaller than the right-hand side of Eq. (6) for any 0.33 ≤ b/a ≤ 1. Therefore, the upper equilibrium position cannot be attained. Single lower equilibrium positions for heavy spheroids are shifted closer to the bottom wall, and the gaps between the spheres and oblate spheroids are two times greater than for neutrally-buoyant particles. The predictions obtained from Eq. (6) (dashed curve) are in good agreement with the results of direct simulations of particle motion.

V. CONCLUSION
We studied inertial focusing of oblate spheroidal particles in channel flow at moderate Reynolds numbers Re = 11 − 44 using lattice Boltzmann simulations. We found that all spheroids reorient to perform log-rolling motion. They focus to equilibrium positions which depend only on their equatorial radius a, and not on the polar radius b. This is in agreement with experimental observations for particles of various shapes 10 , but has not been directly confirmed by simulations. Based on the measurements of the lift force on logrolling spheroids, we proposed a scaling formula for the lift force Eq. (5), analogous to the widely used scaling for spheres Eq. (2). This scaling is valid throughout the channel except for narrow near-wall regions, for the whole range of studied Reynolds numbers and aspect ratios. Equations (5,6) allow us to predict equilibrium positions of spheroids, both neutrally buoyant and non-neutrally buoyant, by using the lift coefficient c l for the spheres. Therefore, our scaling formula can be used for practical estimates to develop methods of inertial shape-based separation. Our scaling holds for b/a > 0.3 and Re ≤ 44, at large Reynolds numbers equilibrium positions or low aspect ratios could depend on both particle radii. It would be important therefore to estimate the range of validity of Eq. (5) in terms of the particle's aspect ratio and channel Reynolds number.
The natural generalization of this work would be to study inertial migration of prolate particles and to develop some time-averaged analogue of Eq. (5). Indeed, experimental evidence suggests that the largest diameter of the particle is the key parameter, defining its equilibrium position 11 , while its rotation depends on the aspect ratio. However, the problem for prolate spheroids is much more complex, because particles in an equilibrium state tumble in the flow plane and therefore the lift force is unsteady.