Design and Simulation of Pinched Flow Fractionation for Blood Cell

Document Sample
Design and Simulation of Pinched Flow Fractionation for Blood Cell Powered By Docstoc
					                            Ba.Sc. Thesis




Design and Simulation of Pinched Flow Fractionation for
                    Blood Cell Separation


                      Karsten B. Andersen, s042515
                        Simon Levinsen, s042570




   Supervisors: Winnie Edith Svendsen, Maria Dimaki & Fridolin Okkels

            MIC – Department of Micro and Nanotechnology
                   Technical University of Denmark


                            June 22, 2007
ii
Abstract

This thesis contributes to the effort of making a Chromosome Total Analysis System
(CTAS) in the NaBiS group at MIC - Department of Micro- and Nanotechnology, Tech-
nincal University of Denmark (DTU). The system requires that white blood cells are
separated from the remaining blood particles. We have developed a theoretical model for
blood cell separation using a Pinched Flow Fractionation (PFF) system. The model is
based on basic microfluidic theory and supported by numerical simulations.
    Separation of particles by PFF has experimentally been demonstrated by Takagi et
al. [3] , Yuushi et al. [4] and Masumi et al.[9]. But the theory presented in these
articles was not able to predict the particle behavior. Therefore a new model concerning
fundamental predictions of various critical parameters such as the alignment demand, and
the separation position in a PFF system are presented in this thesis.
    In the last part of this thesis, a design guide for the fabrication of a separation device
using PPF is included. The guide contains vital optimal design parameters and informa-
tion about how the critical separation radius can be manipulated. The guide can be read
independently of the rest of the thesis.
    Parallel with this project an experimental PFF system was fabricated by the NaBiS
group. But at the date for handing in of this thesis experiment has not yet been performed
on the design.




                                             iii
iv   ABSTRACT
     e
Resum´

Denne afhandling bidrager til forsøget med at lave et Total Kromosom Analyse System
                     a
i NaBiS gruppen p˚ MIC - Institut for Mikro- og Nanoteknologi, Danmarks Tekniske
Universitet (DTU). Systemet kræver, at hvide blodceller adskilles fra resten af blodpar-
tiklerne. Vi har udviklet en teoretisk model for blodcelleadskillelse med brug af Pinched
                                                           a
Flow Fractionation (PFF) system. Modellen er baseret p˚ grundlæggende mikrofluid teori
og understøttet af numeriske simuleringer.
    Adskillelse af partikler med PFF er blevet eksperimentelt demonstreret af Takagi et
al. [3], Yuushi et al [4] og Masumi et al. [9]. Men den teori, som præsenteres i disse
afhandlinger, er ikke i stand til at forudsige partiklernes opførsel. Derfor er en ny model
vedrørende grundlæggende forudsigelser af adskillige kritiske parametre, s˚ adan som align-
ment krav, og seperationspositionen i et PFF system præsenteret i denne afhandling.
    Den sidste del af afhandlingen indeholder en design vejledning til fabrikationen af en
adskillelseschip med brug af PFF. Vejledningen indeholder afgørende optimale designpara-
metre og information om, hvordan den kritiske adskillelsesradius kan manipuleres. Denne
vejledning kan læses uafhængigt af resten af afhandlingen.
    Parallelt med dette projekt blev et eksperimentalt PFF system fabrikeret af NaBiS
                  a
gruppen. Men p˚ den dato, hvor denne afhandling er afleveret, er eksperimentet endnu
                     a
ikke blevet udført p˚ designet.




                                            v
vi   RESUME
Preface

This Bachelor thesis is submitted in order to meet the requirements for obtaining a bachelor
of science degree (Ba Sc.) at the Technical University of Denmark. The projekt has been
made in the Nano Bio Integrated Systems (NaBIS) group at Department of Micro- and
Nanotechnology (MIC), DTU.
    The work has been carried out under supervision of Winnie Edith Svendsen, Maria
Dimaki and Fridolin Okkels in the period from Feburary 2007 to June 2007.


Kongens Lyngby, June 22nd 2007




         Simon Levinsen                                Karsten Brandt Andersen
            s042570                                            s042515


                    MIC - Department of Micro- and Nanotecnology
                           Technical University of Denmark
                               DTU - Building 345 east
                              DK-2800 Kongens Lyngby
                                      Denmark




                                            vii
viii   PREFACE
Contents

List of figures                                                                                                                                          xiv

List of tables                                                                                                                                           xv

List of symbols                                                                                                                                         xvii

1 Introduction                                                                                                                                            1
  1.1 Basic Properties of Pinched Flow Fractionation Systems . . . . . . . . . . .                                                                        1
  1.2 Outline of the Thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .                                                                 3

2 Governing Equations                                                                                                                                     5
  2.1 Equation of Continuity .      .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .     5
  2.2 Navier–Stokes Equation        .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .     6
  2.3 Diffusion Equation . . .       .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .     7
  2.4 Hagen–Poiseuille Law .        .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .     7
  2.5 Boundary Conditions . .       .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .     8

3 Nature of Blood                                                                          9
  3.1 Red Blood Cells . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
  3.2 White Blood Cells . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
  3.3 Buffer Fluid . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10

4 Fluid Motion                                                                                                                                           11
  4.1 Flow Profile . . . . . . . . . . . . .                 .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .    11
  4.2 Reynolds Number . . . . . . . . . .                   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .    14
  4.3 Hydraulic Resistance . . . . . . . .                  .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .    15
  4.4 Finite Channel Height Dependency                      .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .    16
  4.5 Streamlines . . . . . . . . . . . . .                 .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .    16

5 Particle Motion                                                                                                                                        17
  5.1 Drag Force . . . . . . . . . . . . . .                    .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .    17
  5.2 Other Forces Acting on the Particle                       .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .    20
  5.3 Drag Number . . . . . . . . . . . . .                     .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .    23
  5.4 Diffusion Constants . . . . . . . . . .                    .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .    24

                                                        ix
x                                                                                                                                CONTENTS

    5.5   Motion of the Particles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24

6 Analytical Calculations                                                                                                                            25
  6.1 Lateral Penetration Distance .             .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   25
  6.2 Borderline Position . . . . . . .          .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   27
  6.3 Diffusion Process . . . . . . . .           .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   27
  6.4 Theoretical Expectation Values             .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   32
  6.5 Particle Aligning . . . . . . . .          .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   33

7 Numerical Simulation Method                                                                                                                        35
  7.1 Finite Element Method . . . . . . . . . . . .                              .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   35
  7.2 Mesh Structure . . . . . . . . . . . . . . . . .                           .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   36
  7.3 Conclusion . . . . . . . . . . . . . . . . . . .                           .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   39
  7.4 Quasi 3D Simulation Method . . . . . . . . .                               .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   39
      7.4.1 Determination of the Damping Factor                                  .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   40
      7.4.2 Implementing the Quasi 3D Model . .                                  .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   42

8 Fluid Interaction Simulation                                                                                                                       43
  8.1 Model Definition . . . . . . . . . . . . . . .                          .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   43
      8.1.1 Characterization of a Good Result .                              .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   45
  8.2 Physical Expectation . . . . . . . . . . . . .                         .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   45
  8.3 Simulation Method . . . . . . . . . . . . . .                          .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   46
  8.4 Results and Discussion . . . . . . . . . . . .                         .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   47
  8.5 Conclusion and Optimal Design Parameters                               .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   49

9 Separation Mechanism                                                                                                                               53
  9.1 Identical Outlet Channel . . . . . . . . . . .                         .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   54
  9.2 Non Identical Channels . . . . . . . . . . .                           .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   57
      9.2.1 Undesired Pressure Variation . . . .                             .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   58
      9.2.2 Separation by Nonidentical Channels                              .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   60
  9.3 Conclusion and Optimal Design Parameters                               .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   63
  9.4 Comparison with other Models . . . . . . .                             .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   64

10 Consideration Regarding the Final Design                                                                                                          67
   10.1 Motion of Particles Around Corners . . . .                           .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   67
   10.2 Narrowing of the Pinched Segment . . . . .                           .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   69
   10.3 Channel Height . . . . . . . . . . . . . . . .                       .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   70
   10.4 Separating White- from Red Blood Cells . .                           .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   71

11 Designing a Pinched Flow                                                                                                                          73
   11.1 Definition of Parameters . .      .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   73
   11.2 Aligning the Particles . . .     .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   75
   11.3 Separating the Particles . .     .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   76
   11.4 Overall Design Parameters .      .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   77
   11.5 Use of the Design Guide . .      .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   78
CONTENTS                                                                                    xi

12 Conclusion                                                                              81

A 2D Script File                                                                           87

B Q3D Script File                                                                          95

C Separation Mechanism Simulation                                                         101
  C.1 Simulation Definition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101
  C.2 Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102

D Maple Worksheet for Calculating the Hydraulic Resistance                               105
xii   CONTENTS
List of Figures

 1.1   Overview of a PFF system . . . . . . . . . . . . . . . . . . . . . . . . . . . .                  2
 1.2   Channel profile of the SU8 stamp . . . . . . . . . . . . . . . . . . . . . . . .                   3
 1.3   The pinched segment of the final PFF system made in PDMS . . . . . . . .                           3

 3.1   Red blood cell geometry . . . . . . . . . . . . . . . . . . . . . . . . . . . . .                 9

 4.1   Geometry used to determine the parallel plate flow profile . . . . . . . . . . 11
 4.2   Geometry used to determine the flow profile for an rectangular cross section 13
 4.3   Velocity profile in a channel with a rectangular cross section . . . . . . . . . 14

 5.1   Fixed sphere in a steady state flow . . . . . . . . . . . . . . . . . . . . . . . 18
 5.2   Moving sphere in a steady flow . . . . . . . . . . . . . . . . . . . . . . . . . 20

 6.1   Geometry used in analytical calculations . . . . . . . . . . . . .    .   .   .   .   .   .   .   25
 6.2   Borderline position as function of the flow rate ratio . . . . . .     .   .   .   .   .   .   .   28
 6.3   Model for the description of diffusion . . . . . . . . . . . . . . .   .   .   .   .   .   .   .   29
 6.4   Cross section of the analyticaly determied concentration profile       .   .   .   .   .   .   .   31
 6.5   Theoretically determined FWHM as a function of x position . .         .   .   .   .   .   .   .   33

 7.1   FEM domain for a Gaussian-shaped channel . . . . .        . . . . . . . . .       . . . .         36
 7.2   FEM - Basis function . . . . . . . . . . . . . . . . .    . . . . . . . . .       . . . .         36
 7.3   Profile for the mesh structure . . . . . . . . . . . . .   . . . . . . . . .       . . . .         37
 7.4   Intersection width as of function Usf, Mesh A . . .       . . . . . . . . .       . . . .         38
 7.5   Intersection width as of function Usf, Mesh B . . . .     . . . . . . . . .       . . . .         38
 7.6   Intersection width as of function Usf, Mesh C . . . .     . . . . . . . . .       . . . .         38
 7.7   Example of wriggles due to a course mesh . . . . . .      . . . . . . . . .       . . . .         39
 7.8   Schematic illustration of the model for determination     of the damping          factor          40

 8.1   Simplified 2D geometry used in the quasi 3D simulations . . . .        .   .   .   .   .   .   .   44
 8.2   Schematic illustration of some of the quality parameters . . . .      .   .   .   .   .   .   .   45
 8.3   FWHM as a function of the flow rate ratio . . . . . . . . . . . .      .   .   .   .   .   .   .   47
 8.4   FWHM as a function of the flow rate ratio with a course mesh           .   .   .   .   .   .   .   48
 8.5   FWHM as a function of the flow rate ratio with a fin mesh . .           .   .   .   .   .   .   .   48
 8.6   Concentration profile with a flow rate ratio of 0.03 . . . . . . .      .   .   .   .   .   .   .   49
 8.7   Concentration profile with a flow rate ratio of 0.10 . . . . . . .      .   .   .   .   .   .   .   49

                                           xiii
xiv                                                                                  LIST OF FIGURES

      8.8    Corner concentration for various channel widths . . . . . . . . . . . . . . . . 50
      8.9    Intersection width as a function of flow rate ratio for various widths . . . . 50

      9.1    Presentation of the channel number convention . . . .       .   .   .   .   .   .   .   .   .   .   .   .   54
      9.2    Schematic illustration of the separation position . . . .   .   .   .   .   .   .   .   .   .   .   .   .   54
      9.3    Separation position variation due to pressure variation     .   .   .   .   .   .   .   .   .   .   .   .   59
      9.4    Illustration of the three channel length design . . . . .   .   .   .   .   .   .   .   .   .   .   .   .   61
      9.5    Illustration of the Geometry used in [3] . . . . . . . .    .   .   .   .   .   .   .   .   .   .   .   .   65

      10.1 Streamlines around a corner in 3D . . . . . . . . . . . . . . . . . . . . . . . 68
      10.2 Schematic illustration of a narrowing . . . . . . . . . . . . . . . . . . . . . . 69
      10.3 Cross section of a butterfly diffusion profile in 3D . . . . . . . . . . . . . . . 71

      11.1   Basic overview of a pinched flow device . . . . . . . . .        . . . .         .   .   .   .   .   .   .   74
      11.2   Alignment of two particles . . . . . . . . . . . . . . . . .    . . . .         .   .   .   .   .   .   .   76
      11.3   Optimization plot of inlet channel width versus flow rate        ratio .         .   .   .   .   .   .   .   77
      11.4   Geometry used as concrete example . . . . . . . . . . .         . . . .         .   .   .   .   .   .   .   78
      11.5   Determination of the intersection width ratio . . . . . .       . . . .         .   .   .   .   .   .   .   79

      C.1 Geometry 1 for particle separation simulation . . . . . . . . . . . . . . . . . 102
      C.2 Streamlines simulated with Comsol . . . . . . . . . . . . . . . . . . . . . . . 103
List of Tables

 7.1   Parameters to control the mesh . . . . . . . . . . . . . . . . . . . . . . . . . 37
 7.2   Details of the tree mesh proposals . . . . . . . . . . . . . . . . . . . . . . . 38

 8.1   Parameters determining the quality of the result . . . . . . . . . . . . . . . 46

 9.1   Influence on the separation position due to pressure drops . .     . . . . .   .   .   .   59
 9.2   Change in separation position do to pressure drop in different     channels    .   .   .   60
 9.3   Result of theory and simulation for separation position . . . .   . . . . .   .   .   .   63
 9.4   Separation position using article theory and new theory . . .     . . . . .   .   .   .   66




                                          xv
xvi   LIST OF TABLES
List of symbols

   Symbol   Description                    Unit
   ρ        Mass density                   kg /m3
   x        Position vector                m
   v        Fluid velocity vector          m /s1
   u        Particle velocity vector       m /s1
   n        Surface normal vector
   F        Force vector                   N
   f        Force density                  N/m3
   v        Fluid velocity                 m /s1
   u        Particle velocity              m /s1
   r        Particle radius                m
   m        Mass                           Kg
   c        Concentration                  m−3
   σ        Cauchy stress tensor           N /m2
   κB       Boltzmanns constant            1.38 × 10−23 J/K
   p        Pressure                       N /m2
   Rhyd     Hydraulic resistance           Pa s/m3
   η        Viscosity of fluids             Pa s
   Q        Volume flow rate                m3 /s
   Usf      Scaling factor
   L        channel Length                 m
   h        Channel height                 m
   w        Channel width                  m
   w′       Intersection width             m
   T        Temperature                    K
   α        Quasi 3D damping factor        Pa s/m2
   R        Angle                          radian
   R′       Intersection angle             radian
   V0       Characteristic velocity        m/s
   L0       Characteristic length          m
   τ        Characteristic time            s
   N0       Characteristic concentration   m−3
   D        Diffusion constant              m2 /s
   Ω        Rotation velocity              s−1

                              xvii
xviii   LIST OF SYMBOLS
Chapter 1

Introduction

In the treatment of various diseases a chromosome analysis is necessary in order to make
a precise diagnostic. Therefore the chromosomes need to be isolated. Today this is done
by the use of different machines and is very time consuming. In the Nano Bio Integrated
Systems (NaBIS) group at MIC - Department of Micro and Nanotechnology, Technical
University of Denmark (DTU), work is initiated on the development of a device, which
should be able to isolate the chromosomes from a blood sample with the use of a microflu-
idic device. Since the white blood cells are the only chromosome containing particles in
the human blood, these need to be separated from the rest of the blood particles. This
thesis is a theoretical treatment of this separation process, which is build on the concept
which is as known Pinched Flow Fractionation.
    The work on this thesis has been made in close collaboration with another bachelor
group, which is working on the actual fabrication and testing of the device. There will
throughout the thesis be references to the work made by this group.


1.1    Basic Properties of Pinched Flow Fractionation Systems
The concept of particle separation with the use of ”Pinched Flow Fractionation” (PFF)
in micro fluidic systems has been demonstrated by Takagi et al. [3] , Yuushi et al. [4]
and Masumi et al. [9]. An advantage of PFF is that it is possible to separate particles
according to their sizes continuously.
    Generally the system consists of a inlet channel, a main channel, a pinched segment,
a separation chamber and various outlet channels as illustrated in Fig. 1.1. The particles,
which are suspended in a fluid, are introduced via the inlet channel. The main channel
must have a larger flow rate than the inlet channel. This should insure that the particles
in the pinched segment, which is a short narrow channel, are aligned along the upper
sidewall. The vital thing in the pinched segment is the particle alignment. To fabricate
this device various clean-room processes are necessary. The present fabrication method
used by the NaBiS group at MIC, DTU is outlined below.

   • A photomask containing the entire system is processed.

                                            1
      2                                                           CHAPTER 1. INTRODUCTION

                                           Pinched segment             Outlet channels
                    Inlet channel


Introduction of particles




                     Main channel       Separation chamber

      Figure 1.1: This is the simplified geometry of a PFF system. Particles suspended in liquid is
      introduced into the system via the inlet channel. When entering the pinched segment they will be
      aligned to the upper sidewall due to the higher flow rate from the main channel. In the separation
      segment the particles are separated according to their size and thereby divided into various outlet
      channels.



          • By photolithography and curing a negative of the system structures is transferred
            to the photoresist SU8.

          • The SU8 is then used as a stamp. The stamp is covered with Polydimethylsiloxane
            (PDMS) The PDMS is cured and a negative of the stamp structures (positive of the
            system structures) is thereby transferred to the PDMS.

          • The PDMS is finally bonded with a pyrex wafer and is ready to be connected with
            inlet and outlet hoses.

      When this processing method is used, the channels of the systems will have nearly rec-
      tangular cross sectional shape. In Fig. 1.2 a cross section of the channel stamp is shown.
      From this figure it is seen that the cross section of the channels will be nearly rectangular.
      The cross section of the channels will throughout this thesis be assumed to be rectangular.
      The wriggles at the right side of the stamp is only an imaging effect. It is due to, that the
      PDMS is sticky and the needle that is measuring the height gets stuck and thereby does
      not move smooth over the surface.
      In Fig. 1.3 an image of the pinched segment of an actual device is shown. The final sepa-
      ration device, has to be connected with hoses. This is done by pressing a needle through
      the PDMS and afterwards replacing it with a tube.
1.2. OUTLINE OF THE THESIS                                                                                3


Cross depth [nm]




                                                                                     50 µm


                            Cross channel direction[µm]

Figure 1.2:       This is the profile of the SU8               Figure 1.3: This is a part the final separation
stamp used in the production of the channels, for             device made in PDMS. This is a picture of the
a 40µm wide channel. The width of this stamp                  pinched segment. The thin channel is the inlet
is not 40µm. This is because the measurement                  channel and the other the main channel. This
has not been perpendicular to the channel direc-              picture is made by the practical project [10].
tion. The measurement is made by the practical
project N. O. Christiansen and M. G. Hansen
[10]. It is seen that the stamp has a nearly rec-
tangular shape.


1.2                    Outline of the Thesis
In this section the outline of the thesis is shortly summarized.
                   • Chapter 2 : ”Governing Equations”. The basic equation used throughout the thesis
                     is presented.
                   • Chapter 3 : ”Nature of Blood”. A presentation of the characteristic properties of
                     both the blood cells and the biological buffer.
                   • Chapter 4 : ”Fluid Motion”. A derivation of various velocity profiles that are used
                     throughout this thesis.
                   • Chapter 5 : ”Particle Motion”. A consideration of forces that influence on the particle
                     motion in microfluidic channels.
                   • Chapter 6 : ”Analytical Calculations”. In this chapter some theoretical analyses have
                     been carried out and various theoretical expectations are presented. Some design
                     rules that should be used in the design process is also determined.
                   • Chapter 7 : ”Numerical Simulation Method”. The commercial simulation program
                     Comsol is introduced. The so called quasi 3D model is also presented.
                   • Chapter 8 : ”Fluid Interaction Simulation”. In this chapter simulations on the align-
                     ment of particles in the pinched segment have been performed. This leads to the
                     determination of some optimal design parameters.
4                                                      CHAPTER 1. INTRODUCTION

    • Chapter 9 : ”Separation Mechanism”. The actual separation mechanism in the device
      is treated. This is done both by analytical calculations and numerical simulations.
      Optimal design parameters for the separation to work are also presented.

    • Chapter 10 : ”Considerations Regarding the final design”. In this chapter various
      parameters influencing on the final behavior of the device are presented, and the
      optimal sizes of these are determined. In this chapter the optimal design parameters
      for the separation of blood cells is also presented.

    • Chapter 11 : ”Designing a Pinched Flow”. In this chapter the optimal design para-
      meters, which have been presented in the above chapters, are summarized. And a
      number of general design rules are presented.

    • Chapter 12 : ”Conclusion”. In this chapter we have made a conclusion on the thesis.

Chapter 11 ”Designing a Pinched Flow” is included in the thesis as a design guide, which
can be used if one has to design a pinched flow, and does not have any idea on how to do
it. It is made in such a way that one should be able to understand it, even though this is
the only chapter one reads. The chapter can thereby be used in the future work so that
one quickly gets an overview of the optimal design parameters.
Chapter 2

Governing Equations

In this chapter the basic governing equations in microfluidic are presented. These equations
are necessary for the understanding the fluid motion.
    When the expression ”fluid particle” is used later in this thesis, it is not to be under-
stood as an actual particle. Instead it should be understood as a small volume of fluid,
which contains many molecules. Therefore the equation of motion does not describe the
motion of each individual atom, but rather the motion of such a fluid particle. In the
following section four fundamental equations will be described. These are the equation
of continuity, the Navier–Stokes equation, the diffusion equation and finally an equation
relating pressure drop to the total flowrate, which is known as the Hagen–Poiseuille law.



2.1    Equation of Continuity
The continuity equation follows directly from the conservation of mass. The total mass of
a given volume V0 in terms of its density is given by:


                                      m0 =         ρdV                                 (2.1)
                                              V0


The amount of fluid leaving this volume through a small part of the surface ds in an unit
time, is ρvds where v is the fluid velocity. The total amount of fluid leaving this volume
will be given by the following surface integral:


                                      ∆m0 =        ρvds                                (2.2)


The decrease of mass in an unit time can also be expressed via the derivative of the total
mass with respect to time.
                                            ∂
                                  ∆m0 = −         ρdV                                (2.3)
                                            ∂t V0

                                             5
6                                                     CHAPTER 2. GOVERNING EQUATIONS

