Parabolic velocity profile causes shape-selective drift of inertial ellipsoids

Abstract Understanding particle drift in suspension flows is of the highest importance in numerous engineering applications where particles need to be separated and filtered out from the suspending fluid. Commonly known drift mechanisms such as the Magnus force, Saffman force and Segré–Silberberg effect all arise only due to inertia of the fluid, with similar effects on all non-spherical particle shapes. In this work, we present a new shape-selective lateral drift mechanism, arising from particle inertia rather than fluid inertia, for ellipsoidal particles in a parabolic velocity profile. We show that the new drift is caused by an intermittent tumbling rotational motion in the local shear flow together with translational inertia of the particle, while rotational inertia is negligible. We find that the drift is maximal when particle inertial forces are of approximately the same order of magnitude as viscous forces, and that both extremely light and extremely heavy particles have negligible drift. Furthermore, since tumbling motion is not a stable rotational state for inertial oblate spheroids (nor for spheres), this new drift only applies to prolate spheroids or tri-axial ellipsoids. Finally, the drift is compared with the effect of gravity acting in the directions parallel and normal to the flow. The new drift mechanism is stronger than gravitational effects as long as gravity is less than a critical value. The critical gravity is highest (i.e. the new drift mechanism dominates over gravitationally induced drift mechanisms) when gravity acts parallel to the flow and the particles are small.


Introduction
There are numerous practical examples where the understanding of the motion of heavy particles suspended in lighter fluid is important. Typically, this will apply to any solid particle in air, commonly called aerosols. In the atmosphere, aerosols scatter light and thus also affect the global radiation budget (Holländer 1993). These particles also serve as condensation nuclei for cloud formation and their motion inside the clouds is important for understanding rain initiation (Balkovsky, Falkovich & Fouxon 2001;Falkovich, Fouxon & Stepanov 2002) as well as snow crystal growth (Gavze, Pinsky & Khain 2012).
Closer to Earth, aerosols are typically associated with vehicle and industrial emissions causing severe health problems (Morawska & Zhang 2002). Also asbestos fibres in isolation materials in buildings can be directly linked to formation of lung cancer (Miserocchi et al. 2008). When controlling and separating these particles from the suspending fluid, as well as understanding the deposition in the airways, it is important to find how these particles behave in channel flows. Especially important is to understand any physical process causing lateral migration of particles as this will determine if particles, for example, are concentrated in the middle of a channel or towards the walls.
Suspension flows are commonly modelled with spherical particles, due to the many models available to calculate the motion (Crowe et al. 2012). For spherical particles it is known that inertia of the surrounding fluid can cause three types of lateral motion. Firstly, if the particle is not moving with the same velocity as the surrounding fluid, i.e. there is a slip velocity, and the particle is rotating in a constant flow without gradients, there is a lateral force on the particle called a Magnus force (Crowe et al. 2012). Secondly, if the particle has a slip velocity and is rotating in a linear shear flow, there is a lateral force on the particle called a Saffman force (Crowe et al. 2012). Thirdly, if the particle is suspended in a regular pipe flow, the parabolic velocity profile causes itself a migration of particles away from the centreline. As a particle move closer to the wall, the lateral migration force is balanced by pressure that is built up between the particle and the wall, causing the particle to find an equilibrium radial position in the channel. This effect is referred to as the Segré-Silberberg effect (Segré & Silberberg 1961). All these effects are, however, only present when fluid inertia is relevant. If fluid inertia is negligible, which typically is true for aerosols, there is no lateral drift of spherical particles neither caused by a slip velocity nor a quadratic velocity profile.
When it comes to modelling non-spherical particles, it is common to study ellipsoids, and particularly ellipsoids with rotational symmetry called spheroids. In the absence of fluid inertia, Lamb (1932) provided analytical expressions of the force on an ellipsoid in a constant flow given the orientation of the particle and Jeffery (1922) provided analytical expressions of the torque on an ellipsoid in a linear flow field. The force on a suspended inertial ellipsoid can thus be determined by the instantaneous slip velocity and the torque can be determined by local velocity gradients. The results by Jeffery (1922) have, however, mainly been used to study particles that are not affected by particle inertia and therefore will assume a rotation that gives zero torque. Lundell & Carlsson (2010) coupled the torques by Jeffery (1922) to the equations of motion of the particle, and thus the rotational motion of an inertial particle could be simulated assuming that the surrounding flow still had zero inertia. This approach has been used in several studies using Lagrangian particle tracking (LPT) methods in e.g. turbulent channel flow with included particle inertia (e.g. Zhao et al. 2015).
In order to analyse the particle motion, it is convenient to make a Galilean transformation to a system following the particle (see figure 1). The forces on the particle (and thus its motion) can be investigated based on the flow relative to the particle in this system.

