Parabolic velocity profile causes drift of inertial prolate spheroids -- but gravity is stronger

Motion of elongated particles in shear is studied. In applications where particles are much heavier than the carrying fluid, e.g. aerosols, the influence of particle inertia dominates the particle dynamics. Assuming that the particle only experiences a local linear velocity profile, its rotational and translational motion are independent. However, we show that quadratic terms of the local velocity profile combined with particle inertia cause a lateral drift of prolate spheroidal particles. We find that this drift is maximal when particle inertial forces are of the same order of magnitude as viscous forces, and that both extremely light and extremely heavy particles have negligible drift. In the non-inertial case, the particle rotates according to the local linear velocity profile, with each instantaneous orientation corresponding to a velocity that gives zero force on the particle. This results in a translational motion in the flow direction with periodic velocity fluctuations. With added particle inertia, the particle is slow to react to the surrounding fluid motion and the particle will switch between slower and faster rotation compared to the zero-force solution. The final motion that gives zero integrated force over a rotational period, is a motion with a lateral drift. We show that this drift is purely an effect of the non-sphericity of the particle and its translational inertia, while rotational inertia is negligible. Finally, although this inertial drift will contribute to the lateral motion of heavy elongated particles in channel flow, sedimentation due to gravity will dominate in any practical application on earth.


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 et al. 2001;Falkovich et al. 2002) as well as snow crystal growth (Gavze et al. 2012).
Closer to earth, aerosols are typically associated with vehicle and industrial emissions causing severe health problems (Morawska & Zhang 2002). Also asbestos fibers 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 of the spherical particles. 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 either caused by a slip velocity or 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) has 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).
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 would be using a LPT method for determining the particle dynamics that 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 inertia 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 without either particle or 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  Figure 1: 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 CM and the orientation is given by the symmetry axis and the Euler angles and .
when fluid inertia is neglected. The present results thus provide new fundamental knowledge about the migration of aerosols in channel flows. The flow problem is defined in Section 2 and the numerical method used is described in Section 3. The results are presented in Section 4 and discussed in Section 5. Two important aspects: consequences for LPT simulations and the effect of gravity, are investigated in Sections 6 and 7, respectively. Finally the conclusions are summarized in Section 8.