Because both Eq. (2.2) and Eq. (2.3) express the amount fluid leaving a given volume,
these two equation must be equal. Using the Gaussian theorem one gets:
                                       ∂
                                     −            ρdV =        ρvds ⇔
                                       ∂t    V0
                                       ∂
                                     −            ρdV =        ∇ · (ρv)dV ⇔
                                       ∂t    V0           V0
                               ∂ρ
                                  + ∇ · (ρv)dV =0                                       (2.4)
                          V0   ∂t
For this equation to hold true, for any volume V0 , it is necessary that the integrand equals
zero.

                                   ∂ρ
                                       + ∇ · (ρv) = 0                                 (2.5)
                                   ∂t
This equation is known as the continuity equation. It states that the decrease in the total
mass of fluid volume is equal to the amount of fluid that leaves through its surface.


2.2     Navier–Stokes Equation
The equation of motion in microfluidic systems is called the Navier–Stokes equation. It is
derived from Newton’s second law for a fluid particle:
                                             d
                                         m      v=        Fj ⇔
                                             dt
                                                      j

                                            ρDt v =       fj                            (2.6)
                                                      j

The second equation is divided by the volume, and it is therefore the density, ρ, and the
force density, f , that are used in the equation. Normally when using Newton’s second
law particles have well defined positions. In the transformation from a mowing lagrangian
coordinate system, in which the fluid particle is fixed, to a fixed coordinate system in
which the fluid particle moves, one must replace the time derivative with the material
derivative:
                                           ∂
                                     Dt =     + (v · ∇)                             (2.7)
                                          ∂t
When this expression for the material derivative is used, Eq. (2.6) takes the form:
                                      ∂
                                 ρ       v + (v · ∇)v =               fj                (2.8)
                                      ∂t
                                                                  j

By introducing the body forces acting on the fluid particle, the Navier–Stokes equation
for an incompressible fluid is as stated below.
                                ∂
                           ρ       v + (v · ∇)v = −∇p + η∇2 v                           (2.9)
                                ∂t
2.3. DIFFUSION EQUATION                                                                     7

In this equation p is the pressure in the channels. This term therefore represents the
pressure gradient forces, which is the force acting on the fluid particle due to a difference in
pressure around it. η is the fluid viscosity. The second term of the body forces is therefore
the viscous force, which arrises due to frictional forces between a fluid particle and its
neighboring fluid particles. To arrive at this result it is further assumed, that the fluid is
incompressible which simplifies the equation of continuity Eq. (2.5) into ∇ · (ρv) = 0.
    Along with the continuity equation, the Navier–Stokes equation is used to determine
the velocity fields for a fluid system.
    The equation is a nonlinear partial differential equation which for most cases can not
be solved analytically. But for some ideal cases, it simplifies significantly and makes it
possible to determine an analytical solution. The condition for this simplification will be
discussed in Chapter 4.


2.3     Diffusion Equation
In the study of particle motion in fluids it is necessary to consider the diffusion process.
All particles will experience this effect to a certain amount depending on their diffusivity
in the particular fluid.
    From the continuity equation it is possible to determine an equation which describes
the diffusion process. According to L.D. Landau and E. M. Lifshitz [5] the following
equation holds true when the concentrations are small.
                                      ∂c
                                Dt c =    + v · ∇c = D∇2 c                           (2.10)
                                       ∂t
In this equation c is the concentration, and D is the particular diffusivity, which is inde-
pendent of the concentration at low concentrations. This equation is used to predict the
concentration profile.
    A general measure for the system diffusion relative to the system convection time is
      e
the P´clet number. The value of this number gives an idea of the amount of mixing
                             e
in the system. When the P´clet number is large the diffusion will be small relative to
the convection, which is the fluid velocity driven motion. The number is given as the
characteristic diffusion time divided by the characteristic convection time as stated below.
          e                                           e
At high P´clet numbers convection will and at low P´clet numbers diffusion dominates.
                                                       w2
                                 Diffusion time        V0 w
                           Pe =                 = D =
                                                  w                                    (2.11)
                                Convection time        D
                                                       V0


2.4     Hagen–Poiseuille Law
When a constant pressure drop is present over a given straight channel section the flow
rate Q and the pressure drop ∆p will be related through the following equation, which is
known as the Hagen–Poiseuille law H. Bruus [13].
                                                            ∆p
                             ∆p = Rhyd Q       ⇔     Q=                                (2.12)
                                                            Rhyd
8                                            CHAPTER 2. GOVERNING EQUATIONS

where Rhyd is the hydraulic resistance.
    This relation between the hydraulic resistance, the pressure drop and the flow rate will
be used throughout this thesis to determine both flow rates and hydraulic resistances. The
hydraulic resistance depend on a number of different parameters for instance the channel
length, width, hight and fluid viscosity. The actual hydraulic resistance for particular
channels is determined in Chapter 4.


2.5    Boundary Conditions
In order to solve the Navier–Stokes equation, determine the velocity profile and hence the
flow rate, some boundary condition are needed.
    In fluidic systems this boundary condition is the ”no-slip” condition on the side walls
of the channel. This condition implies that the fluid must be at rest relative to the channel
sidewalls. If the fluid were allowed to slide along the sidewalls it would continuously loose
kinetic energy due to friction and eventually terminate the flow. This ”no-slip” boundary
condition will be used throughout this thesis.
Chapter 3

Nature of Blood

Before the actual theoretical treatment is initiated it will be necessary to determine some
characteristic parameters for the blood cells.
    The effect that will be used to separate the red from the white blood cells, exploits the
fact that the average white blood cell is substantially larger than the average red blood
cell. There are some white blood cells that are smaller than the red blood cells, but on
average they are bigger. In the following section a more detailed description of the two
types of cells, along with a description of the biological buffer in which they are suspended,
will be given.



3.1     Red Blood Cells

A red blood cell does not have a standard sphere shape. It have a more disc-shaped
appearance as shown in Fig. 3.1. The radius of 3.6 µm stated in this figure is the average
radius of red blood cells. There are of course both larger and smaller cells. the thickest
part of a red blood cell is on average 2.2 µm according to J. N. Mazunder [1].




                                                                             tr = 2.2 µm

                                           r = 3.6 µm


Figure 3.1: The overall structure of a red blood cell. It is seen that it have a disc shaped
apearence with a thickness of tr = 2.2 µm. Its radius is r = 3.6 µm.



                                             9
10                                                CHAPTER 3. NATURE OF BLOOD

3.2    White Blood Cells
The white blood cells tends to have a sphere like appearance. The average diameter of
these cells is around 10 µm, they can be up to 20 µm [1].
    The average drained density of both white and red blood cells is ρpart = 1125 kg/m3
according to G. Elart et al. [11]. Even though both red and white blood cells are not
solid, we will assume that the forces exerted on them in the device not will be sufficient
to change their shapes.


3.3    Buffer Fluid
The red and withe blood cells which will be going to be separated is suspended in a
standard buffer called Phosphate Buffered Saline (PBS). This buffer is unlike normal blood
newtonian, and therefore it will have a viscosity which is independent of for instance the
local shear of the fluid. The fact that this buffer is newtonian simplifies the calculation
drastically. This buffer consist mainly of saltwater, and it will therefore be assumed
that the density and viscosity of the buffer is around that of water, i.e. a density of
ρ = 1 × 103 kg/m3 and a viscosity of η = 1 × 10−3 Pa · s will be used. On the data sheet
for this buffer the density and viscosity is denoted as not available.
Chapter 4

Fluid Motion

In determination of fluid motion the central equation that has to be solved is the Navier–
Stokes equation. Solving this equation is only possible in some ideal cases. In the following
sections the flow profile for two type of channels will be derived, and the Reynolds number
will be introduced.


4.1     Flow Profile
To determine the actual flow profiles one has to solve the Navier–Stokes equation shown
below. As stated in Chapter 2, it is rather complicated to solve this equation analytically.

                                 ∂
                             ρ      v + (v · ∇)v = −∇p + η∇2 v                          (4.1)
                                 ∂t
When certain ideal condition are met this nonlinear equation simplifies into a linear equa-
tion. When steady state is considered the time derivative will vanish. If the geometry
further more has translation invariance then the nonlinear term vanishes. It is then pos-
sible to solve the equation analytically.
    Some of the analytical calculations will be made in 2D, i.e. with infinitely thick chan-
nels. Therefore it is necessary to derive the flow profile for such channels. In Fig. 4.1 the
geometry used in this derivation is shown. Due to translation invariance and steady state,
the Navier–Stokes equation simplifies substantially. The appropriate boundary condition


                                                   y
              w    p0 + ∆p                             x                     p0



                                               L
Figure 4.1: Illustration of the parallel plate geometry used to determine the 2D flow profile.
Due to the pressure difference the fluid is forced in the positive x direction.

                                             11
12                                                           CHAPTER 4. FLUID MOTION

is as stated below:
                                       ∂2            ∆p
                                           v (z) = −
                                          2 x
                                                                                         (4.2)
                                       ∂z            ηL

                                           vx (0) = 0                                  (4.3a)
                                          vx (w) = 0                                   (4.3b)

Together with the boundary condition stated in Eq. (4.3a) and Eq. (4.3b) which represents
the ”no-slip” condition, this equation is easily solved to yield the velocity profile stated in
Eq. (4.4).
                                             ∆p
                                   vx (z) =      (w − z)z                                (4.4)
                                            2ηL
It is seen that the velocity profile has a parabolic appearance as a function of the distance
to the sidewalls.
     In Chapter 1 it was mentioned that all the channels in the final device will have a
nearly rectangular cross section. It was also argued that the deviation from rectangular
cross section could be neglected. Therefore it is necessary to determine the flow profile for
channels having this particular cross section.
     In the final device the main part of the channels will have a constant cross section and
thereby possess the desired translation invariance. Therefore it is possible to determine
an analytical expression for the flow profile of the fluid in these channels. Steady state
will be assumed in finding the solution to the Navier–Stokes equation.

                                          ∇p = η∇2 v                                     (4.5)

To determine the flow profile one has to solve the partial differential equation along with
the appropriate boundary conditions stated below. The given cross sectional geometry is
shown in Fig. 4.2.

                          ∇p = η∇2 v                                                   (4.6a)
                                                    1
                      vx (y, z) = 0,      on y =   ± w,     z = 0 or z = h             (4.6b)
                                                    2

 When a channel has a constant pressure drop across it the left hand side of Eq. (4.6a) is
easily solved. This is because the pressure distribution in this kind of channel is given by:

                                                     ∆p
                                       p(x) = p0 −      x                                (4.7)
                                                     L

Due to the translation invariance in the x direction this derivative vanishes from the ∇2
term. Since we also know that the flow only has a non zero component in the x direction
Eq. (4.6a) simplifies into:

                                   ∆p      ∂2   ∂2
                               −      =        + 2 vx (y, z)                             (4.8)
                                   ηL      ∂y 2 ∂z
4.1. FLOW PROFILE                                                                                             13

                                                              z

                                                         h




                                                                                     y
                                       −1w
                                        2
                                                                            1
                                                                            2w


Figure 4.2: This is an illustration of the cross sectional geometry which is used in the determi-
nation of the velocity profile.


To solve this problem both sides of the above equation is expanded in a sine series along
the z direction. The left side and the velocity profile can be expressed as:
                                                              ∞
                                       ∆p    ∆p 4                     1        z
                                   −      =−                            sin nπ                             (4.9a)
                                       ηL    ηL π                     n        h
                                                             n, odd
                                               ∞
                                                                           z
                                 vx (y, z) =         fn (y) sin nπ                                         (4.9b)
                                                                           h
                                               n=1

The arguments in the sine functions are due to the boundary conditions. Insertion of these
functions in Eq. (4.8) yields:

                        ∞                                                      ∞
              ∆p 4              1        z            ∂2   ∂2                                      z
            −                     sin nπ   =              + 2                      fn (y) sin nπ       ⇔
              ηL π              n        h            ∂y 2 ∂z                                      h
                       n, odd                                               n=1
                         ∞                           ∞
                ∆p 4            1        z                    ′′           n2 π 2               z
            −                     sin nπ   =                 fn (y) −             fn (y) sin nπ            (4.10)
                ηL π            n        h                                  h2                  h
                       n, odd                      n=1


It is easily seen that in order for this equation to hold true fn (y) must vanish, when
n is even, because only the odd terms in the sum are included on left hand side. To
determine the expansion coefficient fn (y) for odd n, one must solve the following second
order differential equations.

                                4∆p 1    ′′     n2 π 2
                           −          = fn (y) − 2 fn y,                       for n odd                   (4.11)
                                πηL n            h

                                                                                         1
Those equations together with the boundary condition fn (± w) = 0, can be solved to
                                                          2
give:
                                              y
                        4h2 ∆p 1       cosh(nπ )
                fn (y) = 3         1−         h    ,      for n odd           (4.12)
                                              w
                         π ηL n3      cosh(nπ )
                                                                      2h
14                                                          CHAPTER 4. FLUID MOTION




                                                                                 Velocity bar [arbitrary units]
                                              Height[µm]
                                       Width [µm]

Figure 4.3: This is an illustration of the velocity profile for a channel of width 40 µm and height
30 µm. It is further assumed that the length of the channel is 100 µm and there is a pressure drop
of 50Pa across it.



When this is inserted into Eq. (4.9b) the velocity profile for a channel with a rectangular
cross section is found to be:

                                                              y
                             4h2 ∆p
                                       ∞
                                            1     cosh(nπ )                  z
                  vx (y, z) = 3                1−        h          sin nπ                                        (4.13)
                                                         w
                              π ηL          n3    cosh(nπ )                  h
                                      n,odd                  2h



To illustrate how this profile will look in a real channel it is shown in Fig. 4.3 for a channel
with a width of 40 µm and a height of 30 µm. When the flow is a pressure driven steady
state flow in channels, it is also known as a poiseuille flow. This is the type of flow that
will be in the final device.




4.2     Reynolds Number

There is another case in which it is possible to solve the Navier–Stokes equation. In order
to verify this, Eq. (4.1) is made dimensionless by introducing a dimensionless number
called the Reynolds number. This number is introduced with the use of the characteristic
4.3. HYDRAULIC RESISTANCE                                                                15

length, L0 , pressure, P0 , and velocity, V0 .

                               ∂
                            ρ     v + (v · ∇)v      = −∇p + η∇2 v ⇔
                               ∂t
                       V0 ∂        V2                    P0     V0 η
                     ρ        v + 0 (v · ∇)v        =−      ∇p + 2 ∇2 v ⇔
                       T0 ∂ t      L0                    L0     L0
                      ρV0 L0 ∂
                                  v + (v · ∇)v      = −∇p + ∇2 v ⇔
                        η      ∂t
                               ∂
                          Re      v + (v · ∇)v      = −∇p + ∇2 v                     (4.14)
                               ∂t

From this equation it is seen that when the Reynolds number (Re) is much smaller than
unity, the time derivative and the nonlinear term can be neglected. The Reynolds number
is a measure of the size of the inertia term compared to the viscous term. For large
Reynolds numbers the inertia term dominates, and for small Reynolds numbers the viscous
term dominates.
                                     |Inertia term|   ρV0 L0
                              Re =                  =                             (4.15)
                                    |Viscous term|      η
Due to the small length scale and velocities in microfluidic systems this is often well below
unity. When the nonlinear term vanishes the flow is also said to be a laminar flow. A
laminar flow is characterized by the absence of flow line rolls. When channels have more
than one length scale the characteristic length scale, L0 , used in Eq. (4.15) should be the
smallest length in the channels.
   In the device considered in this thesis the Reynolds number will be:

                         103 kg/m3 × 10−5 m/s × 20 × 10−6 m
                  Re =                                      ≈ 2 × 10−4               (4.16)
                                      10−3 Pa s
Since the Reynolds number is much smaller than unity the nonlinear can be neglected
even though translation invariance is not presence.


4.3     Hydraulic Resistance
As described in Chapter 2 there is a connection between the pressure drop and the flow
rate in a long straight channel. They are connected via the socalled hydraulic resistance
in the Hagen–Poiseuille law.
                                      ∆p = Rhyd Q                                  (4.17)

With the use of this equation it is straight forward to determine the structures needed in
order to achieve a desired flowrate. The flowrate of a rectangular channel is derived in
Chapter 9 to be:
                                             ∞
                                h3 w∆p              1 192h           w
                        Q=             1−                    tanh nπ                 (4.18)
                                 12ηL               n5 π 5 w         2h
                                            n,odd
16                                                         CHAPTER 4. FLUID MOTION

When this is compared with Eq. (4.17) it is easily seen that the hydraulic resistance for a
channel having a rectangular shaped cross section is:
                                          ∞                           −1
                             12ηL               1 192h           w
                     Rhyd   = 3   1−              5 π5 w
                                                         tanh nπ                       (4.19)
                              h w               n                2h
                                        n,odd

It is interesting to note that the hydraulic resistance is directly proportional to the length
of the channel. This will be used later.


4.4     Finite Channel Height Dependency
In the final design the height of the channel will be constant. When channels with a finite
height are used instead of infinitely high channels the fluid velocities will be lowered. In
Chapter 7 it is described how this knowledge can be used to make a 2D simulation which
would yield the same results as a 3D model. The damping effect can also be seen in
Fig. 4.3. In this figure the velocity profile is shown and it is obvious that it has a parabolic
appearance in the height direction as well. Therefore the finite channel height acts as a
damping force an the total flow rate.


4.5     Streamlines
A small volume of fluid that flows in a system from a certain point will follow a specific path
through the system. This fluid path is called a streamline. In order for two streamline to
cross each other a specific point in the geometry must have two different velocity gradients.
This is of course not physically possible. Therefore streamlines will newer cross. When
the flow have reached steady state the streamlines will remain constant.
    These streamlines can in certain cases be used to track actual particles through the
system. Requirements for this is the case is discussed in Chapter 5.
Chapter 5

Particle Motion

In this chapter the motion of particles suspended in a fluid will be investigated. In order
to understand the motion of particles in a fluidic system it is essential to understand the
forces acting on these particles.
    The Stokes drag is the force that makes a particle approximately follow the streamlines
of a fluid, therefore it plays a vital role. Besides this force there will be other forces acting
on the particles. These are for instance a lateral migration force and the gravity force.


5.1     Drag Force
When a particle is suspended in a fluid it will tend to follow the fluid motion. This is due
to the Stokes drag which is present when the particle velocity is different from that of the
fluid.
    In the derivation of the forces acting on the particles, it is used that the fluid motion
around a sphere is the same whether or not the sphere is mowing or fixed at a certain
position. Since the forces are identical, a fixed sphere is used in this derivation, because
then the Stokes drag is easier to derive. Furthermore it is assumed that the particle is
placed in a velocity field that have a uniform velocity at infinity, v∞ .
    To derive the Stokes drag the two fundamental microfluidic equations, the Navier–
Stokes equation for a steady state situation with a Re ≪ 1, so that the nonlinear term is
negligible, and the continuity equation are used.

                                         η∇2 v = ∇p                                       (5.1)

                                           ∇·v =0                                         (5.2)


    By defining v as the velocity around the sphere, and v∞ as the fluid velocity at infinity
it follows that ∇ · v∞ = 0. Combining this with the continuity equation Eq. (5.2) yields:

                        ∇ · (v − v∞ ) = ∇ · v − ∇ · v∞ = ∇ · v = 0                        (5.3)

                                              17
18                                                      CHAPTER 5. PARTICLE MOTION


                                                      v


     v∞                                                   R                                v∞



                                                      v


Figure 5.1: this is an illustration of a fixed sphere in a steady state flow. The fluid will have the
velocity v∞ at infinity and v around the sphere.




From vectorial mathematics it is known that the divergence of a rotation field equals zero.
Comparing this with the fact that ∇ · (v − v∞ ) = 0 it is clear that v − v∞ can be written
as the rotation of a vector field:

                                       v − v∞ = ∇ × A                                       (5.4)

The equation of motion for the liquid around the sphere together with the appropriate
boundary condition at surface of the sphere are:

                 ∂v
               ρ     + (v · ∇)v = −∇p + η∇2 v                                              (5.5a)
                  ∂t
               v(r = R) = 0,      No slip condition on the sphere surface                  (5.5b)

Due to the low Reynolds number Eq. (5.5a) can be simplified like Eq. (5.1). To determine
the Stokes drag force on the sphere due to the fluid motion, one has to determine the
rotation of the vector field A and thereby get a relation between v and v∞ . In order to
determine A, a scalar function f (r) is introduced so that A is given by:

                       A = f ′ (r)n × v∞ = ∇f × v∞ ,            f = f (r)                   (5.6)

f (r) and v can be found to be as stated below.

                                c2 v∞ + n(v∞ · n)      3n(v∞ · n) − v∞
                    v = v∞ −                      − c3                                     (5.7a)
                                2        r                    r3
                                                               c3 c2
                                                     f (r) = − + r                         (5.7b)
                                                                r   2

Where in the following n is the sphere surface normal. From the boundary condition
5.1. DRAG FORCE                                                                        19

stated in Eq. (5.5b) the two constants c2 and c3 can be determined.

                             c2 v∞ + n(v∞ · n)      3n(v∞ · n) − v∞
           v(r = R) = v∞ −                     − c3                  =0⇔
                              2       R                   R3
                                  c2    c3                c2    3c3
                          v∞ (1 −    − 3 ) + n(v∞ · n)(−      + 3)=0⇔
                                  2R R                    2R    R
                                      c2   c3              c2    3c3
                                  1−     −     =0 ∧ −          + 3 =0⇔
                                     2R R3                 2R    R
                                                         3            1
                                                     c2 = R ∧ c 3 = R 3              (5.8)
                                                         2            4
Combining this with Eq. (5.7a) and Eq. (5.7b) f and v are finally found to be:

                        3     1 R3
                f (r) = Rr +                                                        (5.9a)
                        4     4 r
                             3R v∞ + n(v∞ · n) R3 3n(v∞ · n) − v∞
                    v = v∞ −                  −                                    (5.9b)
                              4       r         4       r3
Rewriting Eq. (5.9b) to spherical coordinates yield the following expression.

                                                   3R   r3
                              vr = |v∞ | cos θ(1 −    + 3)                        (5.10a)
                                                   2r  2r
                                                   3R   r3
                              vθ = |v∞ | sin θ(1 −    + 3)                        (5.10b)
                                                   4r  4r
When the velocity field around the sphere is known it is possible to determine the pressure
distribution on the surface of the sphere. Starting from the simplified Navier–Stokes
equation Eq. (5.1) the pressure distribution in terms of f and v∞ is found to be:

                                  p = p0 + ηv∞ · ∇∇2 f                              (5.11)

Insertion of Eq. (5.9a) into Eq. (5.11) yields the following expression for the pressure
distribution:
                                           3 v∞ · n
                                  p = p0 − η         R                            (5.12)
                                           2    r2
The drag force on the sphere is then given by:

                             ′                          ′           ′
              F=     (pni − σik nk )ds =   (−p cos θ + σrr cos θ − σrθ sin θ)ds     (5.13)

Here it is used that the vertical forces counteracts each other. The stress tensors are
According to [5] given by:

                                 ′         δvr
                                σrr = 2η                                          (5.14a)
                                            δr
                                 ′         1 δvr   δvθ  vθ
                                σrθ   = η(       +     − )                        (5.14b)
                                           r δθ    δr   r