Particle in vertical channel
x 2 mG mG Flow around particle after galilean transformation Linear decomposition of u f Velocity Shear (velocity gradient) Shear gradient Figure 1. Illustration of the Galilean transformation for a particle moving in a plane Poiseuille flow. After the transformation, the particle experiences a flow that can be decomposed as shown. The effects of the different flows are additive thanks to the linearity of the Stokes equations for very viscous flow.
Furthermore, the decomposition of the flow seen by the particle into four fundamental flows (one of which is the parabolic velocity profile studied in the present work) is illustrated. The superposition principle of creeping (Stokes) flow implies that the total force on the particle is equal to the sum of the forces exerted by each flow alone. Thus, the effects of curvature studied in the present paper can be added to the effects of shear and homogeneous flow to study particles in all possible quadratic flows.
For the sake of discussion, we will now consider a plane laminar Poiseuille flow with a quadratic velocity profile with suspended spheroidal particles and a dilute concentration. If we were using a LPT method for determining the particle dynamics, which uses the gradient of the velocity but no higher derivatives, we would find that the particles are nicely following the straight streamlines and rotate according to the local shear given by the distance from the centreline. The orientational dynamics of the particles will then be fully determined by their behaviour in a simple shear flow. We can add both fluid and particle inertias to the rotation owing to the recent effort in mapping out the orientational behaviour of spheroids in simple shear flow (Rosén et al. 2016). Due to the symmetry of the simple shear flow, there is no lateral drift of the spheroidal particles. Any drift must thus be caused by higher-order derivatives of the local velocity field. Chwang (1975) studied the spheroidal particle with neither particle nor fluid inertia in a quadratic flow and found no lateral drift. It is still likely that inertia of the surrounding fluid will induce a Segré-Silberberg effect also for the spheroidal particles due to the parabolic velocity profile. In this work, we will show that the influence of particle inertia combined with a quadratic velocity profile will cause a lateral drift even when fluid inertia is neglected. The present results thus provide new fundamental knowledge about the migration of aerosols in channel flows. Figure 2. Illustration of the flow problem; a prolate spheroidal particle is suspended in a quadratic velocity profile; the position of the particle is given by the centre of mass x CM and the orientation is given by the symmetry axis s and the Euler angles θ and φ.
The flow problem is defined in § 2 and the numerical method used is described in § 3. The results are presented in § 4 and discussed in § 5. Two important aspects: consequences for LPT simulations and the effect of gravity, are investigated in § § 6 and 7, respectively. Finally the conclusions are summarized in § 8.