Flow problem
A prolate spheroid with major semi-axis is suspended in a quadratic background flow according to figure 1. The unit vectors 1 , 2 and 3 denote the flow direction, the velocity gradient direction and the vorticity direction, respectively. The spatial coordinates are given in dimensional form as * = ( * 1 , * 2 , * 3 ). In this work, we will mainly use the non-dimensional coordinates scaled by the particle major semi-axis, i.e. = * / = ( 1 , 2 , 3 ). The nondimensional coordinates of the particle centre-of-mass is denoted by CM .
The prolate spheroidal particle with major semi-axis and equatorial radius / p , where p is the particle aspect ratio (length/width), is described through 2 1 + 2 p 2 2 + 2 p 2 3 = 1. (2.1) The non-dimensional coordinates 1 , 2 and 3 are scaled by the particle length and refer to the body-fixed system spanned by unit vectors 1 , 2 and 3 . The orientation of the particle is given by the direction of the unit vector along the symmetry axis (note that = 1 ). We also express the orientation using the spherical coordinate angles and such that = ( 1 , 2 , 3 ) = (sin cos , sin sin , cos ) according to fig. 1. The (dimensional) background flow is given by * bg ( ) = (1/2) 2 2 2 1 , where (1/2) is the curvature of the velocity profile. The local shear rate at the particle position thus becomes L ( * CM,2 ) = | * CM,2 |. As a global time scale in this flow problem, we choose We assume that fluid inertia is neglected, i.e. that the local particle Reynolds number Re p,L ( CM,2 ) = f 3 | 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 which are non-dimensionalized using characteristic length , time ( ) −1 and pressure . The boundary conditions of the problem is that the background flow is obtained far away from the particle, i.e.
and that there is no slip of the fluid on the particle surface Γ, i.e.
where is the (non-dimensional) velocity of the particle centre-of-mass located at CM and the (non-dimensional) particle angular velocity is . The motion of the particle is determined by the equations of motion where St trans. and St rot. are the Stokes numbers characterizing the translational inertia and rotational inertia, respectively. We initially set these to the same value, corresponding to a spheroid with uniform density p . Both translational and rotational inertia will thus be set using the global Stokes number St G : (2.7) The parameter Φ = (4/3) −2 p is the non-dimensional volume of the spheroidal particle. The non-dimensional force, torque and inertial tensor are given by , and I, respectively. The equations are non-dimensionalized with characteristic torque 4 and the characteristic inertial tensor element p 5 . Note that in order for the local Re p,L ( CM,2 ) to be negligible, while the local St L ( CM,2 ) = p 3 | CM,2 |/ is relevant, the solid-to-fluid density ratio must fulfill = p / f 1. Of course, this condition causes the particle to sediment in the presence of gravity. Initially in this study, we will neglect gravity effects but we will discuss these effects in Section 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 and torque 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 f is expressed as integrals over the particle surface Γ: (3.1) Figure 2: The spheroidal grid. Due to symmetry it is enough to store precomputed matrices for the /2 grid points along the first half meridian, marked with circles.
Here, the Stokes double layer potential D with density 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 is the Kronecker delta and is the alternating symbol. By letting 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 : 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 = = 0. In this case, (3.4) must instead be solved for the velocities and , for which relations similar to (3.5) 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 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 2. We use the trapezoidal rule with equidistant points in the periodic direction (along a circle of latitude) and an -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.2) is singular when the evaluation point coincides with a point 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 which is close to the boundary, but not on it, the integration kernel of (3.2) 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 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 away from the boundary. This expansion can be used to evaluate D at exactly one point on the boundary, as shown in figure 3. 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 such that D QBX [Γ, ] ( ) = D for each grid point on the surface, where is a vector containing the values of the density in all grid points. The matrices D 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 /2 grid points along the first half meridian, shown in fig. 2. 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.