20                                                    CHAPTER 5. PARTICLE MOTION

When these tensors components and the pressure distribution are evaluated at the surface
of the sphere, they are given as:

                                   σrr = 0                                              (5.15a)
                                           3η
                                    σθ = − v∞ sin θ                                     (5.15b)
                                           2R
                                              3η
                                     p = p0 −    v∞ cos θ                               (5.15c)
                                              2R
When this is inserted into Eq. (5.13) the formular for the Stokes drag on a fixed sphere is
found to be:
                                       F = 6πηRv∞                                  (5.16)
If the particle has a velocity itself for instance u as it is illustrated in Fig. 5.2 the Stokes
drag will be given as:
                                      F = 6πηR(v∞ − u)                                    (5.17)
It is seen that the force is proportional to the velocity difference. Therefore will this force
make the particles follow the streamline motion.




                 v∞                                      u                v∞




Figure 5.2: A mowing particle in a constant flow. It will experience a drag force proportional to
the velocity difference between the sphere and the fluid.




5.2     Other Forces Acting on the Particle
Even when a particle is moving in an unidirectional flow it will experience a lateral force
perpendicular to the streamlines of the flow. This is due to the fact that the velocity profile
in a channel of a finite size isn’t uniform. In a channel having the cross sectional form
described in Chapter 1 the velocity will have a parabolic dependency on the distance to
all of the channel sidewalls. The nonconformity leads to a difference of stress and pressure
on the two sides of a particle suspended in this flow. Therefore the particle will experience
a force acting perpendicular to the fluid streamlines.
    In the previous section we considered the case where the particle was placed in an
infinitely large channel with an uniform velocity profile. The pressure and stress on the
5.2. OTHER FORCES ACTING ON THE PARTICLE                                                      21

two sides of the particle are then equal and there would be no lateral migration. In the
actual device the velocity field will not be uniform, and therefore this physical effect needs
to be considered.
    Since it is vital for the functionality of the pinched flow that the particles are aligned
along one sidewall, the lateral migration could have a big influence on the final device.
    The lateral migration effect has been described via a perturbation theory presented
by Ho and Leal [6]. They investigated the lateral migration principle for particles with
Reynold’s number smaller than unity. Their investigation was focussed on circular chan-
nels, but since the lateral migration force is due to the nonuniform velocity profile, which
is also present in rectangular channels, the description should be equivalent for rectangular
channels as well.
    The starting point of this deviation is the Navier–Stokes- and continuity equation for
the velocity and pressure field v = U − V and q = P − Q respectively. Here U and Q are
the fields that the fluid would have if the particle wasn’t present. Where as V and P are
the fluid fields due to the present of the particle. These equations are together with the
appropriate boundary condition:

                       ∇2 v − ∇q = Re(v · ∇v + v · ∇V + V · ∇v)                          (5.18a)
                             ∇·v =0                                                      (5.18b)
                                 v = Ωs × r − V         on r = 1                         (5.18c)
                                 v=0          on the sidewalls                           (5.18d)
                                 v→0           as r → ∞                                  (5.18e)

The two first equations are the Navier–Stokes equation and the continuity equation respec-
tively. Eq. (5.18c) is the equation of motion for the rotation of the particle and is actually
a way of applying the ”no slip” condition on the surface of the particle. In this equation
Ωs is the rotation velocity of the sphere. Eq. (5.18d) dictates the ”no slip” condition on
the sidewalls of the channel and Eq. (5.18e) states that at infinity in the channel direction
the particles should have the same velocity as the fluid.
    Using these equations together with the assumptions of a Poiseuille flow profile, [6]
derive the force acting on the sphere across the streamlines. From this force they derive
the velocity of a particle in steady state traveling in an unidirectional Poiseuille flow. This
velocity is determined to be:
                              2
                            ρvm w 3
                     uL =         κ 36[(1 − 2s)2 G1 (s) − (1 − 2s)G2 (s)]                 (5.19)
                            Re6πη

In this equation vm is the maximum fluid velocity, κ is defined as the particle radius divided
by the channel width, s is the particle distance from the sidewall divided by the channel
width and G1 (s) and G2 (s) are constants that they define in the paper. For s = 0.19 these
are 0.982 and 0.536 respectively. This s value represents a big white blood cell with radius
of 7.5 µm.
    In order to get a feeling of the size of this lateral migration, let us consider a particular
case with a maximum fluid velocity vm = 10 µm/s. When a Poiseuille flow is assumed the
22                                                           CHAPTER 5. PARTICLE MOTION

                                                                              2
mean velocity will be given, in term of the maximum velocity, as vm . If a 100 µm long
                                                                   3
channel is considered, the time the fluid spends in the channel will then on average be:

                                                  100 × 10−6 m
                                  ttravel =   2                   = 15s                 (5.20)
                                              3   × 10 × 10−6 m/s

Using this time together with the lateral migration velocity stated in Eq. (5.19) it is
straight forward to determine the lateral migration distance.

             dL =uL ttravel                                                             (5.21)
                    103 kg/m3 (10 ×     10−6 m/s)2
                                               × 40 ×         10−6 m
             dL =                                              0.18753                  (5.22)
                            2 × 10 −4 × 6π × 10−3 Pa s

                    × 36[(1 − 2 × 0.19)2 × 0.982 − (1 − 2 × 0.19) × 0.536] × 15s
                ≈ 0.171 µm                                                              (5.23)

This distance is very small compared to other distances in the channel. It is around 1
pro mille of the particle radius. This small lateral migration distance even for the largest
particle is indicating that the effect can be neglected. But in order to be sure let us compare
the lateral migration force to the Stokes drag force derived in the previous section. As an
example it is assumed that the particle velocity is 1% below that of the flow.
    The lateral migration force is according to [6] given by Eq. (5.24).

                      FL = ρu2 r2 κ2 36[(1 − 2s)2 G1 (s) − (1 − 2s)G2 (s)]
                             m                                                          (5.24)

In this equation r is the radius of the particle. Again values for the largest white blood
cells are used in the calculation.

              FL =103 kg/m3 × (10 × 10−6 m/s)2 × (7.5 × 10−6 m)2 × 0.18752
                     × 36[(1 − 2 × 0.19)2 × 0.982 − (1 − 2 × 0.19) × 0.536]
                 ≈ 3.215 × 10−19 N                                                      (5.25)

The Stokes drag on the same particle is, as stated in the previous section, proportional to
the velocity difference between the particle - and the fluid velocity. When the same mean
                          2
fluid velocity as above of × 10−5 m/s is assumed the Stokes drag will be:
                              3

                FStokes = 6πηr(v∞ − u)
                                                                       2
                FStokes = 6π10−3 Pa s 7.5 × 10−6 m × 0.01 ×                × 10−5 m/s
                                                                       3
                         = 9.42 × 10−15 N                                               (5.26)

By comparing the lateral migration force with the Stokes drag force, it is seen that the
lateral migration force is around 0.000341% of the Stokes drag.

                                  FL       3.215 × 10−19 N
                                       =                   = 0.000341%                  (5.27)
                           FStokes          9.42 × 10−15 N
5.3. DRAG NUMBER                                                                          23

Since the lateral migration distance is so small it is a valid simplification to neglect it.
    The blood cells will of course also be influenced to some extend by the gravitational
force. But as described in Chapter 3 will the blood cells have almost the same density as
the buffer in which they are suspended. Therefore the gravitational force acting on the
particles will be almost balanced by the lifting force from the fluid. And even though it
had an effect it will make no difference whether the particle is aligned against the top or
bottom of the channels.


5.3     Drag Number
In order to determine how particles suspended in a fluid moves let us consider Newton’s
second law for these particles. In this section only the Stokes drag force is considered,
Newton’s second law for a particle of mass m reads:
                                           ∂u
                                        m       = Fdrag                                (5.28)
                                            ∂t
                                                                                 4
The particle mass, m, can be rewritten in terms of the particle density (m = 3 πr3 ρ) when
sphere shaped particles is assumed. Then Eq. (5.28) takes the form:
                               4 3        ∂u
                                 πr ρpart      = 6πηa(v∞ − u)                          (5.29)
                               3           ∂t
In this equation v∞ is the fluid velocity at infinity and u is the particle velocity. This
equation can be rearranged to give the following:
                             4 3        ∂u
                               πr ρpart      = 6πηr(v∞ − u) ⇔
                             3          ∂t
                                          2ρpart r2 ∂u
                                                       = v∞ − u                        (5.30)
                                              9η    ∂t
By introducing the characteristic velocity, V0 and length scale, L0 , Eq. (5.30) can be made
dimensionless. Then by introducing the dimensionless drag number Dr the equation
simplifies further.
                           ∂u˜                            2ρpart r2 V0
                        Dr         ˜
                                = v∞ − u, ˜          Dr =                              (5.31)
                            ∂t                              9ηL0
When this drag number is small (≪ 1) the time dependent part can be neglected and
Eq. (5.31) simplifies into:
                                            ˜ ˜
                                            u = v∞                                     (5.32)
This equation tells us that for small drag numbers the particle will follow the streamline
going though its center of mass.
    To determine whether this streamline motion is valid for the red and white blood cells
mowing in the buffer, let us calculate the drag number for the largest of the white blood
cells.
              2 × 1125kg/m3 × (7.5 × 10−6 m)2 1 × 10−5 m/s
        Dr =                                                  ≈ 3.516 × 10−6           (5.33)
                        9 × 10−3 Pa s × 40 × 10−6 m
This is well below unity. It is therefore a valid assumption that a particle will follow the
streamline, which goes through its center of mass.
24                                                    CHAPTER 5. PARTICLE MOTION

5.4     Diffusion Constants
The particles will to some extend be influenced by the diffusion process.
   A critical parameter when considering diffusion of particles suspended in a liquid is
the diffusion constant. This can according to E. Guyon et al. [7] be determined by the
Einstein relation:
                                            kB T
                                       D=                                        (5.34)
                                            6πηr
When the characteristic values for the red and the white blood cells are inserted the
diffusion constants for these cells yields:

                             1.38 × 10−23 J/K × 300K
                    Dred =                           ≈ 8.79 × 10−14 m2 /s               (5.35a)
                               6π10−3 Pa s × 2.5 µm
                             1.38 × 10−23 J/K × 300K
                  Dwhite   =                         ≈ 2.93 × 10−14 m2 /s               (5.35b)
                               6π10−3 Pa s × 7.5 µm

In the following work a diffusion constant of 10−13 will be used for both red and white blood
cell. This is done in order to be sure that the effect from diffusion is not underestimated.
If the actual diffusion constants are smaller the device will only improve its functionality.


5.5     Motion of the Particles
Of the three forces considered here, it was seen that both the lateral migration force and
gravity could be neglected.
    It was further seen in the previous section that when the drag numbers are much
smaller than unity a particle will follow the streamline going through its center of mass.
In the present situation this is the case, since it was seen that even for the largest particles
considered, this drag number was of the order 10−6 . In the rest of the thesis this will
therefore be assumed whenever particle motion is considered.
Chapter 6

Analytical Calculations

The present analytical calculations are made in order to determine the expected result of
the simulations. All these analytical calculations are made on the basis of the geometry
shown in Fig. 6.1. All the analytical calculation in this chapter are made in quasi 3D
mode.
                                                  y                 Outlet channel
               Main channel


                                                                vout
                                 v
                       w                                       Qtotal
                                 Q2
                                                                                          x

                                             Q1                 L
                                                                             Borderline
           Borderline asymptote           Usf v


               Particle enters                          Inlet channel
                                            d

Figure 6.1: The geometry used in the analytical calculations in this chapter. It is assumed that
the velocity in the inlet channel is controlled so that it will be a fraction of the velocity in the
main channel by the scaling factor Usf according to vinlet = Usf v. The red borderline in the figure
represents the boundary between the two fluids. The asymptote of this borderline is also shown.




6.1     Lateral Penetration Distance
In order to verify that the inertia of the cells does not play a role, in other words that the
particles just follow the streamlines going through its center of mass, let us investigate how
a particle from the inlet channel is accelerated in the main channel. It is in this calculation

                                                  25
26                                       CHAPTER 6. ANALYTICAL CALCULATIONS

assumed that the fluid will have a constant velocity of vmean in the main channel and a
velocity of 0.3vmean in the inlet channel in the direction illustrated in Fig. 6.1.
    When a particle enters the main channel from the inlet channel it will be accelerated
in the x direction by the fluid due to the Stokes drag force. In this section it will be
considered how long time the Stokes drag needs to accelerate the particle up to a certain
amount of the fluid velocity. The velocity in the x direction is described by the following
differential equation.
                             4           ∂ux
                               πρpart r3     = FStokes ⇔
                             3            ∂t
                                         ∂ux      9η
                                             =           (vx − ux )                      (6.1)
                                          ∂t   2ρpart r2
This ODE can easily be solved to give the particle velocity as a function of time:
                                                             9η
                             ux (t) = vx + c1 exp −                  t                   (6.2)
                                                         2ρpart r2
To determine the constant c1 it must be used that ux (0) = 0. With the use of this
boundary condition the velocity of the particle as a function of time is found to be:

                                                            9η
                           ux (t) = vx 1 − exp −                     t                   (6.3)
                                                         2ρpart r2

When an exponential function is present one often talks about a characteristic time of τ ,
which is defined so that when t = τ the exponential would be e−1 . The characteristic time
in this case is then given by:
                                         9η
                             exp −                τ = exp(−1) ⇔
                                      2ρpart r2
                                                        2r2 ρpart
                                                   τ=                                    (6.4)
                                                           9η
With the use of this characteristic time it is possible to determine the penetration depth for
the largest white blood cell, which is the one that will experience the largest penetration
depth. In this calculation it is assumed that the velocity in the y direction is 3 µm/s.

              Ly =uy τ ⇔
                     2 r2 ρpart
                 =uy            ⇔
                     9 η
                                 2 (7.5 × 10−6 m)2 × 1125kg/m3
              Ly =3 × 10−6 m/s                                 ≈ 0.042 nm                (6.5)
                                 9           10−3 Pa s
From this result it is seen that the inertia of the cells is completely negligible in the
description of their motion. The actual penetration depth could be even smaller. This is
because it has not been taken into account in this calculation that the velocity in the y
direction drops due to the Stokes drag because the main channel fluid have zero velocity
in this direction.
6.2. BORDERLINE POSITION                                                                  27

6.2     Borderline Position
A central point in the sorting mechanism is that the particles must be aligned up against
one of the sidewalls. This is achieved by the introduction of a higher flow rate in the main
channel. This causes the inlet flow to be squeezed against the sidewall
    In this section the particles diffusion process will be neglected. It will be considered
in the following section.
    When diffusion is neglected the inlet channel flow and the main channel flow will remain
separated. Therefore it is possible to determine the asymptotic borderline position shown
in Fig. 6.1 as a function of the relative flow rate difference between the main channel and
the inlet channel.
    To derive the expected asymptotic borderline position (w′ ) the fact that the fluids will
                                                                        Qinlet
remain separated is used. When this is the case the flow rate ratio             , which is the
                                                                        Qmain
inlet channel flow rate divided by the main channel flow rate, will be constant through
out the outlet channel, which yields:

                                           w′
                               Qinlet     0 dy    vx (y)
                                      =    w             ⇔
                               Qmain      w′ dy   vx (y)
                                           w′
                                          0 dy    K(w − y)y
                                     =     w                ⇔
                                          w′ dy   K(w − y)y
                                           1 ′2      1
                                             w w − w′ 3
                                     =1    2         3                                  (6.6)
                                        w 3 + 1 w′ 3 − 1 w′ 2 w
                                      6       3        2

The solution to this equation is shown in Fig. 6.2. There is no further information in the
expression for the result, and therefore it is not shown here. When diffusion is neglected
the borderline will not change after it have reached equilibrium. Therefore the position
shown in Fig. 6.2 is also the position of it at the end of the outlet channel.


6.3     Diffusion Process
In the previous section the diffusion process was neglected. This is of course not a valid
simplification, so in order to determine the influence due to the diffusion process a simpli-
fied model is used.
     In this chapter the treatment of particles will be done by a concentration profile. In
such a way that this profile can be seen as a probability distribution for the particles.
Through this concentration treatment we will get an understanding of how different pa-
rameters influence on the particle position.
     To get an idea of whether diffusion or convection will dominate, let us determine the
  e
P´clet number for the simplified model showed in Fig. 6.1. As the characteristic velocity,
                                                                     e
V0 , the mean velocity in the outlet channel is used. This gives a P´clet number when the
scaling factor, which is the maximum velocity in the inlet channel divided by the maximum
28                                             CHAPTER 6. ANALYTICAL CALCULATIONS




             Intersection width w′ [µm]




                                                            Qinlet
                                          Flow rate ratio   Qmain
Figure 6.2: This is the analytically determined borderline position between the two flows. It is
seen that it as expected grows for increasing flow rate ratios.



velocity in the main channel, is 0.3 of:

       v(w + Usf d)w 1 × 10−5 m/s × (40 × 10−6 m + 0.3 × 20 × 10−6 m) × 40 × 10−6 m
Pe =                =
         (w + d)D               (40 × 10−6 m + 20 × 10−6 m) × 10−13 m2 /s
                    ≈3067                                                      (6.7)

                       e
This is a rather big P´clet number, the v is the maximum velocity in the main channel.
Therefore it will be expected that convection dominates, and that the fluids will not
manage to mix completely before the end of the channel.
    In order to derive an expression for the concentration profile, the model shown in
Fig. 6.3 is used. In this model the inlet channel is interpreted as a point source with a
concentration described by a delta function. The actual size of the concentration does
not matter, since it is only a factor that scales the result, therefore unity is chosen. The
velocity profile in the main channel is now assumed to be constant. When the velocity
in the main channel is assumed constant the x coordinate can, besides being a position
coordinate, be interpreted as a time scale. This is easily realized if one imagines a moving
coordinate system that travels along with the flow in the main channel. In the cross
sections of this coordinate system the diffusion process will look like the diffusion as
constant total concentration process instead of the constant concentration process. In this
moving coordinate system the delta function describing the concentration profile would
be the boundary condition at x = 0.
6.3. DIFFUSION PROCESS                                                                          29

                                         y




                                             vmean

                                                                                  x

                                         Point source
Figure 6.3: The model used for the derivation of an analytical expression for the diffusion process.
The inlet channel is treated as a point source with a constant concentration described by a delta
function. The fluid in the main channel is assumed to have a constant velocity of vmean .



   The diffusion equation along with the appropriate boundary conditions is, as described
in Chapter 2, given as stated below. Due to steady state the time derivative can be
neglected.
                                     ∂c
                                        + v · ∇c = D∇2 c                                    (6.8a)
                                     ∂t
                                          c(0, y) = N0 δ(y)                                 (6.8b)
In this model the second derivative with respect to the x coordinate can be neglected due
to the moving fluid, since this will make the curvature of the concentration very small in
the x direction. The velocity in the main channel is assumed to be a constant in the x
direction and thereby zero in the y direction. This makes the y derivative vanish on the
left hand side of Eq. (6.8a). These simplifications transform Eq. (6.8a) into:
                                                 ∂c   ∂2c
                                         vmean      =D 2                                     (6.9)
                                                 ∂x   ∂y
In order to solve this PDE the Fourier transform method can be used. According to S. Lea
[8] is the Fourier transform of the concentration profile and the second order derivative
given by:
                                                +∞
                                          1
                              C(x, k) = √           dy c(x, y) exp(−iky)                   (6.10a)
                                          2π −∞
                         ∂ 2 c(x, y)
                       F              = (ik)2 C(x, k)                                      (6.10b)
                             ∂y 2
In this equation C is the Fourier transform of the original concentration. The method is
used to transform the original PDE into an ODE (Eq. (6.11)) which is easier to solve.
                                       ∂c(x, y)    ∂ 2 c(x, y)
                                 vmean          =D             ⇔
                                         ∂x            ∂y 2
                                      ∂C(x, k)
                                vmean           = −Dk 2 C(x, k)                             (6.11)
                                         ∂x
30                                          CHAPTER 6. ANALYTICAL CALCULATIONS

This new ODE is easily solved to give the solution:
                                                        −k 2 D
                                  C(x, k) = C0 exp             x                            (6.12)
                                                        vmean
The constant C0 can be determined with the use of the boundary condition stated in
Eq. (6.8b).
                                           +∞
                                      1
                     C0 = C(0, k) = √          dy c(0, y) exp(−iky)
                                      2π −∞
                                  +∞
                           N0
                        =√           dy δ(y) exp(−iky)
                            2π −∞
                           N0            N0
                        = √ exp(0) = √                                                      (6.13)
                            2π            2π
Then the solution to the original PDE is found using the inverse Fourrier transform.
                                  +∞
                             1
                  c(x, y) = √        dk C(x, k) exp(iky)
                              2π −∞
                                      +∞
                             1 N0                  −k 2 D
                          =√ √           dk exp           x exp(iky)
                              2π 2π −∞             vmean
                                  N0             −vmean y 2
                          =                exp                                              (6.14)
                              2     πDx            4Dx
                                   vmean

It is seen that the concentration profile has the same structure as for a normal diffusion
                                                                               x
process, except that the time parameter in this expression is replaced by          . A plot
                                                                             vmean
of this concentration profile is shown in Fig. 6.4 which is a cross section at x = 80 µm.
In this figure it is seen that the effect due to the diffusion process is that the boundary
between the two flows will be made less abrupt. There will be established a mixed layer
between them. The width of this mixed layer can be interpreted as the Full Width Half
Max (FWHM) of the concentration gradient. In order to be able to control the position
of the particles this FWHM should be as small as possible.
    It is of course possible to get an analytical expression for this FWHM. It will be
determined in the following. The concentration gradient in the y direction is:
                                     −N0 y                    −vmean y 2
                       c′ (x, y) = √
                        y               Dx
                                                  exp                                       (6.15)
                                  4 π( vmean )3/2               4Dx
In order to determine the position of the maximum concentration gradient it is necessary
to differentiate this with respect to the y coordinate again.
                    −N0                    −vmean y 2       N0 y 2             −vmean y 2
     c′′ (x, y) = √
      yy               Dx
                                 exp                    + √    Dx
                                                                         exp
                 4 π( vmean )3/2             4Dx         8 π( vmean )5/2         4Dx

                       −1            y2            N0          −vmean y 2
              =      Dx
                               +    Dx 5/2
                                                   √ exp                                    (6.16)
                   ( vmean )3/2 2( vmean )        4 π            4Dx
6.3. DIFFUSION PROCESS                                                                         31




    Concentration ratio c(y)
                        c(0)




                                                     Channel positiony [µm]

Figure 6.4: Cross section of the concentration profile at x = 80 µm. It is seen that at this x
position the maximum gradient position is a little below 2 µm. This figure is made with a diffusion
constant of 1 × 10−13 m2 /s




The zero point for Eq. (6.16) should then be determined. This is quite trivial, since it is
only necessary to solve the parenthesis equal to zero which yields the value for y stated
below. There is also a negative y solution but since this makes no physical sense it is
neglected.


                                                                 2Dx
                                                  ymax grad =                               (6.17)
                                                                 vmean