Flow problem
A prolate spheroid with major semi-axis l is suspended in a quadratic background flow according to figure 2. The unit vectors e 1 , e 2 and e 3 denote the flow direction, the velocity gradient direction and the vorticity direction, respectively. The spatial coordinates are given in dimensional form as x * = (x * 1 , x * 2 , x * 3 ). In this work, we will mainly use the non-dimensional coordinates scaled by the particle major semi-axis, i.e. x = x * /l = (x 1 , x 2 , x 3 ). The non-dimensional coordinates of the particle centre of mass is denoted by x CM .
The prolate spheroidal particle with major semi-axis l and equatorial radius l/r p , where r p is the particle aspect ratio (length/width), is described through The non-dimensional coordinates x 1 , x 2 and x 3 are scaled by the major semi-axis l and refer to the body-fixed system spanned by unit vectors e 1 , e 2 and e 3 . The orientation of the particle is given by the direction of the unit vector along the symmetry axis s (note that s = e 1 ). We also express the orientation using the spherical coordinate angles θ and φ such that s = (s 1 , s 2 , s 3 ) = (sin θ cos φ, sin θ sin φ, cos θ) according to figure 2. The (dimensional) background flow is given by u * bg (x) = (1/2)γ l 2 x 2 2 e 1 , where (1/2)γ is the curvature of the velocity profile. The local shear rate at the particle position thus becomesγ L (x * CM,2 ) =γ |x * CM,2 |. As a global time scale in this flow problem, we With these spatial and temporal scalings, the non-dimensional form of the background flow becomes u bg (x) = 1 2 x 2 2 e 1 . (2.2) We assume that fluid inertia is neglected, i.e. that the local particle Reynolds number Re p,L (x CM,2 ) = ρ fγ l 3 |x CM,2 |/μ is zero (ρ f is the fluid density, μ is the fluid dynamic viscosity), so that the flow is governed by the incompressible Stokes equations where u is the fluid velocity and p is the pressure. These equations are non-dimensionalized using characteristic length l, time (γ l) −1 and pressure μγ l. The boundary conditions of the problem are that the background flow is obtained far away from the particle, i.e. lim and that there is no slip of the fluid on the particle surface Γ , i.e.
where V is the (non-dimensional) velocity of the particle centre of mass located at x CM and the (non-dimensional) particle angular velocity is ω. The motion of the particle is determined by the non-dimensional equations of motion where F is the non-dimensional force, M is the non-dimensional torque and I is the non-dimensional inertial tensor, which are scaled by the characteristic force μγ l 3 , torque μγ l 4 and inertial tensor element ρ p l 5 , respectively. Here, ρ p is the density of the spheroidal particle, and Φ = (4π/3)r −2 p is the non-dimensional volume of the particle. In (2.6), St trans. and St rot. are the Stokes numbers characterizing the translational and rotational inertia of the particle, respectively. When deriving the non-dimensional equations of motion using the characteristic quantities mentioned above, these numbers will have the same value, namely where St G is called the global Stokes number. Later, in § 5.2, we will permit St trans. and St rot. to have separate values, but initially they will have the same value and be set according to (2.7). Note that in order for the local Reynolds number Re p,L (x CM,2 ) to be negligible, while the local Stokes number St L (x CM,2 ) := ρ pγ l 3 |x CM,2 |/μ is significant, the solid-to-fluid density ratio must fulfil ρ p /ρ f 1. Of course, this condition will cause the particle to sediment in the presence of gravity. Initially in this study, we will neglect gravitational effects, but these will be discussed in § 7.

Boundary integral formulation
For St G > 0, the equations of motion (2.6) are solved using Matlab's ordinary differential equation solver ode113 with relative tolerance 10 −8 . The force F and torque M are computed numerically from the position, orientation and velocity of the particle at every time step using a boundary integral formulation based on Power & Miranda (1987) and Gonzalez (2009). In this formulation, the flow field in the fluid domain D f is expressed as integrals over the particle surface Γ Here, the Stokes double layer potential D with density q is given by (Einstein's summation convention is used here, with all indices ranging over {1, 2, 3}) The completion flow V is needed to represent the force and torque, and is given by where δ ij is the Kronecker delta and ε ijk is the alternating symbol.
By letting x go to Γ in (3.1) and using the no-slip condition together with a jump property of D, one arrives at a boundary integral equation for q which is closed using the relations For details we refer to af Klinteberg & Tornberg (2016). For St G = 0, which corresponds to a massless particle, the solution procedure is different since we know from (2.6) that F = M = 0. In this case, (3.4) must instead be solved for the velocities V and ω, for which relations similar to (3.5a,b) hold (see af Klinteberg & Tornberg 2016).
Hence, the flow field produced by (3.1) is the solution to the Stokes equations (2.3), given that the density q solves (3.4).

