Effect of inertial lift on a spherical particle suspended in flow through a curved duct

Effect of inertial lift on a spherical particle suspended in flow through a curved duct

Brendan Harding1, Yvonne M. Stokes1, Email address for correspondence: brendan.harding@adelaide.edu.au    Andrea L. Bertozzi2
Abstract

We develop a model of the forces on a spherical particle suspended in flow through a curved duct under the assumption that the particle Reynolds number is small. This extends an asymptotic model of inertial lift force previously developed to study inertial migration in straight ducts. Of particular interest is the existence and location of stable equilibria within the cross-sectional plane towards which particles migrates. The Navier–Stokes equations determine the hydrodynamic forces acting on a particle. A leading order model of the forces within the cross-sectional plane is obtained through the use of a rotating coordinate system and a perturbation expansion in the particle Reynolds number of the disturbance flow. We predict the behaviour of neutrally buoyant particles at low flow rates and examine the variation in focusing position with respect to particle size and bend radius, independent of the flow rate. In this regime, the lateral focusing position of particles approximately collapses with respect to a dimensionless parameter dependent on three length scales, specifically the particle radius, duct height, and duct bend radius. Additionally, a trapezoidal shaped cross-section is considered in order to demonstrate how changes in the cross-section design influence the dynamics of particles.

1 Introduction

Inertial lift force influences the motion of particles suspended in fluid flow through a duct. The effect was first demonstrated in the classical experiment of Segre & Silberberg (1961) where particles suspended in flow through a (straight) cylindrical pipe were observed to migrate towards an annulus approximately times the radius of the pipe. This sparked many studies of the simplified scenario of a spherical particle suspended in Poiseuille flow between two plane parallel walls (e.g. Ho & Leal (1974), Schonberg & Hinch (1989) and Asmolov (1999)). Subsequently, particle migration in a circular pipe was investigated at much larger Reynolds numbers (Re) with studies showing that the radius of the focusing annulus grows with increasing Re (Matas et al. 2004, 2009). In recent years inertial lift has been used in the field of microfluidics for the separation of cells by size and has been extensively studied experimentally (Di Carlo 2009; Martel & Toner 2014; Warkiani et al. 2014; Geislinger & Franke 2014). Much of the behaviour, especially in complex geometries, is still only understood at an empirical level.

In the case of fully enclosed non-circular ducts the inertial lift force is generally much more difficult to approach from an analytical perspective. Several experimental studies of straight ducts with square or rectangular cross-sections have shown that particles typically tend to migrate towards one of four (stable) equilibrium positions located a finite distance from the centre of each wall (Abbas et al. 2014; Amini et al. 2014; Di Carlo et al. 2009; Ciftlik et al. 2013). Recent work by Hood et al. (2015) extended the analysis of Schonberg & Hinch (1989) from the case of Poiseuille flow between two walls to Poiseuille flow through square ducts. They develop a model of the inertial lift force at small particle Reynolds number which is calculated via the Lorentz reciprocal theorem and utilises the finite element method (FEM) to approximate terms in the expansion that are analytically intractable. Their results agree well with experimental data and demonstrate that much of the behaviour can be characterised by effects that occur for small particle Reynolds number. Their approach was also applied to rectangular ducts where it is shown that the relative size of the attraction zone for equilibria near the side walls becomes increasingly small as the aspect ratio increases, and even disappears entirely for large enough particles (Hood 2016). Other studies have considered the behaviour at higher Reynolds numbers, in particular Nakagawa et al. (2015) used an immersed boundary method to approximate the full Navier–Stokes equations and demonstrate that particles may additionally focus near the corners at high flow rates in agreement with experiments of Miura et al. (2014). These studies also clarify that particles are less focused at high flow rates suggesting that lower flow rates are more useful for cell separation and sorting.