In order to determine an analytical expression for the FWHM of the concentration gradient
it is necessary to find the value of the y gradient for the y value stated in Eq. (6.17) which
yields:



                                                                                        2
                                                                                 2Dx
                                                    −N0 v2Dxmean
                                                                       −vmean   vmean
                                ′
                               cy (x, ymax grad ) = √    Dx 3/2
                                                                 exp
                                                   4 π( vmean )             4Dx
                                               −N0 vmean
                                              = √        exp(−1/2)                          (6.18)
                                                  8πDx



It is now possible to determine the FWHM using this expression as the maximum value.
32                                       CHAPTER 6. ANALYTICAL CALCULATIONS

First the two half max positions are found:
                                                     1
                                          c′ (x, y) = c′ (x, ymax grad )
                                           y
                                                     2 y
                                                     −N0 vmean
                                          c′ (x, y) = √
                                           y                     exp(−1/2) ⇔
                                                      2 8πDx
                  −N0 y               −vmean y 2       −N0 vmean
                √    Dx
                               exp                 =    √        exp(−1/2) ⇔
               4 π( vmean )3/2          4Dx            2 8πDx
                                      −vmean y 2          Dx
                            y exp                  =            exp(−1/2)           (6.19)
                                        4Dx              2vmean

Eq. (6.19) is the last step that can be done analytically. To solve this last equation a
numerical method must be used. It is solved with the use of the commercial programme
Mathematica 5.1 and the resultant FWHM are:

                                                    Dx
                                     HM1 =2.718                                   (6.20a)
                                                   vmean
                                                    Dx
                                     HM2 =0.451                                   (6.20b)
                                                   vmean
There were also two negative solutions which of course are rejected. By subtracting
Eq. (6.20b) from Eq. (6.20a) the FWHM is found.

                                                               Dx
                         FWHM = HM1 − HM2 = 2.266                                   (6.21)
                                                              vmean
This predicted position of the FWHM is shown as a function of the x position in the outlet
channel in Fig. 6.5. When the values for the actual channels are inserted, the theoretical
expectation value can be determined. It is the FWHM at the end of the outlet channel
that is interesting. At this position the FWHM will be:

                                 1 × 10−13 m2 /s × 100 × 10−6 m
              FWHM = 2.266                                      ≈ 2.27 µm           (6.22)
                                         1 × 10−5 m/s

Therefore the theoretical expectation value for the FWHM is 2.27 µm.


6.4    Theoretical Expectation Values
It is now possible to state some theoretical expectation values prior to the simulation
process. The simulation will be on the same geometry as was shown in Fig. 6.1. Therefore
it should be possible to do direct comparison.
    In the section concerning diffusion it was found that the FWHM of the concentration
gradient at the end of the outlet channel is around 2.3 µm. In the theoretical treatment
of the diffusion process it was seen that the FWHM was independent of the scaling factor,
 6.5. PARTICLE ALIGNING                                                                            33


y position [µm]

                      HM1

                                                                                          FWHM
                            M


                           HM2



                    x position [10−5 m]     Concentration gradient [abitrary units]

 Figure 6.5: On the left the theoretical prediction of the maximum gradient position is shown
 together with the positions of the two half maximum values. It is seen that the upper of the HM’s
 is further away from the maximum position. This is due to the presence of the wall which dictates
 that the derivative of the concentration should be zero here. This non symmetry can also be seen
 on the right part, which is the theoretical prediction of the concentration gradient. It is here seen
 that it is not symmetric around the maximum position. This is due to the wall effect. These plots
 are made with the following constants: vmean = 1 × 10−5 m/s and D = 1 × 10−13 m2 /s



 and thereby also the of the flow rate ratio. This independence could be due to the fact that
 the velocity profile in the treatment of the diffusion process was taken to be a constant,
 which will certainly not be the case in the actual device as well as in the simulations.
    It is further more expected that the intersection width, which is another word for the
 asymptote of the borderline, is given by the borderline position shown in Fig. 6.2.


 6.5              Particle Aligning
 A central concept in the pinched flow is as described earlier the alignment of particles
 towards one of the side walls in the pinched segment.
    According to Yuushi et al. [4] the particles can be assumed to be aligned against the
 channel sidewall when Eq. (6.23) is valid.
                                             Qinlet        2r
                                                         ≤                                     (6.23)
                                          Qmain + Qinlet   w
 Here Qinlet and Qmain is the inlet channel flow rate and the main channel flow rate re-
 spectively. Eq. (6.23) is only valid when the velocity profile of the fluid is assumed to be
 constant, which is not a valid simplification.
     In order to get a more correct alignment demand one must take the parabolic behavior
 into account.
     When the situation is assumed to be in 2D, the demand for alignment is that the inlet
 channel flow rate divided by the total flow rate should be less than the flow rate between
34                                     CHAPTER 6. ANALYTICAL CALCULATIONS

the wall and the intersection width divided by the total flow rate.
                                           2r
                          Qinlet          0 dyK(w − y)y
                                      ≤    w                 ⇔
                       Qmain + Qinlet     0 dyK(w − y)y
                                            1            1
                          Qinlet        K w(2r)2 − (2r)3
                                          2        3
                                      ≤                             ⇔
                       Qmain + Qinlet       1 3
                                         K w − w
                                                  1 3
                                               2       3
                                          1          1
                          Qinlet            w(2r)2 − (2r)3
                                          2          3
                                      ≤          1 3
                                                                 ⇔
                       Qmain + Qinlet              w
                                                 6
                          Qinlet                                     2r
                                      ≤3β 2 − 2β 3      , for β =                    (6.24)
                       Qmain + Qinlet                                w

This inequality is just another way of determining the intersection width (in this equation
denoted 2r). Therefore this inequality express the same as Eq. (6.6).
    To get a feeling of the mistake one makes by using Eq. (6.23) instead of Eq. (6.24),
let’s compare the inlet flow rate ratios demand for aligning of a particle with a radius of
5 µm in a 40 µm wide pinched segment. The original flow rate and the new flow rate
demand is given by:
                  Qinlet      2r 2 × 5 µm
                             ≤ =          = 0.25                                    (6.25a)
               Qmain + Qinlet w    40 µm
                                                2                    3
                  Qinlet         2 × 5 µm              2 × 5 µm
                              ≤3                    −2                    ≈ 0.156   (6.25b)
               Qmain + Qinlet     40 µm                 40 µm

It is seen that the difference due to the new alignment demand is approximately 60%.
If one uses the original alignment inequality one can not be sure that the particles are
aligned towards the channel wall.
Chapter 7

Numerical Simulation Method

In this chapter the simulation method, that has been used in the project, will be described.
A short explanation of the Finite Element Method (FEM) will be given. It will furthermore
be discussed how the mesh structure influences on the results. Finally a quasi 3D method
is derived, so that simulations in 3D can be performed on a 2D model.
    The Navier–Stokes- and the diffusion-convection equation can only be solved analyti-
cally in very simple, highly symmetrical, cases. Therefore it will be necessary to perform
numerical simulation in order to fully analyze the PFF-system. Nonlinear partial differ-
ential equation like the Navier–Stokes equation can solved numerically in many different
ways. The method used, which is called the Finite Element Method (FEM), divides the
analyzed structure into discreet elements and transforms the continuous equation into a
matrix problem. These problems can be solved with the use of different softwares. In this
projet all solutions are determined with the use of the simulation program: Comsol 3.3,
which is specialized in solving partial differential equations. In script mode the simulation
was preformed by Comsol with Matlab, using Matlab version 7.0.4.365(R14) Service Pack
2. If nothing else is mentioned the simulations are made on a standard PC.


7.1     Finite Element Method
When dealing with fluid dynamics one often wants to determine the velocity field of the
flow. This field can be represented by the vector field v(r, t) which is the solution to a
general partial differential equation.

                               Lv(r, t) = f (r, t),   for r ∈ Ω                          (7.1)

In this equation L is a differential operator describing the physical properties of the system,
f is the source or force term and Ω is the domain in which the problem should be solved.
The boundary condition on ∂ Ω are typically given by either Diriclet or Neumann condition.
    The problem is discretized and it is then possible to determine an approximate solution
which satisfies the boundary condition. An example of a division of a 2D Gaussian-shaped
channel into finite elements is shown in Fig. 7.1. Each element is then associated with a
basis function φj , which can be a variety of shapes depending on the situation. In Fig. 7.2

                                               35
36                                 CHAPTER 7. NUMERICAL SIMULATION METHOD

a linear basis function within an element is illustrated. Because of these function it is
possible to construct the solution to Eq. (7.1) on discrete form:
                                                   N
                                      v(r, t) =          vi φi (r, t)                             (7.2)
                                                   i=1

In order to determine the solution one needs to determine the vi ’s in this equation for all
elements N.




Figure 7.1:      Illustration of the Finite Ele-         Figure 7.2:     This is an illustration of how
ment Method. The Gaussian-shaped domain Ω                a single element which is described by a linear
is divided into a mesh of finite sized elements.          basis function meets the boundary condition of
Within this mesh one of the linearly interpolat-         unity in the raised corner. The figure comes
ing basis functions φj is shown. This figure is           from [2]
from L. H. Olesen [2]



7.2     Mesh Structure
When the FEM method is used in simulations, there will be some limitations on the pre-
cision of the simulated results do to the finite element size. This can result in nonphysical
solutions and must therefore be considered. This section will comment on the different
results that can arise do to a chance in element size.
     In our simulation, the mesh structure has been defined so that the important parts
of the PFF-system have a finer mesh structure than the rest of the model. In Comsol
it is possible to control the mesh structure by various parameters. The parameters used
in defining the mesh in this project are listed in Table 7.1. When a model is defined in
Comsol, each corner, boundary and subdomain are indexed. Therefore it is possible to
assign the parameters described in Table 7.1 to specific parts of the model. In the fluid
interaction simulation described in Chapter 8 it was chosen to optimize the mesh around
point 4 and 6 and at boundary 1 and 2, see Fig. 7.3. At the intersection between the two
7.2. MESH STRUCTURE                                                                           37

                hmax:        Maximum element size in the overall structure
                hgradvtx:    Element growth rate from a specific point
                hgradedg:    Element growth rate from a specific boundary
                hmaxvtx:     Maximum element size at a point
                hmaxedg:     Maximum element size at a boundary


Table 7.1: Above are listed the parameters, which are generally used in Comsol, to optimize the
mesh structure. These can be defined, so that certain important parts of the structure will have a
finer mesh.




Figure 7.3: This is an illustration of the mesh stucture made with the use of the parameters
described in Table 7.1. At point 4 and 6 the maximum element size are respectively 3 µm and
0.5 µm and 1 µm at boundary 1 and 2. The overall maximum element size is 6µm. The mesh is
optimized so that the critical parts of the structure have a finer elements.


fluids from the inlet and the main channel the velocity and concentration profile can be
rather complex. Therefore the mesh structure is made finer at point 4 and 6. Since the
interesting parameter in the simulation is how the fluids moves down the outlet channel
the element size along boundary 1 is also decreased. This finer mesh is indeed necessary
due to the small diffusion constant which results in a sharp concentration drop. A sharp
drop needs fine mesh in order to be described right.
    The mesh structures cannot be infinitely small, since this would result in infinitely
many elements, which of course is not solvable. When a finite number of elements is used
it can for instance result in a stepwise curve where the real result should be a smooth curve.
One needs to make a compromise such that the problem is solvable and the simulation
yields somewhat valid physical results. To determine this mesh related compromise, three
different mesh structures was investigated, their details are shown in Table 7.2.
    In order to test the effect of using the tree different mesh structures they were sub-
jected to the same simulation. The results for the position of the maximum concentration
gradient as function of the scaling factor can be seen at Fig. 7.4, Fig. 7.5 and Fig. 7.6.
 38                                                                         CHAPTER 7. NUMERICAL SIMULATION METHOD

                                                                                                                   Mesh A               Mesh B                      Mesh C
                                                         Max. element size                                         10µm                 6µm                         3µm
                                              Max. element size at point 4                                         6µm                  3µm                         1µm
                                              Max. element size at point 6                                         6µm                  0.5µm                       0.5µm
                                      Maximum element size at a boundary 2                                         6µm                  1µm                         0.5µm
                                      Maximum element size at a boundary 1                                         6µm                  1µm                         0.5µm
                                                      Number of elements                                           537                  6422                        16430

 Table 7.2: Details of the 3 mesh that were investigated in order to find a proper mesh for the
 simulations. In all mesh structures an element growth rate of 1.1 were used.




                                                                                                                              Intersection width [µm]
Intersection width [µm]




                                                                Intersection width [µm]


                          14                                                       14                                                              14


                          12                                                       12                                                              12


                          10                                                       10                                                              10


                          8                                                               8                                                             8


                          6                                                               6                                                             6


                          4                                                               4                                                             4
                           0    0.1     0.2   0.3   0.4   0.5                              0    0.1   0.2   0.3   0.4   0.5                              0    0.1   0.2   0.3   0.4   0.5
                               Scaling factor Usf                                              Scaling factor Usf                                            Scaling factor Usf
 Figure 7.4: The graph shows                                    Figure 7.5:     For the same                                  Figure 7.6: Similarly, when
 the position of the maximum                                    graph, using mesh B it is seen                                using mesh C, a comparison
 contraction gradient at bound-                                 that this simulation is more                                  with mesh A and B yields that
 ary 1 as of function Usf. This                                 smooth than the one made whit                                 this simulation gives the best
 simulation is made with mesh A                                 mesh A. The solution time was                                 result, but the simulation time
 and it is seen that the function                               785s.                                                         was also as long as 1958s
 increases in steps. The simula-
 tion time with this simulation
 was 206s.




 In these figures it is evident that the finer mesh structures that were used, the smoother
 curve one gets.

     Another effect that can arise when a course mesh is used is wriggles at position where
 the physical solution is a constant. An example of this can be seen in Fig. 7.7. These
 wriggles arises when for instance a parabolic basis function is chosen. On the boundary
 between mesh elements it is demanded that the derivative is continues. When large el-
 ements are used the concentration step from zero to unity can be over one element and
 therefore the element after this must have a large gradient at the boundary. In order for
 it to get down to unity it will also have a big gradient at its other boundary. But not as
 big as the starting gradient. The next element then again must have a large gradient for
 the derivative to be continues and so on.
7.3. CONCLUSION                                                                                           39

7.3     Conclusion
In the above section the obvious result that the finer mesh, the better result, was seen. But
since the time it takes to perform a simulation increases with finer mesh a good simulation
will be a simulation that fulfill the two following conditions:

  1. The meshstructure should be fine enough to yield a proper result.

  2. The mesh structures must not be so fine that the simulation time increases unin-
     tended.

Even if these requirements are maintained, the simulated results still have to be analyzed.
Nonphysical results can occur if the meshstructure is insufficient at the critical positions.
   It was seen that mesh A was not fine enough and the fact that mesh C’s simulation
time was 3.2 times longer than mesh B’s, it was decided that mesh B will be used as the
standard mesh in the coming simulations.


7.4     Quasi 3D Simulation Method
Even rather simple simulations in 3D requires a huge number of elements. The number of
elements in the mesh for a 3D model quickly exceeds 40000.
    When channels with a finite height instead of infinitely thick channels are used the
velocity profile will be limited. As described in Chapter 4, the finite height of the channels
acts as an artificial damping force. This leads to the idea that a 3D model with a constant



                                                  1.6
                Concentration [arbitrary units]




                                                                                        d = 20 µm
                                                  1.4

                                                  1.2

                                                   1

                                                  0.8

                                                  0.6

                                                  0.4

                                                  0.2
                                                     0   0.1       0.2       0.3          0.4       0.5
                                                                               Qinlet
                                                               Flow rate ratio Qmain

Figure 7.7: Illustration of the wriggle effect, that arises due to a to course mesh. It arises
because the derivative on the element boundary should be continues.
40                               CHAPTER 7. NUMERICAL SIMULATION METHOD

height can be modeled in a 2D simulation when a damping force, which originates from
the finite channel height, is included in the Navier–Stokes equation. This results in the
modified Navier–Stokes equation:

                       ρ(∂t v + (v · ∇)v) = −∇p + η∇2 v + α(h)v                        (7.3)

In this equation α(h) is the damping factor due to the finite channel height. When α(h)
has been determined, a 3D model with a constant channel height can be simulated in a
2D model.

7.4.1   Determination of the Damping Factor
To determine the damping factor α(h) from Eq. (7.3) one must examine and compare the
results from a setup of parallel plates which is considered both perpendicular to the plates
and along the plate direction. The geometry used is shown in Fig. 7.8. Between the two
plates a fluid is assumed to flow along the x direction. At first the situation is considered
along the y axis. This result in a parallel plate flow. In Chapter 4 the velocity profile for
this geometry was found to be:
                                              ∆p
                                   vx (z) =       (h − z)z                             (7.4)
                                              2Lη
    The model is now considered in the xy-plane and the z dependency is neglected which
results in a 2D model. On this model the Navier–Stokes equation with the artificial
damping term is applied to determine the velocity profile. This equation simplifies even
further than the previous equation because the z dependence is neglected. Eq. (7.3) then
simplifies into:
                                      α(h)v = ∇p                                    (7.5)
This equation is easily solved and yields the following velocity profile.
                                                  ∆p
                                        vx = −                                         (7.6)
                                                 Lα(h)



                 z    y
                     x

                          p0 + ∆p                            p0
                h                                                 b

                                    L
Figure 7.8: This is the geometry used to derive the quasi 3D model. Between two infinitely
large parallel plates a fluid is assumed to flow in the x direction.
7.4. QUASI 3D SIMULATION METHOD                                                         41

An expression for the damping factor α can be found by demanding that the flowrate out
of the model shown in Fig. 7.8 should be identical regardless of which method that is used
to determine the velocity profile. Therefore let us determine the two flowrates. First for
the situation considered along the y axes.
                                        b            h
                             Q1 =           dy           dzvx (z)
                                    0            0
                                        b            h
                                                              ∆p
                                =           dy           dz       (h − z)z
                                    0            0            2Lη
                                  b∆p 1 2 1 3                       h
                                =        hz − z                         ⇔
                                  2Lη 2      3                      0
                                   b∆p 3
                             Q1 =      h                                              (7.7)
                                  12Lη

The flowrate for the situation considered in the xy plane is determined in the same way:
                                        b            h
                             Q2 =           dy           dzvx
                                    0            0
                                        b            h
                                                                 ∆p
                                =           dy           dz −         ⇔
                                    0            0              Lα(h)
                                     b∆p
                             Q2 = −       h                                           (7.8)
                                    Lα(h)

The damping factor α(h) is determined by demanding that the two flow rates should be
identical.

                                     Q1 = Q2 ⇔
                                 b∆p 3       b∆p
                                      h =−       h⇔
                                 12Lη      Lα(h)
                                           12η
                                   α(h) = − 2                                         (7.9)
                                            h
When this damping factor is included in the Navier–Stokes equation as an artificial damp-
ing force, as stated below, a 3D simulation can when the system has a constant height, be
replaced by a less computer capacity demanding 2D simulation.

                                                                             12η
                       ρ(∂t v + (v · ∇)v) = −∇p + η∇2 v −                        v   (7.10)
                                                                              h2
The simulations preformed in this project will all be made by a model with constant
height. The results form these simulations should be used to determine optimal design
parameters for a three dimensional device, and since it was seen by test simulations that
the 2D simulations and the quasi 3D simulations yielded approximately the same result
the quasi 3D model will be used throughout this thesis. The two script used in this
comparison can be seen in Appendix A and Appendix B.
42                             CHAPTER 7. NUMERICAL SIMULATION METHOD

7.4.2   Implementing the Quasi 3D Model
To implement the quasi 3D model, which were described in Chapter 7, in comsol 3.3, the
damping force term α(h)v must be included in the Navier–Stokes equation. This is done
by introducing an artificial damping force in the subdomain settings for the structure.
    When a channel height of 30 µm is used, the hight dependent damping constant α is
13.6 × 106 kg/ms. This damping term will limit the velocity profile and thereby make the
fluid less easy flowing.
Chapter 8

Fluid Interaction Simulation

The simulations in this chapter will concentrate mainly on the interaction process of fluids
in the pinched segment. Due to the diffusion process the interaction process will cause the
transition from one fluid to the other to be less abrupt. Due to the small length scales the
fluids will not mix completely, but to some extend remain more or less separated.
    In the actual device it is the same fluid, that flows in both main- and inlet channel,
the only exception is that the flow from the inlet channel contains particles. In these
simulations the cells will be represented by a concentration profile, which as described in
Chapter 6 can be interpreted as a probability distribution of finding the particle at the
exact position.


8.1    Model Definition
The simulations in this chapter are made in quasi 3D mode on a 2D model. The influence
of a finite channel height is as described in Chapter 7 included as an artificial damping
term in the Navier–Stokes equation. A channel height of 30 µm is used in the damping
term. A schematic illustration of the geometry used in the simulation is shown in Fig. 8.1.
In this geometry it is assumed that the particle containing flow is introduced via the inlet
channel of width d. A flow without particles is introduced in the main channel of width
w. The two flows enter the pinched segment and leave again through the outlet channel,
which in this model acts as the pinched segment. The outlet channel is made 100 µm
long. This particular length was chosen because earlier experiments have documented the
sorting effect with a pinched segment of this length, see [3] and [9]. The simulation will
not contain actual particles. Instead a concentration profile will be used to illustrate the
particle distribution. In order to meet the requirement that particles are introduced via
the inlet channel the initial concentration conditions are unity in the inlet channel and
zero in the main channel.
    One of the parameters, which will be investigated, is how the relationship between the
two entrance velocities influence on the final result. It is only the relationship between the
                                                                             e
velocities in the two channels that influence on the result, as long as the P´clet number is
kept well above unity.

                                            43
44                                        CHAPTER 8. FLUID INTERACTION SIMULATION

                                                      y               Outlet channel
                    Main channel


                                                                  vout
                                     v
                          w                                      Qtotal
                                     Q2
                                                                                            x

                                                 Q1               L
                                                                               Borderline
                Borderline asymptote         Usf v


                   Particle enters                        Inlet channel
                                                d

Figure 8.1: The geometry used in the simulations. It is seen that it is almost identical with the
geometry used in the analytical calculations. This is done so that the results are comparable. The
arrows indicate the velocities in their respective channels. Thereby they also indicate the flow rate
direction in the channels, in this figure illustrated by the Q’s. The line denoted ”Borderline” is
the line illustrating the boundary between the two fluids.



   The velocity in the inlet channel is defined from the velocity in the main channel
according to:
                                             vinlet = Usf v                                     (8.1)

In this equation Usf is the scaling factor between the two velocities. By ensuring that this
factor is below unity the inlet channel velocity will always be smaller than the velocity of
the main channel flow. This is important since the purpose of the main channel flow is to
ensure the alignment of the particles. Alignment is achieved by a larger flow ratio from
the main channel, which thereby will force the inlet channel flow up against the wall.
                                    e
    During the simulations the P´clet number was kept constant. This is done because
              e
the higher P´clet number the better result, so in order to be able to compare the results
                                                               e
for different scaling factors and inlet channel widths, the P´clet number must be kept
constant. This is achieved by varying the main channel velocity v according to:

                                                P e · D(w + d)
                                           v=                                                   (8.2)
                                                 w(w + Usf d)

    The velocity fields in finite sized channels are, as mentioned in Chapter 4, not constant.
Therefore the velocities, stated in Eq. (8.1), are the maximum velocities in the respective
channels.
    The only varying geometric size is the width of the inlet channel. This is done to