Discretization and quadrature by expansion
The boundary integral equation (3.4) is discretized using the Nyström method (Atkinson 1997, ch. 4), which enforces the equation at the grid points of the discretized particle surface, shown in figure 3. We use the trapezoidal rule with equidistant points in the periodic direction (along a circle of latitude) and an n θ -point Gauss-Legendre quadrature rule in the non-periodic direction (along a meridian). This choice gives us spectral accuracy for smooth and well-resolved integrands on the particle surface (af Klinteberg & Tornberg 2016).
Applying the quadrature rule directly to the integral equation (3.4) is problematic for two reasons. Firstly, the integration kernel of the double layer potential (3.2a,b) is singular when the evaluation point x coincides with a point y on the boundary, and can therefore not be handled by a quadrature rule for smooth functions. Moreover, when (3.1) is evaluated in a point x which is close to the boundary, but not on it, the integration kernel of (3.2a,b) is sharply peaked, which causes a significant loss of accuracy close to the particle surface. We use a recent method called quadrature by expansion (QBX) to treat both of Figure 3. The spheroidal grid. Due to symmetry it is enough to store precomputed matrices for the n θ /2 grid points along the first half-meridian, marked with circles. these problems. The idea behind this method is described briefly below; for a detailed description, we refer to af Klinteberg & Tornberg (2016). QBX is based on the observation that the double layer potential D is a smooth function away from the boundary. To avoid the problems close to the boundary, we can create a local expansion of D around a point c away from the boundary. This expansion can be used to evaluate D at exactly one point on the boundary, as shown in figure 4. To solve the integral equation we thus need one expansion centre for every grid point. Since D has different limits from the interior and the exterior, we use both an inner and an outer expansion centre and compute the potential D QBX as an average to get the value on Γ .
There exist matrices D j such that D QBX [Γ, q](x j ) = D j Q for each grid point x j on the surface, where Q is a vector containing the values of the density q in all grid points. The matrices D j depend only on the geometry of the spheroid and can be precomputed; the result can be seen as a regular but target-specific quadrature rule. Due to rotational and mirror symmetry, it is sufficient to store matrices for the n θ /2 grid points along the first half meridian, shown in figure 3. The precomputation allows the method to be both fast and accurate, since precomputation can be done once with high accuracy and the result is reused in every time step.

Validation
The method has been validated against test cases for both inertial and massless particles, in linear shear flow and the quadratic background flow considered here. The analytical solutions in these cases are provided by Jeffery (1922) and Chwang (1975); the validation is further described in Bagge (2016). Using the parameters given in Appendix A, the relative errors are below 10 −6 in all test cases. x 3 x 1

Results
We start by considering a prolate spheroid of r p = 3 and St G = 50 initialized at rest at position x CM = (0, 1, 0) with a slightly oblique orientation. The trajectories of the particle translation and orientation are illustrated in figure 5. The orientation of the particle is clearly drifting towards a rotation around its minor axis, an intermittent rotation that we call tumbling. This is the same type of behaviour seen in a simple shear flow by Lundell & Carlsson (2010). What is more striking is what is happening to the translational motion of the particle. After an initial transient the particle ends up laterally drifting towards regions of higher shear.
To quantify the observed drift, we consider a particle of r p = 3 that is initialized at position x CM = (0, 1, 0) with velocity V = u bg (x CM ), aligned in the flow-gradient direction, i.e. at s = (0, 1, 0) with zero angular velocity. Note that at the given initial position, the local and global Stokes numbers are the same, i.e. St G = St L (x CM,2 = 1). The particle is then free to both rotate and translate up until t = 80. This is repeated for a number of different St G in the range St G ∈ [0, 100]. The resulting trajectory of the particle centre of mass is shown in figure 6. While the spheroid with no inertia (St G = 0) does not show any lateral motion at all, as soon as there is particle inertia, the particle starts to migrate in positive x 2 -direction, i.e. towards higher shear. Note, that the trajectories are plotted as function of the global time scale even though the local relevant time scale is changing (due to increasing local shear) with x CM,2 . However, this effect is negligible since the particle has not moved significantly in the x 2 -direction at t = 80. For all St G studied here, the particle assumes its final state (a periodic rotation) quite quickly after an initial transient while the particle is accelerating to the surrounding flow. The final average lateral drift velocity V drift = V 2 for each St G is seen to be constant, and is estimated by fitting a linear function to x CM,2 (t) between times t = 20 and 80 as seen by the dashed lines in figure 6. From the figure it is also evident that there is a critical Stokes number St c ≈ 30 such that V drift has a maximum value when Considering now a particle with different aspect ratios r p = 1, 2, 3, 4 at St G = 50 in figure 7(a), we find that a spherical particle, as expected, has no lateral drift. As soon as the particle gets a prolate shape, it gets a constant drift velocity, which at St G = 50 is increasing with r p . In figure 7(b), the final drift velocity V drift is plotted as function of St G for the different r p and we find that both the maximum V drift,max and the Stokes number St c where this occurs are indeed increasing with r p .
As previously mentioned, if we would release the particle at different heights x CM,2 in the flow, the particle would experience different shear rates and the relevant time scale would change. In this case the local St L = |x CM,2 |St G depends on lateral position in the flow, which in turn determines the drift velocity. To demonstrate this, we released the particle of r p = 3 at two more locations x CM,2 = 1/2 and x CM,2 = 2. We can see in figure 8(a) that V drift is only dependent on St L as all the curves collapse when the velocity is plotted vs this parameter. Consequently, there will always be a particle position in the channel x CM,2 = x CM,2,c where we have a maximum drift velocity, which is where If we e.g. start at |x CM,2 | < 1, we will always eventually reach this position since the particle is drifting towards regions of higher shear. Eventually, the drift will vanish as St L → ∞, which is demonstrated in figure 8(b) where the results of V drift at very high St L are illustrated.