This paper extends the work of Hood et al. (2015) to the case of curved ducts. The focusing behaviour of curved and spiral ducts has been largely confined to experimental studies which have shown it to be an effective mechanism to separate particles/cells by size (Martel & Toner 2012, 2013; Ramachandraiah et al. 2014; Russom et al. 2009; Nivedita et al. 2017). There have been some studies utilising direct numerical simulation, see for example Liu et al. (2016) and Ookawara et al. (2006), in which the inertial lift effect is often added as an external force on the particle via an over-simplified model. Curved ducts are well-known to develop a secondary flow consisting of two counter-rotating vortices within the duct cross-section commonly referred to as Dean flow in reference to the early studies of Dean (1927) and Dean & Hurst (1959). The experimental literature demonstrates that size based separation is caused by an interplay between the inertial lift force and the drag forces coming from the secondary fluid flow. Here we carefully analyse the Navier–Stokes equations to develop a model of the forces, acting on a particle, within the cross-sectional plane. The model is then applied to the special case of low flow rate to gain some insight into how the inertial lift force and secondary flow drag interact and lead to focusing positions that differ with respect to particle size. Motivated by the experiments of Guan et al. (2013) and Warkiani et al. (2014) we also consider a curved duct having a trapezoidal cross-section to better understand how the shape of a duct influences the separation of different size particles.

Microfluidic devices used for particle separation often have a spiral design in order to accommodate the required focusing length relative to the bend radius. Although the bend radius in such devices is continuously changing, it has been demonstrated that the slowly changing curvature in such devices has negligible impact on the fluid flow in the sense that it is reasonable to treat the bend radius as being locally constant at each point within the spiral (Harding & Stokes 2018). Therefore, ducts with constant bend radius are considered in this study and, by smoothly interpolating results obtained for several different bend radii (using the same cross-section), we are able to approximate the particle dynamics at any location within a spiral duct.

A leading order model of the relevant cross-sectional forces is developed in Section 2. There are several key steps in the case of a curved duct. The first is the introduction of a rotating reference frame in which the fluid and particle motion is approximately steady. This allows us to neglect time dependence of the cross-sectional forces on the particle, but also introduces inertia and Coriolis forces into the Navier–Stokes equations and force model. The second step is the introduction of disturbance flow variables followed by a non-dimensionalisation of the problem and a perturbation expansion with respect to the particle Reynolds number. This is similar in spirit to the approach of (Hood et al. 2015) for a straight duct. The Lorentz reciprocal theorem is used to re-express first order forces, which include the inertial lift force, as a volume integral depending on the leading order approximation of the disturbance flow. Third, we break up the background flow velocity into its axial and secondary components and use this to further expand the force and torque experienced by the particle into distinct parts. Lastly, a model for the trajectory of a particle is developed based on its terminal velocity resulting from the forces driving its motion within the plane of the duct cross-section.

In Section 3 we apply the general model from the preceding section to the specific case of particle behaviour in the limit of low flow rate. This is an interesting case to consider for two reasons. Firstly it allows us to consider a simplified approximation of the background fluid flow and eliminate dependence on flow velocity so that we may focus on the effect of bend radius and particle size on the focusing behaviour. In particular, by utilising an expansion of the background flow developed in (Harding 2018a) we demonstrate that it is sufficient to take a leading order approximation of the axial and secondary components of the flow. It then becomes clear that the interaction between the secondary flow drag and inertial lift force converges towards a limiting behaviour (for a fixed particle size and duct dimensions) that is different from that obtained for straight ducts. Second, it allows us to validate the model by demonstrating that the focusing behaviour observed in straight ducts is recovered for large bend radii. This provides new insights into how the focusing behaviour in curved ducts differs from that in straight ducts.

In Section 4 we present and discuss the focusing behaviour resulting from our model for ducts having square, rectangular and trapezoidal cross-sections. Curved square ducts are not be particularly good at focusing or separating particles but illustrate some interesting dynamics which we comment on briefly. On the other hand, rectangular and trapezoidal ducts are quite effective at focusing particles and the focusing position depends on the size of the particle and bend radius of the duct. Furthermore, the behaviour approximately collapses to a single curve when plotted with respect to a dimensionless parameter which describes the relative size of the inertial lift and secondary drag forces. These results provide insights into recently published experimental results.

2 Governing equations and perturbation expansion

Here we derive the equations for our leading order force model.

2.1 Governing equations in the lab frame of reference

Figure 1: Curved duct with rectangular cross-section containing a spherical particle located at . The enlarged view of the cross-section around the particle illustrates the origin of the local coordinates at the centre of the duct. The bend radius is with respect to the centre-line of the duct and is of modest size for illustration.

