Feedback control of inertial microfluidics using axial control forces [original]

Lab on a Chip

PAPER

Cite this: Lab Chip,2014,14,2115

Received 3rd February 2014,

Accepted 24th March 2014

DOI: 10.1039/c4lc00145a

www.rsc.org/loc

Feedback control of inertial microfluidics using

axial control forces†

Christopher Prohm*and Holger Stark

Inertial microfluidics is a promising tool for many lab-on-a-chip applications. Particles in channel flows

with Reynolds numbers above one undergo cross-streamline migration to a discrete set of equilibrium

positions in square and rectangular channel cross sections. This effect has been used extensively for particle

sorting and the analysis of particle properties. Using the lattice Boltzmann method, we determined the

equilibrium positions in square and rectangular cross sections and classify their types of stability for different

Reynolds numbers, particle sizes, and channel aspect ratios. Our findings thereby help to design microfluidic

channels for particle sorting. Furthermore, we demonstrated how an axial control force, which slows

down the particles and shifts the stable equilibrium position towards the channel center. Ultimately, the

particles then stay on the centerline for forces exceeding the threshold value. This effect is sensitive to the

particle size and channel Reynolds number and therefore suggests an efficient method for particle separation.

In combination with a hysteretic feedback scheme, we can even increase the particle throughput.

1. Introduction

In recent years, a number of devices using fluid inertia in

microfluidic setups have been proposed for applications such

as particle steering and sorting or for the whole range of flow

cytometric tasks in biomedical applications. They include cell

counting, cell sorting, and mechanical phenotyping.

1–4

These

devices rely on cross-streamline migration of solute particles

subjected to fluid flow where fluid inertia cannot be neglected

as is commonly done in microfluidics. In this article we have

demonstrated how control forces along the channel axis influ-

ence inertial cross-streamline migration and how feedback

control using axial forces enhances particle throughput.

Segré and Silberberg, who investigated colloidal particles

in circular channels, were the first to attribute cross-streamline

migration to fluid inertia.

They observed that flowing particles

gathered on a circular annulus about halfway between

the channel center and the wall. This effect is connected to an

inertial lift force in the radial direction. It becomes zero right

on the annulus which marks degenerate stable equilibrium

positions in the circular cross section. For microfluidic appli-

cations, channels with a rectangular cross section are used

since they can be fabricated more easily. The reduced symmetry

qualitatively changes the lift force profile and only a discrete

set of equilibrium positions remain.

In square channels,

they are typically found halfway between the channel center and

the centers of the channel walls.

In numerical studies, migra-

tion to positions on the diagonals are also observed.

8,9

In rect-

angular channels, the number of equilibrium positions is

further reduced to two when the aspect ratio strongly deviates

from one.

1,10

The particles all gather in front of the long channel

walls. The exact equilibrium positions are of special importance,

as they ultimately determine how devices function based on

inertial microfluidics.

1–3

Inertial lift forces that drive particles away from the chan-

nel center are caused by the non-zero curvature or the para-

bolic shape of the Poiseuille flow profile.

6,11

Only close to the

channel walls, wall-induced lift forces push particles towards

the center. In channels with a rectangular cross section, the

curvature of the flow profile is strongly modified. Along the

short main axis, the flow profile remains approximately para-

bolic, while along the long main axis it almost assumes the

shape of a plug flow with strongly reduced curvature in the

center when the cross section is strongly elongated.

The

large difference in curvature along the two main axes

modifies the lift force profiles in both directions.

We will

investigate them in more detail in this article.

The method of matched asymptotic expansion allows an

analytic treatment of inertia-induced migration and calcula-

tion of the lift force profiles.

13,14

As the method requires the

particle radius to be much smaller than the channel diame-

ter, it is hardly applicable to microfluidic particle flow, where

this assumption is often violated. Here, numerical approaches

provide further insight. Previous studies in three dimensions

Lab Chip,2014,14, 2115–2123 | 2115This journal is © The Royal Society of Chemistry 2014

Institute of Theoretical Physics, Technische Universität Berlin, Hardenbergstr. 36,

10623 Berlin, Germany. E-mail: Christopher.Prohm@TU-Berlin.de

†Electronic supplementary information (ESI) available: Includes a summary of

the lattice Boltzmann method and the implementation of the Inamuro

immersed boundary method. See DOI: 10.1039/c4lc00145a

Published on 01 April 2014. Downloaded by TU Berlin - Universitaetsbibl on 24/02/2016 13:55:06.

View Article Online

View Journal

| View Issue

2116 |Lab Chip,2014,14, 2115–2123 This journal is © The Royal Society of Chemistry 2014

have used the lattice Boltzmann method,

the finite element

method,

or multi-particle collision dynamics.

Using additional control methods such as optical lattices

or optimal control

can increase the efficiency of microfluidic

devices. In an attractive experiment, Kim and Yoo demonstrated

amethodtofocusparticlestothechannelcenter.

They

applied an electric field along the channel axis to slow down

the particles relative to the external Poiseuille flow, which

induces a Saffmann force towards the channel center.

The

experiments were performed at Reynolds number Re ≈0.05.