Discussion
The previous section has shown the results that particle inertia is sufficient to cause a lateral drift of a particle in a quadratic velocity profile. The source of the resulting translational and rotational motion will be discussed in this section.  Jeffery (1922) and Lamb (1932) present expressions for the force F 0 and torque M 0 on a translating and rotating spheroidal particle in a quiescent fluid, while Chwang (1975) derived expressions for the force F Chw and torque M Chw on a non-translating and non-rotating spheroidal particle in a paraboloidal velocity profile. One important thing to note is that the results by Chwang (1975) were derived for a paraboloidal profile u Chw bg (x) = (1/2)(x 2 2 + x 2 3 )e 1 . By adjusting these results, we can express the force F par and torque M par valid for the parabolic background velocity profile considered here u bg (x) = (1/2)x 2 2 e 1 . This derivation is given in Appendix B (it turns out that the torque is the same in both cases, M par = M Chw ).

Without inertia
Through the principle of superposition, the total force and torque on a moving and rotating spheroidal particle in a parabolic velocity profile in the absence of fluid inertia can be found by combining the adjusted results by Chwang (1975) with the expressions of Jeffery (1922) and Lamb (1932), i.e. F = F 0 + F par and M = M 0 + M par . If we want to study the free motion of the spheroid with no particle inertia (St trans. = St rot. = 0 in (2.6)), the resulting force F and torque M are set to zero.

Rotation without inertia
Let us now consider a particle that is oriented in the flow-gradient plane (θ = π/2). The torque on a spheroidal particle rotating in the flow-gradient plane with angular velocityφ in a quiescent fluid is given in non-dimensional form by Jeffery (1922) as where The torque on a stationary particle with orientation φ and location x CM in a parabolic velocity profile is given through the (non-dimensional) expression by Chwang (1975) as x CM,2 e 2 , (5.3) where E = 1 − r −2 p is the eccentricity of the spheroid. The total torque on the particle rotating arbitrarily in a parabolic velocity profile is given by M = M 0 + M par . The resulting expression is actually exactly equivalent to the expression by Jeffery (1922) in a linear shear flow. The conclusion drawn already by Chwang (1975) is thus that the quadratic terms in the background flow do not affect the particle rotation. Consequently, the particle without rotational inertia (St rot. = 0), which rotates according to the solution of M = 0, will just rotate according to the local shear rate in an intermittent tumbling motion described by Jeffery (1922) aṡ The rotational period is given by

Translation
The force on a spheroidal particle moving with instantaneous velocity V and orientation φ in a quiescent fluid is given by Lamb (1932) as where the tensor K in the body-fixed coordinate system (equivalent to φ = 0) is given by . (5.8) If φ / = 0, the tensor will simply be transformed according to The force on a stationary particle with orientation φ (between s and u CM ) and location x CM in a parabolic velocity profile is given through the adjusted expression by Chwang (1975) The total force on the particle arbitrarily translating and rotating in a parabolic velocity profile is given by F = F 0 + F par . If the particle is oriented with an angle φ and free to translate in the absence of translational inertia (St trans. = 0), the velocity that satisfies F = 0 is given by (5.10) The particle has a velocity in the flow direction, which depends on the instantaneous angle φ. If the particle is free to rotate without rotational inertia, it will perform an intermittent tumbling motion according to (5.4). Consequently, this thus results in an intermittent translational motion in the flow direction.