determine which relationship between the widths of the inlet and the main channel that
yields the best result. Since it is only the relative difference between the two channel
widths that influences on the result, it is sufficient to vary this parameter instead of both
widths.
8.2. PHYSICAL EXPECTATION                                                                       45

                                             y              Outlet channel
             Main channel


                                                           vout
                            v
                                                          Qtotal
                            Q2
                                                                                    x

                                        Q1       Corner concentration
                                     Usf v                              Intersection width


                  Inlet channel

Figure 8.2: A schematic illustration of some of the parameters that determine the quality of the
result. The corner concentration point is shown. This concentration should be as close to unity as
possible. It is also shown where the intersection width is determined. This distance is a measure
for the width of the particle containing fluid at the end of the outlet channel. It should therefore
be as small as possible, in order to ensure the alignment.



8.1.1    Characterization of a Good Result
Before simulations are initiated it is necessary to define the characteristics of a good result.
On Fig. 8.2 a illustration of the parameters, that characterize the good result is shown.
At the end of the outlet channel the velocity profile should have evened out. So that it
looks like it consisted of only one fluid. A good result is further characterized by a sharp
concentration drop at the end of the outlet channel. The optimal concentration profile
at this position would be a step function. This profile would appear if the fluids did not
mix. The steepness of the concentration drop can be determined by investigation of the
concentration derivative profile versus the y coordinate. As earlier described the FWHM
of this graph will be a good approximation for the width of the mixed layer.
    The fluid in the main channel will have a velocity to the right and therefore the corner
concentration, shown in Fig. 8.2, is a good indicator for the amount of mixing inside the
inlet channel. A good result is therefore characterized by a high corner concentration.
    The last quality factor that will be considered is the intersection width. In the simula-
tions it will be defined as the position of the maximum concentration derivative at the end
of the outlet channel. When it is less than the diameter of the smallest particle we will
expect that all particles are aligned along the sidewall. All this characteristic parameters
determining the quality of the result are summarized in Table 8.1.


8.2     Physical Expectation
We expect that the inlet channel width is determined by the flow ratio relationship between
the main and the inlet channel. Therefore it is expected that the intersection width will
increase when the flow rate ratios increases.
46                                 CHAPTER 8. FLUID INTERACTION SIMULATION

          Parameter          Symbol                      Explanation
      Intersection width       w′           Width of the particle containing flow
                    ∂c
        FWHM of ∂y           FWHM        Width of the mixed layer between the flows
     Corner concentration      —        Concentration at the corner of the inlet channel

Table 8.1: These are the parameters which will be used when determining the quality of the
result. The intersection width and the FWHM should be as small as possible and the corner
concentration should be as large as possible.



               e
    When the P´clet number is kept constant a valid expectation would be that the FWHM
was independent of the scaling factor and thereby the flow rate ratio. This is because the
  e
P´clet number is, as described in Chapter 6, a measure for the convection time relative to
the diffusion time.
    We will expect that the corner concentration grows for increasing flow rate ratios, since
this increment will force the mixing process away from the inlet channel. Furthermore we
will expect that as the width of the inlet channel increases so will this concentration,
because this increment will also force the mixing process away from the inlet channel.


8.3     Simulation Method
The simulations were made with the commercial simulation program ”Comsol3.3” which
is described in Chapter 7. The two equations that were solved was the quasi 3D Navier–
Stokes equation introduced in Chapter 7 and the convection diffusion equation described
in Chapter 2.
    The simulations were made in script mode. This script can be seen in Appendix B.
In order to test whether the above expectations are valid, a number of design parameters
were changed during the simulations.
    The first step in each simulation is to calculate the velocity constants for that particular
scaling factor and inlet channel width. The problem is then solved and for instance the
FWHM is plotted as a function of the flow rate ratio.
    In order to be able to compare the result from the simulations with the result which
was determined in Chapter 6, the scaling factors must be rewritten so that it is a scaling
factor between the flow rates rather than between the velocities:
                                             d/2      h
                                  Qinlet     −d/2 dy 0 dzUsf v
                            Qsf =        =    w/2       h
                                                                 ⇔
                                  Qmain
                                              −w/2 dy 0 dzv
                                     d/2      h
                                     −d/2 dy 0 dzv
                            Qsf =Usf w/2      h
                                                                                         (8.3)
                                     −w/2 dy 0 dzv

From this it is possible to determine a scalar, that the scaling factor should be multiplied
with, in order to rewrite the scaling factor to concern the flow rates.
8.4. RESULTS AND DISCUSSION                                                                                      47


                        4.8

                        4.6             Wall limitation                             d = 20 µm, h = 30 µm

                                                                                    d = 40 µm, h = 30 µm

                        4.4                                                         d = 20 µm, h = 10 µm

            FWHM [µm]                                                               d = 40 µm, h = 10 µm

                        4.2                    Velocity limitation

                         4

                        3.8

                        3.6

                        3.4

                        3.2

                         3

                        2.8
                              0   0.1          0.2        0.3        0.4   0.5        0.6          0.7     0.8
                                                                           Qinlet
                                                       Flow rate ratio     Qmain

Figure 8.3: In this figure the FWHM is shown as a function of the flow rate ratio. It is seen
that for both small and large flow rate ratios the FWHM decreases. From this plot it is also seen
that the effects are controlled mainly by the flow rate ratios, since the graphs are simulated for
different width and heights of the cahnnels. The dashed lines are included to illustrate the effect
of the two limitations.



8.4     Results and Discussion
             e
When the P´clet number is kept constant it is expected that the FWHM would be constant
and independent of the Flow rate ratio. But the simulation showed that this was not the
case. On Fig. 8.4 the FWHM is shown as a function of the flow rate ratio. The simulations
have been performed with various heights and thereby with various damping factors.
    From Fig. 8.3 it is seen that the FWHM of the concentration gradient is controlled
mainly by the flow rate ratio. This can be seen because that all the graphs, regardless of
the parameters used coincide with each other when they are plotted as a function of the
                                                                 e
scaling factor. From this figure it is seen that even though the P´clet number was constant
the FWHM of the concentration gradient, which as stated earlier can be interpreted as
the width of the mixing layer between the two fluids varied between 3 µm and 4.4 µm.
    The FWHM has a maximum of approximately 4.4 µm at a flow rate ratio of around
0.05. In order to understand this maximum it must be realized that two different limiting
effects influence on the FWHM behavior as illustrated in Fig. 8.3. For small flow rate
ratios (< 0.05) the FWHM drops as the Flow rate ratio decreases. When the flow rate
ratio is small so will the width of the concentration profile be. When this is small so is
the FWHM of its gradient, since the concentration gradient at the channel sidewall must
48                                                      CHAPTER 8. FLUID INTERACTION SIMULATION

                 4.8                                                           4.6


                 4.6
                                                                               4.4




                                                                   FWHM [µm]
     FWHM [µm]

                 4.4
                                                                               4.2

                 4.2
                                                                                4
                  4

                                                                               3.8
                 3.8


                                                                               3.6
                 3.6


                 3.4                                                           3.4
                       0   0.05   0.1   0.15      0.2     0.25                       0   0.05   0.1   0.15     0.2    0.25
                                               Qinlet                                                        Qinlet
                           Flow rate ratio     Qmain
                                                                                         Flow rate ratio     Qmain

Figure 8.4: The FWHM at the end of the                           Figure 8.5: This is the FWHM when around
outlet channel as function of the Flow rate ra-                  35000 elements are used. When this graph is
tio. It is seen that there is a maximum FWHM                     compared to Fig. 8.4 it is seen that the small
for a flow rate ratio of around 0.05. This par-                   variation is simulation effect due to the finite
ticular behavior of the FWHM is due to two                       element size. They will therefore not be com-
different limiting effect on the FWHM. These                       mented further.
are explained elsewhere



be zero as described in Chapter 6. The other limiting factor is the non uniformity of the
velocity field. The velocity in the channels will increase towards the center due to the
”no-slip” condition. The actual diffusion process takes place at the boundary between
the two flows. The longer into the channel, this boundary is, the shorter time the fluid
at the boundary is in the channel. When bigger flow rate ratios are used the boundary
position tends towards the center of the channel and thereby limits the diffusion time. A
shorter diffusion time will of course limit the FWHM. This maximum is therefore just the
maximum between these two limiting effects. Therefore a flow rate ratio of 0.05 must be
avoided.
    The wriggles on the right side of the FWHM maximum seen in Fig. 8.3 are not a
physical result. They arise due to the finite mesh size. This is clear when one considers
Fig. 8.4 and Fig. 8.5, in which the FWHM for a model with 30 µm high channels and an
inlet channel width of 20 µm is shown when a course and a fine mesh is used respectively.

    Another effect that needs to be commented arises if one compares the concentration
plots shown in Fig. 8.6 and Fig. 8.7, which are made with a flow rate ratio of 0.01 and
0.02 respectively. On these figures it is seen that when small flow rate ratios are used
the mixing process of the two fluids actually takes place inside the inlet channel, and the
concentration plot show a kind of wave structure. This makes the velocity profile more
complex at the intersection between the two fluids and will therefore increase the particle
diffusion complexity.
    As earlier noted is the corner concentration a parameter illustrating the amount of
mixing in the inlet channel. Therefore it will be investigated how this concentration
depends on the flow rate ratio for two different widths. This can be seen in Fig. 8.8. From
8.5. CONCLUSION AND OPTIMAL DESIGN PARAMETERS                                                                                                        49




                                         Concentration [arbitrary units]




                                                                                                                   Concentration [arbitrary units]
      [10−5 m]




                                                                              [10−5 m]
                       [10−4 m]
                                                                                              [10−4 m]


Figure 8.6: This is the concentration profile                               Figure 8.7: In this figure a flow rate ratio
when a flow rate ratio of 0.01 is used. It is seen                          of 0.02 is used. Here the mixing process only
that the mixing of the fluids actually initiates                            takes place in the early beginning of the inlet
far into the inlet channel due to the small in-                            channel far away from the corner concentration
let velocity relative to the velocity in the main                          measurement point in the inlet channel. It is
channel.                                                                   also seen that the fluids still remain more or
                                                                           less separated at the end of the outlet channel.



this figure it is seen that the corner concentration follows the expected results. It increases
smooth and rapidly as the flow rate ratio increases. It is further seen that around a flow
rate ratio of 0.02 the concentration reaches a constant of unity. This is due to the fact
that the earlier discussed wave structure only exist for flow rate ratios below this value.
    In order to determine whether the theoretically determined intersection widths were
valid, it has been shown in Fig. 8.9 along with a number of simulated results as a function of
the flow rate ratio. From this figure it is seen that the theoretically determined intersection
width is in good agreement with the simulated results regardless of the parameters used.
Therefore the theoretically determined alignment demand will be valid as long as the
diffusion constants of the particles are below 1 × 10−13 m2 /s. The step like behavior is due
to the finite mesh size. This effect is commented in Chapter 7.
    From Fig. 8.3 it is seen that in order to minimize the FWHM either large or small Flow
rate ratios must be used. It was determined from the concentration plot that in order to
avoid mixing in the inlet channel flow rate ratios above 0.02 must be used. Therefore the
interesting part of Fig. 8.3 is the velocity limited region. It is seen that if larger intersection
widths can be accepted it is possible to limit the FWHM by using flow rate ratios above
0.02.
    It is possible to manipulate this intersection widths by changing the width of the
pinched segment. Thereby it is possible to use flow rate ratios above 0.02.


8.5              Conclusion and Optimal Design Parameters
It was found in this chapter that the behaviour of the fluid was mainly governed by the
flow rate ratio, since that the graphs shown in this chapter coincide when they are plotted
50                                                                                                               CHAPTER 8. FLUID INTERACTION SIMULATION




                                                     1.1
            Corner concentration [arbitrary units]
                                                      1


                                                     0.9                                                                                                         d = 20 µm, h = 30 µm

                                                                                                                                                                 d = 40 µm, h = 30 µm

                                                                                                                                                                 d = 30 µm, h = 10 µm
                                                     0.8
                                                                                                                                                                 d = 40 µm, h = 10 µm



                                                     0.7


                                                     0.6


                                                     0.5


                                                     0.4



                                                           0                                               0.1            0.2          0.3                           0.4                0.5
                                                                                                                                        Qinlet
                                                                                                                   Flow rate      ratio Qmain

Figure 8.8: The corner concentration as a function of the flow rate ratio for various channel
widths and heights. It is seen that it approaches unity for flow rate ratios above 0.02.




                                                                                                 0.4
                                                               Intersection width w′ [10−5 m]




                                                                                                0.35


                                                                                                 0.3


                                                                                                0.25


                                                                                                 0.2                                         d = 30 µm, h = 30 µm

                                                                                                                                             d = 40 µm, h = 30 µm

                                                                                                                                              d = 30 µm, h = 10 µm
                                                                                                0.15
                                                                                                                                              d = 40 µm, h = 10 µm

                                                                                                                                             Analytical result
                                                                                                 0.1


                                                                                                0.05


                                                                                                  0
                                                                                                       0         0.1        0.2      0.3              0.4                0.5
                                                                                                                                           Qinlet
                                                                                                                       Flow rate ratio     Qmain

Figure 8.9: The intersection width as a function of the Flow rate ratio for two inlet channel
widths. It is seen that it behaves as expected and increases for for increasing flow rate ratios. It
is the finite mesh size that is reason for the step like behavior.
8.5. CONCLUSION AND OPTIMAL DESIGN PARAMETERS                                          51

as a function of the flow rate ratio.
                                                                              e
     In the previous section in the examination of the FWHM for constant P´clet numbers,
it was seen that for small Flow rate ratios a wave structure appears at the boundary
between the two flows. Therefore one should use flow rate ratios above 0.02. But since the
intersection width increases for increasing flow rate ratios, one must make a compromise.
     In order to get the smallest FWHM it would be preferred to use large flow rate ratios.
They must at least be above 0.05, in order to avoid the FWHM maximum.
     If the outlet channel could be narrowed larger Flow rate ratios could be used. Remem-
ber that at elevated Flow rate ratios the FWHM actually decreased due to the velocity
limitation. This could also be interesting since the larger channels that can be used the
easier the fabrication process becomes and less expensive production method can be used.
     The manipulation of the intersection width by a narrowing of the pinched segment is
commented in Chapter 10
     In the final design a width of 40 µm and a height of 30 µm will be used. From Fig. 8.9
it is seen that in order to insure alignment with this pinched segment width a flow rate
ratio of 0.1 must be used.
52   CHAPTER 8. FLUID INTERACTION SIMULATION
Chapter 9

Separation Mechanism

The most important part of the final device is of course the actual separation part. The
flow in the pinched segment now contains both the inlet flow, in which the particles are
suspended, and the main flow. This flow is lead into a semicircle shaped chamber, that
contains a number of outlet channels. These outlet channels can either be identical or
nonidentical. In the following it is now assumed that the particles are perfectly aligned
up against the upper channel wall in the pinched segment. The geometry is rotated 180◦
from the previous chapter.
    The channel number convention is presented in Fig. 9.1. A term, which is called
separation position, will be used in the following. It is defined as, the distance from the
top sidewall to the streamline that terminates at the separation between channel 1 and 2,
as illustrated in Fig. 9.2. This position will act as a separation radius, so that particles
with radius’ below this value will end up in channel 1 and particles with a radius above
this value will end up in one of the remaining channels. In order to sort between red and
white blood cells a separation position around 4 µm is desired. The reason for this is that
the average radius of a red blood cell is 3.6 µm. But according to [3] is the behavior of
a particle dominated by its minimum length. The minimum length of a red blood cell is
its thickness of 2.2 µm. Therefore smaller separation position could be used if one was
interested in isolating the red blood cell. But in this device however it is the white blood
cell that should be isolated. Therefore it is preferred that some of the white blood cells
end up in channel 1 rather than if some of the red blood cells end up in channel 2.
    The Hagen–Poiseuille law, which was introduced in Chapter 2, will be used to deter-
mine the flow rates for the individual channels. Furthermore a width of 40 µm and a
height of 30 µm, for the pinched segment, will be used in the calculations.
    The theoretical treatment in this chapter will be made in full three dimension, which
ensures that the theory holds true in the actual device. Throughout this chapter the
theoretically determined separation positions will be supported by simulations made with
the quasi 3D model. These simulations are made on the model presented in Fig. 9.1. For
detailed information we reefer to Appendix C.
    In the following it will be assumed that the semicircle contains six outlet channels.
Identical channels means that both the size, shape and thereby their hydraulic resistance

                                            53
54                                               CHAPTER 9. SEPARATION MECHANISM



                 1

                           2



                               3



                               4
                                                             w   w’

                       5

                 6




Figure 9.1:     This is an illustration of the     Figure 9.2: To determine the separation po-
channel number convention. In the following        sition, in the figure denoted as w′ , a streamline
the top channel will be referred to as channel     is backtraced from the separation point between
1. Furthermore will the fluid that is predicted     channel 1 and 2. The separation position is de-
to end up in channel 1, be refereed to as the      fined as the distance from the upper side wall in
channel 1 flow.                                     the pinched segment to the separation stream-
                                                   line in the pinched segment. Particles with ra-
                                                   dius below this value will end up in channel 1.



is identical, as it is seen at geometry A in Fig. 9.4.
     In Chapter 5 it was derived that the particle motion is governed by the streamline
motion, in which a particle just follows the streamline going through its center of mass.
This particle motion will be used in the calculations through out this chapter. Therefore
it is of vital importance to investigate the separation of the fluid due to various geometry
changes. This will be the focus in this chapter.


9.1     Identical Outlet Channel
The simplest fluid motion would arise if the ”no slip” condition on the sidewalls could be
neglected. If this was the case the fluid would have a constant velocity throughout the
channel. A very simple relation between the streamline position in the pinched segment
and its position in the semicircle would then hold true. The two ratios stated below would
then be identical.
                                          w′    R′
                                             =                                        (9.1)
                                          w     R
In this equation w′ is the separation position, R represents the full angle of the semicircle
and R′ is the angle to the fluid streamline. These parameters is illustrated in Fig. 9.2.
    In an actual device the negligence of the ”no slip” condition is not a valid simplification.
Therefore Eq. (9.1) will not hold true. It neglects the fact that the flow rate per channel
9.1. IDENTICAL OUTLET CHANNEL                                                               55

width at the center of the channel is larger than it is at the side of the channel. When six
outlet channels are used they will receive one sixth of the total flow rate each.


    It will be assumed that the fluid will act as a creeping flow. A creeping flow means that
it will ”creep” around corners. In other words no rolls is excepted. The limited height of
the device contributes to the absence of these rolls, because for these rolls to arise a rather
easy flowing fluid is necessary, and this is exactly what the limited height of the channel
counteracts. With the assumption of creeping flow are the separation positions between
the six channels controlled by the flow profile in the pinched segment. In such a way that
the upper sixth of the flow rate will end up in channel 1. As described in Chapter 4 is the
flow profile for a channel with a rectangular shaped cross section given by:




                                                                    y
                            4h2 ∆p
                                      ∞
                                           1     cosh nπ
                                                         h                           z
                 vx (y, z) = 3                1−                            sin nπ       (9.2)
                             π ηL          n3    cosh nπ
                                                         w                           h
                                     n,odd
                                                                    2h




When the velocity profile is known, it is possible to determine the flow rate for a given
channel flow. The flow rate for the channel 1 flow is given by:




                                      −w/2+w′            h
                              Q=                dy           dz vx (y, z)                (9.3)
                                     −w/2            0




The separation position will be defined as the integration boundary that yields one sixth
of the total flow rate.


   The total flow rate is determined using Eq. (9.3) with w′ = w so that the integral is
56                                                          CHAPTER 9. SEPARATION MECHANISM

over the entire channel.




                   w/2              h
      Qtotal =            dy            dz vx (y, z)
                  −w/2          0
                                                                               y
                   w/2              h
                                          4h2 ∆p
                                                        ∞
                                                             1     cosh nπ
                                                                           h                   z
             =            dy            dz 3                    1−                  sin nπ
                  −w/2          0          π ηL              n3    cosh nπ
                                                                           w                   h
                                                       n,odd
                                                                              2h
                                                                        y                           h
              4h2 ∆p           w/2
                                         1
                                             ∞ cosh                 nπ           h
                                                                        h                 z
             = 3                dy          1−                                −    cos nπ
               π ηL        −w/2          n3    cosh                 nπ
                                                                       w        nπ        h
                                   n,odd                                                            0
                                                                       2h
                                                                        y
              4h2 ∆p
                           ∞
                                    1        w/2            cosh    nπ       2h
                                                                        h
             = 3                                   dy 1 −
               π ηL                 n3      −w/2            cosh    nπ
                                                                       w     nπ
                          n,odd
                                                                       2h
                                                              y       w/2
              8h3 ∆p
                           ∞
                                1      h sinh nπ
                                                 h
             = 4                   y−
               π ηL             n4    nπ cosh nπ
                                                  w
                          n,odd                                       −w/2
                                                               2h
                           ∞                                                               ∞
              8h3 ∆p                1     2h         w                                            1    π4
             = 4                       w−    tanh nπ                         using that              =
               π ηL                 n4    nπ         2h                                           n4   96
                          n,odd                                                           n,odd
                                           ∞
                 8h3 ∆p wπ 4                         2h                     w
             =               −                               tanh (2n − 1)π
                  π 4 ηL 96                      (2n − 1)5 π                2h
                                            1

                 8h3 ∆p    wπ 4
             ≈                          − 0.46535w                                                          (9.4)
                 π 4 ηL        96




In this calculation it was exploited that a sine and a hyperbolic cosine have uniform
convergence. Therefore it is allowed to change the order of summation and integration. In
the last step the summation is evaluated numerically, since it is not possible to determine
the last term in the sum analytically.


   As stated earlier is the separation position given as the integration boundary w′ in
Eq. (9.3), which result in a flow rate that is one sixth of the total flow rate. It is found
that a separation position of w′ = 9.88 µm results in this ratio. To verify this lets determine
9.2. NON IDENTICAL CHANNELS                                                                                    57

the channel 1 flow rate for this separation position.
                −w/2+w′                  h
       Q1 =                 dy               dz vx (y, z)
               −w/2                  0
                                                                                         y
                −w/2+w′                  h
                                               4h2 ∆p
                                                             ∞
                                                                  1     cosh nπ
                                                                                h                        z
          =                 dy               dz 3                    1−                         sin nπ
               −w/2                  0          π ηL              n3    cosh nπ
                                                                                w                        h
                                                            n,odd
                                                                                        2h
                                                                               y
           4h2 ∆p
                        ∞
                                1            −w/2+w′                cosh   nπ          2h
                                                                               h
          = 3                                          dy 1 −
            π ηL                n3       −w/2                      cosh    nπ
                                                                              w        nπ
                      n,odd
                                                                              2h
                                                                 nπ                             nπ
           8h3 ∆p
                        ∞
                                1        h sinh                     (−w/2 + w′ )       − sinh      (−w/2)
                                                                  h                              h
          = 4                      w′ −
            π ηL                n4      nπ                                 cosh nπ
                                                                                       w
                      n,odd
                                                                                       2h
              8h3 ∆p π 4 w′
          ≈                 − 1.59125w                            , for w′ = 9.88µm                          (9.5)
               π 4 ηL 96