We will take up this idea and study, at moderate Reynolds

numbers, how the inertial lift force profile changes under an

axial control force.

A more sophisticated method to operate a system is feed-

back control where the control action depends on the current

state. It is widely used in engineering and everyday life.

microfluidic systems, optical tweezers combined with feed-

back control provide a strategy to measure microscopic forces

in polymers and molecular motors.

21–23

In lab-on-a-chip devices,

several strategies are suggested for sorting particles. They all

monitor particle flow directly and use the recorded signal to

implement feedback-controlled optical manipulation.

24–26

will apply a simple form of feedback control to keep the parti-

cles in the channel center.

In this paper, we use the lattice Boltzmann method to

investigate several aspects of inertial microfluidics. We study

in detail the equilibrium particle positions in microfluidic

channels with square and rectangular cross sections and cat-

egorize their types of stability. In particular, we show how,

for channels with sufficiently elongated cross sections, colloi-

dal particles are constrained to move in a plane. We also

show how the inertial lift force profile is manipulated by

applying an axial control force such that the stable equilib-

rium position gradually moves to the channel center. The

effect strongly depends on the particle size and therefore can

be applied for particle sorting. Finally, using the axial force,

we implement hysteretic feedback control to keep the particle

close to the channel center and demonstrate how this enhances

particle throughput compared to the case of constant forcing.

In the conclusions, we refer to the experiments of Kim and

Yoo

as potential experimental approaches to realize the axial

feedback control.

The article is organized as follows: in sect. 2, we introduce

the microfluidic geometry, explain details of the lattice Boltzmann

implementation, and shortly introduce Langevin dynamics

simulations; our results on equilibrium positions and lift-

force profiles in square and rectangular channels are reported

insect.3;wedemonstratetheinfluence of axial control forces

on the lift force profile in sect. 4; we combine it with feedback

control in sect. 5; and we finish with the conclusions in sect. 6.

2 Methods

In this section, we first introduce the microfluidic system.

We then shortly discuss the lattice Boltzmann method and

refer to the details in the appendix. We introduce the procedure

used to determine the inertial lift forces, and finally the

Langevin dynamics for our feedback-control scheme.

2.1 Microfluidic system

We investigated a microfluidic channel with a rectangular

cross section of height 2h, width 2w, and length Las illus-

trated in Fig. 1. We chose the coordinate system such that

the zaxis coincides with the channel axis and the xand yaxes

define the horizontal and vertical directions in the cross sec-

tion, respectively. The channel center corresponds to x=y=0.

The channel was filled by a Newtonian fluid with density ρ

and kinematic viscosity νand a pressure driven Poiseuille flow

was applied.

The maximum flow velocity u

at the channel

center determines the Reynolds number Re = 2wu

/ν. The

implementation of the Poiseuille flow within the lattice

Boltzmann method will be discussed in the next section.

Inside the channel, we placed a neutrally buoyant colloid

with radius a. It follows the streamlines of the applied

Poiseuille flow with an axial velocity v

close the external

Poiseuille flow velocity. Due to the fluid inertia the colloidal

particle experiences a lateral lift force f

lift

, which leads to

cross-streamline migration. In sects. 4 and 5, we also applied

an additional axial control force f

ctl

to the colloidal particle.

We used periodic boundary conditions along the axial direc-

tion and a channel length L=20ato ensure that the periodic

colloidal images do not interact with each other and thereby

do not influence our results.

2.2 The lattice Boltzmann method

We used the lattice Boltzmann method (LBM) to solve the

Navier–Stokes equations of a Newtonian fluid.

27,28

LBM employs

an ensemble of point particles that perform alternating steps of

free streaming and collisions. The particles are constrained to

move on a cubic lattice with a lattice spacing Δx. This restricts

the particle velocities to a discrete set of vectors 

cisuch that

after each streaming step with duration Δtthe new particle

positions again lie on the lattice. In LBM, one describes the

Fig. 1 A schematic of the microfluidic channel (a) and the xz plane at

y= 0 (b). Further explanations are given in the main text. It is sufficient

to only determine the inertial lift force f

lift

in the red quadrant due to

the symmetry of the rectangular cross section.

Lab on a ChipPaper

Published on 01 April 2014. Downloaded by TU Berlin - Universitaetsbibl on 24/02/2016 13:55:06.

View Article Online

Lab Chip, 2014, 14, 2115–2123 | 2117This journal is © The Royal Society of Chemistry 2014

number of particles at lattice point 

xwith velocity 

ciby the

distribution function fxt

i()

. The first two moments of this

distribution function give the hydrodynamic variables: number

density



() ()



xt f xt

 

and velocity   

uxt cf xt

() ,







Further details of our implementation are found in sect. S1 of

the ESI.†

We implemented the pressure driven Poiseuille flow by

imposing a constant body force 

gon the fluid such that the

fluid velocity 

uxt()used to calculate the equilibrium distri-

bution f

ieq

(see ESI†) is replaced by

 

uug 



(1)

We confirmed as shown in Fig. 2 that this procedure does

indeed reproduce the analytically known Poiseuille flow profile.