Only rotational inertia
Although not easily achievable in practice, let us now consider the case where St rot. > 0 but translational inertia is neglected St trans. = 0. Since we know from the expressions above that the translational motion is only in the flow direction in the absence of translational inertia and the rotational motion is only dependent on the motion in the gradient direction, the translation will not influence the particle rotation. With added rotational inertia, the particle will thus behave exactly according to Lundell & Carlsson (2010), where the instantaneous torque during the final rotational motion typically is non-zero but fulfils T 0 Mdt = 0 for the rotational period T. At extremely large St rot. → ∞, the particle will rotate with a constant angular velocityφ = −0.5|x CM,2 | and period T H = 4π/(|x CM,2 |) as the particle inertial forces overcome viscous forces from the surrounding fluid. A critical rotational Stokes number St 0.5 was introduced by Lundell & Carlsson (2010) to describe the transition. This number is defined as the St rot. where the tumbling period is T = (T J + T H )/2. With corrections given by Nilsen & Andersson (2013), this number can be found through Lundell & Carlsson (2010) also found that the particle initially oriented out of the flow-gradient plane always still drifted towards the tumbling motion due to the centrifugal forces on the particle. The rate of this orbit drift was seen to be close to maximum when St rot. = St 0.5 . Interestingly, the critical rotational Stokes number St 0.5 for maximum orbit drift seems to scale similarly with aspect ratio as the critical translational Stokes number St c for the lateral drift as illustrated in figure 9. This indicates that the maximum lateral drift also arises in the transition when inertial forces overcome viscous damping. The fact that the instantaneous torque on the particle is exactly represented by the solutions by Jeffery (1922) is also confirmed by the Stokes flow simulations in the present work. In figure 10, we plot the torque M 3 on the particle as a function of orientation φ and angular velocityφ. Furthermore, we superimpose the trajectories of the numerical simulations from our results at higher St L = St rot. = St trans. . Here, we can clearly see how the particle travels through regions of both positive and negative torque during a rotation and approaching constant angular velocityφ = −0.5|x CM,2 | as St L increases. We find in the numerical simulations also that the translational inertia has no effect on the rotation, except for the fact that the lateral drift changes the local time scale (through the local shear rate) and thus also the local Stokes number St L . -8 -π/2 -π/4 π/4 π/2 0 φ φ Figure 10. Torque M 3 given on a particle (r p = 3) in a quadratic flow at x CM,2 = 1 with no translational velocity in a tumbling orbit (s 3 = cos θ = 0) as function of orientation φ and angular velocityφ; the solid, dashed and dotted lines show the superimposed path of a particle that is free to translate at St L = 0, 50 and 100, respectively.

Adding translational inertia
With no translational inertia, the particle moves with oscillating V 1 -velocity (V 2 = 0) corresponding to F = 0. The oscillation period corresponds to the oscillation period of φ.
With added rotational inertia (still no translational inertia) as we observed previously, there will be a different oscillation period, but V 1 is still given by the F = 0 solution, which does not imply any lateral drift (V 2 = 0). With added translational inertia (no rotational inertia) the particle will be slow to react to the forces and the particle will typically experience non-zero forces but eventually leading to a translational motion that fulfils T 0 F dt = 0. This solution actually has V 2 / = 0 and the particle will drift laterally. With added rotational inertia, the oscillation period of the velocity will decrease as the oscillation period of φ will decrease. The fact that the instantaneous force on the particle is exactly represented by the analytical expressions of F = F 0 + F par is confirmed in the present Stokes flow simulations. In figure 11, we illustrate the instantaneous forces π/2 -π/4 π/4 0 -π/2 π/2 -π/4 π/4 0 φ φ Figure 11. Forces given on a particle (r p = 3) with fixed orientations of a tumbling orbit (s 3 = cos θ = 0) as functions of orientation φ and streamwise velocity V 1 ; the solid, dashed and dotted lines show the superimposed path of a particle that is free to translate at St L = 0, 50 and 100, respectively: (a) force component  parallel (F 1 ) and normal (F 2 ) to the flow, for a given orientation φ and streamwise velocity difference V 1 − u bg,1 (x CM,2 ). Superimposed, we again plot the trajectories of the numerical Stokes flow simulations at higher St L = St rot. = St trans. . Here, we can again observe how the particle travels through regions of positive and negative forces, and that the velocity oscillations will decrease in amplitude and approach a constant velocity as St L is increased.
To really conclude what is causing the drift, additional simulations were performed where St trans. was varied independently from St rot. . The resulting drift velocities V drift after the initial transient are summarized in table 1. It is found that V drift is mainly dependent on the translational inertia of the particle. Even though the oscillation period changes, the effect of St rot. at a constant St trans. is almost negligible.
It is also quite clear what will happen with oblate spheroids. Since the inertial oblate spheroid will drift towards a rotation around its symmetry axis, there will be no 'jerk' in the translational motion, the force will be the same regardless of rotational phase and the particle will assume the velocity corresponding to zero force. There will thus not be any lateral motion for oblate particles.