To verify whether this separation position will bound one sixth of the total flow rate.
Channel Eq. (9.5) is divided by Eq. (9.4).

                      8h3 ∆p     π 4 w′
                       π 4 ηL     96          − 1.59125w
           Q1                                                                      1
                 =                                                ≈ 0.166664 ≈              , for w = 40µm   (9.6)
          Qtotal                                                                   6
                      8h3 ∆p     wπ 4
                       π 4 ηL     96         − 0.46535w


Because the infinite sum is replaced by a final sum the result is not exactly one sixth of
the total flow. But when more and more terms of the sum is included the result tends
toward one sixth. So if the pressure were the same at each outlet channel the separation
position would be w′ = 9.88 µm. This separation position would cause all particles with
diameters below 19.76 µm to end up in channel 1. Therefore six identical channels cannot
be used to separate red from white blood cell, since even the largest white blood cells have
radius below this value.
   In the simulation the separation position for six identical outlet channels were found
to be 9.39 µm. This is in good agrement with the theoretically determined separation
position of 9.88 µm.


9.2    Non Identical Channels
Since the earlier determined separation position is more than twice the desired separation
position of 4 µm, it is necessary to modify the design. This could be achieved by allowing
the channels to be nonidentical. This could either be a difference in the hydraulic resis-
tance across the channels, or a difference in the pressure drop across them. Before the
58                                            CHAPTER 9. SEPARATION MECHANISM

determinination of the optimal design parameters, the influence due to unwanted pres-
sure differences on the separation position will be investigated. This is done in order to
determine the amount of precision that is needed when controlling the outlet pressure.


9.2.1    Undesired Pressure Variation
In this section it is investigated how a pressure variation or rise in one of the outlet channel
influence on the separation position. It will further be considered whether it influences on
the result if the pressure variation exists in different channels.
    It is expected that if the pressure at some of the lower channels is raised, so that the
pressure drop across this channel decreases, the separation position would increase. This
is due to the Hagen–Poiseuille law which can be rearranged to yield the flow rate:

                                                ∆p
                                          Q=                                              (9.7)
                                                Rhyd

From this equation it is easily seen that when the pressure drop decreases, so does the
flow rate for the given channel. This decrement makes the flow rates in the other channels
increase, and therefore the separation position will increase according to this increment in
flow rate. We will expect that the flow rate variation will depend linearly with the flow
                                     1
rate, with the linearity coefficient Rhyd .
    In order to verify whether this expectation is valid. Some simulations have been per-
formed for various pressure differences at channel 6, for detailed information see Appen-
dix C. The results from these simulations is presented in Table 9.1, and plotted in Fig. 9.3.
In these simulations it is also investigated whether the hydraulic resistance of the channels
influence on the results. This is done by varying the height, and thereby the damping
factor, of a section in the outlet channels. From the results presented in Table 9.1 and in
Fig. 9.3 it is easily seen that thinner channels, and thereby larger hydraulic resistances,
limits the influence from the pressure variation. It is also seen that when the pressure
difference is zero the three simulation yields the same separation position. Therefore it
will be preferred, if it is possible, to use thin channels. In order to investigate whether
it changes the result if the unwanted pressure difference were at channel 4 rather than at
channel 6, some simulations are made, where it is investigated how pressure variances in
channel 4 changes the result. The results from these simulations are presented along with
the equivalent results from the channel 6 simulations in Table 9.2. From these results it is
seen that when the pressure variation are in channel 4 the influence on the result is more
significant. This is due to the fact that when there are pressure variations in the channels,
the streamline will not be as ideal as predicted. This could be due to the fact that there
are more buffer separation positions in the first example than in the second.
    But the important result from these simulations is that the separation position varia-
tion due to pressure variation can be limited. In the final design the hydraulic resistance
will be even larger than for the 3 µm high outlet channels due to the length of the channels.
So this variation will be of less importance.
9.2. NON IDENTICAL CHANNELS                                                                                    59



  Pres. diff. [Pa]                           Sep. pos. for 30µm [µm]   S.p. for 10µm [µm]   S.p. for 3µm [µm]
        -50                                           4.62                    5.50                7.90
        -40                                           5.80                    6.35                8.27
        -30                                           6.77                    7.20                8.56
        -20                                           7.77                    8.03                8.85
        -10                                           8.64                    8.75                9.12
         -5                                           9.03                    9.07                9.24
          0                                           9.39                    9.39                9.39
          5                                           9.67                    9.62                9.45
         10                                          10.15                    9.97                9.67
         20                                          10.79                   10.61                 9.92
         30                                          11.48                   11.19                10.22
         40                                          12.26                   11.85                10.50
         50                                          13.17                   12.58                10.79

Table 9.1: In this table the influence on the separation position due to pressure variations in
channel 6 is presented. The heights stated in the top line is the height of a damping section in all
the channels. A variation of this height changes the hydraulic resistance of the channels. It is seen
that the thinner the channels are, the smaller is the influence from the pressure variation.
                 Separation position [µm]




                                                                                       h = 30 µm
                                                                                       h = 10 µm
                                                                                        h = 3 µm




                                                          Pressure difference [Pa]
Figure 9.3: This is the separation position variation due to pressure variations differences in
channel 6 for various damping heights. It is seen that lower channels limits the separation position
variation.
60                                           CHAPTER 9. SEPARATION MECHANISM

 Pres. diff. [Pa]   Sep. pos. variation in channel 6 [µm]     S.p. variation in channel 4 [µm]
       -40                          5.80                                    4.31
       -30                          6.77                                    5.94
       -20                          7.77                                    7.2
       -10                          8.64                                    8.4
        -5                          9.03                                   8.93
         0                          9.25                                    9.4
         5                          9.67                                    9.86
        10                         10.15                                   10.34
        20                         10.79                                   11.26
        30                         11.48                                   12.12
        40                         12.26                                   13.14
        50                         13.17                                   14.08

Table 9.2: These simulations are made in order to investigate whether the unwanted pressure
difference in channel 6 or 4 changes the the result. It is seen that when the unwanted pressure
difference is in channel 4 the effect on the result is more significant.




9.2.2    Separation by Nonidentical Channels


Since it was not possible to separate the blood cells with identical channels, nonidentical
channels must be used, if we want to use 6 outlet channels.

   On Fig. 9.3 it is seen that when a lower pressure is present in channel 6, the separation
position is shifted upwards. This indicates that it is possible to manipulate the separation
position, so even though six symmetrically placed channels are used, a separation between
red- and white blood cells could be possible. From Eq. (9.7) it is seen that instead of
varying the pressure drop, one could vary the hydraulic resistance of the outlet channels
and still get the same flow rate. From this variation one could achieve the separation
position which is necessary for blood cell separation.

    In order to limit the channel 1 flow rate, it is seen form Eq. (9.7), that one would have
to increase the hydraulic resistance of this channel relative to the other. To determine
how large the difference in hydraulic resistance must be, it is necessary to reconsider the
flow rate in the pinched segment. Since the particle separation must take place between
channel 1 and 2 let us assume that these channels will have the same hydraulic resistance,
as it is seen in the B geometry Fig. 9.4. The velocity profile in the pinched segment is
independent of the arrangements of the outlet channels. The velocity profile stated in
Eq. (9.2), will be used to determine the channel 1 flow rate that will result in a separation
position around 4 µm. To determine this flow rate Eq. (9.3) is evaluated with w′ = 4 µm.
9.2. NON IDENTICAL CHANNELS                                                                     61




             A                                    B                                     C




Figure 9.4: This is an illustration of the three geometries that are used in this chapter. In A all
the outlet channels have the same length and thereby the same hydraulic resistance. In B channel
1 and 2 are seven times longer than the rest. And in C it is only channel 1 that is seven times
longer than the rest.



Eq. (9.3) is calculated analog to Eq. (9.5), and the flow rate yields:
                                                −w/2+4 µm            h
                                 Q4µm =                     dy           dz vx (y, z)
                                               −w/2              0

                                              8h3 ∆p π 4 4 µm
                                          ≈                   − 0.0831w                      (9.8)
                                               π 4 ηL    96

This flow rate is then divided by the total flow rate stated in Eq. (9.4) to give the channel
one flow rate:

                       8h3 ∆p     π 4 4 µm
                        π 4 ηL        96      − 0.0831w
            Q4µm
                   =                                        ≈ 0.033395 , for w = 40µm ⇔
            Qtotal
                        8h3 ∆p     π4 w
                         π 4 ηL     96    − 0.46535w

             Q4µm ≈ 0.033395 Qtotal                                                          (9.9)

Since we would like channel 1 and 2 to be identical, channel 2 will also receive 0.033395 of
the total flow rate. The rest of the outlet channels needs to split the remaining flow rate
if they are assumed to be identical. This results in a flow rate for each of the remaining
channels of:
                              1 − 2 × 0.033395
                   Qother =                       Qtotal ≈ 0.2333 Qtotal              (9.10)
                                      4
62                                           CHAPTER 9. SEPARATION MECHANISM

Since the flow rate of each channel is determined by its hydraulic resistance, it is possible
to determine what the ratio of resistance between channel 1 and 2 and the 4 lower channels
must be.
                            Qother 0.2333 Qtotal
                                  =              ≈ 6.997 ⇔
                            Q4µm 0.0333 Qtotal
                            Qother ≈ 6.997 Q4µm ⇔
                             ∆p             ∆p
                                   ≈ 6.997      ⇔
                            Rother         R4µm
                             R4µm ≈ 6.997 Rother                                      (9.11)

The hydraulic resistance of channel 1 and 2 should therefore be approximately seven times
the resistance of the remaining channels. In Chapter 4 the hydraulic resistance of a channel
having a rectangular cross section was found to be:

                                         ∞                            −1
                            12ηL                1 192h           w
                    Rhyd   = 3   1−                      tanh nπ                      (9.12)
                             h w                n5 π 5 w         2h
                                        n,odd


In this equation it is seen that the hydraulic resistance is directly proportional to the
channel length. This can be used to manipulate the hydraulic resistances for the individual
channels. A way of insuring that the hydraulic resistance of channel 1 and 2 is seven times
larger than the resistance of the others is therefore to make them seven times as long.
    The separation position between channel 2 and 3 must be determined, in order to make
sure that the white blood cells end up in channel 2. This is determined in the same way
as the original separation position. The separation position between channel 2 and 3 was
found to be at 5.85 µm. This will cause particles who’s diameter is between 8 µm and
11.7 µm to end up in channel 2.
    In the simulation a separation position between channel 1 and 2 of 3.82 µm were found,
and the separation position between channel 2 and 3 were determined to be at 5.59 µm.
Good agrement between theory and simulations are seen.
    It is seen that with this rather small separation position between channel 2 and 3 the
main part of the white blood cells will end up in channel 3 to 6 instead of in channel 2.
    Since the separation position between channel 1 and 2 still have to be around 4 µm it
is necessary to find another way to increase the separation position between channel 2 and
3. This can be done by allowing the hydraulic resistance of channel 2 to be lower than
that of channel 1. Therefore let us assume that the hydraulic resistance of channel 2 is
as low as that of channel 3 to 6. This lower hydraulic resistance will cause the channel 1
flow rate ratio to decrease, because when the hydraulic resistance of channel 2 is decreased
the flow rate for this channel will increase at the expend of all the other channels. When
channel 2 to 6 is seven times shorter than channel 1. The channel 1 flow rates ratio of the
total flow rate will given by:

                                                                            1
            Q1 + 5 × Qother = Q1 + 5 × 7Q1 = Qtotal         ⇔     Q1 =        Q       (9.13)
                                                                           36 total
9.3. CONCLUSION AND OPTIMAL DESIGN PARAMETERS                                                 63

                                        Geometry C              Geometry B
                                   Simulation Theory        Simulation Theory
              Channel 1 and 2        3.47µm     3.6µm          3.82     4.0µm
              Channel 2 and 3       11.54µm    11.8µm          5.59     5.9µm

Table 9.3: This table shows the result for both theory and simulation for particle separation
position between channel 1 and channel 2 and between channel 2 and 3. At geometry C channel
1 has a hydraulic resistance 7 times higher than that of channel 2 to 6. At geometry B channel 1
and 2 have a hydraulic resistance 7 times higher than that of channel 3 to 6. It is seen that the
theory and simulation match within a maximum of 5%.


If this design was used the separation position between channel 1 and channel 2 would
decrease to 3.62 µm. A device with this separation position would still be able to separate
white from red blood cell since the red blood cell have an average radius of 3.6 µm. It
will still work even though this separation position is close to the radius of a red blood
cell, because, as stated earlier in this chapter, is the motion of a particle dominated by
its smallest dimension, which for the red blood cell is its thickness of 2.2 µm. If this
design was used the separation position between channel 2 and 3 would be at 11.76 µm
and thereby it is insured that even the largest white blood cells end up in channel 2.
    When theses parameters are used in a simulation the separation position between
channel 1 and 2 is found to be at 3.47 µm, and 11.54 µm between channel 2 and 3, which
is in good agreement with the theoretically found value. For a schematic overview of the
simulated- and theoretically determined separation positions see Table 9.3.


9.3     Conclusion and Optimal Design Parameters
In Chapter 5 it was derived that the drag number was well below unity, for even the largest
particles, which indicates that a particle just will follow the streamline going through its
center of mass. If creeping flow and the exact alignment of the particles towards the upper
wall further more is assumed the expected particle position can be predicted with the use
of the above presented method.
    In the previous section it was found that when channel 2 to 6 were made seven times
shorter than channel 1 it was possible to separate the red from the white blood cells.
When these length were used the separation position between channel 1 and 2 were found
to be at 3.62 µm and the separation position between channel 2 and 3 were 11.76 µm.
With these two separation positions the red blood cells will end up in channel 1 and the
white blood cells in channel 2.
    It was found that when channels with relative high hydraulic resistances were used,
the pressure could vary more without changing the result substantially. But since the
height of the channels already are 30 µm, they cannot be made thinner due to the particle
sizes. Instead relative narrow channels can be used in order to get this damping effect.
Therefore long and narrow outlet channels are preferred.
    To make the separation mechanism as flexible as possible, the channels can be designed,
so that all channel will have the same length. There could then be placed a drain at
64                                           CHAPTER 9. SEPARATION MECHANISM

one seventh of the channel length in the six lower channels. In this way the separation
position depending on the numbers of open drains could vary between the lower limit,
where all channels except channel 1 were short, of 3.62 µm, and upper limit where the
hydraulic resistance of all the channels are identical, which results in a separation position
of 9.88 µm. If this design is used it is therefore possible to vary the separation position
between 3.62 µm and 9.88 µm.
    The method, which was presented above, used to determine the relative hydraulic
resistances can be automated, so that by dictating the numbers of outlet channels one
wants to use along with the desired separation positions the relative hydraulic resistances
can be determined.
    In Appendix D a maple worksheet is presented, that is capable from these information
to calculate the necessary relative hydraulic resistances between the channels.
    This worksheet is based on the method presented in this chapter. By defining the
height and width of the pinched segment the total flow rate is determined. Then the
channel 1 flow rate determined by the given separation position is calculated. In the
worksheet there are two ways of calculating the hydraulic resistances. In the first simpler
method the user can control the number of outlet channels and the separation position
between channel 1 and 2. The channel 1 flow rate is determined using Eq. (9.3). Then
is the channel 1 flow rate ratio of the total flow rate determined, and the remaining part
of the total flow is divided equally between the remaining channels. The relative flow
rate difference is found by dividing the flow rate ratio of the remaining channels with the
channel 1 flow rate ratio.
    In the advanced calculations one is able to control as many of the separation positions
as one wants. The procedure is the same only is the flow rate for each channel now
determined by:
                                        ′
                                  −w/2+wn             h
                          Qn =               dy           dz vx (y, z)dzdy             (9.14)
                                       ′
                                 −w/2+wn−1        0

If one chooses not to control all the separation positions then the remaining flow rate is
divided equally between the remaining channels as before. The resistances relative to the
hydraulic resistance of channel 1 is determined by dividing the individual channel 1 flow
rate ratios with the flow rate ratio of channel 1.


9.4     Comparison with other Models
In this section the model presented in this chapter will be applied to the results obtained
by Takagi et al. in [3]. This is done in order to investigate whether the method presented
in this thesis, is capable of predicting the separation position when non identical channels
are used.
    In [3] the geometry of the outlet channels, which is shown in Fig. 9.5, is defined in
such a way that the hydraulic resistance of the A channels is 48 times larger than that of
                                                      1
channel B. Therefore each A channel will receive 60 of the total flow and channel B will
receive 4 of it. In the article theory the separation positions were given by:
         5
9.4. COMPARISON WITH OTHER MODELS                                                           65