We placed a colloid in the Poiseuille flow and studied its

position 

r, velocity 

v, and angular velocity 



. We com-

bined the colloid to the fluid using the Inamuro Immersed

Boundary (IB) method

with “five iterations”. For reference,

we presented a short summary of our implementation in

sect. S2 of the ESI.†

We used the palabos LB code

to implement the LB algo-

rithm. We modified the immersed boundary (IB) algorithm

to correctly account for the periodic boundary conditions

along the channel axis and implemented the colloid dynam-

ics according to eqn (S9) in the ESI.†

Along the channel width, we used a total of 101 lattice

sites including the boundaries. We implemented a cubic sim-

ulation grid and chose the number of lattice sites in the

other two directions accordingly. We chose the maximum

flow velocity in the channel such that the Mach number sat-

isfies Ma = u

max

≤0.1. Finally, we adjusted the kinematic

viscosity by the relaxation time τand thereby fixed the

desired Reynolds number Re = u

max

w/v. When τ>1, we

readjusted the Mach number such that τ= 1, as it has been

shown that the accuracy of the combined LBM-IB methods

greatly decreases for relaxation times larger than one.

32,33

2.3 Determining inertial lift forces

To determine inertial lift forces from LB simulations, we

constrained the colloid to a fixed lateral position by simply

disregarding any colloid motion in the cross-sectional plane.

However, we did not impose any constraints on the axial and

rotational motions.

To speed up our simulations, we initialized the system

with the analytical solution of the rectangular Poiseuille flow

and gave the colloid an initial axial velocity v

= 0.8u

, where

is the flow velocity at the channel center. Going through

transient dynamics, the system relaxes rapidly into a unique

steady state within the first 1000 time steps. We continued the

time-evolution up to the vortex diffusion time T= 0.5w

/vand

determined the inertial lift force by averaging the colloidal force

fluid



Ffrom eqn (S7) in the ESI†over the last 2000 time steps

of the simulations. We demonstrated in previous studies

15,17

that this procedure does indeed reproduce correct lift-force

profiles.

2.4 Langevin dynamics simulations

As demonstrated below, axial control forces influence the

inertial lift-force profiles which we determined in the LB sim-

ulations. We then used these profiles in Langevin dynamics

simulations of the colloidal motion to investigate the poten-

tial benefit of feedback control using axial control forces.

We restricted ourselves to channels with an aspect ratio

w/h= 1/3, which ensures that the colloidal dynamics essen-

tially takes place in the xz plane as discussed in sect. 3.2. As

we will discuss in sect. 4, the inertial lift force f

lift

and the

axial velocity v

depend on the applied axial control force f

ctl

We also included thermal noise to exploit the stability of the

fix points of the colloidal motion under feedback control. Fol-

lowing our work in ref. 17, we only included thermal noise

along the lateral direction, as the axial velocities are much

larger than the lateral ones. Then, the Langevin equations of

motion in the lateral and axial directions are given by



dlift ctl

txfxf t()() (2)

dctl

tzvxf

() (3)

where the white noise force has zero mean, 〈η(t)〉= 0, and its

variance obeys the fluctuation–dissipation theorem, 〈η(t)η(t′)〉=

Tξδ(t−t′). We solved the Langevin equations using the con-

ventional Euler scheme.

The parameters were chosen for a

channel with width 2w=20μm and at temperature T= 300 K.

3. Inertial lift forces for different

channel geometries

In channels with circular cross sections, inertial lift forces

drive colloids to a circular annulus with a radius of about

Fig. 2 Velocity profiles of the Poiseuille flow in the rectangular

channel plotted at x= 0 along the yaxis (the long cross sectional axis)

for different aspect ratios w/h. The solid lines show the analytical form

of the profiles,

while the symbols show results of the LB simulations.

Lab on a Chip Paper

Published on 01 April 2014. Downloaded by TU Berlin - Universitaetsbibl on 24/02/2016 13:55:06.

View Article Online

2118 |Lab Chip,2014,14, 2115–2123 This journal is © The Royal Society of Chemistry 2014

half the channel radius. The axial symmetry is reduced in chan-

nels with square or rectangular cross sections and instead of an

annulus, particles accumulate at a discrete set of stable equi-

librium positions.

In addition, the system also shows unstable

equilibrium positions, where the lift force also vanishes but

particles migrate away from them upon a small disturbance.

In the following two sections, we will investigate the loca-

tion and the stability of the equilibrium positions for different

particle sizes, Reynolds numbers, and channel geometries.

We will discuss in detail how we can tailor colloidal motion

by varying the aspect ratio of the channel cross section. Due

to symmetry, we can restrict our discussion to the upper right

quadrant shown in Fig. 1. In Fig. 3, we show the forces acting

on a particle with radius a/w= 0.4 and the resulting trajecto-

ries at Reynolds number Re = 10 for different channel cross

sections. We will discuss the relevant features first for a chan-

nel with square cross sections and then for general rectangu-

lar cross sections.

3.1 Square channels

In a channel with a square cross section and at Reynolds

number Re = 10, a particle with radius a/w= 0.4 experiences