Results
We start by considering a prolate spheroid of p = 3 and St G = 50 initialized at rest at position CM = (0, 1, 0) with a slightly oblique orientation. The trajectories of the particle translation and orientation are illustrated in figure 4. 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 p = 3 that is initialized at position   CM = (0, 1, 0) with velocity = bg ( CM ), aligned in the flow gradient direction, i.e. at = (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 ( CM,2 = 1). The particle is then free to both rotate and translate up until = 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 5. 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 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 CM,2 . However, this effect is negligible since the particle has not moved significantly in the 2 -direction at = 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 drift = 2 for each St G is seen to be constant, and is estimated by fitting a linear function to CM,2  figure 6a, 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 p . In fig. 6b, the final drift velocity drift is plotted as function of St G for the different p and we find that both the maximum drift,max and the Stokes number St c where this occurs are indeed increasing with p .
As previously mentioned, if we would release the particle at different heights 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 = | 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 p = 3 at two more locations CM,2 = 1/2 and CM,2 = 2. We can see in figure 7a that drift is only dependent on St L as all the curves collapse when the velocity is plotted versus this parameter. Consequently, there will always be a particle position in the channel CM,2 = CM,2,c where we have a maximum drift velocity, which is where St L ( CM,2,c ) = St c . If we e.g. start at | 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 fig. 7b where the results of 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 0 and torque 0 on a translating and rotating spheroidal particle in a quiescent fluid, while Chwang (1975) derived expressions for the force Chw and torque Chw on a non-translating and non-rotating spheroidal particle in a parabolic velocity profile. One important thing to note is that the results by Chwang (1975) were derived for a paraboloidal profile Chw bg ( ) = (1/2) ( 2 2 + 2 3 ) 1 . By adjusting these results, we can express the force par and torque par valid for the here considered parabolic background velocity profile bg ( ) = (1/2) 2 2 1 . This derivation is given in appendix B.

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. = 0 + par and = 0 + 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 and torque 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 The torque on a stationary particle with orientation and location CM in a parabolic velocity profile is given through the (non-dimensional) adjusted expression by Chwang (1975)  where = √︃ 1 − −2 p is the eccentricity of the spheroid. The total torque on the particle rotating arbitrarily in a parabolic velocity profile is given by = 0 + 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 = 0, will just rotate according to the local shear rate in an intermittent tumbling motion described by Jeffery (1922) as ( , CM ) = | CM,2 | 2 p sin 2 + cos 2 2 p + 1 .

Translation
The force on a spheroidal particle moving with instantaneous velocity 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 K ( ) = R ( )K (0)R −1 ( ) with the rotation matrix R. The force on a stationary particle with orientation (between and CM ) and location CM in a parabolic velocity profile is given through the adjusted expression by Chwang ( The total force on the particle arbitrarily translating and rotating in a parabolic velocity profile is given by = 0 + 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 = 0 is given by ( , CM ) = 1 6 3 2 CM,2 + 1 − 2 cos 2 1 . (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.

Adding inertia 5.2.1. 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 fulfilling ∫ 0 d = 0 for the rotational period . At extremely large St rot. → ∞, the particle will rotate with a constant angular velocity = −0.5| CM,2 | and period = 4 /(| CM,2 |) as the particle inertial forces overcome viscous forces from the surrounding fluid. A critical rotational Stokes number Lundell & Carlsson (2010) also found that the particle initially oriented out of the flowgradient 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 8. 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 9, we plot the torque 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| 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  Figure 9: Torque given on a particle ( p = 3) in a quadratic flow at CM,2 = 1 with no translational velocity in a tumbling orbit ( 3 = cos = 0) as function of orientation and angular velocity ; the solid, dashed and dotted lines shows the superimposed path of a particle that is free to translate at St L = 0, 50 and 100, respectively.
fact that the lateral drift changes the local time scale (through the local shear rate) and thus also the local Stokes number St L .

Adding translational inertia
With no translational inertia, the particle moves with oscillating 1 -velocity ( 2 = 0) corresponding to = 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 1 is still given by the = 0 solution, which does not imply any lateral drift ( 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 fulfills ∫ 0 d = 0. This solution actually has 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 = 0 + par is confirmed in the present Stokes flow simulations. In figure 10, we illustrate the instantaneous streamwise and lateral forces 1 and 2 for a given orientation and streamwise velocity difference 1 − bg,1 ( 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 drift after the initial transient are summarized in table 1. It is found that 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 analyzed 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 Lagrangian particle tracking (LPT) methods. In these schemes, the force on the particle is calculated by knowing the (dimensional) velocity difference Δ * = * CM − * 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 1 and gradient direction 2 with constant curvature of the velocity profile given by . The particle experiences the velocity Δ * and will thus have dimensional force contribution according to Lamb (1932) 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 CM,2 = 0 in (5.9). The dimensional result thus becomes * curv. = 16 2 1 − 2 cos 2 3 cos 1 2 1 − sin 2 2 = 16 2 1 − 2 cos 2 3 K ( ) 1 , (6.2) where 1 = (−2 + (1 + 2 ) log[(1 + )/(1 − )])/ 3 and 2 = (2 + (3 2 − 1) log[(1 + )/(1 − )])/ 3 are parameters determined by the particle geometry. Furthermore, all geometry dependent parameters K , , 1 and 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 | | 2 /|Δ * | 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.

Comparison with gravitational forces
In this work we have observed that prolate spheroidal particles in a parabolic flow gives rise to a lateral drift. This is however only true as long as fluid inertia is neglected while particle inertia is dominating, i.e St L Re p,L and Re p,L 1, which means that the density ratio = p / f must fulfill 1 (7.1) and This new mechanism then seems to be relevant for solid particles suspended in air. However, as soon as p ≠ f , there is sedimentation in the presence of gravity. We must therefore carefully examine when the dimensional drift velocity found here * drift = drift 2 ( drift is the dimensionless velocity found in fig. 6  with the correction term ( p , ) being given by where is the angle between gravity and particle orientation and the functions 1 and 2 are defined as For argument's sake, let us consider a slender fiber and that we need to fulfill = p 10.
(7.11) Given that the solid particle is suspended in air somewhere on earth, we can assume that = 10 −5 Pa s, f = 10 0 kg/m 3 , p = 10 3 kg/m 3 and = 10 1 m/s 2 . This gives us a condition on the velocity profile curvature of 10 10 m −1 s −1 . (7.12) How much is this? If we superimpose a constant flow in opposite direction, we will end up with a parabolic velocity profile similar to that in a planar channel flow. Given a plane Poiseuille flow, can be expressed in terms of the pressure gradient (7.13) which means that in order to achieve the previous inequality, we need a pressure gradient of d * d * 10 5 Pa/m, (7.14) i.e. a pressure drop of 1 bar per meter of the channel. In order to keep such a flow laminar given a channel of height , the channel Reynolds number Re must fulfill Re = f 3 10 3 , (7.15a) i.e. 10 −4 m. (7.15b) For the particle to fit inside such a channel, it must be 10 −4 m. Given our assumption of the particle and fluid density, the creeping flow assumption now reads i.e. 10 3 m 3 10 3 × 10 10 3 10 −5 , (7.16c) i.e. 10 −5/3 m, (7.16d) i.e. the particle must be smaller than 2 mm, which is less strict than the previous assumption of the particle size. One might argue that the effects presented in this work then is dominating for slender fibers with the size of tenths of microns suspended in air. However, in any practical application, having these microfibers suspended in a channel of sub-millimeter size with a velocity around * ≈ 100 m/s, lacks true relevance. Furthermore, we have not began to comment on other effects that might arise due to these higher pressures and small length dimensions of the particle, such as compressibility, temperatures, Brownian motion and the break-down of the continuum assumption of the fluid. These arguments quickly diverge into absurdity, and it is safe to assume that the inertial drift due to the quadratic velocity profile will never be dominating over sedimentation in any earthly application. Even if it is not dominating, the inertial migration described in this work should be regarded as a correction factor to a sedimentation process in a flow where fluid inertia of the flow around the particles is negligible but particle inertia is dominating due to a high solid-to-fluid density ratio.

Conclusions
In this work we have discussed the influence of a parabolic velocity profile on the dynamics of an inertial prolate spheroidal particle both in terms of translation and rotation. 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 the 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 can 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 is overcoming 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 p . However, since the intermittent tumbling motion is essential to have this type of drift, inertial oblate particles and spheres ( 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 is 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 Lagrangian particle tracking (LPT) methods, the force on the particle is usually calculated using only the instantaneous velocity difference Δ * 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 | | 2 /|Δ * | 1 must hold for this contribution to be negligble. = 16 = p /ℎ = 0.6 = 30 = 30 Tolerance for gmres = 10 −12 Tolerance for ode113 = 10 −8 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 analysis will however break down due to sedimentation if the particles are not density matched. We find that the inertial drift described in this work will never be the dominating cause for migration of heavy particles in a channel flow on earth, but will rather enter as a correction factor to a sedimentation process. that differ between par Chw and par are highlighted blue in (B 1); cf. (5.9). The torque on the particle is the same in the parabolic flow bg ( ) and the paraboloidal flow Chw bg ( ). Chwang (1975) found that a particle with St trans. = 0 moving in the paraboloidal flow Chw bg ( ) will have the velocity Chw ( , CM ) = 1 6 3 2 CM,2 + 2 − 2 − 2 cos 2 1 .

(B 2)
In the parabolic flow bg the velocity is instead given by (5.10). Note that the difference Chw − = 1 6 (1 − 2 ) 1 is independent of and CM .