Figure 9.5: The geometry used by Takagi et al. in [3]. This is the geometry that the new
method of determining the separation positions presented in this chapter will be applied to, in
order to predict the separation positions correctly. The figure is from [3].


                                                        4
                                       ′     w × (1 −   5
                                      wN   =                N                           (9.15)
                                                 nB
In this equation nB is the number of outlet A channels. We would like to compare the
theoretically separation positions from the article with the separation positions which are
derived when the method presented in this thesis is used. In this method the separation
positions is given as the w′ that makes the following equation true.
                                                        ′
                                 4          h     −w/2+wN
                           1−              0 dz −w/2      dy vx (y, z)
                                 5N   =                                                 (9.16)
                            nB                 h    w/2
                                              0 dz −w/2 dy vx (y, z)

In Table 9.4 the separation positions between all the channels determined by both models
are presented. In [3] they expect that particles with a radius of 2.5 µm would exit through
channel A7. But that does not match their own theory, according to this the particles
should exit through channel A8, see Table 9.4.
    From Table 9.4 it is seen that with the use of the new method we expect that particles
with a radius of 2.5 µm end up at the boundary between channel A3 and A4, but on the
A3 side of the boundary.
    When they introduced particles into the system, particles with a diameter of 5 µm
ends up in channel A4 and very close to the boundary of A3. From this it is seen that
when the new method is used, it is possible to predict the particle separation position with
good precision.
    In [3] they introduced a ”scaling parameter” by measuring the flow rates in the individ-
ual outlet channels, with the use of small particles. they stated, from this measurement,
that the hydraulic resistance of the drain only were 18.9 times larger than that of the A
channels. When this hydraulic resistance was used in their theory the theoretical expec-
tation became correct.
66                                               CHAPTER 9. SEPARATION MECHANISM

 Separation position between         Sep. pos. (article theory)[µm]        Sep. pos. new theory [µm]
     Channel A1 and A2                            0.33                                1.43
     Channel A2 and A3                            0.67                                2.06
     Channel A3 and A4                            1.00                                2.57
     Channel A4 and A5                            1.33                                3.00
     Channel A5 and A6                            1.67                                3.39
     Channel A6 and A7                            2.00                                3.76
     Channel A7 and A8                            2.33                                4.10
     Channel A8 and A9                            2.67                                4.42
    Channel A9 and A10                            3.00                                4.73
    Channel A10 and A11                           3.33                                5.03
    Channel A11 and A12                           3.67                                5.32
     Channel A12 and B                            4.00                                5.60

Table 9.4: This is the analytically determined separation positions both with the use of the
article theory and with the model presented in this chapter. In collom two it is the separation
positions that arise when the article theory is used. In the right collom it is the separation positions
that arises when the new method is used. It is seen that these are substantially larger.



    If the theory presented in the articles which is equivalent to Eq. (9.1) were applied to
the geometries in the previous sections. The resulting expected separation position would
be off from the simulations by 29 %, 67 % and 69 % for design A, B and C respectively.
The method used in this thesis were off by 5.24 %, 4.7 % and 4.3 % respectively.
    Therefore we are a bit puzzled with the use of this correction factor, which makes
theory and experiment fit perfectly. Especially when the theory presented in this thesis is
able better to predict the separation position, without the use of this scaling parameter.
Chapter 10

Consideration Regarding the Final
Design

In this chapter some consideration regarding the final device , such as the particle motion
around corners, narrowing of the pinched thickness and influence from from the height of
the channels will be presented, and the optimal design parameters for blood cell separation
will be presented.


10.1     Motion of Particles Around Corners
In the calculations and simulations so far particles have been treated by a concentration
profile or streamline motion. This treatment is not sufficient for all purposes. For instance
does the way of treating the particles not respect the fact that the particles have a finite
size.
    An effect that needs to be investigated is how the fluid flow pasts a corner in 3D. It
will be preferred if the streamline behavior didn’t depend on their respective z position
in the channel. Because otherwise this could change the results found in the quasi 3D
simulations. In Fig. 10.1 a result from a full 3D simulation showing streamline behavior
around a corner is shown. From this figure it seen that it is only the streamline that is
closest to the walls that shows some abnormalities compared to the others. But because
of the finite particle sizes the streamlines closest to the wall will be of no importance.
    Due to this smooth flow of the streamlines around the corners in 3D, we will expect
that it will not change the result substantially if the simulation were performed in 3D.
Therefore the parameters determined in the quasi 3D simulation could be used in the
actual device.
    In the final design the sharpest corner that the particles have to pass by is a 90◦ corner.
Therefore it is the streamlines around this corner that is investigated, because when the
corner becomes less sharp, the streamline will tend to be more smooth.
    Another effect of the particles passing past a corner that must be taking into account.
Is whether the streamlines at outlet of the pinched segment comes closer to the corner
than their separation position in the pinched segment. This must be investigated because

                                             67
68             CHAPTER 10. CONSIDERATION REGARDING THE FINAL DESIGN




Figure 10.1: Fluid streamline for the flow around a corner in 3D. It is seen that except for
the streamlines closest to the wall, their behavior is independent on their internal position in the
channel.
10.2. NARROWING OF THE PINCHED SEGMENT                                                        69




                                                                          w0

                       w1                                             ′
                             ′                                       w0
                            w1

Figure 10.2: This is a schematic illustration of a narrowing process. Even during this narrowing
process streamlines will still not cross each other. Therefore the narrowing process in a channel
will contribute to the particle alignment towards the channel wall




if this is the case the effective separation radius would increase since particles cannot move
closer to the wall than their radius.
     In order to investigate the effect from this corner motion a simulation have been per-
formed on an actual separation example. From this it is seen that when a separation
position of 4.78 µm is found in the pinched segment, the streamline actually flow past the
corner at a distance of 3.52 µm This is a simulation of an extreme situation. From this
it seen that the actual separation position will be larger than the one determined in the
calculations, but as long as somewhat realistic circumstances is used it will not be more
than 1 µm bigger.
     The angle between the inlet- and main channel is of less importance for the functionality
of the device. Since it is seen in the simulation that the velocity profile quickly established
it will not make any difference on the result whether the angle between the inlet- and
main channel were 45◦ or 90◦. Therefore it will not change the result when the inlet and
main channels in the final device is arranged in another way than in the calculations and
simulations.



10.2      Narrowing of the Pinched Segment

In Chapter 6 and Chapter 8 the optimal design parameters regarding the alignment of
particles were determined. In Chapter 8 it was seen that the position of the maximum
gradient, and thereby the intersection width for the main part of the flow rate ratio was
larger than the diameter of a red blood cell. When this is the case we cannot be sure that
the red blood cells are aligned against the wall.
    Therefore in order to be able to use elevated flow rate ratios it will be necessary to
manipulate the width of the pinched segment.
    In Chapter 5 it was found that a particle follows the streamline going through its
center of mass. In this frame it is straight forward to determine the narrowing effect. Let
us consider the narrowing segment which is shown in Fig. 10.2. The narrowing effect is
derived in the below calculations. It is assumed that the flow rate below the intersection
70             CHAPTER 10. CONSIDERATION REGARDING THE FINAL DESIGN

streamline is constant, which is equivalent with the demand that streamlines cannot cross.
                                         Q1     Q0
                                              =       ⇔
                                        Qtotal Qtotal
                         w1′                           w′ 0
                         0     dy k1 (w1 − y)y        0     dy k0 (w0 − y)y
                         w1                    =       w0                     ⇔
                         0     dy k1 (w1 − y)y         0 dy k0 (w0 − y)y
                               1           1         1            1
                                 w w′ 2 − w′ 3         w w′ 2 − w′ 3
                               2 1 1       3 1   =   2 0 0        3 0 ⇔
                                 1 3     1 3           1 3      31
                                  w − w1                w    − w0
                                 2 1     3             2 0       3
                                                     w′ 1
                                           w1 =           · w0                             (10.1)
                                                     w′ 0
From this equation it is easy to determine which pinch segment width that should be used,
in order to decrease the intersection width to a given value. If for instance one wants to
have a intersection width of 4 µm and the intersection width was 8 µm in a 40 µm wide
channel the outlet channel should be narrowed to a width of 20 µm.
    During this narrowing the particles which already are aligned towards wall will expe-
rience a pressure up against the wall, which could cause them to be squeezed. But we will
expect that even though the blood cells are not solid as described in Chapter 1 they will
be able to withstand the pressure and remain more or less spherical in shape.


10.3      Channel Height
Another parameter that needs to be considered before designing the actual device is how
high the channels should be. The lower limit of the height is that it should be larger than
the diameter of the largest particles. In order to have a certain safety margin channel
depth below 30 µm will not be considered when separating blood cells. A difference in
channel depth will influence on the flow profile. In Chapter 4 it was described that when
thinner channels are used, the velocity profile will be limited, in contrast to the velocity
profile for infinitely thick channels.
    In this device it is preferred to use rather low channels since it is of vital importance, for
the functionality of the device, that it is possible to control the position of the particles
with great precision. In higher channels the vertical position of the particles is harder
to control. In lower channels the vertical particle position is not able to change simply
because there is no place for the particles to migrate. Therefore a channel hight of around
20 µm to 30 µm is preferred in the final device. Since the largest white blood cell can be
up to 20 µm, 30 µm high channels will be used.
    When working with diffusion in channel with a finite thickness the so called butterfly
effect needs to be considered. This effect and its influence on the diffusion profile is
described by A. E. Kamhols et al. [12].
    When diffusion takes place in channels with a limited height the velocity profile due to
the hight of the channels will influence on the concentration profile. The actual diffusion
process will not change, it is just the circumstances around it which leads to the non
constant diffusion profile. The vertical variation of the velocity leads to the fact that
10.4. SEPARATING WHITE- FROM RED BLOOD CELLS                                                                                           71




                                                                                                     Concentration [arbitrary units]
                                                                        velocity [arbitrary units]
Figure 10.3: This figure shows the concentration profile in a channel having a rectangular cross
section due to the butterfly effect. Due to the higher fluid velocity at the center of the channel the
concentration profile becomes hour-glass shaped.


the fluid at the center of the channel will have a higher velocity than the fluid at the
top- and bottom walls. This velocity difference will result in a change in the equilibrium
concentration profile. As the particle via diffusion moves towards the center of the channel
the higher velocity at this position will force the particles to spend less time in the channel,
this will again reduce the effective diffusion time. Due to this effect the diffusion profile
will be wider at the bottom and the top of the channel, where the velocity difference
towards the horizontal center of the channel will be less than the difference at the center
of height. This will then result in a hourglass chapped concentration profile, see Fig. 10.3
According to [12] this effects minimizes for low channels. This effect therefore implies as
well that channels with heights of 20 µm to 30 µm should be used in the final device.


10.4      Separating White- from Red Blood Cells
From all the above presented optimal parameters it is now possible to state the optimal
design parameter for PFF system that is made in order to sort between white- and red
blood cells.
    Since the largest white blood cells can be as large as op to 20 µm the channels cannot
be made lower or narrower than this hight. In order to have a certain safety margin, so
that clogging of the device can be prevented a constant channel height of 30 µm must be
used throughout the device, since lower channel are preferred.
    In Chapter 6 an alignment demand for the flow rates was derived. From this inequality
it was found that in order to insure alignment of red blood cells in a pinched segment of
72            CHAPTER 10. CONSIDERATION REGARDING THE FINAL DESIGN

width 40 µm a flow rate ratio of 0.1 must be used.
    In order to separate the particles it is a question of temper about the number op
channels that should be used. But if we decide on the use of six outlet channels it was
found that it was possible to sort between the cells so that the red- and white blood cells
leaves through channel 1 and 2 respectively. This could be done by making channel 1
seven times longer than channel 2 to 6. Thereby it will also have a hydraulic resistance
that is seven times larger than the others. In the final device a simpler geometry can be
used, since there no need for changing the separation position, when the optimal value has
been determined. For the simplest case only 2 outlet channels can be used, a channel for
red blood cells and one for the rest of the particles. But since the white blood cells after
the pinched segment needs to be captured it will be preferred that the flow rate, in which
the white blood cells is suspended in after the pinched segment, not is to large. Therefore
three outlet channels is preferred, one for red- and the smallest of the white blood cells,
one for ´the largest white blood cells and a drain. As an example consider a 3 outlet
system with a pinched segment of width 30 µm, which requires a flow rate ratio of 0.15
for alignment, and separation positions of 4 µm and 6 µm on must design the channels so
that the hydraulic resistance of channel 2 is 1.1 times larger than that of channel and the
hydraulic resistance of the drain channel should be 16.7 times larger than that of channel
1.
    In the final design the outlet channels should all be made so that they are rather long
and narrow in order to avoid unwanted effects due to pressure variations.
Chapter 11

Designing a Pinched Flow

In this chapter the optimal design parameters, that have been determined throughout the
thesis in both the theoretical calculation and the numerical simulations, will be presented.
    This chapter is made so that it can be understood by a person who has not read the
rest of the thesis. Therefore new terms introduced in this section will be explained.
    First parameters concerning the alignment of the particles in the pinched segment
will be commented. Then parameters concerning the actual sorting mechanism will be
presented. At the end of the chapter a short example of the use of this design guide is
included.


11.1     Definition of Parameters
Some of the terms used in this chapter are presented in Fig. 11.1. The basic geometry
of a Pinched Flow Fractionation (PFF) system is, that two flows are introduced into a
pinched segment, which is a short narrow channel. The particles are introduced into one
of the channels and this is therefore named the inlet channel. The other channel is called
the main channel. In the pinched segment the flow from the inlet channel must be forced
up against one of the sidewalls by the larger flow coming from the main channel. This
insures that the particles are aligned towards the channel wall.
    After the pinched segment, which should be made rather short for instance 100 µm,
the two flows are let into a semicircle shaped chamber that contains a number of outlets.
When the particles are aligned in the pinched segment they will be sorted into these
channels according to their size. The pinched segment must be short, because otherwise
diffusion would prevent the alignment of the particles. Aligned particles are particles that
are in contact with the channel sidewalls.
    The separation position is the most important parameter in the device. It is defined
as the distance from the channel side wall, unto which the particles are aligned, to the
position of the streamline that terminates at the boundary between two channels. If
nothing else is mentioned it is the separation position between channel 1 and 2 that is
used. A streamline is defined as the path a small volume of fluid will take though the
system. Particles with a radius below the separation position between channel 1 and 2

                                            73
74                                       CHAPTER 11. DESIGNING A PINCHED FLOW




                                    Channel 1

                                                                       Channel 2
                Particle alignment

       Inlet channel
                       Q
                        in
                          le
                               tt

     Main channel
                               n
                             ai
                            m
                           Q




                                                                      Separation chamber
                 Pinched segment



Figure 11.1: This is the basic geometry of a pinched flow. With the use of different velocities
in the inlet- and the main channel, the particles are aligned up against the upper side wall in the
pinched segment. In the separation chamber they are sorted between the outlet channels according
to their sizes. In this figure separation between red and white blood cells is illustrated.
11.2. ALIGNING THE PARTICLES                                                               75

will end up in channel 1.


11.2     Aligning the Particles
In Chapter 6 it was derived, that in order to insure particle alignment towards the channel
sidewall in the pinched segment, one must meet the requirement stated in the inequality
below.
                          Qinlet        1       1                   2r
                                      ≤ β2 − β3,           for β =                    (11.1)
                      Qmain + Qinlet    2       3                   w
Here Qinlet and Qmain are the flow rates from the inlet channel- and main channel respec-
tively, r is the radius of the smallest particle that must be aligned and w is the width of
the pinched segment.
    If the flow rates cannot be changed enough to fulfill the the inequality stated Eq. (11.1),
one must use another method in order to insure particle alignment. This method is to
narrow the pinched segment, because when w is decreased, β increases which also makes
the right hand side of Eq. (11.1) increase, since it is not physically possible for β to be
above unity. Thereby particle alignment can still be insured.
    In the inequality presented in Eq. (11.1) diffusion has been neglected. When this is not
a valid simplification more complicated methods to determine the flow rate ratio required
for particle alignment must be used.
    The requirement stated in Eq. (11.1) is valid as long as the diffusion constant of the
smallest particles is below 1 × 10−13 m2 /s, which covers particles with a radius above
0.2 µm suspended in the Phosphate Buffered Saline buffer.
    In Fig. 11.2 an illustration of the alignment of two particles are shown. The particle
denoted A has a radius above the separation position value and will therefore not end up
in channel 1. But the particle denoted B will end up in channel 1 because its radius is
below the separation position value.
    Due to diffusion the two flowrates will not remain separated. The Full Width Half
Max (FWHM) of the concentration gradient is a measure for the width of the mixed layer
between the two flows. From Chapter 8 it was found that in order to get the smallest
FWHM rather large flow rate ratios, which is the inlet channel flow rate divided by the
main channel flow rate, should be used. this not necessarily the optimal design parameter,
since the intersection width (w′ ), which is the width of the particle containing flow at the
end of the pinched segment, increases with growing flow rate ratios.
    When the diffusion constant exceeds 1 × 10−13 m2 /s diffusion needs to be taken into
account. One must then make simulations on the simplified model, which is presented in
Chapter 8. But as long as a diffusion constant below 1×10−13 m2 /s are used Eq. (11.1) can
be used to determine the necessary flow rate ratios. Instead of calculating the necessary
flow rate ratios from Eq. (11.1) one can look at Fig. 11.3, which is a plot of the intersection
                                                                                   2r
width ratio, that must be used for alignment, of the pinched segment width              versus
                                                                                    w
the total flow rate ratio.
                                          2r
    When one has decided on a ratio of , in other words if one has decided which width
                                           w
he or she wants to make the pinched segment, it is possible from this figure to determine
76                                       CHAPTER 11. DESIGNING A PINCHED FLOW




                             w w′       B
                                               A




Figure 11.2: This is an illustration of two particles that are aligned towards the upper channel
wall. The streamline representing the separation position is also shown. It is seen that particle A
has a radius above the separation position will end up in channel 2. Particle B has a radius below
it will end up in channel 1.



the maximum flow rate ratio that can be used when one still wants to be sure that the
particles are aligned. Of course it is allowed to use smaller flow rate ratios than this value,
since this would just insure the alignment even further. But flow rate ratios below 0.02 is
not recommended since this small flow rate ratios will force the diffusion process far into
the inlet channel, which will increase the velocity complexity. The position of the particles
then becommes more difficult to control.
    If it is not possible to achieve the alignment of particles with flow rate ratios above
0.02 one must choose a narrower pinched segment in order to get a larger intersection
width ratio which would allow larger flow rate ratios to insure the alignment.


11.3      Separating the Particles
In Chapter 9 it was seen that even when the number of outlet channels has been chosen,
it is possible to manipulate the separation position by varying the hydraulic resistance of
the outlet channels.
     In order to determine the relative difference in hydraulic resistance that should be used
in the design of the device, some computational methods must be used. To help with these
computations we refer to the maple worksheet which is shown in Appendix D.
     In this worksheets one has two options. Either one can use the standard calculations
or the advanced calculations. The standard calculations are capable, when the number of
outlet channels have been chosen, to determine the hydraulic resistances, that results in
the desired separation position between channel 1 and 2. In the advanced calculation it is
possible to control as many separation positions one wants to.
11.4. OVERALL DESIGN PARAMETERS                                                               77




         2r
         w
         Intersection width ratio




                                                      Qinlet
                                    Flow rate ratio   Qmain

Figure 11.3: This figure is used to determine the flow rate ratios between the inlet- and main
flow rate, that will insure the alignment of particles with radius r in a pinched segment of width
w.



    A way of changing the hydraulic resistances of the outlet channels is to change their
length scales. This is because the hydraulic resistance is directly proportional to the
length of the channels. So by doubling the length of the channel one doubles the hydraulic
resistance.
    A way of making the device more flexible when designing the outlet channels is to
make all the channels the same length as a starting point. Then they can be made shorter
by introducing a drain at the appropriate positions down the outlet channels. Then one
during the experiment can change which outlet of the channel one wants to use and thereby
manipulate with the hydraulic resistance of the channel. This variation will change the
separation position as well.


11.4          Overall Design Parameters
In Chapter 9 it was found that in order to minimize the effects due to unwanted pressure
variation the channels should have as high a hydraulic resistance as possible.
    This increase in hydraulic resistance can be achieved by changing the geometry of
the channels, that is to make them thinner, narrower and longer. So when designing a
78                                      CHAPTER 11. DESIGNING A PINCHED FLOW

             20 µm
                               Inlet channel




Pinched segment width

                                               100 µm


                                     Main channel

                      40 µm
Figure 11.4: This is the basic geometry that is used as an example of the use of the design guide
in practise. The width of the pinched segment will be determined in the calculations.



separation mechanism keep in mind that long narrow channels limit the effects due to
unwanted pressure variations.
    In a pinched flow one would like the flow to be a creeping flow in the separation chamber
in order to avoid rolls that would interfere with the sorting mechanism. A creeping flow
is a flow that creeps around corners and thereby has no rolls.
    A way of insuring that the flow will be almost like creeping flow is to use rather thin
channels throughout the device. This is because the limited height will act as a damping
factor on the fluid velocity, due to the ”no-slip” boundary condition on the top and bottom
of the channels. This damping factor will make the fluid less easy flowing, and thereby
limit the rolls in for instance the separation chamber. Therefore use as thin channels as
possible, without clogging the device.


11.5      Use of the Design Guide
In this section the above presented design guide will be used on a specific example.
    In this example we would like to separate particles with radius of 2 µm and particles
with a radius of 4 µm. Therefore we would like to achieve a separation position of 3 µm.
We would like to achieve this with an inlet channel width of 20 µm and a main channel
width of 40 µm. Due to instrumental limitation on the flow rates, a flow rate ratio below
0.15 cannot be used. In the actual separation mechanism we would like to use six outlet
channels, and with only the top channel having a change in hydraulic resistance. The
geometry is shown in Fig. 11.4.
    We would like to use a pinched segment with a width of 40 µm. If this width is used the
                           2 × 2 µm
intersection ratio will be          = 0.1. When this intersection ratio is used in Fig. 11.5,
                            40 µm
11.5. USE OF THE DESIGN GUIDE                                                              79




        w′
        wIntersection width ratio




                                                       Q1
                                    Flow rate ratio    Q2
Figure 11.5: This is a normalized plot of the separation position ratio as function of the flow
rate ratio. It can be used to find the required ratios to design a optimal system.



the red dashed line, it is seen that a very low flow rate ratio must be used in order to insure
particle alignment. But since this is not possible the pinched segment must be narrowed.
From Fig. 11.5 it is seen that the optimization also can be used in the opposite direction.
From this it is obvious that an intersection width ratio of 0.23 should be used when the
flow rate ratio is 0.15, the blue dashed line. Therefore particle alignment can be achieved
                                                    3 µm
by narrowing the pinched segment to a width of            = 13 µm.
                                                     0.23
    With the use of these parameters the particles are now aligned towards the upper
wall in the pinched segment. We are now ready to focus on the actual separation of the
particles. As stated above we would like to use six channels, and only the top channel
should have a difference in hydraulic resistance relative to the others. To calculate the
precise hydraulic resistances, that should be used to achieve this, we have used the maple
worksheet that has been included in Appendix D. It is the simple calculation that has
been used. And the numbers that are inserted in the appendix are the numbers used to
determine the hydraulic resistances for this particular example. It was found that the
separation position of 3 µm can be achieved by making the hydraulic resistance of channel
1, 1.15 times larger than that of channel 2 to 6.
    Since minimal channel height is preferred we have chosen a channel height of 10 µm
in the pinched segment so that clogging is still prevented.
80                                   CHAPTER 11. DESIGNING A PINCHED FLOW

    To sum up: To achieve the separation between particles with a radius of 2 µm and
4 µm a inlet channel width of 20 µm, a main channel width of 40 µm and a pinched
segment width of 13 µm can be used along with a flow rate ratio of 0.15. The separation
position of 3 µm can be achieved by six outlet channels with channel 1 having a hydraulic
resistance that is 1.15 times larger than that of channel 2 to 6.
Chapter 12

Conclusion

In this thesis we have focused on the development of a new method for the operation of a
Pinched Flow Fractionation (PFF) device. Through this development the optimal design
parameters regarding the separation of human blood cells have been determined. This
theory and optimal parameters have been derived using both analytical calculations and
numerical simulations.
    So far, in the present literature on PFF, there has not been presented a valid theory
for various critical parameters in the pinched flow, such as particle alignment demand and
particle separation position. We have derived a theoretical model that predict some of
the critical parameters. The theoretical demand for the alignment of particles have been
modified, so that it takes the ”no-slip” boundary condition on the channel sidewalls into
account. It was found that the fluid behavior was controlled mainly by the relative flow
rates between the inlet channels.
    Until now the published models have not been able to predict the particle separation
positions within the pinched segment. We have in this thesis presented a model that is
able to predict the separation positions in the pinched segment even when nonidentical
channels are used. This model have been supported by the prediction of the separation
position found in [3], which was not possible by the original theory. The developed model
showed also good agreement with numerical simulations.
    Main focus of this thesis was the determination of optimal design parameters for a
Pinched Flow Fractionation system that should separate red- from white blood cells. These
optimal design parameters are presented along with a description on how the design could
be changed in order to separate between other particle sizes in a design guide. The
intention with the design guide is that any person who wants to fabricate a PFF system
can read it and from this derive the parameters needed to get the separation position he
or she wants.
    Unfortunately the device has not yet been tested by the practical part of the project.
Therefore we cannot comment on the results.




                                           81
82   CHAPTER 12. CONCLUSION
Acknowledgements

At first we would like to give a special thanks to our supervisors: Winnie Edith Svend-
sen, Maria Dimaki & Fridolin Okkels. Especially Fridolin for the cosy weekly meetings
throughout this project. We would also like to thank the rest of the NaBIS group at MIC,
for inspiring group meetings and fact finding assistance.
                                           u
    We will as well like to thank Camilla B¨low Christiansen for proof-reading of the thesis.
    At last we would like to thank Nikolaj Ormstrup Christiansen and Martin Græsvænge
Hansen, which have performed the practical part of the Pinched Flow Fractionation
project, for great and inspiring co-operation throughout the period.




                                             83
84   CHAPTER 12. CONCLUSION
Bibliography

[1] Biofluid Mechanics
    J. N. Mazunder,
    1. edition 1992

[2] Computational Fluid Dynamics in Microfuidic Systems
    L. H. Olesen,
    Master Thesis, Mikroelektronik Centret (MIC) Technical University of Denmark
    31 July 2003

[3] Continuos particle separation in a microchannel having asymmetrically arranged
    multiple branches
    J. Takagi, M. Yamada, M. Yasuda & M. Seki,
    Lab Chip, 2005, 5, 778-784

[4] Continuous separation of particles using a microfluidic device equipped with flow rate
    control valves
    Y. Sai, M. Yamada, M. Yasuda & M. Seki,
    Juonal of Chromatrography A, 1127 (2006) 214-220

[5] Fluid Mechanics
    L.D. Landau & E. M. Lifshitz,
    Second edition, 2003

[6] Inertial migration of rigid spheres in two-dimensional unidirectional flows
    B. P. Ho & L. G. Leal,
    J. Fluid. Mech. (1974), vol. 65, part 2, pp. 365-400

[7] Mathematics for physicists
    S. M. Lea,
    first edition, 2004

[8] Physical hydrodynamics
    E. Guyon, J. P. Hulin, L. Petit and C. D. Mitescu
    First edition, 2001




                                         85
86                                                                BIBLIOGRAPHY

 [9] Pinched Flow Fractionation: Continuous Size Separation of Particles Utilizing a
     Laminar Flow Profile in a Pinched Microchannel
     M. Yamada, M. Nakashima & M. Seki,
     Anal. Chem. 2004, 76, 5465-5471

[10] Pinched Flow Fractionation for Blood Cell Separation
     M. G. Hansen & N. O. Christiansen
     Bachelor Thesis, Department of Micro- and Nanotechnology (MIC) Technical Uni-
     versity of Denmark
     To be handed in at the 25. of June 2007

[11] The Physics Factbook
     G. Elart et al.
     Web page: http://hypertextbook.com/facts/2004/MichaelShmukler.shtml

[12] Theoretical Analysis of Molecular Diffusion in Pressure-Driven Laminar Flow in
     Microfluidic Channels
     A. E. Kamholz & P. Yagar,
     Biophysical journal, January 2001, vol. 80, 155-160

[13] Theoretical microfluidics
     H. Bruus,
     Second edition, fall 2005
Appendix A

2D Script File

%%%%%%%%%%%%%%%%%%%%%%
%     Model 2D      %
% %%%%%%%%%%%%%%%%%%%%
%
%Simulation information:
%
%Mesh B, 5906 elements for a inlet channel width of 20E-6 m
%Peclet number - PE = 3000
%When simulation for different widths velocities are changed according
%to the width of the channel and the scaling factor

%COMSOL Multiphysics Model M-file
% Generated by COMSOL 3.3 (COMSOL 3.3.0.405, $Date: 2007/06/21 08:05:46 $)

clear all; clc; flclear fem

% COMSOL version
clear vrsn vrsn.name = ’COMSOL 3.3’; vrsn.ext = ’’; vrsn.major = 0;
vrsn.build = 405; vrsn.rcs = ’$Name: $’; vrsn.date = ’$Date:
2007/06/08 07:26:36 $’; fem.version = vrsn;

%Number of widths used in the simulation
nr=6;
%Starting width
ns=2;
%Distance between widths
t=10e-6;
%Number of Scaling factors
nrUSf=100;
%Distance between the scaling factors

                                    87
88                                           APPENDIX A. 2D SCRIPT FILE

tUSf=0.005;
% Pe number
Pe=3000;
%Diffusivity for the cells
D = 1e-13;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The for loop runs 1 time at each width
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for w = ns:1:nr
    width=t*w; % The Width of the side channel is defined
    % Geometry. The base point of the inle channel is defined so the oulet
    % channel will have a length of 100e-6 m regardless of the inlet
    % channel width.
    g1=rect2(’180e-6’,’40e-6’,’base’,’corner’,’pos’,{’-110e-6’,’0’},’rot’,’0’);
    g2=rect2(width,’40e-6’,’base’,’corner’,’pos’,{’-10e-6’,’-40e-6’},’rot’,’0’);

     % Analyzed geometry
     clear s
     s.objs={g1,g2};
     s.name={’R1’,’R2’};
     s.tags={’g1’,’g2’};

     fem.draw=struct(’s’,s);
     fem.geom=geomcsg(fem);

   % Initialize mesh B
fem.mesh=meshinit(fem, ...
                   ’hmax’,[6e-6], ...
                  ’hgrad’,1.1, ...
                  ’hmaxvtx’,[4,3e-6,6,0.5e-6], ...
                  ’hmaxedg’,[2,1e-6,6,1e-6]);



     % Descriptions
     clear descr
     descr.const= {’D’,’Diffusion constant’,’Us’,’Side flow’,’Ui’,’Inflow’};
     fem.descr = descr;

     % (Default values are not included)

     % Application mode 1 (Navier-Stokes model)
     clear appl
     appl.mode.class = ’FlNavierStokes’;
                                                      89

appl.gporder = {4,2};
appl.cporder = {2,1};
appl.assignsuffix = ’_ns’;
clear prop
prop.analysis=’static’;
appl.prop = prop;
clear bnd
bnd.u0 = {0,0,0,’-4*Ui*s*(1-s)’,0};
bnd.type = {’noslip’,’neutral’,’strout’,’uv’,’uv’};
bnd.v0 = {0,0,0,0,’4*Us*s*(1-s)’};
bnd.ind = [3,1,1,1,5,2,1,1,4];
appl.bnd = bnd;
clear equ
equ.rho = 1e3;
equ.eta = 1e-3;
equ.cporder = {{1;1;2}};
equ.gporder = {{1;1;2}};
equ.ind = [1,1];
appl.equ = equ;
fem.appl{1} = appl;

% Application mode 2 (Convection Diffusion model)
clear appl
appl.mode.class = ’FlConvDiff’;
appl.assignsuffix = ’_cd’;
clear prop
prop.analysis=’static’;
clear weakconstr
weakconstr.value = ’off’;
weakconstr.dim = {’lm3’};
prop.weakconstr = weakconstr;
appl.prop = prop;
clear bnd
bnd.type = {’N0’,’cont’,’Nc’,’C’,’C’};
bnd.c0 = {0,0,0,0,1};
bnd.ind = [3,1,1,1,5,2,1,1,4];
appl.bnd = bnd;
clear equ
equ.D = ’D’;
equ.v = ’v’;
equ.u = ’u’;
equ.ind = [1,1];
appl.equ = equ;
fem.appl{2} = appl;
90                                            APPENDIX A. 2D SCRIPT FILE

     fem.frame = {’ref’};
     fem.border = 1;
     clear units;
     units.basesystem = ’SI’;
     fem.units = units;

     % Multiphysics
     fem=multiphysics(fem);

     % Extend mesh
     fem.xmesh=meshextend(fem);

     %Defining empty vectors
     maxi = zeros(1,2);
     usfvec = zeros(1,2);
     xposvec = zeros(1,2);
     FWHMvec = zeros(1,2);

     % This loops runs one time for each scaling factor
     for i=1:1:nrUSf
         Usf = tUSf*i;
        % Claculate Velocity in main channel
         Ui=(Pe*D*(40e-6+width))/(40e-6*(40e-6+Usf*width));
         usfvec(1,i)=Usf;

        % Definition of constants
         fem.const = {’D’,D, ...
             ’Ui’,Ui, ...
             ’Us’,Usf*Ui};



         disp(sprintf(’Ui: %11.3e   Pe: %9.3f Usf: %9.3f’,Ui,Pe,Usf))

         % Solve problem (Navier-Stokes)
         fem.sol=femstatic(fem, ...
             ’solcomp’,{’u’,’p’,’v’}, ...
             ’outcomp’,{’u’,’p’,’v’});

         % Solve problem (Convection diffusion)
         fem.sol=femstatic(fem, ...
             ’init’,fem.sol, ...
             ’solcomp’,{’c’}, ...
             ’outcomp’,{’c’,’u’,’p’,’v’});
                                                                  91

figure(1)
subplot(2,1,1)
% Plot solution
postplot(fem, ...
    ’tridata’,{’U_ns’,’cont’,’internal’,’unit’,’m/s’}, ...
    ’trimap’,’jet(1024)’, ...
    ’title’,’Surface: Velocity field [m/s]’);
axis equal tight
% Plot solution
subplot(2,1,2)
postplot(fem, ...
    ’tridata’,{’c’,’cont’,’internal’,’unit’,’mol/m^3’}, ...
    ’trimap’,’jet(1024)’, ...
    ’title’,’Surface: Concentration, c [mol/m^3]’);
axis equal tight;

% k set to 1 for messuring the FVHM at the end of the channel
k=1;
% tvarsnit set to -100e-6
tvarsnit(k)=-100e-6*k;
% np is the number of division of the y coordinate in the cross
% sections
np = 100;
px =-100e-6*k*ones(1,np); % prefactor changed to -100e-6
py = linspace(0,40e-6,np);
%%%%%%%%%
% Finding maximum contration gradient at main channel end
%%%%%%%%%
carr(:,k)=postinterp(fem,’c’,[px;py]);
dcarr(:,k) = postinterp(fem,’grad_c_cd’,[px;py]);
[maximum,ypos] = max(dcarr(:,k));
yposvec(i,k)=ypos*40e-8;
%%%%%%%%
% Finding FWHM of contration gradient
%%%%%%%%
maxi(i,k)=maximum;
dcarri=dcarr(1:ypos,k);
pyi=py(1:ypos);
FWHM1 = interp1(dcarri,pyi,1/2*maximum);
dcarrb = dcarr(:,k);
pyb = py;
dcarrb(1:ypos)=[];
pyb(1:ypos)=[];
FWHM2 = interp1(dcarrb,pyb,1/2*maximum);
92                                             APPENDIX A. 2D SCRIPT FILE

           FWHM = FWHM2-FWHM1;
           FWHMvec(i,k)=FWHM;
           % barr saves information for width,
           % usf, maximum c gradient position (Ypos)and FWHM at the end of the
           % channel
           barrw(i,w)=w*4;
           barrUsf(i,w)=Usf;
           barrFWHM(i,w)=FWHM;
           barrYpos(i,w)=yposvec(i,k);
           cont(i,w)=postinterp(fem,’c’,[-1e-5;0]);
     end
      w
end
%%%%%%%%%%%%%%%%%%%%%%%%%
% Plotting functions    %
%%%%%%%%%%%%%%%%%%%%%%%%%
%figure; plot(py,carr(:,k));
%figure; plot(py,dcarr(:,k));
%figure; plot(usfvec,maxi(:,k));
%figure; plot(usfvec,FWHMvec(:,k));
figure(2) hold off
%Plotting as FWHM function of Usf and different widths
for i=ns:1:nr
    plot(barrUsf(:,i),barrFWHM(:,i),plotcolor(i));
    hold on
    le = [’M’,int2str(i),’=’,’’’’,’Width ’,int2str(i*t*1e6),’’’’];
    eval(le);
end legend(M2,M3,M4,M5,M6); xlabel(’USF’,’FontSize’,15);
ylabel(’FWHM’,’FontSize’,15); title(’FWHMvsUSF’,’FontSize’,15);
set(gca,’FontSize’,15) hold off figure(3)
%Plotting Ypos function as of Usf and different widths
hold on for i=ns:nr
    plot(barrUsf(:,i),barrYpos(:,i),plotcolor(i));
    le = [’M’,int2str(i),’=’,’’’’,’Width’,int2str(i*t*1e6),’’’’];
    eval(le);
end legend(M2,M3,M4,M5,M6); title(’PosofMax’,’FontSize’,15);
xlabel(’USF’,’FontSize’,15); ylabel(’Posc’,’FontSize’,15);
set(gca,’FontSize’,15) shg hold off figure(4)
%plotting of contration in [-1e-5,0] as funktion of wide and usf
Usfliste=linspace(tUSf,tUSf*nrUSf,nrUSf); hold on for i=ns:nr
    plot(Usfliste,cont(:,i),plotcolor(i));
    le = [’M’,int2str(i),’=’,’’’’,’Width’,int2str(i*t*1e6),’’’’];
    eval(le);
end
                                                                      93

legend(M3)%,M3,M4,M5)%M8,M9,M10);
title(’Contration’,’FontSize’,15); xlabel(’USF’,’FontSize’,15);
ylabel(’contration’,’FontSize’,15); set(gca,’FontSize’,15) shg hold
off