the inertial force profile shown in Fig. 3(a). The gray lines

indicate possible trajectories followed by the particles, which

are free to migrate. Stable and unstable equilibrium positions

are also indicated. We observed that the migration roughly

occurs in two steps. From the channel center and the chan-

nel walls, strong radial forces drive the particle onto an

almost circular annulus at about r≈0.4w. Since the forces

are strong, migration occurs very rapidly. Then, the particle

slowly migrates along the annulus to its equilibrium position

here situated on the diagonal direction.

Together with the channel center there are in total nine

equilibrium positions or fix points in the channel cross sec-

tion. Four of them are indicated in Fig. 3(a). The channel

center is always unstable and the particle migrates away from

it. There are four fix points along the diagonal axes and four

along the main axes (x,ydirections) of the channel cross

section. We plotted their distances from the center versus

the colloid radius for several Re as shown in Fig. 4 and also

indicated their stability. The fix points along the diagonals are

always positioned further away from the channel center as

there is more space for the particle. Consistent with the previ-

ous results,

6,9,15

we observed how both types of equilibrium

positions move closer towards the channel center with increas-

ing particle size and decreasing Reynolds number. Most impor-

tantly, small particles at high Reynolds numbers have their

stable equilibrium positions on the main axes, while larger

particles at lower Reynolds number move to the equilibrium

positions on the diagonals. This is a new result compared to

previous treatments.

6,8,9

In the literature, equilibrium positions in square channels

have been reported along the main axes,

along the diagonals

for large deformable drops

or on both axes.

In contrast to

ref. 9 we observed that the particles move either to the diago-

nal equilibrium positions or to the fix points on the main

axes but the equilibrium positions are never stable at the

same time as illustrated in Fig. 4. It has been demonstrated

in spiral channels with a trapezoidal cross section

that such

Fig. 3 Inertial lift forces (black arrows) of a colloidal particle with

radius a/w= 0.4 in a pressure driven flow at Re = 10. The forces are

plotted in the upper right quadrant of the microchannel cross section

for aspect ratios w/h= 1 (a), w/h= 1/2 (b) and w/h= 1/3 (c). The gray

lines indicate the trajectories of the colloidal particle as it experiences

the lift forces. Lift forces larger than 0.35ρν

are not shown, as their

positions are indicated by squares. Also indicated are the stable and

unstable equilibrium positions. With “saddle”we denote the

equilibrium positions unstable only along one direction.

Fig. 4 Equilibrium positions d

(distance from the center) in a square

channel are plotted versus the colloid radius for different Re.

Equilibrium positions exist along the main axis in the x,ydirections

(solid lines) and along the diagonal (dashed line). Closed circles indicate

stable equilibrium positions, whereas open circles are unstable.

Lab on a ChipPaper

Published on 01 April 2014. Downloaded by TU Berlin - Universitaetsbibl on 24/02/2016 13:55:06.

View Article Online

Lab Chip,2014,14, 2115–2123 | 2119This journal is © The Royal Society of Chemistry 2014

a sudden change in stability can be used to efficiently sort

particles by size. We noted that while the particle size is fixed by

the specific system under investigation, the Reynolds number

remains a free parameter and can be used to tune the stability

of the equilibrium positions. Stable equilibrium positions on

the diagonals were not observed in the experiments and finite

element simulations performed by Di Carlo et al.

While we

are not aware of any obvious reason how to resolve this discrep-

ancy, we have noted that for channel aspect ratios different

from unity a subtle change from stable diagonal equilibrium

positions towards axial positions occurs as discussed in the

next section. Furthermore, diagonal positions are only stable

for sufficiently large particles.

3.2 Rectangular channels

In the experiments, channels typically with rectangular cross

sections were used since the number of stable equilibrium

positions reduces to two situated on the short main axis.

observed the same behavior in the force profiles in Fig. 3 for

large colloids, while for channels with a square cross section

a particle migrates to its stable position on the diagonal

[Fig. 3(a)], this fix point vanishes with a decreasing aspect

ratio w/hand the stable equilibrium position switches to the

short main axis along the xdirection [Fig. 3(b)]. Further

decreasing the aspect ratio w/h, the saddle fix point on the

yaxis vanishes completely and moves to the center at x=y=0,

where it keeps its stability along the yaxis [Fig. 3(c)]. This has

an important consequence (already exploited by us

) that the

colloid is constrained to the xz plane at y= 0 and its dynamics

becomes two-dimensional. In contrast, in Fig. 3(b), a particle

starting close to the centerline moves out of the y= 0 plane on

its way to the stable equilibrium position at x≈0.4. This is

consistent with the simulations performed by Gossett et al.,

where a similar behavior was observed. We will now elaborate

in more detail on these observations.

We first plotted the lift force along the yaxis, i.e.,atx= 0 for

several aspect ratios w/has shown in Fig. 5. Due to symmetry,

the lift force always points along the ydirection. For the qua-

dratic cross section, w/h= 1, zero lift forces indicate the unsta-

ble fix point in the center (x=y= 0) and the saddle fix point

