Modeling and simulation of particle dynamics in microfluidic channels
Hampus Söderqvist (hasa0038@student.umu.se)
Abstract
Microfluidics is the science of fluid flows in the microscopic domain and has a wide range of applications. Examples in microbiology include us- ing microfluidic devices to study cell adhesion, or to separate cancer cells from healthy cells in a blood sample. Microchannels are also used in stud- ies of colloidal particles, for example. Cells and other types of particles move as expected in the direction of the flow, but they also tend to move away both from the walls and the center in a microchannel. The forces involved in this phenomenon are caused by the velocity gradient of the fluid and a disturbance in the flow field as a particle moves close to a wall.
The former drives particles out towards the walls and the latter gives rise to an opposite force. An example of this effect is that small, spherical particles tend to settle at 0.6 channel radii from the center of cylindrical microchannels at low Reynolds numbers.
The ability to predict the behavior of particles in microfluidics is impor- tant during development and fabrication of microchannels. However, no general analytical models exist for particle migration in channels, which can be a problem if one wants to design a microfluidic channel for a specific purpose. In the present work we show how computational fluid dynamics software can be used for developing a particle dynamics model. We first simulate the interaction between a single particle and the surrounding fluid flow, and then base a force model on these simulations and apply the model to a system of several particles, which is a computationally efficient simulation procedure. We use this method to show the migration behav- ior of spherical particles in several channel types, including channels with square and circular cross sections for which we verify our results against the literature. We observe how particles of different sizes can be sepa- rated in a cylindrical channel and that the separation can be enhanced in a curved channel. In addition to the simulations, a channel prototype is fabricated with a 3D printer and is used to verify the correctness of a flow simulation.
Sammanfattning
Mikrofluidik är läran om flöden i mikroskopisk skala och tillämpas in- om många områden. Mikrofluidikutrustning kan exempelvis användas in- om mikrobiologi för att studera cellers vidhäftningsförmåga, eller för att separera cancerceller från friska celler i ett blodprov. Ett annat använd- ningsområde är vid studier av kolloidala partiklar. Celler och andra sorters partiklar rör sig som väntat i flödesriktningen, men de tenderar även att röra sig bort både från väggarna och från mitten i mikrokanaler. Kraf- terna bakom detta fenomen orsakas dels av skjuvtöjningen i flödet, dels av störningen i flödesfältet kring partiklar när de rör sig nära en vägg.
Skjuvtöjningen driver partiklar ut mot väggarna medan störningen i flö- desfältet ger upphov till en motriktad kraft. Ett exempel på denna effekt är att små, sfäriska partiklar tenderar att stabilisera sig kring 0,6 kanal- radier från mitten av cylindriska mikrokanaler vid låga Reynoldstal.
Det är viktigt att kunna förutsäga hur partiklar beter sig i flöden när mikrokanaler ska tillverkas. Det saknas dock generella analytiska model- ler för partikeldynamik i kanaler, vilket kan vara ett hinder om man vill formge en mikrofluidikkanal för ett särskilt ändamål. I det här arbetet visar vi hur beräkningsströmningsmekanik kan användas för att utveckla en modell för partikeldynamik i mikrofluidik. Vi simulerar först växelver- kan mellan en enskild partikel och den omgivande vätskan, och baserat på dessa simuleringar tar vi sedan fram en kraftmodell som används för att simulera ett system av flera partiklar, vilket är en beräkningsmässigt effektiv simuleringsprocedur. Vi använder denna metod för att kartlägga hur sfäriska partiklar beter sig i flera olika kanaltyper, inklusive kanaler med kvadratiska och cirkulära tvärsnitt, som vi jämför med litteraturen.
Vi observerar hur partiklar av olika storlekar kan separeras i en cylindrisk kanal och att separationen förstärks i en krökt kanal. Utöver simulering- arna tillverkar vi en kanalprototyp med en 3D-skrivare och använder den för att verifiera giltigheten i en flödessimulering.
Contents
1 Introduction 1
1.1 Background . . . . 1
1.2 Purpose of this project . . . . 1
1.2.1 Goals . . . . 2
2 Theory 2 2.1 Navier-Stokes equations . . . . 2
2.2 Reynolds number . . . . 2
2.3 Flow profiles in microchannels . . . . 3
2.3.1 Shear flow . . . . 3
2.3.2 Secondary flow . . . . 4
2.4 Equations of motion for a particle . . . . 4
2.5 Forces on particles in microchannels . . . . 4
2.5.1 Drag force . . . . 5
2.5.2 Wall lift force . . . . 5
2.5.3 Slip-shear lift force . . . . 6
2.5.4 Shear-gradient lift force . . . . 7
2.6 Focusing and separation of particles . . . . 8
3 Method 9 3.1 Computational fluid dynamics . . . 10
3.2 Direct numerical simulations . . . 10
3.2.1 Simulation geometry . . . 11
3.2.2 Physics simulation settings . . . 12
3.2.3 Computational mesh . . . 13
3.2.4 Selecting particle positions . . . 14
3.2.5 Collecting lift force data . . . 15
3.3 Particle trajectory simulations . . . 16
3.4 Channel fabrication and experimental evaluation . . . 17
4 Results 17 4.1 Direct numerical simulations . . . 17
4.1.1 Cylindrical and square channels . . . 17
4.1.2 Rectangular channel . . . 22
4.1.3 Semicircular and semi-elliptical channels . . . 23
4.2 Particle trajectory simulations . . . 25
4.2.1 Cylindrical and square channels . . . 25
4.2.2 Rectangular channel . . . 26
4.2.3 Semicircular and semi-elliptical channels . . . 27
4.2.4 Cylindrical channel with expanding outlet . . . 28
4.2.5 Spiral channels . . . 29
4.3 Experimental verification . . . 31
5 Discussion 32 5.1 Direct numerical simulations . . . 32
5.2 Particle trajectory simulations . . . 33
5.3 Channel fabrication and experimental verification . . . 34
6 Conclusions 35
7 References 37
1 Introduction
1.1 Background
Microfluidics is the science of fluid dynamics in the microscopic domain. It has been studied at least since the early 1800s when L. J. Poiseuille investigated the properties of flows of human blood in narrow tubes [1]. Microfluidics is relevant in many fields of research today, an example being that microchannels can be used for high-throughput sorting of cancer cells and healthy blood cells in a blood sample [2]. In addition, using microchannels for bacterial adhesion studies [3] can reveal how piliated pathogenic bacteria bind to host cells [4,5,6].
These applications are based on the phenomenon that particles that initially are randomly distributed in a narrow channel tend to focus along certain trajectories in the flow, even if no external forces act on the particles, such as electric or magnetic forces. This phenomenon was examined in the famous work by Segre and Silberberg [7, 8] and has since then been thoroughly investigated, both analytically and numerically.
In order to determine the forces that cause the particles to focus along their equilibrium trajectories, computational fluid dynamics (CFD) methods can be used [9,10]. Simulations of this kind are called direct numerical simulations and are computationally expensive, so one typically wants to minimize the number of such simulations and use their results in a more computationally efficient way.
Many studies have been published concerning particle dynamics in straight chan- nels with circular and rectangular cross sections, and attempts have been made at capturing the lift forces in mathematical expressions [9]. These expressions typically involve regression coefficients that depend on the velocity of the fluid and the sizes of the specific channel and particles, and are not particularly use- ful for general flow fields and channel geometries. Particle separation can be achieved in many types of channels, for example channels which expand and con- tract in the flow direction, zigzag channels or spiral channels [11]. Most channels in the literature have rectangular cross sections since these are the easiest chan- nels to fabricate with photolithographic techniques [11, 12, 13]. However, with the recent development of soft lithographic methods such as 3D printing, the channels are not limited to these shapes, which motivates the study of particle dynamics in other channel types.
1.2 Purpose of this project
The general purpose of this project is to gain an understanding of particle dy- namics in microfluidics and to apply this theory in simulations, since simulations can give valuable information when developing and optimizing microfluidic de- vices for particle sorting.
We want to investigate how particle migration in a microchannel depends on the shape and size of the channel and the velocity of the fluid. We aim to draw on previous studies and investigate channel geometries that are not commonly
1.2.1 Goals
The following list outlines the concrete goals of this project.
1. Devise procedures both for lateral force computations and particle trajec- tory simulations.
2. Replicate results from published studies.
3. Design and optimize simulated microchannels with respect to their ability to separate particles of different species.
4. Print the microchannels with a 3D printer and perform experiments to verify the simulations.
2 Theory
2.1 Navier-Stokes equations
We want to study particle dynamics in fluid flows, and in order to do this we must know some things about fluid mechanics. We begin with the Navier-Stokes equations, which describe the motion of an incompressible fluid with viscosity µand density ρ as
∂u
∂t + (u · ∇)u = −1
ρ∇p +µ
ρ∇2u + f , (1)
∇ · u = 0, (2)
where u is the fluid velocity, p is the pressure and f is the acceleration due to external forces [14]. We assume the fluid to be Newtonian having constant viscosity and density.
Fluids cannot permeate the solid walls that are considered in this work, and they cannot have a velocity parallel to the wall at the fluid-wall interface, so we set the no-slip condition u = 0 as a boundary condition on the walls.
2.2 Reynolds number
The flow Reynolds number, Re, is a dimensionless quantity important in fluid mechanics. It characterizes the behavior of a fluid flow and is defined as
Re = ρU l
µ , (3)
where U is the mean flow velocity and l is the characteristic length, which in the context of microchannels is the hydraulic diameter DH [14], defined as
DH=4Ac Pc
, (4)
where Ac and Pc are the area and perimeter of the channel’s cross section, respectively.
2.3 Flow profiles in microchannels
The motion of particles in a fluid flow depends on the flow profile, by which we mean the velocity distribution in the fluid. The behavior of a fluid flow can be divided into two main categories: the laminar regime and the turbulent regime. A flow is called laminar if the fluid moves in parallel layers. A flow with a more chaotic behavior of the pressure and velocity fields is called turbulent.
Laminar flows in cylindrical tubes correspond to Re < 2100 and turbulent to Re > 4000 [14]. The present work is only concerned with flows well within the laminar regime at Reynolds numbers no higher than 200.
A laminar flow in a straight, cylindrical microchannel with no-slip condition on the wall is called a Poiseuille flow and has a parabolic velocity profile [14]:
u(r) = umax
1 − r2
R2
, (5)
where umax is the flow speed at the center of the channel cross section, R is the channel radius and r is the radial distance from the center. A channel of length Land pressure difference ∆p between inlet and outlet has umax= R2∆p/(4µL).
2.3.1 Shear flow
A shear flow is a flow in which different layers of the fluid move at different speeds, see figure 1. The Pouiseuille flow is an example of a shear flow. Shear flows are extremely important for the dynamics of particles in microfluidic chan- nels, which is further discussed in section 2.5.
2.3.2 Secondary flow
An undisturbed, laminar flow through a straight channel has no velocity com- ponent perpendicular to the direction of the channel. In contrast, a fluid flowing through a curved channel experiences a centrifugal pressure gradient which to- gether with the wall friction acts on the fluid near the walls [15]. This causes a secondary flow, so called Dean vortices, to arise in the fluid. These vortices are perpendicular to the primary flow and are important since they affect the motion of particles in the channel.
2.4 Equations of motion for a particle
In order to compute the dynamics of a particle in a flow, we need to know how the flow described by equation (1) and (2) affects the particle. In the absence of external forces, the motion of a particle is caused by the pressure and shear stress of the fluid. Shear stress is a physical quantity of force per unit area:
τ = F
A, (6)
where F is a force acting tangentially on a surface with area A. Pressure acts normally on a surface and is thus a form of normal stress.
The translational motion Up of a particle is governed by Newton’s second law on the following form [9]:
mpdUp dt =
Z
S
(−p1 + τ )n ds, (7)
where mp is the mass of the particle, p is the pressure, 1 is the unit tensor, τ describes the shear stress acting on the particle and S is the surface of the sphere.
The rotational motion Ωp of the particle is described by the corresponding angular moment equation [9]:
d
dt(Ip· Ωp) = Z
S
(x − xcm) × [(−p1 + τ )n] ds, (8)
where Ip is the moment of inertia tensor of the particle and xcm is its center of mass. In the present work we only consider particles of spherical shape with uniform density, for which the moment of inertia is Ip= 2mprp/5, where rp is the radius of the particle [16].
2.5 Forces on particles in microchannels
As we discuss in section 2.4, the motion of a particle in a fluid flow is caused by the stress on the particle. In the context of microfluidics, the forces on a particle
flow, called drag force, and forces perpendicular to the primary flow, called lateral forces or lift forces. The lateral forces are the cause of the equilibrium trajectories discussed in the introduction. The most important lateral forces are the wall-induced force, the slip-shear induced force and the shear-gradient induced force. These forces are discussed in the following sections.
2.5.1 Drag force
An object moving through a fluid is subjected to a force known as drag force, which acts in the direction of the fluid’s motion relative to the object. A well- known model for the drag force on a spherical object at low Reynolds numbers is Stokes’ law [14]:
FD= 6πµrpus, (9)
where rp is the radius of the sphere and us is the velocity of the fluid relative to the sphere, also called slip velocity.
The drag force on an object in a general flow can be described in terms of the dimensionless drag coefficient:
CD≡ FD 1 2ρu2sAp
, (10)
where Apis called the characteristic area of the object and typically is defined as the object’s projected area in the flow direction [14]. Stokes’ law corresponds to the drag coefficient
CD= 24 Res
, (11)
where Res is a Reynolds number related to the slip velocity and is defined as
Res= ρusd
µ , (12)
for a sphere with diameter d. We are mainly interested in small particles for which we assume that the slip velocity is small and that Stokes’ law applies.
2.5.2 Wall lift force
The wall-induced lift force is a lateral force that is significant on particles that are close to the wall of a channel. The interaction between the particle and the flow field in this region induces a force on the particle directed away from the wall. If the particle is small relative to the channel cross section, we may assume that the particle’s motion only is affected by one wall at any given
When a particle moves straight towards the wall, the effect of the wall is to decelerate the particle. We are interested in the behavior when the primary motion of the particle is parallel to the channel wall, as is the case when a fluid flows through the channel. In this case the disturbance in the flow field in the vicinity of the particle gives rise to a force on the particle with a significant component perpendicular to the wall. The magnitude of the wall lift force on a spherical particle in a channel with square cross section is
FW= CWρu2maxd6
H4, (13)
where CW is a coefficient that depends on the particle position in the channel cross section and the Reynolds number, d is the particle diameter and H is the height of the channel cross section [11]. This force results in a lateral migration speed ulat of the particle which near the wall can be approximated to the first order as [12]
ulat= 3
64Resedused, (14)
where used is the sphere’s velocity component parallel to the main flow, called sedimentation speed, and Resed = ρusedd/µ is a Reynolds number related to the sedimentation speed. The right hand side of equation (14) is approximately constant, which means that the particle moves away from the wall at a constant speed near the wall, but ulatquickly tends to zero further away from the wall [12].
2.5.3 Slip-shear lift force
The flow in a channel is typically sheared near a wall as a consequence of the no-slip boundary condition. The particle’s interaction with the wall makes the particle lag behind the surrounding flow, which causes a different relative fluid speed on opposite sides of the particle, see figure 2. This speed difference induces a pressure difference across the particle, and a corresponding force on the particle in the direction towards higher relative flow speed, perpendicular to the main flow. Depending on the flow profile, this force can have a significant component parallel to the wall in the channel cross section, as indicated by arrows in figure 3.
The magnitude of the slip-shear lift force on a spherical particle in a cylindrical channel is [17]
FSS= Kusr2pγ1/2
ν1/2, (15)
where γ is the shear rate, ν = µ/ρ is the kinematic viscosity and K = 81.2 is a numerical constant.
Figure 2: A spherical particle lagging in a flow with constant shear. The arrows indicate the undisturbed flow profile relative to the center of the particle.
Figure 3: A cylindrical channel with arrows indicating the direction parallel to the wall in a cross section of the channel.
2.5.4 Shear-gradient lift force
A particle may experience different fluid speeds on opposite sides even if it does not lag behind the flow. This is the case for a flow with a curved profile, for example a Poiseuille flow, see figure 4a. Like in the case with the slip-shear lift force, this results in a force on the moving particle in the direction towards higher relative flow speed, which typically is the direction towards the nearest wall. The magnitude of the shear-gradient lift force on a spherical particle in a channel with square cross section is
FSG= CSGρu2maxd3
H, (16)
where CSGis a coefficient that depends on the particle position and flow Reynolds number [11]. If the relative flow speed is linear, as in figure 4b, each side of the particle experiences the same flow speed and no lift force is achieved.
(a) (b)
Figure 4: A spherical particle in a flow. The arrows indicate the undisturbed flow profile relative to the center of the moving particle. (a) The relative flow speed is different on opposite sides of the sphere in a flow with varying shear, which induces a lift force on the sphere. (b) The relative flow speed is the same on opposite sides of the sphere in a flow with constant shear.
2.6 Focusing and separation of particles
The lift forces discussed in the previous section cause particles to move slowly in the cross section of the channel, so called lateral migration. The slip-shear lift force is relatively weak compared with the wall and shear-gradient lift forces [12, 13] and is mainly responsible for slow migration in the direction parallel to the channel walls, as indicated in figure 3. Neglecting the slip-shear induced force, we expect the shear-gradient force and the wall force to cancel each other at some points in the channel cross section since their directions are opposite. This is true in particular for cylindrical channels, where all lift forces are perpendicular to the walls due to symmetry, so the slip-shear force does not make the particles migrate along the walls. The equilibrium positions thus form a ring at a certain radius from the center of the channel [8]. Figure 5 illustrates the balance between the wall lift force and the shear-gradient lift force.
Figure 5: Schematic of wall-induced and shear-gradient induced lift forces on particles in a channel. The wall-induced lift dominates near the channel wall and the shear- gradient induced lift dominates near the center of the channel.
The equilibrium positions of particles in the channel cross section depend both on the size of the particle and the size of the channel. An important dimension- less parameter in this context is the blockage ratio κ, which we define as
κ = d
H, (17)
where d is the diameter of the particle and H is the height of the channel cross section. All channel cross sections considered in this work have a reflective symmetry across the vertical axis passing through the center of the cross section, so H is more specifically the channel height along this axis.
The blockage ratio can affect a particle’s equilibrium position and is the param- eter we consider in this work. Particles that focus along different trajectories can be collected in separate outlets placed near their respective trajectories.
This is the basic idea behind particle sorting. The separation distance between particles of different species can be improved by expanding the channel after the particles have reached their equilibrium trajectories, since the expansion of the channel will increase the distance between streamlines in the flow.
Another parameter important for channels with rectangular cross section is the aspect ratio, AR, which is defined as
AR = W
H, (18)
where W is the width of the cross section. The aspect ratio can affect the stability of equilibrium positions [11, 18].
3 Method
In this work, we devise a procedure for the development and evaluation of microfluidic channels. As a first step we specify the parameters for the problem, which are the shape and size of the channel and the particle, as well as the flow rate through the channel. In the next step we simulate the particle dynamics in the channel. This step is divided into two parts for computational efficiency.
The last step of the procedure is to fabricate the channel using a 3D printer and perform experiments in order to verify the correctness of the simulations. The procedure is summarized in the flowchart in figure 6.
Figure 6: The procedure for developing and evaluating microfluidic channels.
3.1 Computational fluid dynamics
We are interested in determining the trajectories of all particles in a channel flow, since this gives us valuable information about how to optimize the design of particle separation devices. We determine the trajectories by simulating the system using computational fluid dynamics tools.
Simulating the dynamics of a particle by including the particle when solving the Navier-Stokes equations, see equation (1), is extremely computationally expen- sive. The strategy in this work is to perform these heavy computations, called direct numerical simulations, for just one particle at numerous positions in the cross section of a straight channel, which gives the lift force on the particle. The lift force data from these computations then serve as a basis for an interpolation model for the lift force, which can be applied to a system of multiple particles in a channel. This channel must have the same cross-sectional shape as in the direct numerical simulations, but can have a more general shape in the flow direction.
The computational tool we use for this is Comsol Multiphysics 5.2a, which is a software package for finite element method (FEM) simulations for many types of partial differential equations, including the Navier-Stokes equations and other equations that are relevant in fluid dynamics. A simulation using FEM means that the geometric domain of the problem is discretized as a grid of smaller elements, called mesh elements, and that a solution to the computational problem is found by searching for functions that satisfy the relevant equations in the mesh.
All simulations in this project are performed with an Intel Core i7-3770 CPU and 24 GB RAM. All simulations use a fully coupled solver and P2 + P1 discretiza- tion, which is a FEM setting in Comsol Multiphysics specifying that physical quantities are approximated with piecewise quadratic functions instead of the default which uses piecewise linear functions. The fluid viscosity and density are those of water, and the fluid is assumed to be incompressible, see equation (2). The laminar flow interface is used since the flow Reynolds numbers are low.
3.2 Direct numerical simulations
In order to simulate the motion of a particle, we need to know the force acting on the particle. We find this force by including the particle, represented by a sphere, in the simulation geometry when solving the Navier-Stokes equations, which means that the interaction between the particle and the fluid is taken into account. Such a simulation is called a direct numerical simulation (DNS), and stands in contrast to a simulation where the fluid flow is computed by solving the Navier-Stokes equations without the particle, and where the particle motion is computed with pre-defined force models as a post-processing step. In the latter case the particles are assumed not to affect the flow.
3.2.1 Simulation geometry
The DNS in the present work is conducted for a spherical particle in a straight channel that is oriented with the flow direction along the x-axis. The y-axis is chosen to point upwards and the z-axis to point rightwards, as seen from the channel inlet. The origin lies on the inlet, with the exact position depending on the channel’s cross-sectional geometry. The geometries investigated in this work are channels with cross sections of square, rectangular, circular, semicircular and semi-elliptical shapes. See figure 7 for a schematic of these channels and their respective coordinate systems.
Figure 7: A schematic diagram of channels with square, rectangular, circular, semi- circular and semi-elliptical cross sections. The crosses on the inlets mark the origins of the coordinate systems.
All geometries are built in Comsol Multiphysics’ CAD environment using the general procedure presented in the list below. The CAD procedure is designed to give a geometry that is easy to discretize in the finite element method. Figure 8 illustrates a cylindrical channel at four instances in the procedure.
(a) (b)
(c) (d)
Figure 8: Building a cylindrical channel. (a) The geometry after step 2. (b) The geometry after step 4. (c) The geometry after step 6. (d) The geometry after step 8.
is essentially a line segment and acts as a "skeleton" for the channel.
3. Sweep the inlet geometry along the polygon. This steps completes the channel.
4. Add a solid sphere in the cross section halfway through the channel and take the difference between the channel and the sphere. This steps results in the domain which the fluid occupies.
5. Enclose the sphere with a cylinder with its axis parallel to the channel axis, and construct two work planes coplanar with the ends of the cylinder.
The cylinder is only used to partition the geometry and is not part of the physics.
6. Construct a new work plane coplanar with the inlet and project the cylin- der onto the work plane with a circle.
7. Partition the channel with the cylinder from step 5 and the circle from step 6.
8. Partition the channel with the work planes from step 5.
The partitioning of the geometry enables more meshing options later on. The sizes of the various parts of the geometry differ between simulations but the channel is generally around 10–20 times longer than its height, and the par- titioning cylinder is generally around 5 times wider and 10 times longer than the sphere diameter. The cylindrical partition will be referred to simply as the partition in the rest of this text, and the part of the geometry between the work planes built in step 5 will be referred to as the middle segment.
3.2.2 Physics simulation settings
The physics settings are specified once the geometry has been constructed. The flow through the channel is specified by applying a laminar inflow on the in- let with a mean velocity corresponding to the Reynolds number chosen for the simulation, see equation (3), and by setting zero pressure on the outlet. This gives a fully developed laminar flow throughout the channel. The laminar in- flow feature is available in the computational fluid dynamics module in Comsol Multiphysics.
In order to compute the lift force that the flow field induces on a particle placed at a certain position in the cross section, we allow the particle to move freely along the channel and to rotate around its center, but we restrict it to a fixed position in the cross section. Letting the particle reach translational and rota- tional equilibrium under these conditions results in a net lateral force on the particle. This means that the particle would migrate laterally if we did not restrict its position in the cross section. This lateral force is the main interest in the DNS performed in the present work.
The particle has reached its stationary configuration when the x-component of equation (7) and all components of equation (8) are zero. These constraints are specified and solved with a global ODE solver in Comsol Multiphysics.
To simplify the simulation, the boundary conditions are specified as seen from a reference frame following the motion of the sphere’s center. We set a moving wall condition on the channel walls, with velocity −Ux, where Uxis the velocity of the particle as seen from the channel’s frame of reference. Since our frame of reference follows the center of the sphere, we only need to specify the rotational velocity of the sphere. This is done with a sliding wall condition where the velocity of the sphere’s surface is specified in terms of its angular frequencies around the x-, y- and z-axes.
Once the simulation has been performed, the lift force on the particle is obtained by integrating the sum of the shear stress and the pressure across the surface of the sphere.
3.2.3 Computational mesh
The procedure for generating the computational mesh is divided into two main steps: one for the partition and one for the rest of the geometry.
The part of the geometry that is outside the cylindrical partition is meshed by first generating a two-dimensional mesh on the inlet and then using Comsol Multiphysics’ sweep function to sweep this mesh throughout the geometry, which results in prism-shaped mesh elements. These elements can be quite elongated without affecting the solution much, since the fluid velocity varies slowly in the direction of the flow in this region. Such elongated elements reduce the total number of elements and thus the computational load.
In order to fill the inside of the partition with standard tetrahedral mesh ele- ments, Comsol’s software requires that we first split the elements on the partition surface into triangles. A two-dimensional mesh is then generated on the surface of the sphere, and free tetrahedra are generated in the bulk between the sphere and the partition surface. The mesh is finalized by adding boundary layers on the sphere and the channel wall.
The mesh resolution required depends on a number of factors, such as the ge- ometry of the channel and the flow Reynolds number. To verify that the mesh used in this study is sufficiently good, a mesh study is performed on a cylindrical channel with Re = 200 and blockage ratio κ = 0.1 by comparing the lift force on a particle positioned at distance 0.3 channel radii from the center for different mesh resolutions. Figure 9 illustrates an example of such a mesh.
Figure 9: The mesh of a channel with circular cross section. Some faces are made invisible so the interior of the mesh can be displayed, including the cylindrical partition and the sphere. The inset is a zoom-in on the region around the sphere.
As the mesh is refined, the lift force converges to slightly above 1.1 nN, see table 1. This is in agreement with the work of Zhu [19]. Table 2 demonstrates how the lift force is affected if certain regions of the mesh are made coarser.
The values in these tables indicate that the lift force mainly is sensitive to the mesh resolution on the sphere, which is expected since the lift force is obtained by integrating the stress over the surface of the sphere.
Table 1: The lift force for five meshes with increasing number of elements (NoE) in regions of the geometry that are defined in section 3.2.1. DoF stands for degrees of freedom.
Mesh Sphere NoE Middle NoE Inlet NoE DoF Lift force (nN)
1 152 3343 68 44852 1.3628
2 252 5814 206 100365 1.1724
3 344 11868 264 148464 1.1099
4 464 20253 312 209554 1.1269
5 616 28893 358 255887 1.1195
Table 2: Mesh 5 has been modified to be coarser on the sphere (mesh A), in the middle segment (mesh B) and on the inlet (mesh C). Mesh D is the same as mesh 5 in table 1 and is included as a reference. Note that a fine mesh on the sphere causes a fine mesh in the rest of the middle segment.
Mesh Sphere NoE Middle NoE Inlet NoE DoF Lift force (nN)
A 152 12668 358 179031 1.5901
B 616 24991 358 226571 1.1210
C 616 22473 104 148400 1.1187
D 616 28893 358 255887 1.1195
3.2.4 Selecting particle positions
We want to know the lateral force on the sphere at any position in the channel cross section. The DNS is therefore performed for many sphere positions dis- tributed throughout the cross section. Assuming that the force field is smooth, a
the DNS points, giving us an interpolation formula for the lateral force to be used in the particle trajectory simulation.
A simple way of selecting sphere positions for the DNS would be to generate a rectangular grid with the same height and width as the channel, and to perform one simulation for each point in the grid. This may work well if the channel has a rectangular cross section, but more care has to be taken for other geometries.
Unless the channel has a rectangular cross section, some grid points will lie outside the channel, or so close to the channel wall that the spheres intersect with the wall. Such input is unnecessary and could potentially cause errors in Comsol Multiphysics.
In order to generate appropriate sphere positions, a polygon of smaller size but the same shape as the channel cross section is used to restrict which positions that are used. Figure 10 illustrates this for the right half of a channel with a semi-elliptical cross section. The inpoly function in Matlab is used to deter- mine which positions are inside the inner polygon. The resulting positions are visualized in Matlab and inspected by eye. The size of the inner polygon is adjusted if needed.
Figure 10: Particles with their midpoints inside the inner polygon (dashed line) are completely within the channel (black line) and are used in the DNS. These particles are colored blue in the figure. The other positions are discarded.
3.2.5 Collecting lift force data
One DNS is performed per particle position, as described in the previous section.
Comsol Multiphysics can be instructed to perform all simulations in sequence, called a batch sweep, and to output the lift force for each particle position in a text file. Finding mesh settings that work well for all particle positions can be difficult, so the solver typically fails to find a solution for some particle positions. If this happens, the mesh settings must be tweaked and a new batch sweep performed for these specific particle positions. Sometimes this needs to be repeated multiple times before all simulations have succeeded. This results in one text file per batch sweep, and Matlab is then used for collecting the data into a single text file.
3.3 Particle trajectory simulations
The particle trajectory simulations are performed by taking the results from the DNS in the form of interpolation models of the lift force and using them to sim- ulate particle dynamics in any channel with the same cross-sectional geometry and flow Reynolds number as in the DNS. The advantages of this approach over a multi-particle DNS are the vast increase in computational efficiency and the possibility for curved channels. As an example, a spiral channel with a circular cross section is illustrated in figure 11.
Figure 11: Spiral channel with circular cross section in Comsol Multiphysics. The diameter of the channel is 50 µm and the diameter of the spiral is around 1 cm.
These simulations are performed with the particle tracing interface in Comsol Multiphysics. This interface uses virtual particles, which means that the par- ticles are represented as points with translational but no rotational degrees of freedom. The lift force data computed in the DNS is imported into Comsol Multiphysics as an interpolation function. This interpolation can be regarded as representing a continuous force field, and it can be applied to the particles as a custom-made force, in addition to the drag force which is a ready-made feature in the interface.
The trajectories are simulated by taking the undisturbed flow field from a sta- tionary solution and then adding particles at the inlet and computing the tra- jectories based on the undisturbed flow field. This means that the particles do not affect the flow field. It is also assumed in these simulations that the particles are so far apart that there is no interaction between the particles. Their density is set to 1050 kg/m3, corresponding to the density of a bacterium cell, which is so close to the density of water that the buoyancy and gravitational force on the particles can be assumed to cancel.
The mesh sweeping function is well suited for these simulations since the geom- etry typically is very long compared to its width, which implies slow changes in flow speed in the flow direction.
3.4 Channel fabrication and experimental evaluation
We fabricate a channel using a soft lithographic method. The channel geometry in Comsol Multiphysics is exported as an STL file, which is imported into a program called Cura that is used to control an Ultimaker 2+ 3D printer. The geometry is printed as a strip of PVA (polyvinyl alcohol) which then is embedded in a mold of PDMS (polydimethylsiloxane) [20]. PVA is soluble in water and can be rinsed out of the mold, leaving a cavity with the shape of the channel.
A flow through the channel is achieved by connecting a water pump to the channel inlet. In the present work we print a zigzag channel with rectangular cross section and compare the flow field with a simulation.
Experiments are performed using an Olympus IX70 inverted microscope nor- mally used for optical tweezers assays, but modified to also run microfluidic experiments [21]. This setup is optimized for measurements requiring a stable environment [22, 23, 24].
4 Results
This section contains results from the direct numerical simulations, the parti- cle tracing simulations and the channel fabrication for a number of different channels.
4.1 Direct numerical simulations
To show the lift forces on the particles, direct numerical simulations are per- formed as described in section 3.2. We begin by presenting the results for the cylindrical and square channels, for which comparisons are made with the lit- erature. We then present the corresponding results for the wide rectangular channel, the semicircular channel and the semi-elliptical channel.
4.1.1 Cylindrical and square channels
The lift force on a spherical particle with diameter d = 5 µm in a cylindrical channel with diameter D = 50 µm, corresponding to blockage ratio κ = 0.1, is presented in figure 12. The undisturbed flow fields are Poiseuille flows, see section 2.3, and are presented in figure 13.
A positive lift force in figure 12 is directed outwards in the radial direction, and causes the particle to migrate outwards in the cross section. Conversely, a negative lift force is directed inwards and moves the particle towards the center of the cross section. A point where the lift force changes sign from positive to negative in the outwards direction is an equilibrium position, since the lift force pushes the particle towards this point from either side. This point is around r/R ≈ 0.7in figure 12, where R = D/2 is the channel radius, which means that particles in this flow will migrate laterally until they form a ring-shaped layer
Figure 12: Lift force in the radial direction of a cylindrical channel as a function of the radial distance r from the center of the channel, for a sphere with diameter d = 5 µm and a channel with diameter D = 50 µm. The horizontal line is zero. The dashed and dotted vertical lines indicate the equilibrium positions for Re = 100 and Re = 200, respectively, based on Zhu [19].
(a) (b)
Figure 13: The x-component of the undisturbed flow field in the cylindrical channel at (a) Re = 100 and (b) Re = 200.
with results in Zhu [19], as is indicated by the vertical lines in figure 12 which represent the stable equilibrium positions obtained in that work.
Figure 14 presents the pressure and viscous stress in the z-direction on a particle with center at (y, z) = (0, 0.4R), at Reynolds number 200. These patters are also observed by Zhu [19]. The total lift force on the particle is 0.91 nN in this case, of which the contribution from the shear stress is 0.20 nN and the rest is due to the pressure gradient across the sphere’s surface. It is seen in figure 14 that the pressure gradient is negative in the z-direction, which is consistent with the positive lift force at r/R = 0.4 in figure 12.
(a) (b)
Figure 14: The stress in the z-direction on a spherical particle with diameter d = 5 µm positioned at (y, z) = (0, 0.4R) in a cylindrical channel with diameter D = 50 µm at Re = 200. (a) The shear stress. (b) The pressure relative to the atmospheric pressure.
Figure 15 presents the lift force on spherical particles of different sizes, corre- sponding to κ = 0.1 and κ = 0.2, at Re = 100 in the same channel as in figure 12.
The larger particle’s equilibrium positions in the channel cross section are closer to the center, and the corresponding phenomenon is observed by Liu et al. [9]
in rectangular channels. This is a very important result as it demonstrates that different particle species can have different equilibrium trajectories in a flow, which opens up the possibility to sort different particles.
The results are also in general agreement with the famous work of Segre and Sil- berberg [7], who observed equilibrium positions at r/R ≈ 0.6 for lower Reynolds numbers (Re < 30). The simulations are thus verified.
Note that the lift force is zero at the center, but this is not a stable equilibrium position since the force is directed outwards everywhere in the neighborhood of this point.
Figure 15: Lift force on spherical particles with κ = 0.1 and κ = 0.2, respectively, in a cylindrical channel with Re = 100. The horizontal line is zero.
The cylindrical channel is the only channel that is symmetric in the radial direction, so all other results will be displayed as vector fields in the channel cross section.
The lift force on a particle in a channel with square cross section is presented in figure 16. The channel has side length H = 50 µm and the particle has diameter d = 5 µm, corresponding to the blockage ratio κ = 0.1. The different kinds of lift forces discussed in section 2.5 are apparent in this figure, with the outwards-directed shear-gradient lift force near the center, the strong inwards- directed wall lift force near the walls and the relatively weak slip-shear force which causes slow migration parallel to the walls between these regions. The flow field for the square channel is presented in figure 17.
Figure 16: Lift force on a spherical particle in the top right quadrant of a square channel. Region 1 (cyan) and region 2 (light blue) are where the shear-gradient and wall-induced lift forces dominate, respectively. The slip-shear lift force is significant in the white region, as can be seen in the inset. The blue circle depicts the sphere and is included for size reference. An arrow corresponds to the lift force on a sphere with its center at the base of the arrow. Note that arrows overlap near the walls.
Figure 17: The x-component of the undisturbed flow field in the square channel.
Figure 18 illustrates streamlines corresponding to the vector field in figure 16, with red circles indicating the equilibrium positions. This is in agreement with di Carlo et al. [10] who also observe four discrete equilibrium positions in a square channel, albeit for slightly different Reynolds number and blockage ratio.
Figure 18: Streamlines for the vector field in figure 16. The red circles indicate equilibrium positions.
4.1.2 Rectangular channel
The lift force on a particle in a rectangular channel with height H = 50 µm and aspect ratio AR = 4 is presented in figure 19, for Re = 100 and blockage ratio κ = 0.1. The lift force is most significant in the y-direction, which is expected since the shear is low in the z-direction, see figure 20. This indicates that particles will reach a stable position in the y-direction quickly, and possibly migrate slowly in the z-direction towards discrete equilibrium positions.
Figure 19: Lift force on a spherical particle in the top right quadrant of a rectangular channel.
Figure 20: The x-component of the undisturbed flow field in the rectangular channel.
4.1.3 Semicircular and semi-elliptical channels
Figure 21 illustrates the lift force on a particle in a semicircular channel with height H = 50 µm, for Re = 100 and κ = 0.1. The overlaid streamlines in figure 22 show that there is no radial force symmetry along the arc, but rather that there are discrete equilibrium positions. There is also an equilibrium position near the middle of the bottom wall, which is expected based on the results for the square channel, see figure 16. The flow field for the semicircular channel is presented in figure 23.
Figure 21: Lift force on a particle in the right half of a channel with semicircular cross section.
Figure 22: Lift force on a particle in the right half of a channel with semicircular cross section. The streamlines give an indication of where the equilibrium positions are.
Figure 23: The x-component of the undisturbed flow field in the semicircular channel.
Figure 24 illustrates the lift force in a semi-elliptical channel with semi-minor axis H = 50 µm and semi-major axis W = 3H = 150 µm, for Re = 50 and κ = 0.1. The lift force is generally weak in the z-direction, which is expected since the shear is small in that direction, see figure 25.
Figure 24: Lift force on a spherical particle in the right half of a channel with semi-elliptical cross section.
Figure 25: The x-component of the undisturbed flow field in the semi-elliptical channel.
4.2 Particle trajectory simulations
To visualize how the particles are focused due to the lift forces, we import the data from the DNS into the particle tracing interface in Comsol Multiphysics and compute the particle trajectories. For simplicity, we display the trajectories by showing where the particles cross selected cross sections in the channels.
4.2.1 Cylindrical and square channels
Particle trajectories are simulated in a 1 cm long cylindrical channel with diam- eter 50 µm, for both Re = 100 and Re = 200, for particles with κ = 0.1 and lift forces corresponding to the DNS values in figure 12. The results are presented in figure 26. The channel length required for focusing the particles along their equilibrium trajectories is about the same for both Reynolds number, but the time required is shorter for Re = 200 since the particles travel the distance faster in this flow due to the higher flow rate. This may be an important design aspect if large quantities of particles must be focused or sorted.
(a)
(b)
Figure 26: Particle focusing in a cylindrical channel at (a) Re = 100 and (b) Re =
Figure 27 demonstrates the separation of particles with κ = 0.1 and κ = 0.2, respectively, in a cylindrical channel with diameter 50 µm at Re = 100, for lift forces corresponding to the DNS values in figure 15. The particles focus along circles, where the inner circle corresponds to the larger particles. Particle separation is thus achieved.
Figure 27: Separation of particles with κ = 0.1 and κ = 0.2, respectively. The arrow indicates the flow direction. The inner circle corresponds to the larger particles. The channel height and length are not to scale, and the particles are also not to scale.
Figure 28 presents the focusing of particles in a square channel with side length H = 50µm, with lift forces corresponding to the data in figure 16. The particles approach the walls quickly and migrate slowly along the walls towards their distinct equilibrium trajectories. The channel length is 2 cm and this is the length needed for complete focusing.
Figure 28: Particle focusing in a square channel. The red dots represent the particles and the arrow indicates the flow direction. The channel height and length are not to scale.
4.2.2 Rectangular channel
Figure 29 illustrates particle focusing in a straight, 1 m long rectangular channel with aspect ratio AR = 4, corresponding to the data in figure 19. Most particles spread out quickly in the vertical direction and migrate slowly along the top and bottom walls. Some particles stay near the left and right walls.
Figure 29: Particle focusing in a rectangular channel. The red dots represent the particles and the arrow indicates the flow direction. The channel height and length are not to scale.
4.2.3 Semicircular and semi-elliptical channels
Figure 30 illustrates particle focusing in a straight, 20 cm long semicircular channel, with lift forces corresponding to the data in figure 21. The particles quickly migrate out from the center and in from the walls, and then slowly migrate towards their equilibrium positions.
Figure 30: Particle focusing in a semicircular channel. The red dots represent the particles and the arrow indicates the flow direction. The channel height and length are not to scale.
The particle focusing in a straight, 1 m long semi-elliptical channel is illustrated in figure 31, corresponding to figure 24. Note the fast focusing in the vertical direction and very slow migration in the horizontal direction.
Figure 31: Particle focusing in a semi-elliptical channel. The red dots represent the particles and the arrow indicates the flow direction. The channel height and length are not to scale.
4.2.4 Cylindrical channel with expanding outlet
An expanding outlet can be used for increasing the separation distance of dif- ferent particle species into a more practical scale. Figure 32 demonstrates a simulation of particles with κ = 0.1 and κ = 0.2, respectively, in a cylindrical channel with diameter 50 µm, and with a conic expansion at the outlet. The particles form two circles, with the inner circle corresponding to the larger par- ticles. A schematic of this is presented in figure 33, with a channel added for collecting particles of one species (dashed lines).
Figure 32: The particles are focused in a narrow channel, with minor separation.
The expanding outlet amplifies the separation. The arrow indicates the flow direction.
The channel height and length are not to scale.
Figure 33: Particles of different sizes are focused in a narrow channel, with minor separation. The expanding outlet amplifies the separation. The dashed lines indicate the possibility of adding a channel that collects the larger particles.
4.2.5 Spiral channels
Curved channels give rise to secondary flows, as discussed in section 2.3.2. Fig- ure 34 illustrates the equilibrium positions at the outlet of a spiral channel with square cross section of height H = 50 µm, for Re = 100 and particles with κ = 0.1, corresponding to the data in figure 16. The particles form rings due to the Dean vortices, which are presented in figure 35.
Figure 34: Particle focusing in a spiral channel with square cross section.
Figure 35: Dean vortices in a spiral channel with square cross section. The left wall is facing outwards from the center of the spiral.
Figure 36 illustrates the corresponding equilibrium positions in a channel with circular cross section of diameter 50 µm, for two particle species with κ = 0.1 and κ = 0.2, as in figure 15. The Dean vortices, see figure 37, cause the smaller particles to form rings, while the larger particles have a distinct equilibrium position. This separation is significantly larger than that achieved in a straight channel, see figure 27.
Figure 36: Separation of particles of different sizes in a spiral channel with circular cross section. The larger particles have a distinct equilibrium position. The particles are not to scale.
Figure 37: Dean vortices in a spiral channel with circular cross section. The left side is facing outwards from the center of the spiral.
4.3 Experimental verification
Figure 38 demonstrates a flow in a section of a 3D printed zigzag channel with rectangular cross section and two inlets: one for red-colored water and the other for blue-colored water. The boundary between the red and blue water is maintained throughout the channel since the flow is laminar, so no mixing occurs between streamlines, and the time scale is too low for diffusion to be significant. The print is defective in the sense that there is a slight indentation along the top wall, which is due to the fabrication method where two strips of PVA are molded together. The channel cross section is approximately 1 mm wide and 0.2 mm high. See figure 39 for a schematic illustration of the channel cross section. Figure 40 illustrates streamlines from a simulation in a channel with similar shape, with and without indentation. The indentation forces the water to take a slightly more diagonal path through the square region shown in the middle of figure 40b. A similar effect is seen in figure 38. This experiment verifies our simulation of the flow.
Figure 38: Colored water in a zigzag channel as seen from above. The arrow indicates the flow direction.
Figure 39: The channel as seen in figure 38.