Additional remark
The present results of the numerical Stokes flow simulations have demonstrated that the flow problem can be completely analysed analytically by integrating the equations of motion (2.6) using the force and torque expressions by Jeffery (1922) and Lamb (1932) and the modified expressions by Chwang (1975) in (5.1)-(5.3) and (5.6)-(5.9). The QBX method, however, can be utilized in the future for more complicated flow problems, e.g. including complex geometries, where analytical solutions are difficult to obtain.

Consequences for simulations with Lagrangian particles
The fact that there is an additional force that arises from the second spatial derivative of the velocity also has consequences for LPT methods. In these schemes, the force on the particle is calculated by knowing the (dimensional) velocity difference ΔU * = u * CM − V * between the particle and the undisturbed fluid. Similarly, the torque is found only through the orientation and the local velocity gradients, i.e. the first spatial derivatives of the velocity. The second spatial derivative was seen here to only affect the translation and not the rotation. So how can we evaluate if this contribution to the force is relevant?
Consider a particle in a Lagrangian frame centred on the particle with flow direction e 1 and gradient direction e 2 with constant curvature of the velocity profile given byγ . The particle experiences the velocity ΔU * and will thus have dimensional force contribution according to Lamb (1932) as Since the Lagrangian frame is always centred on the particle, the force contribution due to the second spatial derivative of the velocity will be equivalent to setting x CM,2 = 0 in (5.9). The dimensional result thus becomes )/E 3 are parameters determined by the particle geometry. Furthermore, all geometry dependent parameters K , E, A 1 and A 2 are O(1) for moderate aspect ratios. In order to neglect the contribution of the velocity curvature in an LPT simulation when calculating the force, the condition |γ |l 2 /|ΔU * | 1 should thus hold. As a final remark, it is likely that there is a similar correction to the rotation by taking into account the third spatial derivative of the velocity.

Effects of gravity
In order to determine to what extent the drift mechanism scrutinized above will be relevant in a physically realizable system on Earth, it must be compared with the effect of gravity. For inertial particles in Stokes flow, gravity can induce a sideways drift already under the assumption of linear shear (without curvature), shown by Broday et al. (1998). This drift appears since oblique particles sediment sideways and intermediate inertia gives non-symmetric angular distributions. In figure 12, the effect of gravity is investigated. In (a) and (b), the particle motion is shown at different levels of non-dimensional gravitational force where G is the physical (dimensional) gravitational acceleration, and the other parameters are as introduced in § 2. The particle motion is shown for aspect ratio r p = 3 and Stokes number St G = 30 (at which the curvature-induced drift is approximately maximal for this aspect ratio). From the particle motions, the sideway drift can be determined as a function of g, as shown in (c,d). Due to the curvature-induced drift, the net drift is not 0 for g = 0. In fact, for aspect ratio r p = 3, the curvature-induced drift is seen to dominate for |g| < 0.014 with gravity normal to the flow direction and |g| < 0.72 for gravity parallel to the flow direction. The critical gravities for other aspect ratios are shown in  figure 12(c,d) with gravity acting parallel and normal to the flow direction, respectively. Column 5 shows the particle half-length corresponding to g = g crit, , ρ p = 1000 kg m −3 and G = 10 m s −2 in (7.5). Column 6 shows the corresponding Reynolds number due to sedimentation, from (7.6).
due to sedimentation must be small for the creeping condition to be satisfied. The sedimentation velocity U sed is estimated from the vertical force balance between buoyancy and gravity (cf. (6.1)), assuming K to be of O(10) (this has been verified numerically) . The lines are r p = 3 (blue) and r p = 10 (red).