at y≈0.42, which is unstable in the xdirection. As the chan-

nel cross section elongates along the ydirection with decreas-

ing w/h, the lift force driving the particle away from the

channel center becomes weaker and the saddle fix point shifts

towards the channel wall. Below the width w≈0.45h, the lift

force close to the center becomes negative and the unstable

fix point at y= 0 splits into a saddle fix point (now stable in

the ydirection) and an additional unstable fix point. While

we did not show this situation in Fig. 3, it qualitatively looks

the same as in the complete force profile in Fig. 9(a) for

smaller colloids. This is the onset, where the channel center

becomes stable against motion along the yaxis and the col-

loid is constrained to the xz plane at y= 0. Further decreas-

ing the aspect ratio w/h, the unstable and saddle fix points at

y≠0 merge and vanish completely. Only the saddle fix point

at x=y= 0 remains as illustrated in Fig. 3(c). Finally, we

noted that at w/h≈0.75 the stable equilibrium position in

the channel cross section switches from the diagonal to the

xaxis, which is not observed in Fig. 5.

We summarized the situation in the bifurcation diagram

in Fig. 6, where we have plotted the equilibrium positions on

the yaxis versus the aspect ratio w/h. At sufficiently small w/h

only the saddle fix point at y= 0 exists [Fig. 3(c)]. With

increasing w/h, subcritical pitchfork bifurcation occurs. A

second fix point appears which splits into the saddle and

unstable fix points [Fig. 9(a)]. The latter ultimately merges

with the fix point at y= 0 which becomes unstable [Fig. 3(b)].

This resulting situation is illustrated in the left inset for the

whole cross section and with the stable fix point on the

xaxis. The stable fix point moves to the diagonal at w/h≥0.75.

The regime of the subcritical pitchfork bifurcation is much

more pronounced for smaller particles as illustrated in the

Fig. 5 Lift force f

along the y-direction at x= 0 is plotted versus the

yposition of Re = 10 and a/w= 0.4. The different colors correspond to

different aspect ratios w/h.

Fig. 6 Equilibrium positions along the yaxis are plotted versus the

channel aspect ratio w/hof Re = 10 and colloid radius a/w= 0.4.

The equilibrium positions are categorized as stable, saddle, and

unstable, using also their stability with respect to the xdirection. The

blue insets show typical equilibrium positions in the channel cross

section for w/h<1 (left) and for w/h>1.33 (right).

Lab on a Chip Paper

Published on 01 April 2014. Downloaded by TU Berlin - Universitaetsbibl on 24/02/2016 13:55:06.

View Article Online

2120 |Lab Chip,2014,14, 2115–2123 This journal is © The Royal Society of Chemistry 2014

inset of Fig. 7 of a colloid radius a/w= 0.2. The lift force

profiles in this regime (see Fig. 7) even show that the subcrit-

ical transition to a single saddle fix point does not occur for

small aspect ratios. This might be due to the fact that below

w/h= 0.4 the flow profiles close to the wall are the same.

Looking back to Fig. 6, at w/h>1 (now the short main axis

points along the ydirection), the saddle fix point at non-

zero y

first remains. It becomes stable at an aspect ratio

w/h≈1.33 when the fix point on the diagonal vanishes. The

situation is sketched in the right inset.

The basic features of the inertial lift force can be explained

by considering the unperturbed flow field. In particular, the

lift force depends on the curvature of the flow field.

6,11

Along

the long channel axis (ydirection), we observed the flow veloc-

ity shown in Fig. 2. As the channel height hincreases, the flow

profile in the center flattens considerably and the curvature

strongly increases. The differences in the flow profiles of

w/h= 1 and 0.75 are small, which corresponds to the modest

decreaseinthestrengthoftheliftforceshowninFig.5.At

w/h= 0.5 the pronounced flattening of the flow profile sets in

which marks the occurrence of the subcritical bifurcation

and the strong changes in the lift force profiles shown in

Fig. 5. Finally, we observed that the strength of the lift force

in the x-direction also becomes weaker with decreasing

aspect ratio w/hbut the overall characteristics of the profile

(two fix points) remain the same. Especially for aspect ratios

w/h≤0.5, the lift force the in the x-direction does hardly

change.

For many microfluidic applications, such as cytometry

or particle separation,

it is advantageous to ensure that parti-

cles are constrained to move in a plane so that they can easily

be monitored in the focal plane of a microscope. Here we

observed that as soon as the center position becomes a saddle

point, particles will not leave the center plane any more since

they always experience a force driving them back towards the

plane. In Fig. 8, we plotted the necessary aspect ratio to con-

strain particles to the center plane. It becomes smaller with

decreasing particle size and Reynolds number. For the parti-

cle sizes investigated here it is sufficient to choose w/h<0.4 to

ensure that the system is effectively two-dimensional. In the

experiments, it was observed that with increasing Reynolds

number some particles left the center plane,

which we can-

not explain by our results. In Fig. 2a of ref. 36 particles are in

close distance to each other. We therefore suspect that hydro-

dynamic interactions between particles, which are not included

in our simulations, could be relevant for this effect.

4. Axial control of lift forces

