VIEWS: 26 PAGES: 126 POSTED ON: 8/11/2011 Public Domain
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 eﬀort 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 microﬂuidic 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 mikroﬂuid 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 aﬂeveret, 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 ﬁgures 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 Diﬀusion 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 Buﬀer Fluid . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 4 Fluid Motion 11 4.1 Flow Proﬁle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 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 Diﬀusion 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 Diﬀusion 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 Deﬁnition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 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 Deﬁnition 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 Deﬁnition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 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 proﬁle of the SU8 stamp . . . . . . . . . . . . . . . . . . . . . . . . 3 1.3 The pinched segment of the ﬁnal PFF system made in PDMS . . . . . . . . 3 3.1 Red blood cell geometry . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 4.1 Geometry used to determine the parallel plate ﬂow proﬁle . . . . . . . . . . 11 4.2 Geometry used to determine the ﬂow proﬁle for an rectangular cross section 13 4.3 Velocity proﬁle in a channel with a rectangular cross section . . . . . . . . . 14 5.1 Fixed sphere in a steady state ﬂow . . . . . . . . . . . . . . . . . . . . . . . 18 5.2 Moving sphere in a steady ﬂow . . . . . . . . . . . . . . . . . . . . . . . . . 20 6.1 Geometry used in analytical calculations . . . . . . . . . . . . . . . . . . . . 25 6.2 Borderline position as function of the ﬂow rate ratio . . . . . . . . . . . . . 28 6.3 Model for the description of diﬀusion . . . . . . . . . . . . . . . . . . . . . . 29 6.4 Cross section of the analyticaly determied concentration proﬁle . . . . . . . 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 Proﬁle 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 Simpliﬁed 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 ﬂow rate ratio . . . . . . . . . . . . . . . . . . . 47 8.4 FWHM as a function of the ﬂow rate ratio with a course mesh . . . . . . . 48 8.5 FWHM as a function of the ﬂow rate ratio with a ﬁn mesh . . . . . . . . . 48 8.6 Concentration proﬁle with a ﬂow rate ratio of 0.03 . . . . . . . . . . . . . . 49 8.7 Concentration proﬁle with a ﬂow 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 ﬂow 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 butterﬂy diﬀusion proﬁle in 3D . . . . . . . . . . . . . . . 71 11.1 Basic overview of a pinched ﬂow device . . . . . . . . . . . . . . . . . . . . 74 11.2 Alignment of two particles . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 11.3 Optimization plot of inlet channel width versus ﬂow 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 Inﬂuence on the separation position due to pressure drops . . . . . . . . . . 59 9.2 Change in separation position do to pressure drop in diﬀerent 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 ﬂuids Pa s Q Volume ﬂow 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 Diﬀusion 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 diﬀerent 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 microﬂu- 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 ﬂuidic 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 ﬂuid, are introduced via the inlet channel. The main channel must have a larger ﬂow 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 simpliﬁed 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 ﬂow 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 ﬁnally 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 ﬁgure 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 eﬀect. 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 ﬁnal 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 proﬁle of the SU8 Figure 1.3: This is a part the ﬁnal 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 buﬀer. • Chapter 4 : ”Fluid Motion”. A derivation of various velocity proﬁles that are used throughout this thesis. • Chapter 5 : ”Particle Motion”. A consideration of forces that inﬂuence on the particle motion in microﬂuidic 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 ﬁnal design”. In this chapter various parameters inﬂuencing on the ﬁnal 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 ﬂow, 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 microﬂuidic are presented. These equations are necessary for the understanding the ﬂuid motion. When the expression ”ﬂuid 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 ﬂuid, 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 ﬂuid particle. In the following section four fundamental equations will be described. These are the equation of continuity, the Navier–Stokes equation, the diﬀusion equation and ﬁnally an equation relating pressure drop to the total ﬂowrate, 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 ﬂuid leaving this volume through a small part of the surface ds in an unit time, is ρvds where v is the ﬂuid velocity. The total amount of ﬂuid 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 ﬂuid 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 ﬂuid volume is equal to the amount of ﬂuid that leaves through its surface. 2.2 Navier–Stokes Equation The equation of motion in microﬂuidic systems is called the Navier–Stokes equation. It is derived from Newton’s second law for a ﬂuid 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 deﬁned positions. In the transformation from a mowing lagrangian coordinate system, in which the ﬂuid particle is ﬁxed, to a ﬁxed coordinate system in which the ﬂuid 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 ﬂuid particle, the Navier–Stokes equation for an incompressible ﬂuid 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 ﬂuid particle due to a diﬀerence in pressure around it. η is the ﬂuid viscosity. The second term of the body forces is therefore the viscous force, which arrises due to frictional forces between a ﬂuid particle and its neighboring ﬂuid particles. To arrive at this result it is further assumed, that the ﬂuid is incompressible which simpliﬁes 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 ﬁelds for a ﬂuid system. The equation is a nonlinear partial diﬀerential equation which for most cases can not be solved analytically. But for some ideal cases, it simpliﬁes signiﬁcantly and makes it possible to determine an analytical solution. The condition for this simpliﬁcation will be discussed in Chapter 4. 2.3 Diﬀusion Equation In the study of particle motion in ﬂuids it is necessary to consider the diﬀusion process. All particles will experience this eﬀect to a certain amount depending on their diﬀusivity in the particular ﬂuid. From the continuity equation it is possible to determine an equation which describes the diﬀusion 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 diﬀusivity, which is inde- pendent of the concentration at low concentrations. This equation is used to predict the concentration proﬁle. A general measure for the system diﬀusion 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 diﬀusion will be small relative to the convection, which is the ﬂuid velocity driven motion. The number is given as the characteristic diﬀusion 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 diﬀusion dominates. w2 Diﬀusion 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 ﬂow 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 ﬂow rate will be used throughout this thesis to determine both ﬂow rates and hydraulic resistances. The hydraulic resistance depend on a number of diﬀerent parameters for instance the channel length, width, hight and ﬂuid 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 proﬁle and hence the ﬂow rate, some boundary condition are needed. In ﬂuidic systems this boundary condition is the ”no-slip” condition on the side walls of the channel. This condition implies that the ﬂuid must be at rest relative to the channel sidewalls. If the ﬂuid were allowed to slide along the sidewalls it would continuously loose kinetic energy due to friction and eventually terminate the ﬂow. 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 eﬀect 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 buﬀer 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 ﬁgure 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 suﬃcient to change their shapes. 3.3 Buﬀer Fluid The red and withe blood cells which will be going to be separated is suspended in a standard buﬀer called Phosphate Buﬀered Saline (PBS). This buﬀer is unlike normal blood newtonian, and therefore it will have a viscosity which is independent of for instance the local shear of the ﬂuid. The fact that this buﬀer is newtonian simpliﬁes the calculation drastically. This buﬀer consist mainly of saltwater, and it will therefore be assumed that the density and viscosity of the buﬀer 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 buﬀer the density and viscosity is denoted as not available. Chapter 4 Fluid Motion In determination of ﬂuid 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 ﬂow proﬁle for two type of channels will be derived, and the Reynolds number will be introduced. 4.1 Flow Proﬁle To determine the actual ﬂow proﬁles 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 simpliﬁes 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 inﬁnitely thick chan- nels. Therefore it is necessary to derive the ﬂow proﬁle 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 simpliﬁes 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 ﬂow proﬁle. Due to the pressure diﬀerence the ﬂuid 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 proﬁle stated in Eq. (4.4). ∆p vx (z) = (w − z)z (4.4) 2ηL It is seen that the velocity proﬁle 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 ﬁnal 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 ﬂow proﬁle for channels having this particular cross section. In the ﬁnal 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 ﬂow proﬁle of the ﬂuid in these channels. Steady state will be assumed in ﬁnding the solution to the Navier–Stokes equation. ∇p = η∇2 v (4.5) To determine the ﬂow proﬁle one has to solve the partial diﬀerential 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 ﬂow only has a non zero component in the x direction Eq. (4.6a) simpliﬁes 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 proﬁle. 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 proﬁle 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 coeﬃcient fn (y) for odd n, one must solve the following second order diﬀerential 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 proﬁle 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 proﬁle 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 proﬁle 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 ﬂow is a pressure driven steady state ﬂow in channels, it is also known as a poiseuille ﬂow. This is the type of ﬂow that will be in the ﬁnal 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 microﬂuidic systems this is often well below unity. When the nonlinear term vanishes the ﬂow is also said to be a laminar ﬂow. A laminar ﬂow is characterized by the absence of ﬂow 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 ﬂow 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 ﬂowrate. The ﬂowrate 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 ﬁnal design the height of the channel will be constant. When channels with a ﬁnite height are used instead of inﬁnitely high channels the ﬂuid 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 eﬀect can also be seen in Fig. 4.3. In this ﬁgure the velocity proﬁle is shown and it is obvious that it has a parabolic appearance in the height direction as well. Therefore the ﬁnite channel height acts as a damping force an the total ﬂow rate. 4.5 Streamlines A small volume of ﬂuid that ﬂows in a system from a certain point will follow a speciﬁc path through the system. This ﬂuid path is called a streamline. In order for two streamline to cross each other a speciﬁc point in the geometry must have two diﬀerent velocity gradients. This is of course not physically possible. Therefore streamlines will newer cross. When the ﬂow 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 ﬂuid will be investigated. In order to understand the motion of particles in a ﬂuidic 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 ﬂuid, 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 ﬂuid it will tend to follow the ﬂuid motion. This is due to the Stokes drag which is present when the particle velocity is diﬀerent from that of the ﬂuid. In the derivation of the forces acting on the particles, it is used that the ﬂuid motion around a sphere is the same whether or not the sphere is mowing or ﬁxed at a certain position. Since the forces are identical, a ﬁxed 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 ﬁeld that have a uniform velocity at inﬁnity, v∞ . To derive the Stokes drag the two fundamental microﬂuidic 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 deﬁning v as the velocity around the sphere, and v∞ as the ﬂuid velocity at inﬁnity 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 ﬁxed sphere in a steady state ﬂow. The ﬂuid will have the velocity v∞ at inﬁnity and v around the sphere. From vectorial mathematics it is known that the divergence of a rotation ﬁeld 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 ﬁeld: 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 simpliﬁed like Eq. (5.1). To determine the Stokes drag force on the sphere due to the ﬂuid motion, one has to determine the rotation of the vector ﬁeld 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 ﬁnally 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 ﬁeld around the sphere is known it is possible to determine the pressure distribution on the surface of the sphere. Starting from the simpliﬁed 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 ﬁxed 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 diﬀerence. Therefore will this force make the particles follow the streamline motion. v∞ u v∞ Figure 5.2: A mowing particle in a constant ﬂow. It will experience a drag force proportional to the velocity diﬀerence between the sphere and the ﬂuid. 5.2 Other Forces Acting on the Particle Even when a particle is moving in an unidirectional ﬂow it will experience a lateral force perpendicular to the streamlines of the ﬂow. This is due to the fact that the velocity proﬁle in a channel of a ﬁnite 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 diﬀerence of stress and pressure on the two sides of a particle suspended in this ﬂow. Therefore the particle will experience a force acting perpendicular to the ﬂuid streamlines. In the previous section we considered the case where the particle was placed in an inﬁnitely large channel with an uniform velocity proﬁle. 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 ﬁeld will not be uniform, and therefore this physical eﬀect needs to be considered. Since it is vital for the functionality of the pinched ﬂow that the particles are aligned along one sidewall, the lateral migration could have a big inﬂuence on the ﬁnal device. The lateral migration eﬀect 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 proﬁle, 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 ﬁeld v = U − V and q = P − Q respectively. Here U and Q are the ﬁelds that the ﬂuid would have if the particle wasn’t present. Where as V and P are the ﬂuid ﬁelds 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 ﬁrst 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 inﬁnity in the channel direction the particles should have the same velocity as the ﬂuid. Using these equations together with the assumptions of a Poiseuille ﬂow proﬁle, [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 ﬂow. 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 ﬂuid velocity, κ is deﬁned 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 deﬁne 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 ﬂuid velocity vm = 10 µm/s. When a Poiseuille ﬂow 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 ﬂuid 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 eﬀect 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 ﬂow. 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 diﬀerence between the particle - and the ﬂuid velocity. When the same mean 2 ﬂuid 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 simpliﬁcation to neglect it. The blood cells will of course also be inﬂuenced to some extend by the gravitational force. But as described in Chapter 3 will the blood cells have almost the same density as the buﬀer in which they are suspended. Therefore the gravitational force acting on the particles will be almost balanced by the lifting force from the ﬂuid. And even though it had an eﬀect it will make no diﬀerence 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 ﬂuid 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 ﬂuid velocity at inﬁnity 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 simpliﬁes 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) simpliﬁes 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 buﬀer, 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 Diﬀusion Constants The particles will to some extend be inﬂuenced by the diﬀusion process. A critical parameter when considering diﬀusion of particles suspended in a liquid is the diﬀusion 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 diﬀusion 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 diﬀusion constant of 10−13 will be used for both red and white blood cell. This is done in order to be sure that the eﬀect from diﬀusion is not underestimated. If the actual diﬀusion 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 ﬁgure represents the boundary between the two ﬂuids. 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 ﬂuid 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 ﬂuid 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 ﬂuid velocity. The velocity in the x direction is described by the following diﬀerential 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 deﬁned 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 ﬂuid 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 ﬂow rate in the main channel. This causes the inlet ﬂow to be squeezed against the sidewall In this section the particles diﬀusion process will be neglected. It will be considered in the following section. When diﬀusion is neglected the inlet channel ﬂow and the main channel ﬂow will remain separated. Therefore it is possible to determine the asymptotic borderline position shown in Fig. 6.1 as a function of the relative ﬂow rate diﬀerence between the main channel and the inlet channel. To derive the expected asymptotic borderline position (w′ ) the fact that the ﬂuids will Qinlet remain separated is used. When this is the case the ﬂow rate ratio , which is the Qmain inlet channel ﬂow rate divided by the main channel ﬂow 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 diﬀusion 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 Diﬀusion Process In the previous section the diﬀusion process was neglected. This is of course not a valid simpliﬁcation, so in order to determine the inﬂuence due to the diﬀusion process a simpli- ﬁed model is used. In this chapter the treatment of particles will be done by a concentration proﬁle. In such a way that this proﬁle can be seen as a probability distribution for the particles. Through this concentration treatment we will get an understanding of how diﬀerent pa- rameters inﬂuence on the particle position. To get an idea of whether diﬀusion or convection will dominate, let us determine the e P´clet number for the simpliﬁed 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 ﬂows. It is seen that it as expected grows for increasing ﬂow 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 ﬂuids will not manage to mix completely before the end of the channel. In order to derive an expression for the concentration proﬁle, 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 proﬁle 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 ﬂow in the main channel. In the cross sections of this coordinate system the diﬀusion process will look like the diﬀusion as constant total concentration process instead of the constant concentration process. In this moving coordinate system the delta function describing the concentration proﬁle 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 diﬀusion process. The inlet channel is treated as a point source with a constant concentration described by a delta function. The ﬂuid in the main channel is assumed to have a constant velocity of vmean . The diﬀusion 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 ﬂuid, 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 simpliﬁcations 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 proﬁle 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 proﬁle has the same structure as for a normal diﬀusion x process, except that the time parameter in this expression is replaced by . A plot vmean of this concentration proﬁle is shown in Fig. 6.4 which is a cross section at x = 80 µm. In this ﬁgure it is seen that the eﬀect due to the diﬀusion process is that the boundary between the two ﬂows 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 diﬀerentiate 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 proﬁle at x = 80 µm. It is seen that at this x position the maximum gradient position is a little below 2 µm. This ﬁgure is made with a diﬀusion 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 ﬁnd 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 diﬀusion 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 diﬀusion 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 eﬀect. 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 ﬂow rate ratio. This independence could be due to the fact that the velocity proﬁle in the treatment of the diﬀusion 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 ﬂow 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 ﬂow rate and the main channel ﬂow rate re- spectively. Eq. (6.23) is only valid when the velocity proﬁle of the ﬂuid is assumed to be constant, which is not a valid simpliﬁcation. 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 ﬂow rate divided by the total ﬂow rate should be less than the ﬂow rate between 34 CHAPTER 6. ANALYTICAL CALCULATIONS the wall and the intersection width divided by the total ﬂow 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 ﬂow rate ratios demand for aligning of a particle with a radius of 5 µm in a 40 µm wide pinched segment. The original ﬂow rate and the new ﬂow 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 diﬀerence 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 inﬂuences 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 diﬀusion-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 diﬀer- ential equation like the Navier–Stokes equation can solved numerically in many diﬀerent 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 diﬀerent softwares. In this projet all solutions are determined with the use of the simulation program: Comsol 3.3, which is specialized in solving partial diﬀerential 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 ﬂuid dynamics one often wants to determine the velocity ﬁeld of the ﬂow. This ﬁeld can be represented by the vector ﬁeld v(r, t) which is the solution to a general partial diﬀerential equation. Lv(r, t) = f (r, t), for r ∈ Ω (7.1) In this equation L is a diﬀerential 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 satisﬁes the boundary condition. An example of a division of a 2D Gaussian-shaped channel into ﬁnite 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 ﬁnite sized elements. basis function meets the boundary condition of Within this mesh one of the linearly interpolat- unity in the raised corner. The ﬁgure comes ing basis functions φj is shown. This ﬁgure 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 ﬁnite element size. This can result in nonphysical solutions and must therefore be considered. This section will comment on the diﬀerent results that can arise do to a chance in element size. In our simulation, the mesh structure has been deﬁned so that the important parts of the PFF-system have a ﬁner 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 deﬁning the mesh in this project are listed in Table 7.1. When a model is deﬁned in Comsol, each corner, boundary and subdomain are indexed. Therefore it is possible to assign the parameters described in Table 7.1 to speciﬁc parts of the model. In the ﬂuid 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 speciﬁc point hgradedg: Element growth rate from a speciﬁc 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 deﬁned, so that certain important parts of the structure will have a ﬁner 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 ﬁner elements. ﬂuids from the inlet and the main channel the velocity and concentration proﬁle can be rather complex. Therefore the mesh structure is made ﬁner at point 4 and 6. Since the interesting parameter in the simulation is how the ﬂuids moves down the outlet channel the element size along boundary 1 is also decreased. This ﬁner mesh is indeed necessary due to the small diﬀusion constant which results in a sharp concentration drop. A sharp drop needs ﬁne mesh in order to be described right. The mesh structures cannot be inﬁnitely small, since this would result in inﬁnitely many elements, which of course is not solvable. When a ﬁnite 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 diﬀerent mesh structures was investigated, their details are shown in Table 7.2. In order to test the eﬀect of using the tree diﬀerent 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 ﬁnd 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 ﬁgures it is evident that the ﬁner mesh structures that were used, the smoother curve one gets. Another eﬀect 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 ﬁner mesh, the better result, was seen. But since the time it takes to perform a simulation increases with ﬁner mesh a good simulation will be a simulation that fulﬁll the two following conditions: 1. The meshstructure should be ﬁne enough to yield a proper result. 2. The mesh structures must not be so ﬁne 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 insuﬃcient at the critical positions. It was seen that mesh A was not ﬁne 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 ﬁnite height instead of inﬁnitely thick channels are used the velocity proﬁle will be limited. As described in Chapter 4, the ﬁnite height of the channels acts as an artiﬁcial 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 eﬀect, 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 ﬁnite channel height, is included in the Navier–Stokes equation. This results in the modiﬁed 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 ﬁnite 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 ﬂuid is assumed to ﬂow along the x direction. At ﬁrst the situation is considered along the y axis. This result in a parallel plate ﬂow. In Chapter 4 the velocity proﬁle 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 artiﬁcial damping term is applied to determine the velocity proﬁle. This equation simpliﬁes even further than the previous equation because the z dependence is neglected. Eq. (7.3) then simpliﬁes into: α(h)v = ∇p (7.5) This equation is easily solved and yields the following velocity proﬁle. ∆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 inﬁnitely large parallel plates a ﬂuid is assumed to ﬂow in the x direction. 7.4. QUASI 3D SIMULATION METHOD 41 An expression for the damping factor α can be found by demanding that the ﬂowrate out of the model shown in Fig. 7.8 should be identical regardless of which method that is used to determine the velocity proﬁle. Therefore let us determine the two ﬂowrates. 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 ﬂowrate 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 ﬂow 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 artiﬁcial 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 artiﬁcial 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 proﬁle and thereby make the ﬂuid less easy ﬂowing. Chapter 8 Fluid Interaction Simulation The simulations in this chapter will concentrate mainly on the interaction process of ﬂuids in the pinched segment. Due to the diﬀusion process the interaction process will cause the transition from one ﬂuid to the other to be less abrupt. Due to the small length scales the ﬂuids will not mix completely, but to some extend remain more or less separated. In the actual device it is the same ﬂuid, that ﬂows in both main- and inlet channel, the only exception is that the ﬂow from the inlet channel contains particles. In these simulations the cells will be represented by a concentration proﬁle, which as described in Chapter 6 can be interpreted as a probability distribution of ﬁnding the particle at the exact position. 8.1 Model Deﬁnition The simulations in this chapter are made in quasi 3D mode on a 2D model. The inﬂuence of a ﬁnite channel height is as described in Chapter 7 included as an artiﬁcial 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 ﬂow is introduced via the inlet channel of width d. A ﬂow without particles is introduced in the main channel of width w. The two ﬂows 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 eﬀect with a pinched segment of this length, see [3] and [9]. The simulation will not contain actual particles. Instead a concentration proﬁle 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 inﬂuence on the ﬁnal result. It is only the relationship between the e velocities in the two channels that inﬂuence 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 ﬂow rate direction in the channels, in this ﬁgure illustrated by the Q’s. The line denoted ”Borderline” is the line illustrating the boundary between the two ﬂuids. The velocity in the inlet channel is deﬁned 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 ﬂow. This is important since the purpose of the main channel ﬂow is to ensure the alignment of the particles. Alignment is achieved by a larger ﬂow ratio from the main channel, which thereby will force the inlet channel ﬂow 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 diﬀerent 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 ﬁelds in ﬁnite 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 diﬀerence between the two channel widths that inﬂuences on the result, it is suﬃcient 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 ﬂuid 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 deﬁne 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 proﬁle should have evened out. So that it looks like it consisted of only one ﬂuid. A good result is further characterized by a sharp concentration drop at the end of the outlet channel. The optimal concentration proﬁle at this position would be a step function. This proﬁle would appear if the ﬂuids did not mix. The steepness of the concentration drop can be determined by investigation of the concentration derivative proﬁle 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 ﬂuid 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 deﬁned 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 ﬂow ratio relationship between the main and the inlet channel. Therefore it is expected that the intersection width will increase when the ﬂow rate ratios increases. 46 CHAPTER 8. FLUID INTERACTION SIMULATION Parameter Symbol Explanation Intersection width w′ Width of the particle containing ﬂow ∂c FWHM of ∂y FWHM Width of the mixed layer between the ﬂows 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 ﬂow 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 diﬀusion time. We will expect that the corner concentration grows for increasing ﬂow 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 diﬀusion 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 ﬁrst 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 ﬂow 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 ﬂow 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 ﬂow 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 ﬁgure the FWHM is shown as a function of the ﬂow rate ratio. It is seen that for both small and large ﬂow rate ratios the FWHM decreases. From this plot it is also seen that the eﬀects are controlled mainly by the ﬂow rate ratios, since the graphs are simulated for diﬀerent width and heights of the cahnnels. The dashed lines are included to illustrate the eﬀect 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 ﬂow 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 ﬂow 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 ﬁgure 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 ﬂuids varied between 3 µm and 4.4 µm. The FWHM has a maximum of approximately 4.4 µm at a ﬂow rate ratio of around 0.05. In order to understand this maximum it must be realized that two diﬀerent limiting eﬀects inﬂuence on the FWHM behavior as illustrated in Fig. 8.3. For small ﬂow rate ratios (< 0.05) the FWHM drops as the Flow rate ratio decreases. When the ﬂow rate ratio is small so will the width of the concentration proﬁle 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 ﬂow rate ratio of around 0.05. This par- variation is simulation eﬀect due to the ﬁnite ticular behavior of the FWHM is due to two element size. They will therefore not be com- diﬀerent limiting eﬀect 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 ﬁeld. The velocity in the channels will increase towards the center due to the ”no-slip” condition. The actual diﬀusion process takes place at the boundary between the two ﬂows. The longer into the channel, this boundary is, the shorter time the ﬂuid at the boundary is in the channel. When bigger ﬂow rate ratios are used the boundary position tends towards the center of the channel and thereby limits the diﬀusion time. A shorter diﬀusion time will of course limit the FWHM. This maximum is therefore just the maximum between these two limiting eﬀects. Therefore a ﬂow 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 ﬁnite 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 ﬁne mesh is used respectively. Another eﬀect 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 ﬂow rate ratio of 0.01 and 0.02 respectively. On these ﬁgures it is seen that when small ﬂow rate ratios are used the mixing process of the two ﬂuids actually takes place inside the inlet channel, and the concentration plot show a kind of wave structure. This makes the velocity proﬁle more complex at the intersection between the two ﬂuids and will therefore increase the particle diﬀusion 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 ﬂow rate ratio for two diﬀerent 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 proﬁle Figure 8.7: In this ﬁgure a ﬂow rate ratio when a ﬂow 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 ﬂuids 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 ﬂuids still remain more or less separated at the end of the outlet channel. this ﬁgure it is seen that the corner concentration follows the expected results. It increases smooth and rapidly as the ﬂow rate ratio increases. It is further seen that around a ﬂow 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 ﬂow 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 ﬂow rate ratio. From this ﬁgure 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 diﬀusion constants of the particles are below 1 × 10−13 m2 /s. The step like behavior is due to the ﬁnite mesh size. This eﬀect 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 ﬂow 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 ﬂow 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 ﬂow rate ratios above 0.02. 8.5 Conclusion and Optimal Design Parameters It was found in this chapter that the behaviour of the ﬂuid was mainly governed by the ﬂow 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 ﬂow rate ratio for various channel widths and heights. It is seen that it approaches unity for ﬂow 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 ﬂow rate ratios. It is the ﬁnite mesh size that is reason for the step like behavior. 8.5. CONCLUSION AND OPTIMAL DESIGN PARAMETERS 51 as a function of the ﬂow 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 ﬂows. Therefore one should use ﬂow rate ratios above 0.02. But since the intersection width increases for increasing ﬂow rate ratios, one must make a compromise. In order to get the smallest FWHM it would be preferred to use large ﬂow 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 ﬁnal 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 ﬂow rate ratio of 0.1 must be used. 52 CHAPTER 8. FLUID INTERACTION SIMULATION Chapter 9 Separation Mechanism The most important part of the ﬁnal device is of course the actual separation part. The ﬂow in the pinched segment now contains both the inlet ﬂow, in which the particles are suspended, and the main ﬂow. This ﬂow 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 deﬁned 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 ﬂow 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 ﬁgure 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 ﬂuid that is predicted channel 1 and 2. The separation position is de- to end up in channel 1, be refereed to as the ﬁned as the distance from the upper side wall in channel 1 ﬂow. 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 ﬂuid due to various geometry changes. This will be the focus in this chapter. 9.1 Identical Outlet Channel The simplest ﬂuid motion would arise if the ”no slip” condition on the sidewalls could be neglected. If this was the case the ﬂuid 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 ﬂuid streamline. These parameters is illustrated in Fig. 9.2. In an actual device the negligence of the ”no slip” condition is not a valid simpliﬁcation. Therefore Eq. (9.1) will not hold true. It neglects the fact that the ﬂow 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 ﬂow rate each. It will be assumed that the ﬂuid will act as a creeping ﬂow. A creeping ﬂow 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 ﬂowing ﬂuid is necessary, and this is exactly what the limited height of the channel counteracts. With the assumption of creeping ﬂow are the separation positions between the six channels controlled by the ﬂow proﬁle in the pinched segment. In such a way that the upper sixth of the ﬂow rate will end up in channel 1. As described in Chapter 4 is the ﬂow proﬁle 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 proﬁle is known, it is possible to determine the ﬂow rate for a given channel ﬂow. The ﬂow rate for the channel 1 ﬂow is given by: −w/2+w′ h Q= dy dz vx (y, z) (9.3) −w/2 0 The separation position will be deﬁned as the integration boundary that yields one sixth of the total ﬂow rate. The total ﬂow 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 ﬂow rate that is one sixth of the total ﬂow 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 ﬂow 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 ﬂow 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 inﬁnite sum is replaced by a ﬁnal sum the result is not exactly one sixth of the total ﬂow. 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 diﬀerence in the hydraulic resis- tance across the channels, or a diﬀerence in the pressure drop across them. Before the 58 CHAPTER 9. SEPARATION MECHANISM determinination of the optimal design parameters, the inﬂuence due to unwanted pres- sure diﬀerences 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 inﬂuence on the separation position. It will further be considered whether it inﬂuences on the result if the pressure variation exists in diﬀerent 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 ﬂow rate: ∆p Q= (9.7) Rhyd From this equation it is easily seen that when the pressure drop decreases, so does the ﬂow rate for the given channel. This decrement makes the ﬂow rates in the other channels increase, and therefore the separation position will increase according to this increment in ﬂow rate. We will expect that the ﬂow rate variation will depend linearly with the ﬂow 1 rate, with the linearity coeﬃcient Rhyd . In order to verify whether this expectation is valid. Some simulations have been per- formed for various pressure diﬀerences 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 inﬂuence 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 inﬂuence from the pressure variation. It is also seen that when the pressure diﬀerence 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 diﬀerence 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 inﬂuence on the result is more signiﬁcant. 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 buﬀer separation positions in the ﬁrst 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 ﬁnal 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. diﬀ. [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 inﬂuence 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 inﬂuence from the pressure variation. Separation position [µm] h = 30 µm h = 10 µm h = 3 µm Pressure diﬀerence [Pa] Figure 9.3: This is the separation position variation due to pressure variations diﬀerences in channel 6 for various damping heights. It is seen that lower channels limits the separation position variation. 60 CHAPTER 9. SEPARATION MECHANISM Pres. diﬀ. [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 diﬀerence in channel 6 or 4 changes the the result. It is seen that when the unwanted pressure diﬀerence is in channel 4 the eﬀect on the result is more signiﬁcant. 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 ﬂow rate. From this variation one could achieve the separation position which is necessary for blood cell separation. In order to limit the channel 1 ﬂow 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 diﬀerence in hydraulic resistance must be, it is necessary to reconsider the ﬂow 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 proﬁle in the pinched segment is independent of the arrangements of the outlet channels. The velocity proﬁle stated in Eq. (9.2), will be used to determine the channel 1 ﬂow rate that will result in a separation position around 4 µm. To determine this ﬂow 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 ﬂow 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 ﬂow rate is then divided by the total ﬂow rate stated in Eq. (9.4) to give the channel one ﬂow 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 ﬂow rate. The rest of the outlet channels needs to split the remaining ﬂow rate if they are assumed to be identical. This results in a ﬂow 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 ﬂow 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 ﬁnd 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 ﬂow rate ratio to decrease, because when the hydraulic resistance of channel 2 is decreased the ﬂow 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 ﬂow rates ratio of the total ﬂow 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 ﬂow 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 eﬀect. Therefore long and narrow outlet channels are preferred. To make the separation mechanism as ﬂexible 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 deﬁning the height and width of the pinched segment the total ﬂow rate is determined. Then the channel 1 ﬂow rate determined by the given separation position is calculated. In the worksheet there are two ways of calculating the hydraulic resistances. In the ﬁrst simpler method the user can control the number of outlet channels and the separation position between channel 1 and 2. The channel 1 ﬂow rate is determined using Eq. (9.3). Then is the channel 1 ﬂow rate ratio of the total ﬂow rate determined, and the remaining part of the total ﬂow is divided equally between the remaining channels. The relative ﬂow rate diﬀerence is found by dividing the ﬂow rate ratio of the remaining channels with the channel 1 ﬂow 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 ﬂow 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 ﬂow 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 ﬂow rate ratios with the ﬂow 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 deﬁned 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 ﬂow 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 ﬁgure 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 ﬂow 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 oﬀ from the simulations by 29 %, 67 % and 69 % for design A, B and C respectively. The method used in this thesis were oﬀ 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 ﬁt 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 ﬁnal device , such as the particle motion around corners, narrowing of the pinched thickness and inﬂuence 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 proﬁle or streamline motion. This treatment is not suﬃcient for all purposes. For instance does the way of treating the particles not respect the fact that the particles have a ﬁnite size. An eﬀect that needs to be investigated is how the ﬂuid ﬂow 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 ﬁgure 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 ﬁnite particle sizes the streamlines closest to the wall will be of no importance. Due to this smooth ﬂow 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 ﬁnal 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 eﬀect 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 ﬂow 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 eﬀective separation radius would increase since particles cannot move closer to the wall than their radius. In order to investigate the eﬀect 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 ﬂow 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 proﬁle quickly established it will not make any diﬀerence 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 ﬁnal 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 ﬂow 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 ﬂow 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 eﬀect. Let us consider the narrowing segment which is shown in Fig. 10.2. The narrowing eﬀect is derived in the below calculations. It is assumed that the ﬂow 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 diﬀerence in channel depth will inﬂuence on the ﬂow proﬁle. In Chapter 4 it was described that when thinner channels are used, the velocity proﬁle will be limited, in contrast to the velocity proﬁle for inﬁnitely 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 ﬁnal device. Since the largest white blood cell can be up to 20 µm, 30 µm high channels will be used. When working with diﬀusion in channel with a ﬁnite thickness the so called butterﬂy eﬀect needs to be considered. This eﬀect and its inﬂuence on the diﬀusion proﬁle is described by A. E. Kamhols et al. [12]. When diﬀusion takes place in channels with a limited height the velocity proﬁle due to the hight of the channels will inﬂuence on the concentration proﬁle. The actual diﬀusion process will not change, it is just the circumstances around it which leads to the non constant diﬀusion proﬁle. 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 ﬁgure shows the concentration proﬁle in a channel having a rectangular cross section due to the butterﬂy eﬀect. Due to the higher ﬂuid velocity at the center of the channel the concentration proﬁle becomes hour-glass shaped. the ﬂuid at the center of the channel will have a higher velocity than the ﬂuid at the top- and bottom walls. This velocity diﬀerence will result in a change in the equilibrium concentration proﬁle. As the particle via diﬀusion 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 eﬀective diﬀusion time. Due to this eﬀect the diﬀusion proﬁle will be wider at the bottom and the top of the channel, where the velocity diﬀerence towards the horizontal center of the channel will be less than the diﬀerence at the center of height. This will then result in a hourglass chapped concentration proﬁle, see Fig. 10.3 According to [12] this eﬀects minimizes for low channels. This eﬀect therefore implies as well that channels with heights of 20 µm to 30 µm should be used in the ﬁnal 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 ﬂow 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 ﬂow 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 ﬁnal 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 ﬂow 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 ﬂow 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 ﬁnal design the outlet channels should all be made so that they are rather long and narrow in order to avoid unwanted eﬀects 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 Deﬁnition 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 ﬂows 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 ﬂow from the inlet channel must be forced up against one of the sidewalls by the larger ﬂow 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 ﬂows 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 diﬀusion 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 deﬁned 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 deﬁned as the path a small volume of ﬂuid 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 ﬂow. With the use of diﬀerent 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 ﬁgure 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 ﬂow 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 ﬂow rates cannot be changed enough to fulﬁll 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) diﬀusion has been neglected. When this is not a valid simpliﬁcation more complicated methods to determine the ﬂow rate ratio required for particle alignment must be used. The requirement stated in Eq. (11.1) is valid as long as the diﬀusion 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 Buﬀered Saline buﬀer. 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 diﬀusion the two ﬂowrates 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 ﬂows. From Chapter 8 it was found that in order to get the smallest FWHM rather large ﬂow rate ratios, which is the inlet channel ﬂow rate divided by the main channel ﬂow rate, should be used. this not necessarily the optimal design parameter, since the intersection width (w′ ), which is the width of the particle containing ﬂow at the end of the pinched segment, increases with growing ﬂow rate ratios. When the diﬀusion constant exceeds 1 × 10−13 m2 /s diﬀusion needs to be taken into account. One must then make simulations on the simpliﬁed model, which is presented in Chapter 8. But as long as a diﬀusion constant below 1×10−13 m2 /s are used Eq. (11.1) can be used to determine the necessary ﬂow rate ratios. Instead of calculating the necessary ﬂow 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 ﬂow 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 ﬁgure 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 ﬂow 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 ﬂow rate ratios than this value, since this would just insure the alignment even further. But ﬂow rate ratios below 0.02 is not recommended since this small ﬂow rate ratios will force the diﬀusion process far into the inlet channel, which will increase the velocity complexity. The position of the particles then becommes more diﬃcult to control. If it is not possible to achieve the alignment of particles with ﬂow 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 ﬂow 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 diﬀerence 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 ﬁgure is used to determine the ﬂow rate ratios between the inlet- and main ﬂow 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 ﬂexible 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 eﬀects 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 eﬀects due to unwanted pressure variations. In a pinched ﬂow one would like the ﬂow to be a creeping ﬂow in the separation chamber in order to avoid rolls that would interfere with the sorting mechanism. A creeping ﬂow is a ﬂow that creeps around corners and thereby has no rolls. A way of insuring that the ﬂow will be almost like creeping ﬂow is to use rather thin channels throughout the device. This is because the limited height will act as a damping factor on the ﬂuid velocity, due to the ”no-slip” boundary condition on the top and bottom of the channels. This damping factor will make the ﬂuid less easy ﬂowing, 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 speciﬁc 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 ﬂow rates, a ﬂow 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 ﬂow rate ratio. It can be used to ﬁnd the required ratios to design a optimal system. the red dashed line, it is seen that a very low ﬂow 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 ﬂow 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 diﬀerence 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 ﬂow 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 ﬂow, 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 modiﬁed, so that it takes the ”no-slip” boundary condition on the channel sidewalls into account. It was found that the ﬂuid behavior was controlled mainly by the relative ﬂow 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 ﬁrst 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 ﬁnding 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] Bioﬂuid 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 microﬂuidic device equipped with ﬂow 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 ﬂows B. P. Ho & L. G. Leal, J. Fluid. Mech. (1974), vol. 65, part 2, pp. 365-400 [7] Mathematics for physicists S. M. Lea, ﬁrst 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 Proﬁle 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 Diﬀusion in Pressure-Driven Laminar Flow in Microﬂuidic Channels A. E. Kamholz & P. Yagar, Biophysical journal, January 2001, vol. 80, 155-160 [13] Theoretical microﬂuidics 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 Deﬁnition The goal for this simulation was to determine the separation position for three diﬀerent geometries of the outlet channels. The tree diﬀerent 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 diﬀerence 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