Consider a duct which is curved having a constant bend radius and uniform cross-sectional shape, an example of which is depicted in Figure 1. Let denote a Cartesian coordinate system in the lab reference frame with the duct positioned so that the axis is normal to the plane of the bend. Let the bend radius be the distance from the origin to the centre of the duct cross-section at any axial/angular position. We introduce such that are coordinates that describe the location within the duct cross-section relative to its centre and is the angle around the bend measured from the -axis, see Figure 1. Note that we do not consider the flow near the inlet/outlet. With this description is parametrised as

For example, the point on the centre-line is equivalently described by . For a duct having a rectangular cross-section with width and height the walls of the duct are taken to be located at and . For the non-rectangular ducts considered later we take to be the average height and the origin of the cross-section to be such that for the top and bottom walls are located at respectively. Let be a characteristic minimum width of the duct cross-section noting that the ducts considered in this paper all have .

Consider a steady fluid flow, through the duct, satisfying the incompressible Navier–Stokes equations. Letting denote the fluid velocity and denote the pressure, we have

(2.0a)
(2.0b)
(2.0c)

where and denote the fluid density and viscosity respectively (both assumed to be constant), is the (constant) acceleration due to gravity, denotes the interior of the duct and denotes its boundary/surface (i.e. the walls of the duct). The fluid motion is driven by a steady pressure gradient such that can be decomposed into

where denotes the drop in pressure per unit arc-length along the centreline of the duct. Notice that we have assumed gravity does not influence the flow velocity by introducing the hydrostatic pressure to counteract in (2.0a). For simplicity our study is restricted to the specific case where gravity acts in the direction, that is .

The resulting flow, often referred to as Dean flow, has been well-studied for ducts having circular or rectangular cross-sections (Dean & Hurst 1959; Winters 1987; Yanase et al. 1989; Yamamoto et al. 2004). The pressure gradient drives an axial flow from which the inertia of the fluid around the bend then leads to the development of a secondary flow consisting of two counter rotating vortices within the duct cross-section. (At high Dean numbers there can exist multiple solutions, some having several vortices, however such flow conditions are far beyond those of practical interest in the context of microfluidics). Sufficiently far from the inlet/outlet region, the flow velocity components are independent of the angular coordinate (up to their direction) and can be expressed as

(2.0)

where and . The component is the axial flow velocity whereas the components describe the secondary flow. Throughout this paper and are referred to as the background velocity and pressure respectively. In the context of studying the dynamics of a particle suspended in flow through a curved duct will be treated as being known/prescribed. Details of the specific approximation of used to obtain our results is deferred to Section 3 while the treatment throughout the remainder of this section is quite general.

Consider a single spherical particle within the fluid flow through the duct. There are two ways one could analyse the trajectory of a particle within the flow. If the secondary component of the background flow is relatively strong (i.e. such that inertial lift force will be small compared to fluid drag from the secondary flow) then one might view the particle as behaving like a passive tracer (with respect to the background flow velocity ) which is then perturbed by additional drag (due to the nearby walls) and inertial lift forces. Alternatively, if the secondary component of the background flow is moderate or weak (i.e. such that inertial lift force is similar in magnitude to drag from the secondary fluid motion) one could consider the particle to be travelling along a streamline of the axial component of the background flow (i.e. streamlines of ) and that its motion is perturbed by the drag forces from the secondary fluid motion, increased drag from the walls, and inertial lift forces. We consider applications in which the drag forces from the secondary fluid motion is not significantly larger than the inertial lift force and thus take the latter point of view throughout this paper.

Let denote the radius of the spherical particle, denote its density (which we assume to be uniform), and denote the location of its centre at time . It is convenient to parametrise the particle position in cylindrical coordinates as , see Figure 1. The particle follows an axial flow streamline (from which it will be perturbed) and thus is constant, , and where denotes the spin of the particle (i.e. the angular velocity with respect to its centre). The assumption may seem counter-intuitive, i.e. we assume that is fixed only to calculate non-zero forces in these directions such that they cannot remain fixed. However, this is consistent with the approaches in previous studies of inertial lift forces in uni-directional flows, see for example Schonberg & Hinch (1989) and Hood et al. (2015). Furthermore, it can be shown that any additional effects from particle motion in the directions is much smaller than the inertial lift force and can therefore be neglected (given the other assumptions used here). Without loss of generality we can additionally assume that at time such that (or equivalently ) where here denotes the (linear) velocity of the particle.