In their experiments, Kim and Yoo applied an axial electric

field which slows down particles relative to the Poiseuille

flow.

As a result, particles are pushed towards the center-

line. The observed migration can be rationalized with the

Saffman force which is an inherent inertial force.

It acts

perpendicular to a shear flow when particles are slowed down

or sped up relative to the fluid flow. In their experiments,

Kim and Yoo considered the flow with channel Reynolds

numbers Re ≈0.05 well below unity. Our idea is to apply this

concept to moderate Re and manipulate the inertial lift force

using the additional Saffman force. We will show that with

the help of an axial control force, we can modify the inertial

lift force profile such that we can steer a particle to almost

any desired position on the xaxis.

For a particle of radius a= 0.2wwe observed without axial

control the cross sectional force profile shown in Fig. 9(a).

As discussed in the previous section we observed that a parti-

cle is pushed towards the y= 0 plane where it remains con-

fined. When we applied an additional axial control force of

ctl

= 2.5ρν

, the force profile changed drastically [Fig. 9(b)].

In particular, the stable equilibrium position at x/w≈0.46

vanishes and the particle is focussed to the channel center

regardless of its initial position.

Fig. 7 Lift force f

along the ydirection at x= 0 is plotted versus the

yposition of Re = 10 and a/w= 0.2. The different colors correspond to

different aspect ratios w/h. The inset shows the equilibrium positions

along the yaxis plotted versus the channel aspect ratio w/h. The

equilibrium positions are categorized as stable, saddle, and unstable,

using also their stability with respect to the xdirection.

Fig. 8 Aspect ratio w/h

at which the center becomes a saddle fix

point and the particles move in the center plane. w/h

is plotted

versus Re for different particle sizes. The symbols are data points from

the simulations.

Lab on a ChipPaper

Published on 01 April 2014. Downloaded by TU Berlin - Universitaetsbibl on 24/02/2016 13:55:06.

View Article Online

Lab Chip,2014,14, 2115–2123 | 2121This journal is © The Royal Society of Chemistry 2014

We considered in the following a channel of aspect

ratio w/h= 1/3 to ensure that the particle is confined to the

y= 0 plane and focused on the lift force along the xdirection.

In Fig. 10 we plotted the inertial lift force profiles of several

axial control forces f

ctl

. For zero control force (green curve),

the typical lift force profile occurs with the unstable equilib-

rium position at the center and the stable position half way

between the channel center and the wall. When we applied

the axial control force in the flow direction so that the parti-

cle is sped up relative to the flow (negative f

ctl

), the lift force

increases and the stable equilibrium position is pushed fur-

ther towards the channel wall. The stable position moves

closer to the center when we slowed down the particle with a

positive control force f

ctl

acting against the flow. The addi-

tional Saffman force decreases the lift force. Ultimately both

fix points merge and the particle position at the center is sta-

bilized when the inertial force profile becomes completely

negative. We also studied the change in the lift force of the

fixed position xand found that it is nearly linear in the

applied control forces. In Fig. 10 we already observe that

the variation of the inertia lift force with the axial control

force is the strongest close to the wall.

In Fig. 11 we plotted the stable equilibrium position

versus f

ctl

of different particle sizes. Starting from its

uncontrolled value, the equilibrium position continuously

shifts to zero as the control force increases. A negative control

force moves x

towards the wall. In general, we observed that

larger particles require larger control forces in order to move

the equilibrium position. To analyze this effect further, we

investigated the minimum control force f

ctl

needed to steer a

particle towards the channel center for different Reynolds

numbers and particle sizes in the inset of Fig. 11. We observed

that f

ctl

strongly increases both with particle size and

Reynolds number. When fitted to the power law, we obtained

ctl

∝Re

1.02

(a/w)

2.60

. This indicates that we can easily exploit

an axial control force for particle sorting by size. For example,

consider two particle types with sizes a/w= 0.2 and a/w= 0.3.

While the small particle is well focussed to the channel cen-

ter with a control force f

ctl

=3ρv

, the larger particle only

changes its equilibrium position by about 10 %.

5. Feedback control

We have seen already in the previous section that a constant

axial control force allows the manipulation and sorting of

particles of different types. In the following, we will demon-

strate how a simple feedback scheme adds additional control

to the system and, in particular, increases particle throughput.

We will present results where we simulated particle motion

using the Langevin dynamics described in sect. 4.

We used a hysteretic control feedback scheme, which

switches from no control to constant control depending on

Fig. 10 Inertial lift force is plotted versus location xfor different axial

control forces f

ctl

of Re = 10 and a/w= 0.2.

Fig. 11 Equilibrium position x

is plotted versus axial control force f

ctl

for different particle radii a/w. The inset shows the minimum control

force f

ctl

needed to focus a particle to the channel center that is

plotted versus Re for different particle sizes.

Fig. 9 Lift force profile of a/w= 0.2, Re = 10, w/h= 1/3 without the

axial control (a) and with axial control force f

ctl

= 2.5ρν

(b).

Lab on a Chip Paper

Published on 01 April 2014. Downloaded by TU Berlin - Universitaetsbibl on 24/02/2016 13:55:06.