figure(7) hold off
%Plotting as FWHM function of Usf*w and different widths
for i=ns:1:nr
    plot(t*i*barrUsf(:,i),barrFWHM(:,i),plotcolor(i));
    le = [’M’,int2str(i),’=’,’’’’,’Width’,int2str(i*t*1e6),’’’’];
    eval(le);
    hold on
    end
legend(M1,M2,M4,M5,M6); xlabel(’X’,’FontSize’,15);
ylabel(’Y’,’FontSize’,15); set(gca,’FontSize’,15)
title(’FWHMvsWusf’,’FontSize’,15); hold off
%%%%%%%%%%%%%%%%%%%
% Saving function %
%%%%%%%%%%%%%%%%%%%
csvwrite(’barrw’,barrw); csvwrite(’barrFWHM’,barrFWHM);
csvwrite(’barrYpos’,barrYpos); csvwrite(’barrUsf’,barrUsf);
csvwrite(’cont’,cont);
94   APPENDIX A. 2D SCRIPT FILE
Appendix B

Q3D Script File

  %%%%%%%%%%%%%%%%%%%%%%
%     Model Q3D      %
% %%%%%%%%%%%%%%%%%%%%
%
%Simulation information:
%
%Mesh B, 5906 elements for a inlet channel width of 20E-6 m
%Peclet number - PE = 3000, damping force - alpha = 13.3E6 (Pa s)/m^2
%When simulation for different widths velocities are changed according
%to the width of the channel and the scaling factor
%
% COMSOL Multiphysics Model M-file
% Generated by COMSOL 3.3 (COMSOL 3.3.0.405, $Date: 2007/06/21 21:12:36 $)

%
% Clearing all parameters
clear all; clc; clf; flclear fem

% COMSOL version
clear vrsn vrsn.name = ’COMSOL 3.3’; vrsn.ext = ’’; vrsn.major = 0;
vrsn.build = 405; vrsn.rcs = ’$Name: $’; vrsn.date = ’$Date:
2007/06/12 04:00:20 $’; fem.version = vrsn;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%     Simulation parameters %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Number of widths used in the simulation
nr=6;
%Starting width
ns=2;
%Distance between widths

                                    95
96                                          APPENDIX B. Q3D SCRIPT FILE

t=10e-6;
%Number of Scaling factors
nrUSf=100;
%Distance between the scaling factors
tUSf=0.005;
% Pe number
Pe=3000;
%Diffusivity for the cells
D = 1e-13;
% Damping constant
alpha=13.3E6;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The for loop runs 1 time at each width
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for w = ns:1:nr
    % The Width of the side channel is defined
    width=t*w;
    % Geometry. The base point of the inle channel is defined so the oulet
    % channel will have a length of 100e-6 m regardless of the inlet
    % channel width.
    g1=rect2(’180e-6’,’40e-6’,’base’,’corner’,’pos’,{’-110e-6’,’0’},’rot’,’0’);
    g2=rect2(width,’40e-6’,’base’,’corner’,’pos’,{’-10e-6’,’-40e-6’},’rot’,’0’);

     % Analyzed geometry
     clear s s.objs={g1,g2}; s.name={’R1’,’R2’}; s.tags={’g1’,’g2’};

     fem.draw=struct(’s’,s); fem.geom=geomcsg(fem);

   % Iinitialize mesh B
fem.mesh=meshinit(fem, ... ’hmax’,[2e-6], ... ’hgrad’,1.1, ...
’hmaxvtx’,[4,3e-6,6,0.5e-6], ... ’hmaxedg’,[2,1e-6,6,1e-6]);




     % Descriptions
     clear descr descr.const= {’D’,’Diffusion constant’,’Us’,’Side
     flow’,’Ui’,’Inflow’}; fem.descr = descr;

     % (Default values are not included)

     % Application mode 1 (Navier-Stokes model)
     clear appl appl.mode.class = ’FlNavierStokes’; appl.gporder =
     {4,2}; appl.cporder = {2,1}; appl.assignsuffix = ’_ns’; clear
                                                                   97

prop prop.analysis=’static’; appl.prop = prop; clear bnd bnd.u0
= {0,0,0,’-4*Ui*s*(1-s)’,0}; bnd.type =
{’noslip’,’neutral’,’strout’,’uv’,’uv’}; bnd.v0 =
{0,0,0,0,’4*Us*s*(1-s)’}; bnd.ind = [3,1,1,1,5,2,1,1,4];
appl.bnd = bnd; clear equ equ.F_x = ’-alpha*u’; equ.rho = 1e3;
equ.eta = 1e-3; equ.cporder = {{1;1;2}}; equ.F_y = ’-alpha*v’;
equ.gporder = {{1;1;2}}; equ.ind = [1,1]; appl.equ = equ;
fem.appl{1} = appl;

% Application mode 2 (Convection Diffusion model)
clear appl appl.mode.class = ’FlConvDiff’; appl.assignsuffix =
’_cd’; clear prop prop.analysis=’static’; clear weakconstr
weakconstr.value = ’off’; weakconstr.dim = {’lm3’};
prop.weakconstr = weakconstr; appl.prop = prop; clear bnd
bnd.type = {’N0’,’cont’,’Nc’,’C’,’C’}; bnd.c0 = {0,0,0,0,1};
bnd.ind = [3,1,1,1,5,2,1,1,4]; appl.bnd = bnd; clear equ equ.D =
’D’; equ.v = ’v’; equ.u = ’u’; equ.ind = [1,1]; appl.equ = equ;
fem.appl{2} = appl; fem.frame = {’ref’}; fem.border = 1; clear
units; units.basesystem = ’SI’; fem.units = units;

% Multiphysics
fem=multiphysics(fem);

% Extend mesh
fem.xmesh=meshextend(fem);

%Defining empty vectors
maxi = zeros(1,2); usfvec = zeros(1,2); xposvec = zeros(1,2);
FWHMvec = zeros(1,2);

% This loops runs one time for each scaling factor
for i=1:1:nrUSf Usf = tUSf*i;
    %Claculate Velocity in main channel
    Ui=(Pe*D*(40e-6+width))/(40e-6*(40e-6+Usf*width));
    usfvec(1,i)=Usf;

    % Definition of constants
    fem.const = {’D’,D, ... ’Ui’,Ui, ... ’Us’,Usf*Ui,...
    ’alpha’,alpha};



    disp(sprintf(’Ui: %11.3e   Pe: %9.3f Usf: %9.3f’,Ui,Pe,Usf))
98                                      APPENDIX B. Q3D SCRIPT FILE

     % Solve problem (Navier-Stokes)
     fem.sol=femstatic(fem, ... ’solcomp’,{’u’,’p’,’v’}, ...
     ’outcomp’,{’u’,’p’,’v’});

     % Solve problem (Convection diffusion)
     fem.sol=femstatic(fem, ... ’init’,fem.sol, ...
     ’solcomp’,{’c’}, ... ’outcomp’,{’c’,’u’,’p’,’v’}); figure(1)
     subplot(2,1,1)
     % Plot solution
     postplot(fem, ...
     ’tridata’,{’U_ns’,’cont’,’internal’,’unit’,’m/s’}, ...
     ’trimap’,’jet(1024)’, ... ’title’,’Surface: Velocity field
     [m/s]’); axis equal tight
     % Plot solution
     subplot(2,1,2) postplot(fem, ...
     ’tridata’,{’c’,’cont’,’internal’,’unit’,’mol/m^3’}, ...
     ’trimap’,’jet(1024)’, ... ’title’,’Surface: Concentration, c
     [mol/m^3]’); axis equal tight;

     % k set to 1 for messuring the FVHM at the end   of the channel
     k=1;
     % tvarsnit set to -100e-6
     tvarsnit(k)=-100e-6*k;

     % np is the number of division of the y coordinate in the cross
     % sections
     np = 100;
     px =-100e-6*k*ones(1,np); % prefactor changed to -100e-6
     py = linspace(0,40e-6,np);
     %%%%%%%%%
     % Finding maximum contration gradient at main channel end
     %%%%%%%%%
     carr(:,k)=postinterp(fem,’c’,[px;py]); % contration matrix
     dcarr(:,k) = postinterp(fem,’grad_c_cd’,[px;py]); % gradient of contration
     [maximum,ypos] = max(dcarr(:,k)); yposvec(i,k)=ypos*40e-8;
     %%%%%%%%
     % Finding FWHM of contration gradient
     %%%%%%%%
     maxi(i,k)=maximum; dcarri=dcarr(1:ypos,k); pyi=py(1:ypos);
     FWHM1 = interp1(dcarri,pyi,1/2*maximum); dcarrb =
     dcarr(:,k); pyb = py; dcarrb(1:ypos)=[]; pyb(1:ypos)=[];
     FWHM2 = interp1(dcarrb,pyb,1/2*maximum); FWHM = FWHM2-FWHM1;
     FWHMvec(i,k)=FWHM;
     %barr saves information for width,
                                                                        99

        % usf, maximum c gradient position (Ypos)and FWHM at the end of the channel
        barrw(i,w)=w*4; barrUsf(i,w)=Usf; barrFWHM(i,w)=FWHM;
        barrYpos(i,w)=yposvec(i,k);
        cont(i,w)=postinterp(fem,’c’,[-1e-5;0]); end w end
%%%%%%%%%%%%%%%%%%%%%%%%%
% Plotting functions    %
%%%%%%%%%%%%%%%%%%%%%%%%%
%figure; plot(py,carr(:,k));
%figure; plot(py,dcarr(:,k));
%figure; plot(usfvec,maxi(:,k));
%figure; plot(usfvec,FWHMvec(:,k));
figure(2) hold off
%Plotting as FWHM function of Usf and different widths
for i=ns:1:nr plot(barrUsf(:,i),barrFWHM(:,i),plotcolor(i)); hold on
le = [’M’,int2str(i),’=’,’’’’,’Width ’,int2str(i*t*1e6),’’’’];
eval(le); end legend(M2,M3,M4,M5,M6); xlabel(’USF’,’FontSize’,15);
ylabel(’FWHM’,’FontSize’,15); title(’FWHMvsUSF’,’FontSize’,15);
set(gca,’FontSize’,15) hold off figure(3)
%Plotting as Ypos function of Usf and different widths
hold on for i=ns:nr plot(barrUsf(:,i),barrYpos(:,i),plotcolor(i));
le = [’M’,int2str(i),’=’,’’’’,’Width’,int2str(i*t*1e6),’’’’];
eval(le); end legend(M2,M3,M4,M5,M6);
title(’PosofMax’,’FontSize’,15); xlabel(’USF’,’FontSize’,15);
ylabel(’Posc’,’FontSize’,15); set(gca,’FontSize’,15) shg hold off
figure(4)
%plotting of contration in [-1e-5,0] as funktion of wide and usf
Usfliste=linspace(tUSf,tUSf*nrUSf,nrUSf); hold on for i=ns:nr
plot(Usfliste,cont(:,i),plotcolor(i)); le =
[’M’,int2str(i),’=’,’’’’,’Width’,int2str(i*t*1e6),’’’’]; eval(le);
end legend(M2,M3,M4,M5,M6); title(’Contration’,’FontSize’,15);
xlabel(’USF’,’FontSize’,15); ylabel(’contration’,’FontSize’,15);
set(gca,’FontSize’,15) shg hold off

figure(7) hold off
%Plotting as FWHM function of Usf*w and different widths
for i=ns:1:nr plot(t*i*barrUsf(:,i),barrFWHM(:,i),plotcolor(i)); le
= [’M’,int2str(i),’=’,’’’’,’Width’,int2str(i*t*1e6),’’’’]; eval(le);
hold on end legend(M2,M3,M4,M5,M6); xlabel(’X’,’FontSize’,15);
ylabel(’Y’,’FontSize’,15); set(gca,’FontSize’,15)
title(’FWHMvsWusf’,’FontSize’,15); hold off

%%%%%%%%%%%%%%%%%%%
% Saving function %
%%%%%%%%%%%%%%%%%%%
100                                        APPENDIX B. Q3D SCRIPT FILE

csvwrite(’barrw’,barrw); csvwrite(’barrFWHM’,barrFWHM);
csvwrite(’barrYpos’,barrYpos); csvwrite(’barrUsf’,barrUsf);
csvwrite(’cont’,cont);
Appendix C

Separation Mechanism Simulation

In this appendix a short explanation of the separation mechanism simulations used in
Chapter 9 is given.


C.1     Simulation Definition
The goal for this simulation was to determine the separation position for three different
geometries of the outlet channels. The tree different geometries were:

  A Identical outlet channels

  B Outlet channel 1 and 2 have a hydraulic resistance 7 times larger than that of channel
    3 to 6

  C Outlet channel 1 has a hydraulic resistance 7 times larger than that of channel 2 to
    6

A Comsol model of geometry 1 can be seen in Fig. C.1.
    The simulation was made with the quasi 3D method, explained in Chapter 7 and the
increment in hydraulic resistance was achieved by making the appropriate channels 7 times
longer.
    The viscosity was set to η = 10−3 Pa s, and the density was set to ρ = 1000kg/m3 . The
pressure difference between the pinched segment and all the outlet channels were identical.
    The geometries used is stated below:

   • Pinched segment: 40 × 100 µm

   • Separation segment radius: 750µm

   • Normal outlet channel: 388µm × 1000µm

   • Channel height: 30µm

   Mesh for simulation was set as stated below:

                                           101
102                         APPENDIX C. SEPARATION MECHANISM SIMULATION




[m]




                                              [m]

Figure C.1: Here is shown the overall geometry in Comsol 3.3 for particle separation simulations.
In this case all outlet channel are identical.



      • Global maximum element size: 100µm

      • Maximum element size in pinched segment: 1µm

      • Maximum element size at corners between channels: 6µm

The total amount of elements were 24237 for geometry 3.


C.2       Simulation
The system was solved for the incompressible Navier–Stokes equation. To determine the
separation position between channel 1 and 2, a streamline was placed at the border point
between the channels. The line was then backtracked and it was possible to readout
the separation position in the pinched segment from the streamlines position. To see
streamlines for geometry 2 see Fig. C.2.
    Other separation position were determined by moving the starting point of streamlines
to other border points.
C.2. SIMULATION                                                                          103




[m]




                                            [m]
Figure C.2:      This is the Comsol simulation for particle separation. The red lines are the
streamlines corresponding to separation position between channel 1 and 2 and channel 2 and 3.
104   APPENDIX C. SEPARATION MECHANISM SIMULATION
Appendix D

Maple Worksheet for Calculating
the Hydraulic Resistance




               105
106APPENDIX D. MAPLE WORKSHEET FOR CALCULATING THE HYDRAULIC RESISTANCE
107
108APPENDIX D. MAPLE WORKSHEET FOR CALCULATING THE HYDRAULIC RESISTANCE

				
DOCUMENT INFO
Shared By:
Categories:
Tags:
Stats:
views:26
posted:8/11/2011
language:English
pages:126