The presence of the particle alters the fluid flow and we introduce to denote the fluid velocity and pressure in the presence of the particle. The fluid motion is again assumed to be governed by the incompressible Navier–Stokes equations, that is, in general

(2.0a)
(2.0b)
(2.0c)

where is the fluid domain given the presence of the spherical particle and is the surface of the particle. Note that the moving particle results in unsteady fluid flow in the lab frame of reference. In particular, depends on time because is non-stationary. Thus the in (2.0a) cannot be discarded (as it was for the background flow in (2.0a)). Note that the boundary condition (2.0c) has a constant that denotes the spin (or angular velocity) of the particle (as observed in the lab frame). Far from the particle we expect to converge to the background flow .

The particle motion is driven solely by the force and torque exerted on the particle by the surrounding fluid in addition to the gravitational body force. The former are obtained by integrating the fluid stresses over the particle surface:

(2.0a)
(2.0b)

where is the mass of the particle. Note that we take the normal with respect to the fluid domain and thus on the surface of the spherical particle is the normal vector pointing outwards from the centre of the particle.

2.2 Equations of motion in a rotating frame of reference

Although the fluid motion is not steady in the lab frame,we view it as steady relative to the motion of the particle centre (assumed to be moving with constant angular velocity and fixed ). Therefore, in the particle reference frame (noting that the particle may still have a constant spin) we can eliminate the term.

Consider a reference frame rotating about the axis at the (constant) rate . We introduce the coordinate system which we parametrise by such that

The unit vectors in the rotating frame are related to in the lab reference frame by

or conversely

The angular velocity of the rotating reference frame (relative to the lab frame) is

(2.0)

Let be the centre of the particle in the rotating frame of reference. It is always the case that , and . Consequently, the particle velocity in the rotating reference frame is

The particle spin as viewed from the rotating reference frame is reduced by the solid body rotation of the coordinate system and is therefore denoted

Let denote the velocity of a fluid parcel in the rotating frame which is related to the velocity in the lab frame via

Furthermore, the acceleration of a fluid parcel in the two frames are related via

where denotes the material derivative. Letting denote the fluid pressure in the rotating frame (which is such that and thus ), it follows that the equations of motion for the fluid in the rotating reference frame are

(2.0a)
(2.0b)
(2.0c)
(2.0d)

where denotes the fluid domain with respect to the rotating reference frame. The boundary condition (2.0c) arises from the no slip condition on the walls, rotating with angular velocity . Notice that, expressed in this frame, the fluid flow can be considered to be steady since is stationary and each of are constant. Consequently is zero and can be dropped from (2.0a). Observe that it is important for gravity to acts in the direction since otherwise the relative direction of gravity in the rotating frame would be dependent on .

The force and torque on the particle need to be expressed in terms of the rotating reference frame as well. Let denote the force and torque on the particle in the rotating reference frame, respectively. Then one has

(2.0a)
(2.0b)

where and are the mass and moment of inertia of the particle respectively. Note there is no Coriolis like force since the particle velocity is taken to be zero in the rotating reference frame, nor is there an Euler force since is assumed to be constant. Thus, upon re-expressing (2.1) in terms of , and then combining with (2.2) one obtains

(2.0a)
(2.0b)

Note that provides the centripetal force experienced by the particle due to its rotational motion as it travels through the curved duct.

We also need the background flow in the rotating reference frame. Since the background flow is steady in time and independent of angular position then, letting denote pressure and velocity of the background fluid flow respectively, they are related to those in the stationary/lab frame via

(2.0)

We are now equipped to consider the disturbance flow in the rotating reference frame, that is the difference between flow in the presence of the particle and the flow in the absence of the particle.

2.3 Disturbance flow in the rotating frame

In the rotating reference frame we introduce the disturbance velocity and pressure

respectively. Substituting and into (2.2) one obtains

(2.0a)
(2.0b)
(2.0c)
(2.0d)

where we have additionally used the fact that

  • the background flow satisfies (2.1),

  • for any one has ,

  • using (2.0) it can be shown ,

  • for one has .