View Article Online

the lateral particle position with the goal of keeping the parti-

cle close to the channel center. In concrete, we chose a target

interval [−b,b] for the xposition of the particle. We switched

the axial control force to a constant value f

when the particle

is outside the target interval. The modified lift force profile

drives the particle back to the channel center, and we

switched off the control force until the particle leaves the tar-

get interval again. We show the resulting hysteretic control

cycle of f

ctl

in Fig. 12(a), which either acts in the positive or

negative xdirection. The applied control force changes not

only the lift force profile but also the particle velocity along

the channel axis as we demonstrate in Fig. 12(b). When the

control is active, the particle slowed down compared to the

uncontrolled motion.

Fig. 12(c) shows an example of a particle trajectory under

the feedback scheme. The particle starts outside the target

interval [−b,b] and the lift force modified by the axial control

force pulls it towards the channel centerline. As the particle

reaches the centerline, the control is switched off and the

particle is free to evolve. Since the centerline is an unstable

fix point, the particle can leave the target interval and the

control is activated again.

Performing feedback control instead of using a perma-

nently applied control force has the advantage that the effect

of control is reduced. In particular, feedback control gives an

improved particle throughput while keeping the particle close

to the centerline. In Fig. 13 we showed the mean particle

speed (in units of the maximal flow speed u

) as a function

of the target width bfor several control forces f

. For the

uncontrolled motion, v

max

is the particle speed on the center-

line and ν

is the velocity at the equilibrium position. When

the target interval contains the equilibrium position at x= 0.46,

feedback control is not active. The particle stays at the equi-

librium position and moves with ν

. However, for smaller

target widths feedback control sets in. We observed a mean

particle velocity which is almost independent of the target

width and the applied control force f

. Only at b= 0, which

means a permanently applied control force where the particle

always stays in the center, does the mean velocity decrease.

So, the particle throughput is the highest when feedback con-

trol is active.

6. Conclusion

Inertial microfluidics has proven to be particularly useful in

applications such as particle steering and sorting which are

important tasks in biomedical applications. In this paper, we

provided further theoretical insights into inertial microfluidics

using lattice Boltzmann simulations. We put special emphasis

on controlling the particle motion either by designing the

channel geometry or by applying an additional control force.

We first investigated the equilibrium positions in square

and rectangular channels using lift force profiles and catego-

rized them into stable, saddle, and unstable fix points. In

square channels, the stable fix points either sit on the diago-

nals or the main axis of the cross section. This depends on

particle size and Reynolds number and thereby offers the

possibility to sort particles of different sizes. For rectangular

channels we illustrated bifurcation scenarios of fix points sit-

uated on the long main axis. In particular, we showed that

for sufficiently elongated channel cross section particles are

pushed into a plane. Their dynamics becomes two-dimensional

which simplifies the monitoring; and thereby control of parti-

cle motion in experiments.

We then demonstrated how an additional axial control

force allows tuning of the stable equilibrium position, which

moves towards the center with increasing force. Ultimately

the stable position stays on the centerline when the force

exceeds a threshold value which we identified for different

particle sizes and Reynolds numbers. The strong dependence

on these parameters allows separation of particles by size.

Finally, we proposed a hysteretic feedback scheme using the

axial control force to enhance particle throughput compared

to the case when the control force is constantly applied.

Fig. 12 (a) Schematic of the feedback control scheme. The axial

control force f

is switched on when the particle leaves the target