Conclusions
In this work we have investigated a new shape-selective drift mechanism which applies to inertial prolate ellipsoidal particles in a parabolic velocity profile. The results are derived with the assumption that inertia of the fluid flow around the particle is negligible. The work can thus be viewed as a direct continuation of the work by Chwang (1975) to account for particle inertia. It is found that non-sphericity, particle inertia and the quadratic profile contribute to a lateral drift towards regions of higher shear. The particle inertial effects are governed by the local Stokes number St L based on the local shear rate. In the absence of inertia (St L = 0), when the particle is aligned in the flow-gradient plane, it assumes an intermittent tumbling motion where the angular velocity depends on the angle. This intermittent rotation is determined only by the local shear rate and described in detail by Jeffery (1922). This causes the particle in the parabolic profile to have an intermittent streamwise velocity depending on the instantaneous orientation (Chwang 1975).
With higher St L , the particle behaves exactly as expected in a linear shear flow according to the results by Lundell & Carlsson (2010). However, due to its translational inertia, the particle is not quick enough to adapt to the non-zero forces during a particle rotation. The final trajectory of the particle, when the impulse during a rotational period is zero, describes a translational motion with a lateral drift. At very high St L , the particle will approach a rotation with constant angular velocity and also a translation with constant velocity, and the lateral drift vanishes.
If the particle is not initially oriented in the flow-gradient plane it will drift towards such an orbit. Since the particle is seen to have angular dynamics which are identical to a particle in simple shear, we know from Lundell & Carlsson (2010) that the maximum orbit drift occurs at St L ≈ St 0.5 . The critical value of St 0.5 is defined when particle inertia overcomes viscous damping and can be found analytically, see (5.11). We find furthermore that this critical number also predicts when we can expect the maximum lateral drift.
Since the particle tends to migrate towards regions of larger shear, we know that a particle released at a position where St L < St 0.5 will reach both a position where St L ≈ St 0.5 with a maximum drift angle (angle between particle velocity and flow direction) and later have a decreasing drift angle as the particle migrates further in the lateral direction.
The maximum lateral drift velocity was found to be dependent on the particle aspect ratio, with higher drift for larger aspect ratio r p . However, since the intermittent tumbling motion is essential to have this type of drift, inertial oblate particles and spheres (r p ≤ 1) are not affected by this drift since they preferably rotate in a non-intermittent state.
Although the results were found through simulations of the Stokes flow using the numerical QBX method, it is found that the forces and torques on the particle are exactly represented by the analytical expressions by Lamb (1932) and the modified expressions of Chwang (1975) (with modifications presented in Appendix B). This means that the equations of translational and rotational motion can be directly integrated using these expressions to arrive at the same conclusions as in the present work. The strength in the QBX method lies mainly in the fact that we can additionally use the same method for studying other flow problems where analytical solutions are more difficult to obtain, for example having multiple particles and wall-bounded flow through complex geometries.
In LPT methods, the force on the particle is usually calculated using only the instantaneous velocity difference ΔU * between the particle and the undisturbed fluid. Here, it is demonstrated that there is also a contribution to the force from the curvature of the velocity profileγ , and that |γ |l 2 /|ΔU * | 1 must hold for this contribution to be negligible.
Finally, it is known from recent work (e.g. Einarsson et al. 2015a,b) that fluid inertial effects are dominating over particle inertial effects for spheroidal particles of same density as the surrounding fluid. The consequence is that we can only neglect fluid inertia, if the local Stokes number St L is much larger than the local Reynolds number Re p,L 1, i.e. the solid-to-fluid density ratio must be large. In the presence of gravity this will cause the particles to sediment, and we proceeded to investigate under what circumstances the new drift would dominate over gravitational effects. We find that the new drift is more likely to dominate when particles are small and light, and when the flow is vertical, i.e. parallel to the direction of gravity. As an example, on Earth, particles of density 1000 kg m −3 moving through air would need to be smaller than around 10 μm for the new drift to be able to dominate over gravity.

Appendix A. QBX parameters
The parameters used are given in table 3. Here, n φ and n θ are the number of grid points on the spheroid in the periodic and non-periodic direction, respectively (as explained in § 3.2); r/h controls the distance from each expansion centre to the particle surface n φ = 16 n θ = r p n φ r/h = 0.6 κ = 30 p = 30 Tolerance for gmres = 10 −12 Tolerance for ode113 = 10 −8 Table 3. Parameters used for the QBX method.
(h = 2πlr −1 p /n φ ); κ is the upsampling factor of the grid; p denotes the order of the expansions. The Matlab functions gmres and ode113 are used for solving systems of linear equations and ordinary differential equations, respectively. For a more detailed explanation of these parameters, we refer to af Klinteberg & Tornberg (2016).