Similarly the force and torque can be expressed in terms of the disturbance flow as

(2.0a)
(2.0b)

respectively (noting that ). Because are well-defined on the interior of the particle then one may use the divergence theorem and (2.0a) to obtain

(2.0)

recalling that points in to the centre of the particle. Thus the contribution from the background flow to the force is essentially just the inertia of the fluid that would otherwise take up the volume occupied by the particle. For a neutrally buoyant particle (i.e. ) travelling with a velocity close to that of the axial velocity of the background flow we would then expect to be close to the inertia of the displaced fluid so that the difference between the two is negligible. If the particle has lower or higher density than the fluid then the net effect of the two terms is a force directed towards the inside or outside wall respectively. Similarly, for the background contribution to the torque, it may be shown that

(2.0)

At first glance this may appear like a misuse of the divergence theorem, however it indeed holds by first re-arranging the order of the cross and dot products so that the integrand has the form after which the divergence theorem can be applied and then further re-arrangement and use of (2.0a) leads to (2.0). In any case this term is sufficiently small that it can be neglected for the purpose of estimating the inertial lift force to leading order as will be demonstrated upon non-dimensionalising our equations.

2.4 Non-dimensionalisation of the disturbance flow

Let us now non-dimensionalise the variables and equations of motion (working in the rotating reference frame). Spatial variables are non-dimensionalised using the radius of the particle and velocity variables are non-dimensionalised via the characteristic velocity where denotes a characteristic velocity of the background flow , which we take to be the value of its axial component at . This is essentially the same scale as in Hood et al. (2015). The non-dimensionalisation of the remaining variables follow from these via the standard approach for flows in which viscous stresses are dominant. Thus we introduce dimensionless variables denoted by ,

With these scalings, the dimensionless form of (2.3) is

(2.0a)
(2.0b)
(2.0c)
(2.0d)

where and denote the rescaled domains and

(2.0)

is the particle Reynolds number. Writing the duct Reynolds number, , we have , so that the particle Reynolds number can be small even if the duct Reynolds number is not (note that is necessary for the particle to fit in the duct and is typical).

The forces on the particle in the rotating frame must also be non-dimensionalised. Noting that the dimensional force and torque are

we introduce the dimensionless quantities defined by

Noting that , one has from (2.0a)

(2.0)

and similarly from (2.0b)

(2.0)

At this stage, given and an approximation of (for a desired duct shape/size and flow rate), one could solve the non-linear problem (2.4) and then determine the resulting force and torque on the particle via (2.0) and (2.0) respectively. In the absence of a secondary component of the background flow the inertial lift force only becomes the dominant force once the particle has reached ‘terminal’ velocity and spin in relation to its main axial motion. As such, it is typical to determine such that and the axial component of (i.e. ) are negligible through an iterative process. In the simpler case of flow through a straight duct, the remaining components of provide the inertial lift force which leads to drift/migration within the cross-section. However for curved ducts, the secondary component of the background flow results in being influenced by the drag force from the secondary flow motion. A location in the cross-section where the net force from both effects is zero is an equilibrium to/from which a particle may be attracted/repelled depending on the eigenvalues of the Jacobian of the net force (as a function of the particle position within the cross-section, i.e. ).

In the case of a straight duct, the the inertial lift force can be approximated from a perturbation expansion of the disturbance flow for small (Hood et al. 2015). Through this expansion, the non-linear equation (2.4) can be broken up into a sequence of linear Stokes problems that are more easily solved. Furthermore, the resulting decomposition informs the way in which inertial lift influences the particle motion, particularly in our case where the problem is further complicated by the secondary fluid motion. The remainder of this section describes how the approach of Hood et al. (2015) is adopted to the curved duct geometry, and additionally, how separating the background flow into axial and secondary components allows us to separate and identify the terms that influence particle motion within the cross-section.

2.5 Perturbation expansion of the disturbance flow and forces on the particle

We consider a perturbation expansion of the disturbance flow in powers of the particle Reynolds number:

(2.0)

Note that carets are dropped from the perturbation variables since these will always be taken to be dimensionless quantities. Substituting (2.0) into (2.4) and matching powers of , one obtains the zeroth order equations

(2.0a)
(2.0b)
(2.0c)
(2.0d)

and the first order equations

(2.0a)