interval [−b,b|and switched off once the particle reaches the

centerline at x= 0. (b) Particle velocity in the target interval without

and with control. (c) Example of a particle trajectory. The residence

time on the centerline is determined by thermal fluctuations.

Fig. 13 Mean particle velocity 〈v〉in units of the maximum Poiseuille

flow speed u

plotted versus target width b. Different control forces f

within the hysteretic feedback control scheme are used. v

max

is the

particle speed on the centerline and v

is the velocity at the equilibrium

position when no control acts.

Lab on a ChipPaper

Published on 01 April 2014. Downloaded by TU Berlin - Universitaetsbibl on 24/02/2016 13:55:06.

View Article Online

The axial control force can be implemented by applying an

external electric field along the channel axis which induces an

electrophoretic force acting on the particle as demonstrated

by Kim and Yoo.

Since the electric field can easily be tuned,

the method introduced in this article to manipulate the lift

force profile should be realizable in experiments, as well as

to explore the proposed hysteretic feedback scheme.

We plan to extend our work to investigate the collective

colloidal dynamics induced by the hydrodynamic interactions

between particles. Here, the additional axial order such as

particle trains develops.

37,38

Furthermore, it will be challeng-

ing to generalize our control methods to particle suspen-

sions. The theoretical insights developed in this paper and

future work on the collective dynamics will help to generate

novel ideas for devices in biomedical applications based on

inertial microfluidics.

Acknowledgements

We acknowledge support from the Deutsche

Forschungsgemeinschaft in the framework of the Collaborative

Research Center SFB 910.

References

1 S. C. Hur, A. J. Mach and D. Di Carlo, Biomicrofluidics, 2011,

5, 022206.

2 A. J. Mach and D. Di Carlo, Biotechnol. Bioeng., 2010, 107,

302–311.

3 G. Guan, L. Wu, A. A. Bhagat, Z. Li, P. C. Chen, S. Chao,

C. J. Ong and J. Han, Sci. Rep., 2013, 3, 1475.

4 J. S. Dudani, D. R. Gossett, T. Henry and D. Di Carlo,

Lab Chip, 2013, 13, 3728–3734.

5 G. Segré and A. Silberberg, Nature, 1961, 189, 209–210.

6 D. Di Carlo, J. F. Edd, K. J. Humphry, H. A. Stone and

M. Toner, Phys. Rev. Lett., 2009, 102, 094503.

7 D. Di Carlo, Lab Chip, 2009, 9, 3038–3046.

8 Y. Kataoka and T. Inamuro, Philos. Trans. R. Soc., A, 2011,

369, 2528–2536.

9 B. Chun and A. J. C. Ladd, Phys. Fluids, 2006, 18, 031704.

10 J. Zhou and I. Papautsky, Lab Chip, 2013, 1121–1132.

11 J. P. Matas, J. F. Morris and E. Guazzelli, Oil Gas Sci. Technol.,

2004, 59,59–70.

12 H. Bruus, Theoretical microfluidics, Oxford University Press, 2007.

13 E. S. Asmolov, J. Fluid Mech., 1999, 381,63–87.

14 B. P. Ho and L. G. Leal, J. Fluid Mech., 1974, 65, 365–400.

15 C. Prohm, M. Gierlak and H. Stark, Eur. Phys. J. E, 2012, 35,1–10.

16 M. MacDonald, G. Spalding and K. Dholakia, Nature, 2003,

426, 421–424.

17 C. Prohm, F. Tröltzsch and H. Stark, Eur.Phys.J.E, 2013, 36,1–13.

18 W. Y. Kim and J. Y. Yoo, Lab Chip, 2009, 9, 1043–1045.

19 P. G. Saffman, J. Fluid Mech., 1965, 22, 384–400.

20 K. Aström and R. Murray, Feedback Systems: An Introduction

for Scientists and Engineers, Princeton University Press, 2010.

21 A. Jonášand P. Zemánek, Electrophoresis, 2008, 29, 4813–4851.

22 G. J. Wuite, S. B. Smith, M. Young, D. Keller and C. Bustamante,

Nature,2000,404,103–106.

23 M. D. Wang, H. Yin, R. Landick, J. Gelles and S. M. Block,

Biophys. J., 1997, 72, 1335–1346.

24 R. W. Applegate Jr, D. N. Schafer, W. Amir, J. Squier,

T. Vestad, J. Oakey and D. W. Marr, J. Opt. A: Pure Appl. Opt.,

2007, 9, S122.

25 X. Wang, S. Chen, M. Kong, Z. Wang, K. D. Costa, R. A. Li

and D. Sun, Lab Chip, 2011, 11, 3656–3662.

26 M. S. Munson, J. M. Spotts, A. Niemistö, J. Selinummi,

J. G. Kralj, M. L. Salit and A. Ozinsky, Lab Chip, 2010, 10,

2402–2410.

27 B. Dünweg and A. J. Ladd, Advances in Polymer Science,

Springer, Berlin Heidelberg, 2008, pp. 1–78.

28 C. K. Aidun and J. R. Clausen, Annu. Rev. Fluid Mech., 2010,

42, 439–472.

29 X. Shan and H. Chen, Phys. Rev. E: Stat. Phys., Plasmas,

Fluids, Relat. Interdiscip. Top., 1993, 47, 1815–1819.

30 T. Inamuro, Fluid Dynam. Res., 2012, 44, 024001.

31 The Palabos project, 2013, http://www.palabos.org.

32 G. Le and J. Zhang, Phys. Rev. E: Stat., Nonlinear, Soft Matter

Phys., 2009, 79, 026701.

33 T. Krüger, F. Varnik and D. Raabe, Phys. Rev. E: Stat.,

Nonlinear, Soft Matter Phys., 2009, 79, 046704.

34 P. Kloeden and E. Platen, Numerical Solution of Stochastic

Differential Equations, Springer, 2011.

35 D. R. Gossett, H. T. K. Tse, J. S. Dudani, K. Goda, T. A. Woods,

S. W. Graves and D. Di Carlo, Small, 2012, 8, 2757–2764.

36 S. C. Hur, H. T. K. Tse and D. Di Carlo, Lab Chip, 2010, 10,

274–280.

37 W. Lee, H. Amini, H. A. Stone and D. Di Carlo, Proc. Natl.

Acad. Sci. U. S. A., 2010, 107, 22413–22418.

38 K. J. Humphry, P. M. Kulkarni, D. A. Weitz, J. F. Morris and

H. A. Stone, Phys. Fluids, 2010, 22, 081703.

Lab on a Chip Paper

Published on 01 April 2014. Downloaded by TU Berlin - Universitaetsbibl on 24/02/2016 13:55:06.

View Article Online