Passive Mixing in Microchannels - Thesis - 2004

Passive Mixing in Microchannels with Geometric
Variations
by
Hengzi Wang
B.S. (Shanghai Jiao Tong University) 1990

M.S. (Shanghai Jiao Tong University) 1993
A thesis submitted to Industrial Research Institute Swinburne (IRIS) in

fulfillment of the requirements for the degree of Doctor of Philosophy
SWINBURNE UNIVERSITY OF TECHNOLOGY,

HAWTHORN, VICTORIA
AUSTRALIA
January, 2004
Passive Mixing in Microchannels with Geometric
Variations
Declaration
This thesis contains no material which has been accepted for the award of any other de-
gree or diploma at any university and to the best of my knowledge and belief contains no
material previously published or written by another person or persons except where due
reference is made.
Hengzi Wang September 19, 2004

i
Abstract
This research project was part of the microfluidic program in the CRC for Microtechnology,
Australia, during 2000 to 2003. The aim of this research was to investigate the feasibility of
applying geometric variations in a microchannel to create effects other than pure molecular
diffusion to enhance microfluidic mixing. Geometric variations included the shape of a
microchannel, as well as the various obstacle structures inside the microchannel.
Generally, before performing chemical or biological analysis, samples and reagents

need to be mixed together thoroughly. This is particularly important in miniaturized To-
tal Analysis Systems (µTAS), where mixing is critical for the detection stage. In scaling
down dimensions of micro-devices, diffusion becomes an efficient method for achieving
homogenous solutions when the characteristic length of the channels becomes sufficiently
small. In the case of pressure driven flow, it is necessary to use wider microchannels to
ensure fluids can be pumped through the channels and the volume of fluid can provide
sufficient signal intensity for detection. However, a relatively wide microchannel makes
mixing by virtue of pure molecular diffusion a very slow process in a confined volume of
a microfluidic device. Therefore, mixing is a challenge and improved methods need to be
found for microfluidic applications.
In this research, passive mixing using geometric variations in microchannels was stud-
ied due to its advantages over active mixing in terms of simplicity and ease of fabrication.
Because of the nature of laminar flow in a microchannel, the geometric variations were de-
signed to improve lateral convection to increase cross-stream diffusion. Previous research
using this approach was limited, and a detailed research program using computational fluid
dynamic (CFD) solvers, various shapes, sizes and layouts of geometric structures was un-
dertaken for the first time. Experimental measurements, published experimental data and
ii
analytical predictions were used to validate the simulations for selected samples. Mixing
efficiency was evaluated by using mass fraction distributions. It was found that the overall
performance of a micromixer should include the pressure drop in a microdevice, therefore,
a mixing index criterion was formulated in this research to combine the effect of mixing
efficiency and pressure drop. The mixing index was used to determine optimum parameters
for enhanced mixing, as well as establish design guidelines for such devices.
Three types of geometric variations were researched. First, partitioning in channels

was used to divide fluids into mixing zones with different concentrations. Various designs
were investigated, and while these provided many potential solutions to achieving good
mixing, they were difficult to fabricate. Secondly, structures were used to create lateral
convection, or secondary flows. Most of the work in this category used obstacles to disrupt
the flow. It was found that symmetric layouts of obstacles in a channel had little effect
on mixing, whereas, asymmetric arrangements created lateral convection to enhance cross-
stream diffusion and increase mixing. Finally, structures that could create complex 3D
advections were investigated. At high Reynolds numbers (Re = 50), 3D ramping or ob-
stacles generated strong lateral convection. Microchannels with 3D slanted grooves were
also investigated. Mixers with grooved surfaces generated helicity at low Reynolds num-
bers (Re ≤ 5) and provided a promising way to reduce the diffusion path in microchannels
by stretching and folding of fluid streams. Deeper grooves resulted in better mixing effi-
ciency. The 3D helical advection created by the patterned grooves in a microchannel was
studied by using particle tracing algorithms developed in this research to generate streak-
lines and Poincare maps, which were used to evaluate the mixing performance. The results
illustrated that all the types of mixers could provide solutions to microfluidic mixing when
dimensional parameters were optimized.
iii
Acknowledgments
First of all, I want to thank my family. It is a long and difficult period for both my wife,
Hongtao and my son, Thomas. My wife takes care of all the family issues while she was in
Shanghai and many house-works. Thomas always reminds me his existence and I should
accompany with him more. He loves to play with my computer and my car. Except these,
he is a lovely boy and gives me lots of comforts.
I wish to thank my supervisors, Dr. Pio Iovenitti, Professor Erol Harvey and Associate
Professor Syed Masood for their advice and guidance to complete this research project.
All my colleagues have been friendly and helpful. They provided a comfortable envi-
ronment for researchers. We have about 30 colleagues now, and I may not be able to put
all the names here. But I would like to thank the colleagues who gave me advice about my
research project and on personal issues. Andrew Dowling and I have discussed all sorts of
problems and shared the joy of successes - when we have them. Yao Fu, Tony Liu and Irina
Simdikova are lunchtime chatting mates. Sometimes we do discuss research and science!
I consulted with Dr. Rowan Deam on fluid mechanics several times, and I thank him
for helpful advice. I wish to express my gratitude for the helpful discussion about chaotic
mixing with Dr. A. Stroock from Harvard University, and Dr. M. Rudman and Dr. G.
Metcalfe from CSIRO in Australia.
I visited microfluidic diagnostics laboratory in Stanford University, and I wish to thank

Dr. Santiago for spending much time on discussing µPIV systems, and also his students,
Michael Oddy, Klint Rose, Chuanhua Chen and Jung Byoungsok for talking about their
research projects and showing me around their laboratory.
Some papers in the references were not possible to obtain from the local libraries. I
should thank the authors for sending copies of their papers. Just name a few: Professor
David J. Beebe in University of Wisconsin and Dr. Jin-Woo Choi in University of Cincin-
iv
nati.
This research was mostly funded by CRC for Microtechnology, and I wish to thank
the CRC as well.
v
Contents
Abstract i
Acknowledgments iii
List of Figures xi
List of Tables xviii
Nomenclature xix
List of Variables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xix
List of abbreviations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxii
1 Introduction 1
1.1 Background and Significance of the Research . . . . . . . . . . . . . . . . 1
1.2 Statement of the Problem . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.3 Objectives and Scope of the Research . . . . . . . . . . . . . . . . . . . . 3
1.4 Outline of this thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
2 Literature Review 9
2.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
vi
2.2 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.3 Fundamental Theories for Microfluidics . . . . . . . . . . . . . . . . . . . 11
2.3.1 Continuum . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.3.2 Conservation of Momentum . . . . . . . . . . . . . . . . . . . . . 13
2.3.3 Brownian Motion . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.3.4 Taylor-Aris dispersion . . . . . . . . . . . . . . . . . . . . . . . . 15
2.3.5 Chaotic Advection . . . . . . . . . . . . . . . . . . . . . . . . . . 16
2.4 Generation of Fluid Motion in a Microfluidic Device . . . . . . . . . . . . 19
2.4.1 Mechanically Actuated Micropumps . . . . . . . . . . . . . . . . . 20
2.4.2 Electrokinetic Micropumps . . . . . . . . . . . . . . . . . . . . . . 23
2.5 Active Micro-Mixers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
2.5.1 Planar Laminar Bubble Mixer . . . . . . . . . . . . . . . . . . . . 27
2.5.2 Periodic Perturbation Applied to a Main Fluid Stream . . . . . . . 27
2.5.3 Peristaltically Driven Micromixing . . . . . . . . . . . . . . . . . 29
2.5.4 Magnetohydrodynamic (MHD) Micro-mixers . . . . . . . . . . . . 30
2.5.5 Electrohydrodynamic (EHD) Micro-mixers . . . . . . . . . . . . . 32
2.5.6 Magnetic Micro-Stirrer . . . . . . . . . . . . . . . . . . . . . . . . 32
2.6 Passive Microfluidic Mixers . . . . . . . . . . . . . . . . . . . . . . . . . 36
2.6.1 T-type mixers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
2.6.2 Splitting and Recombining . . . . . . . . . . . . . . . . . . . . . . 41
2.6.3 Focusing Fluid Streams . . . . . . . . . . . . . . . . . . . . . . . . 46
2.6.4 Inject Small Fluid Substreams into a Main Stream . . . . . . . . . 47
2.6.5 Passive Microfluidic Chaotic Mixers . . . . . . . . . . . . . . . . . 48
2.7 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
vii
3 Numerical Modelling and Experimental Methodology 57
3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
3.2 Fundamentals of Computational Fluid Dynamics (CFD) . . . . . . . . . . . 59
3.2.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
3.2.2 CFD Solvers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
3.3 CFD Modelling in Microfluidics . . . . . . . . . . . . . . . . . . . . . . . 63
3.3.1 System of Units in Microfluidics . . . . . . . . . . . . . . . . . . . 64
3.3.2 Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . . . . 64
3.3.3 Building Numerical Models for Micro-Flows . . . . . . . . . . . . 65
3.3.4 Simulation Results Analysis . . . . . . . . . . . . . . . . . . . . . 71
3.4 Experimental Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
3.4.1 Excimer Laser Micro-machining . . . . . . . . . . . . . . . . . . . 74
3.4.2 Microfluidic Packaging . . . . . . . . . . . . . . . . . . . . . . . . 75
3.4.3 Microfluidic Mixing Visualization . . . . . . . . . . . . . . . . . . 75
3.5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
4 Understanding Microfluidic Mixing 79
4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
4.2 Flow at Low Reynolds Number . . . . . . . . . . . . . . . . . . . . . . . . 80
4.3 Taylor-Aris Dispersion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
4.3.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
4.3.2 Dispersion of a Solute in Solvent Flowing Slowly Through a Rect-

angular Cross-Section Microchannel . . . . . . . . . . . . . . . . . 86
4.3.3 Applying Taylor-Aris Dispersion to Enhance Mixing . . . . . . . . 89
4.4 Two Important Parameters for Designing a Microfluidic Mixer . . . . . . . 91
4.4.1 Péclet number . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92

viii
4.4.2 Aspect Ratio of a Rectangular Microchannel . . . . . . . . . . . . 95
4.5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98
5 Microfluidic Mixing by Partitioning Flows 101
5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101
5.2 Dividing Fluids into Discrete Concentration Zones . . . . . . . . . . . . . 102
5.3 Dollar Sign Shaped Micromixers . . . . . . . . . . . . . . . . . . . . . . . 104
5.4 3D Ramping Micro-mixers . . . . . . . . . . . . . . . . . . . . . . . . . . 109
5.4.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109
5.4.2 Simulations of 3D Ramping Micro-mixers . . . . . . . . . . . . . . 110
5.4.3 Discussion of Results . . . . . . . . . . . . . . . . . . . . . . . . . 110
5.5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112
6 Optimizing Layout of Cylindrical Obstacles in Microchannels for Enhanced

Mixing 116
6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116
6.2 Numerical Modelling of Mixing in a Microchannel with Cylindrical Ob-

stacles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117
6.3 Results and Discussions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120
6.3.1 Simulation Results . . . . . . . . . . . . . . . . . . . . . . . . . . 121
6.3.2 Experimental Results . . . . . . . . . . . . . . . . . . . . . . . . . 126
6.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128
7 Enhancing Mixing in Microchannels Using Square and Rectangular Obstacles129
7.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129
7.2 Numerical Modelling of Square and Rectangular Obstacles . . . . . . . . . 130

ix
7.3.1 Number of Obstacles . . . . . . . . . . . . . . . . . . . . . . . . . 133
7.3.2 Size of Obstacles . . . . . . . . . . . . . . . . . . . . . . . . . . . 135
7.3.3 Angle of Obstacles . . . . . . . . . . . . . . . . . . . . . . . . . . 139
7.3.4 Obstacle Offsets from Channel Walls . . . . . . . . . . . . . . . . 143
7.3.5 Gap between Obstacles . . . . . . . . . . . . . . . . . . . . . . . . 145
7.3.6 Flow Rate . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 148
7.3.7 Ratio of Height of Obstacles to Depth of Channel . . . . . . . . . . 150
7.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152
8 Numerical Investigation of Mixing in Microchannels with Patterned Grooves 154
8.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 154
8.2 Numerical Setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 155
8.3 A Particle Tracing Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . 158
8.3.1 Point location . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 159
8.3.2 Interpolate velocity . . . . . . . . . . . . . . . . . . . . . . . . . . 160
8.4.1 Flow patterns . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 162
8.4.2 Streaklines . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 163
8.4.3 Poincaré maps . . . . . . . . . . . . . . . . . . . . . . . . . . . . 165
8.4.4 Simulation of Mixing in Microchannels with Patterned Grooves . . 169
8.5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 173
9 Conclusions and Further Research 175
9.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 175
9.2 Contributions of the Research . . . . . . . . . . . . . . . . . . . . . . . . 175
9.2.1 Clarification of Related Theories to Microfluidic Mixing . . . . . . 175

x
9.2.2 Modifications to Serpentine-Shaped Microchannels . . . . . . . . . 178
9.2.3 Design Guidelines for Obstacle Structures . . . . . . . . . . . . . . 179
9.2.4 Design of Microfluidic Mixers with Patterned Grooves . . . . . . . 181
9.3 Limitations of the Research . . . . . . . . . . . . . . . . . . . . . . . . . . 182
9.4 Recommendations for Further Research . . . . . . . . . . . . . . . . . . . 183
9.4.1 Experimental Methods . . . . . . . . . . . . . . . . . . . . . . . . 183
9.4.2 Complex 3D Geometry . . . . . . . . . . . . . . . . . . . . . . . . 184
9.5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 185
References 188
Appendices 200
A List of Publications 200
A.1 Journal Publications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 200
A.2 International Conference Publications . . . . . . . . . . . . . . . . . . . . 201
A.3 Local Publications and Presentations . . . . . . . . . . . . . . . . . . . . . 202
B Selected Fortran Codes Developed to Study Microfluidic Mixing 203
B.1 Fortran Codes for Mesh Conversion . . . . . . . . . . . . . . . . . . . . . 203
B.2 Fortran Codes for Velocity Field in a Rectangular Duct . . . . . . . . . . . 211
B.3 Fortran Codes for Particle Tracing - Main Program . . . . . . . . . . . . . 215

xi
List of Figures
2.1 Statical fluid thermophysical point value . . . . . . . . . . . . . . . . . . . 12
2.2 Stirring of a ‘blob ’of marked fluid by 2D, time-dependent (time series
from t0 ∼ t6), laminar flow . . . . . . . . . . . . . . . . . . . . . . . . . . 17
2.3 Flow rates and pressure heads of micropumps could deliver - a review . . . 20
2.4 Disk piezoelectric reciprocating micropump . . . . . . . . . . . . . . . . . 21
2.5 A typical passive micro check valve . . . . . . . . . . . . . . . . . . . . . 22
2.6 Tesla micropump . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
2.7 2D chaotic advection microfluidic mixer . . . . . . . . . . . . . . . . . . . 28
2.8 Microfluidic mixer by periodical perturbation . . . . . . . . . . . . . . . . 29
2.9 Mechanical travelling waves induced by acoustic streaming . . . . . . . . . 30
2.10 A minute MHD micro-mixer . . . . . . . . . . . . . . . . . . . . . . . . . 31
2.11 (a). Schematic illustration of the active micro-mixer (b). Cross sectional
view of the active micro mixer . . . . . . . . . . . . . . . . . . . . . . . . 33
2.12 Mixing processes in the micro-mixer by applying electrokinetic instability . 34
2.13 Magnetic micro stirrer . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
2.14 Mixing in microchannel by magnetic micro stirrer, (a). 0 rpm, (b). 150
rpm, and (c). 300 rpm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
2.15 Mixing T structures, (a). Parallel configuration; (b). Serial configuration . . 37
2.16 Streamlines of engulfment flow . . . . . . . . . . . . . . . . . . . . . . . . 38

xii
2.17 Double-T micro-reactor . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
2.18 A schematic working principle of T-Sensor . . . . . . . . . . . . . . . . . 40
2.19 Parallel dividing into substreams. . . . . . . . . . . . . . . . . . . . . . . . 42
2.20 In-situ static micro-mixer by subdividing streams . . . . . . . . . . . . . . 43
2.21 Divided flow into substreams (a) vertical arrangement; (b) horizontal ar-
rangement. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
2.22 (a). Channel structures; (b). Schematic drawing of its split-recombine-

principle; (c). Streamlines in a caterpillar mixer. . . . . . . . . . . . . . . . 44
2.23 Danfoss mixer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
2.24 Static Micro-mixer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
2.25 IMM super-focusing micro-mixer . . . . . . . . . . . . . . . . . . . . . . 47
2.26 Hydrodynamic focusing micro-mixer . . . . . . . . . . . . . . . . . . . . . 48
2.27 IMM-Mainz Mixing unit . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
2.28 Injection fluid into a mixing chamber . . . . . . . . . . . . . . . . . . . . . 49
2.29 (a). Mixing Chamber with Micro-plumes (b). Simulation of sub-jets of

reagent into main stream . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
2.30 Serpentine-Shaped Micro-mixers . . . . . . . . . . . . . . . . . . . . . . . 51
2.31 Microfluidic chaotic mixer with grooved surfaces . . . . . . . . . . . . . . 52
3.1 Viscous flow in a rectangular micro-channel . . . . . . . . . . . . . . . . . 67
3.2 A typical meshed fluid volume in a rectangular micro-channel . . . . . . . 68
3.3 Normalized velocity profile u(y, zmid ) in xOy plane . . . . . . . . . . . . . 68
3.4 Normalized velocity profile u(ymid , z) in xOz plane . . . . . . . . . . . . . 69

TM
3.5 Mesh topology, (a). Gambitr 8-Hexahedral mesh; (b). MemCFD 8-
Hexahedral mesh . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
3.6 (a). Structured 8-node Hexahedral mesh by sub-mapping technique; (b).

TM
Irregular 8-node structured mesh by MemCFD . . . . . . . . . . . . . . 70
xiii
3.7 Mesh uniformity study, (a). irregular mesh elements; (b). uniform mesh
elements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
3.8 Mixing of two fluids in a T-channel with grooved surfaces, top: simulation
result, two fluids are represented by blue and red, green colour is the mixed
region, bottom: cited experimental results by Stroock, et al. 2002a, two
fluids are represented by gold and black colours. . . . . . . . . . . . . . . 73
3.9 A schematic illustration of Excimer laser micromachining . . . . . . . . . 74
3.10 Section of Y-channel machined in Polycarbonate by Excimer Laser (SEM

image) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
3.11 Laboratory setup for micro-mixing visualization . . . . . . . . . . . . . . . 76
4.1 Rectangular microchannel . . . . . . . . . . . . . . . . . . . . . . . . . . 81
4.2 Normalized velocity profile in a rectangular microchannel, w/H = 4; u∗ is

normalized velocity, y ∗ = y/w and z ∗ = z/H . . . . . . . . . . . . . . . . . 83
4.3 Normalized velocity profile by the maximum flow velocity plotted in XZ

plane for y = ymid . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
4.4 Normalized velocity profile plotted in XY plane for w/H = 2 (O), w/H =
5 (×) and w/H = 10 (+). . . . . . . . . . . . . . . . . . . . . . . . . . . 85
4.5 Ratio of Taylor dispersivity to pure molecular diffusion coefficient for dif-
ferent Reynolds numbers . . . . . . . . . . . . . . . . . . . . . . . . . . . 88
4.6 Taylor dispersion along a straight channel: numerical simulation for vari-
ous time intervals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90
4.7 Experimental results: the solute dispersed at, (a) time t and (b) time t + δt . 91
4.8 Three microfluidic mixing zones defined by using the Péclet number . . . . 94
4.9 Pressure drop along a rectangular microchannel with the same cross-section
area, for different height to width aspect ratio, at various flow rates q . . . . 97
4.10 Aspect Ratio and flow rate affect the total length to achieve complete mixing 97
5.1 Different concentration gradient separated in the branches, (a). relative

mass fraction distribution in percentage; (b). corresponding mass fraction
in contour plot . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103
xiv
5.2 Mass fraction in a planar $-mixer . . . . . . . . . . . . . . . . . . . . . . . 105
5.3 Streaklines mapped with mass fraction in the planar $-mixer . . . . . . . . 106
5.4 Mass fraction in 3D variation of the planar $-mixer . . . . . . . . . . . . . 107
5.5 Streaklines mapped with mass fraction in the 3D $-mixer . . . . . . . . . . 108
5.6 Two samples of 3D ramping structures made by Excimer Laser microma-

chining. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109
5.7 Computer model of a ramping structure meshed with 8-node Hexahedrons . 110
5.8 A 3D ramping mixer, at Re = 5 . . . . . . . . . . . . . . . . . . . . . . . . 111
5.9 Mass fraction distribution in a 3D ramping mixer, at Re = 50 . . . . . . . . 113
5.10 Streaklines in 3D ramping mixer mapped with mass fraction, at Re = 50 . . 114
6.1 3D model of a T-channel with two arrays of cylindrical obstacles . . . . . . 118
6.2 Layout of square and triangular configuration of obstacle array. (a). con-
figuration No.6. (b). configuration No.7 (c). configuration No.8 . . . . . . . 118
6.3 Meshed fluid volume in a T-type channel with two arrays of cylindrical
obstacles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119
6.4 Relative mass fraction distribution in solution in a microchannel with two

arrays of cylindrical obstacles. . . . . . . . . . . . . . . . . . . . . . . . . 121
6.5 Relative mass fraction distribution in solution in a microchannel with two

arrays of cylindrical obstacles. . . . . . . . . . . . . . . . . . . . . . . . . 122
6.6 Evolution of design with obstacles in the Y-channel for configurations 1 to 8 123
6.7 Finite element simulation results of (a) mixing efficiency, pressure drop
between inlet and outlet versus number of obstacles, and (b) Mixing index
versus number of obstacles(configuration no.1 ∼ no.8). . . . . . . . . . . . 124
6.8 Section of the channel illustrated the convective effect using velocity vector
field . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125
6.9 Channel with cylindrical obstacles (SEM image) . . . . . . . . . . . . . . 126
6.10 Configuration no.1, channel without cylindrical obstacles (food colouring),

Re = 0.4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127
xv
6.11 Configuration no.5, channel with 3 cylindrical obstacles (Ferric Ammo-

nium Sulfate and Ammonium Thiocyanate solutions), Re= 0.4 . . . . . . . 127
6.12 Configuration no.8, channel with 18 cylindrical obstacles (food colouring),

Re = 0.4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127
7.1 T-channel with 4 rectangular obstacles, (a) 3D view; (b) layout. . . . . . . . 131
7.2 Mass fraction distributions, (a) symmetric layout of obstacles, (b) asym-
metric layout of obstacles . . . . . . . . . . . . . . . . . . . . . . . . . . . 132
7.3 Mixing efficiency versus number of square obstacles . . . . . . . . . . . . 133
7.4 Pressure drop versus number of square obstacles . . . . . . . . . . . . . . . 134
7.5 Mixing index, midx , versus number of square obstacles . . . . . . . . . . . 135
7.6 Angle between flow axis and obstacle axis, two and three obstacles . . . . . 136
7.7 Mixing efficiency with two and three obstacles inside a microchannel ver-
sus angle of layout . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137
7.8 Pressure drop versus two and three square obstacles . . . . . . . . . . . . . 137
7.9 Mixing index, midx , versus two and three square obstacles . . . . . . . . . 138
7.10 Mass fraction distribution for two rectangular obstacle layouts with differ-
ent size of obstacles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138
7.11 Mixing efficiency versus size of obstacles . . . . . . . . . . . . . . . . . . 139
7.12 Pressure drop versus size of obstacles . . . . . . . . . . . . . . . . . . . . 140
7.13 Mixing index, midx , versus size of obstacles . . . . . . . . . . . . . . . . . 140
7.14 Mass fraction distribution, various angles of obstacle layouts . . . . . . . . 141
7.15 Mixing efficiency versus angle of rectangular obstacles . . . . . . . . . . . 141
7.16 Pressure drop versus angle of obstacles . . . . . . . . . . . . . . . . . . . . 142
7.17 Mixing index, midx , versus angle of obstacles . . . . . . . . . . . . . . . . 142
7.18 Mass fraction distribution, for different offsets of obstacles from channel
walls, (a). no offset from the center; (b). with offset. . . . . . . . . . . . . . 143
7.19 Mixing efficiency versus offset from the wall . . . . . . . . . . . . . . . . 144

xvi
7.20 Pressure drop versus offset of obstacles from wall . . . . . . . . . . . . . . 144
7.21 Mixing index, midx , versus offset of obstacles from wall . . . . . . . . . . 145
7.22 Mass fraction distribution for different obstacle gaps, (a). 160 µm gap; (b).
400µm gap; . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 146
7.23 Mixing efficiency versus gap between obstacles . . . . . . . . . . . . . . . 146
7.24 Pressure drop versus gap between two obstacles . . . . . . . . . . . . . . . 147
7.25 Mixing index, midx ,versus gap, g, between two obstacles . . . . . . . . . . 147
7.26 Mixing efficiency versus flow rates . . . . . . . . . . . . . . . . . . . . . . 148
7.27 Pressure drop versus Reynolds number . . . . . . . . . . . . . . . . . . . . 149
7.28 Mixing index, midx , versus Reynolds number, Re . . . . . . . . . . . . . . 149
7.29 Mass fraction versus obstacle height to channel depth aspect ratio, h/H,
(a). h/H = 0.2, (b). h/H = 0.9 . . . . . . . . . . . . . . . . . . . . . . . 150
7.30 Mixing efficiency versus height to depth aspect ratio of obstacles . . . . . . 151
7.31 Pressure drop versus obstacle height to channel depth ratio . . . . . . . . . 151
7.32 Mixing index, midx , versus obstacle height to channel depth ratio, h/H . . 152
8.1 Fluid volume in microchannel with grooved bottom surface. L0 is the entry
channel length, Lr is the grooved channel length and Lm is the channel
length to complete the mixing by diffusion . . . . . . . . . . . . . . . . . . 156
8.2 Meshed fluid volume . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 157
8.3 Hexahedral cell decomposed into 5-Tetrahedrons . . . . . . . . . . . . . . 160
8.4 Tetrahedron geometry. Ni , i = 1, 4. are the shape functions . . . . . . . . . 161
8.5 One periodic section of patterned grooves, with section planes labelled
from 0 ∼ 2π, and particles advected from plane 0 to plane 2π . . . . . . . . 163
8.6 Cross section velocity vector planes, 0 ∼ 2π, Re = 5 . . . . . . . . . . . . . 164
8.7 Streaklines of two neighboring fluid particles from two fluid streams re-
spectively, α = 0.20, Re = 5 . . . . . . . . . . . . . . . . . . . . . . . . . 165
8.8 Gap between the two Streaklines in Figure 8.7 . . . . . . . . . . . . . . . . 166

xvii
8.9 Poincaré map, 45o patterned grooves: (a) α = 0.05, Re = 5; (b) α = 0.30,
Re = 5; (c) α = 0.30, Re = 0.01. . . . . . . . . . . . . . . . . . . . . . . . 167
8.10 Length of grooved channel to complete one circulation, for H/w = 0.5 . . 168
8.11 Two fluid streams flow into a T-type channel with patterned grooves, α =
0.10 and Re = 5. (a). The interfacial line is indicated by the bright mixed
region. The helicity was measured by the angle Ω between this interfacial
line and the channel axis, (b). The visualization of mixing . . . . . . . . . . 170
8.12 Interfacial regions when α = 0.30 and Re = 5 . . . . . . . . . . . . . . . 171
8.13 Mean helicity measured by the angle Ω for H/w = 0.5 . . . . . . . . . . . 171
8.14 Helicity at different Reynolds numbers, numerical estimations for α = 0.10 172
B.1 Flow chart to convert gambit neutral mesh to universal mesh . . . . . . . . 205
B.2 Flow chart for viscous flow in a rectangular channel . . . . . . . . . . . . . 212
B.3 Flow chart to compute particle trajectories . . . . . . . . . . . . . . . . . . 216

xviii
List of Tables
3.1 Systems of units in microfluidics . . . . . . . . . . . . . . . . . . . . . . . 65
4.1 Values of k in a rectangular-cross-section channel for various width to

depth aspect ratios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
6.1 Configuration of microchannel with Cylindrical obstacles . . . . . . . . . . 117
6.2 Properties of water and ethanol at 20 ˚ C . . . . . . . . . . . . . . . . . . . 120

xix
Nomenclature
The symbols listed in the chapters follow the convention of international system of
units (SI), SI base units and SI derived units. In practice, some adjustment of the system of
units are described in Chapter 3 for the convenience of microfluidic applications.
List of Variables
Symbol Description SI unit

A Surface of control volume m2
c relative mass fraction
dh Hydraulic diameter m
D Diffusion coefficient m2 s−1
D∗ Virtual (Taylor) dispersivity m2 s−1
E Electrical field V
e Internal energy Nm
F Force N
f Body force N
FD Friction force between liquid and wall N
g Gravitational acceleration m s−2
H Depth of channel m
xx
h Height of obstacles (or depth of groove) m

I Moment of Inertia kg m2
J Capillary force along the channel N
k constant in Equation 4.6
Kn Knudsen number
l Characteristic length m
lm Channel length for a complete mixing m
L Length of the channel m
m Mass kg
mef f mixing efficiency
midx mixing index
N number of molecules in a fluid point
n number of samples
Ni ith shape function
p Pressure Pa
p(t) particle position m
P Power W
Pe Péclet number
Pn Surface force N
q volumetric flow rate m3 s−1
qλ Conductive heat Nm
qR Radiative heat Nm
r radius of channel m
Re Reynolds number
t Time s
tc Convection time s
tD Diffusion time s
T Temperature K
xxi
TCE Thermal Coefficient of Expansion K −1

um Mean velocity in the channel m s−1
v Velocity field m s−1
v1 velocity along x-axis m s−1
v2 velocity along y-axis m s−1
v3 velocity along z-axis m s−1
w width of channel m
Wx Work from outside control volume Nm
x x-coordinate m
y y-coordinate m
z z-coordinate m
α aspect ratio of half groove depth to depth of channel
β aspect ratio of channel depth to channel width
ε permittivity C V −1 m−1
γ Surface tension factor N m−1
λ Mean free path of a molecule m
µ Dynamic viscosity Pa s
ν kinematic viscosity m2 s−1
ρ Density kg m−3
ρE Electric charge density C m−3
τ The lump of fluid in the control volume m3
Ω Angle degree
ψ The electrical potential V
φ Diameter m
ζEO electroosmotic zeta potential V
ζEP electrophoretic zeta potential V
xxii
List of abbreviations
The list of abbreviations used in this thesis is shown below.
Abbreviation Description
CAD Computer Aided Design
CFD Computational Fluid Dynamics
DRIE Deep Reactive Ion Etching
EDL Electrical Double Layer
FEA Finite Element Analysis
FEM Finite Element Method
FVM Finite Volume Method
IC Integrated Circuit
LIGA Lithographic Galvonoformung Abformung
(lithography, electroplating, molding)
MEMS MicroElectroMechanicSystem
MTTF Mean Time To Failure
µTAS Miniaturized total analysis systems
1
Chapter 1
Introduction
1.1 Background and Significance of the Research
This work was conducted in the Industrial Research Institute Swinburne (IRIS) dur-
ing the years from 2000 to 2003. It was part of the Microfluidics project supported by
Cooperative Research Centre (CRC) for Microtechnology in Australia.
Fluidic devices made by MicroElectroMechanicSystem (MEMS) technology can be

termed microfluidics, which have dimensions ranging from microns to a few millime-
ters (Gravesen, Branebjerg & Jensen 1993). Microfluidics research commenced several
decades ago. Early in the mid of 60s, pressure transducers started to utilize microfluidic
techniques. However, it is generally recognized that microfluidics originated at Stanford
University in the development a gas chromatograph system in 1975. Later in 1977, IBM
adopted the idea and incorporated it in the manufacture of inkjet printers (Petersen 1982).
In the years that followed, further research was done but at a quite moderate level . How-
ever, in 1990, Manz introduced the concept of miniaturized total chemical analysis systems
(µTAS) (Manz, Graber & Widmer 1990), which could shrink the whole chemical analysis
2
laboratory to an integrated microfluidic chip that could be held in a human hand. Since
then, the research in microfluidics has dramatically increased, showing its great potential
significance in a wide range of technologies, such as pressure sensors, medical diagnostic
systems, chemical synthesis.
The µTAS concept provides an overarching vision and compelling motivation for de-
velopment of integrated microfluidic systems. However, before fully developed systems
can be put in the market place, sufficient development of all the necessary components to
build such systems needs to be achieved. Among these individual components, the unit
for mixing samples and reagents is critical for the subsequent chemical reactions and de-
tections. The mixing unit has to provide sufficient mixed solution in a confined length
of a micro-device, and before the point where detection is taking place. Mixing in mi-
crofluidic devices can be achieved by active mixing methods and passive mixing methods.
Passive mixing using geometric variations, which includes the shape of a microchannel
and other obstacle structures inside the microchannel, provides a versatile solution to sim-
plify the design and fabrication of microfluidic device, and is totally compatible with the
current microfabrication techniques. Hence, it is a cost-effective mixing method and has
great potential in a practical µTAS. Therefore, as a contribution to building this critical unit
for µTAS, the author presents his research on passive microfluidic mixing with geometric
variations.
1.2 Statement of the Problem
To achieve a homogenous mixed solution of samples and reagents in a microfluidic

device, the dimensions such a device are too small to put reliable mechanical agitation to
stir the fluids. Naturally, the Reynolds numbers of microfluidic devices are low due to their
small dimensions. Hence, viscosity predominates the flows, especially for aqueous solu-
tions of most bio-materials. In such viscous flows, even though one can put mechanical
3
agitation into a microfluidic device, it is not likely to efficiently improve mixing perfor-
mance (Purcell 1977). Hence, mixing relies on the virtue of molecular diffusion.
On the one hand, microchannels are elegant and efficient for mixing fluids in narrow
channels with dimensions up to a few tens of microns. Mixing can be completed by dif-
fusion without any other stirring assistance. This makes microfluidics a very controllable
and predictable technique for rapid mixing. On the other hand, it is difficult to drive flow
through such small channels due to viscosity. Although fluid can be driven electrokineti-
cally, electrokinetically driven flows introduce some other problems, such as surface chem-
istry, gas bubbles created by electrolysis and Joule heating. Pressure driven mechanisms
have commonly been used in microfluidic devices, however, a relatively wider channel is
necessary to reduce the pressure drop. In addition, sometimes a sufficient volume of fluid
is needed for a detection stage, and this also requires a relatively large channel. For mi-
crochannels with dimensions over 100 µm, diffusion is not generally effective within the
confined space of micro-devices, especially for bio-materials with very low diffusion coef-
ficients. These require a long channel for adequate residual time for diffusion to take place.
However, the pressure drop is proportional to the length, and it is difficult to drive the flow
through a long microchannel. For the above reasons, methods have to be found to reduce
the diffusion distance between fluid streams without restricting the flow.
In brief, microfluidic mixing is challenging and attracts the attention of many re-
searchers.
1.3 Objectives and Scope of the Research
The principal objective of this research was to improve the passive mixing efficiency
using geometric configurations in microchannels.
Although many researchers have been attempting to apply active mechanisms to im-
4
prove mixing performance in a microfluidic device, the complexity of active mechanisms

makes it difficult to fabricate and expensive to make. Passive mechanisms have been stud-
ied extensively, however, previous research on passive mixing has concentrated on typically
one of geometric factors, and an overall view of geometric configurations and their effect
on mixing is not clear. Therefore, it is the intention of this research to provide a compre-
hensive research on passive mixing with various geometric configurations.
To achieve this goal, this research systematically investigated the relation between var-
ious geometric configurations and their stirring mechanisms towards microfluidic mixing
by using computational fluid dynamic (CFD) simulations and experimental validations. As
microfluidics is so diverse that a particular design for mixing may be efficient for some
applications, it may not be efficient for other applications. Therefore, an overall systematic
investigation of geometric configurations is necessary to provide useful guidelines for the
design of microfluidic mixing components.
Many micro-fabrication processes create certain surface roughness or structures in-

side microchannel, these surface roughness or structures may have influences on mixing.
Some of these structures are actually obstacles to the flow. Obstacles can easily be micro-
fabricated, and some similar structures had been studied, therefore, this research chose
obstacles or similar structures in microchannels to study passive mixing. The various com-
binations of obstacle structures influence the flow pattern differently, therefore their effects
on mixing can be evaluated.
This research also includes other complementary objectives. These objectives are
necessary for the principal objective to be carried out.
• At the early stage of this research, the theories on microfluidics and microfluidic mix-
ing needed to be identified and clarified. From a certain point of view, microfluidics
is just another terminology for flow in a micro-domain with slow motion and not
necessarily a new technology. Therefore, existing theories for incompressible slow
5
motion flow may help the research on micro-flow and microfluidic mixing.
• Experimental methods needed to be designed to be able to visualize micro-flow and

microfluidic mixing, for the purpose of validation of the chosen theories.
• The major tool for the research is computational fluid dynamics (CFD) package. In
order to provide accurate simulation results, necessary customized coding needed
to be developed to improve the quality of CFD simulations. These include using
external mesh generator to provide uniform mesh and import to the CFD package for
simulation.
• In order to simulate the 3D convection, a set of particle tracing algorithms needed

to be developed to investigate mixing in spatially configured periodic structures in
microchannels. The algorithms can help to reveal the details of 3D convection in
such channels, and provide useful guidelines for the design of this type of micro-
mixer.
The scope of this research was restricted by several aspects. When width of a mi-
crochannel is of the order of a few tens of microns, diffusion mixing is sufficient without
other means to assist. When Re ¿ 1, the velocity in a microchannel is slow and the resid-
ual time is sufficient for diffusion mixing to take place. On the other hand, when Re À 1,
this creates not only high a pressure drop, but also excessive volume of solution for the
subsequent detection. Nevertheless, one advantage of microfluidic devices is the low sam-
ples and reagents consumption, and using high Reynolds number reverses this advantage.
Therefore, high Reynolds numbers are outside of the scope of this research. The research
focused on mixing of two fluids in a microchannel with dimensions larger than 100µm, and
with the Reynolds numbers were typically less than 10, but sufficient enough to produce a
detectable volume of solution.
6
1.4 Outline of this thesis
To achieve the objectives of this research, the very first step was to understand what
had been researched and why there were still problems. This was addressed in Chapter
2 by a review of literature and related theories to microfluidic mixing. Chapter 2 also
demonstrates the complexity and diversity of microfluidic mixing. One outcome from the
literature review was that the periodic disruption to flow potentially favours mixing, be-
cause most of the micromixer designs applied periodic disruptive forms temporally and/or
spatially. However, the periodic disruption is only one of the necessary conditions that
should be satisfied, and might not be always correct, as later addressed in Chapters 6 and
7.
In chapter 3, the methodology to conduct CFD simulations and construct a simple

experimental rig are described. There are various geometric obstacles which may be incor-
porated in a microchannel, however, it is impractical to micro- machine each of the designs
and experimentally test them. Consequently, CFD was applied to extract useful information
for designing passive microfluidic mixers. To validate CFD simulations, selected designs
were micro-fabricated and tested, and compared to the simulation results. In addition, an-
alytical prediction of viscous flow in a rectangular channel was also used as benchmark of
the simulations.
Chapter 4 identified the general issues surrounding microfluidic mixing. In a study

of all the parameter, flow rate was found to be the most important parameters. Flow was
measured by Reynolds and Péclet numbers, and the Péclet number is defined as the ratio of
bulk mass convection to molecular diffusion. Micro-flow can be recognized as flow at very
low Reynolds number, and as time is implicitly defined in flow rate, the residual diffusion
time is also defined. One of the general parameters for microchannels is the height to width
aspect ratio. For predefined flow rate, very high channel aspect ratios reduce the thickness
of fluid layers for fluid to diffuse through. However, high aspect ratio structures need spe-
7
cial microfabrication techniques, such as deep reactive ion etching (DRIE) or LIGA, which
are expensive and not generally available. Therefore, the research concentrated on low
aspect ratio planar micro-fabricated structures, which were much easier for microfabrica-
tion. Moreover, Taylor dispersion was found to have tremendous influence on microfluidic
mixing, and the understanding of Taylor dispersion also provided information about what
microfluidics could, and could not do in an economical sense.
Chapter 5 introduced two novel designs in splitting and recombining fluid streaming
with low pressure drop. From Chapter 6 to 8, obstacles or similar structures inside mi-
crochannels were explored by simulations and experiments. A very important outcome
from these chapters was that asymmetric structures provided more useful disruption of
flows to improve mixing, and that a symmetric layout of structures had little effect on
mixing performance. Chapter 6 applied cylindrical obstacles to investigate influences on
mixing by optimizing the layout and number of obstacles. An extensive study of obsta-
cles was carried out in Chapter 7. Rectangular obstacles were applied to provide necessary
parameters for the study. These parameters included size, gap, angle, obstacle height to
channel depth aspect ratio and position of obstacles. An interesting result was obtained
when the obstacle height to channel depth aspect ratio became negative, which meant the
obstacles actually became grooves in the channel wall. The grooved surfaces provided lo-
cal slip boundary condition, which induced uneven distributed resistance to the flow. The
uneven distributed resistance created lateral convection. A further study of grooved sur-
faces in microchannels was carried out by applying particle tracing algorithms in Chapter
8.
The thesis was concluded with conclusions and recommendations for future research
in Chapter 9. This chapter highlighted a number of important outcomes emerging from the
present study. The investigation of geometric structures to achieve better disruption flows
to create useful convection for mixing provided useful design guidelines for microfluidic
mixing. The extensive utilization of CFD simulations validated selectively by experiments,
8
or published experimental data, was proven to be a rapid and reliable way to optimize the
performance of microfluidic devices.
9
Chapter 2
Literature Review
2.1 Overview
In this chapter, the fundamentals of microfluidics and theories related to mixing are
first reviewed in sections 2.2 and 2.3 to understand the speciality of mixing in micro-scaled
devices. Then, in section 2.4 fluid drive mechanisms are briefly described. The drive
mechanisms determine the flow rate and pressure drop, and are used to confine the research
scope on microfluidic mixing. The existing microfluidic mixing techniques were catego-
rized into active and passive micro-mixers, which were covered in sections 2.5 and 2.6
respectively. Section 2.7 summarizes the findings of the literature review and identifies the
areas of research which form the basis of the present research.
2.2 Background
In fluid dynamics, in order to study the flow pattern transition from laminar to fully
developed turbulent flow, the Reynolds number must exceed the first critical transitional
10
number (Recr1 ∼ 2300) and even the second critical number (Recr2 > 105 ). When the
Reynolds number is larger than the first transitional Reynolds number, the flow turns to tur-
bulent under the disturbance, however, as soon as the disturbance is removed, the flow turns
back to laminar flow. When the Reynolds number is beyond the second critical transitional
number, any disturbance to the flow causes turbulent flow, and remains turbulent even af-
ter the disturbance was removed (Pan 1988). For the majority of microfluidic devices, the
dimensions were well below 1mm, which made the Reynolds numbers well below the crit-
ical transitional thresholds. Therefore, it was generally accepted that flow in microfluidic
devices was laminar, especially when the working fluid was liquid.
Laminar flow was better understood and easier to manipulate than turbulent flow,
therefore, the laminar nature of microfluidics was an advantage. On the other hand, mi-
crofluidic mixing could not obtain any assistance from turbulence. This meant that any
mixing was to be performed by virtue of diffusion, and diffusion was driven by the concen-
tration gradient between two or more fluids, in accordance with Fick’s law (Bird, Stewart
& Lightfoot 1960). In fact, all the mixing processes for miscible fluids are completed by
diffusion eventually. The engineering challenge is to make the diffusion path small enough
that the process can take place rapidly.
In a macroscopic mixing process, mechanical agitations are normally applied to create

turbulent or chaotic advection. Turbulence stretch and break down large fluid streams into
smaller fluid filaments, so that the diffusion path can be reduced to a microscopic scale.
Eventually, the mixing process is completed by diffusion to create a homogenous solution.
However, it is impractical to bring mechanical agitation into a microchannel for two rea-
sons: first, it is difficult to fabricate; and secondly, it is not an effective approach due to the
low Reynolds number in a microchannel (Purcell 1977).
Diffusion in a microchannel up to a few tens of micron width was found to be effi-

cient for most applications. However, when the size of a microchannel was more than a
hundred microns, pure molecular diffusion could be very slow to achieve complete mixing.
11
The length of channel to achieve a complete mixing by virtue of diffusion would be very
long. As the resistance to the flow was proportional to the length of the channel, this also
meant increasing pressure drop and the micropump might fail to drive the fluids through
the channels. This makes it difficult to design microfluidic chips with such long channels,
therefore, it was necessary in the present research to investigate the techniques to improve
mixing within the confined space of microfluidic devices.
In the following sections, the fundamentals of microfluidics and mixing are discussed,
and the existing active and passive microfluidic mixing methods are reviewed.
2.3 Fundamental Theories for Microfluidics
When a fluidic channel is scaled down to micro-scale, the scaling effects make some
physical parameters, such as gravity and inertia, less important than parameters like surface
tension and viscosity. Conventional fluid mechanics still applies to study microfluidics,
however, some assumptions need to be clarified.
2.3.1 Continuum
In fluid mechanics, a fluid is a substance that deforms continuously under shear stress.
Fluids are composed of molecules in motion. A truly detailed picture of flow will pre-
sumably track each individual molecule. However, this was still impractical for current
computational or experimental resources. Instead of tracking individual molecules, the
fluid could typically be treated as a continuum, which could be infinitely divisible and had
a thermophysical property value at each point in space.
In practice, there was a smallest point that the fluid could be divided into. The volume
of this point should contain a sufficient number of molecules, N , that could give a statisti-
12
cally thermophysical value. According to the random process theory, N was approximately
104 to achieve a statistic point value ρ (Figure 2.1). So the continuum assumption holds
true for average scale length L ∼ 1µm for gases and L ∼ 100nm for liquids. To simplify
Figure 2.1: Statical fluid thermophysical point value

(Santiago 2002)
the judgement of continuum of microfluidics, the Knudsen number could be used. Knud-
sen number (Kn) was defined as the ratio of the mean-free-path of the molecules (λ) to the
characteristic length L of the device. It could be stated as,
λ
Kn = (2.1)
L
The flow regions could be divided into four zones in accordance with the Knudsen number.
They were continuum, slip, transition and free-molecular flow. For liquids with low molec-
ular weight, Knudsen numbers are approximately zero. Continuum was valid for such
liquid flow when the channel was over the micron scale (Gad-el Hak 1999). Therefore, the
continuity equation could be written in Equation 2.2,
∇·v =0 (2.2)
where v is velocity field.

13
2.3.2 Conservation of Momentum
The conservation of linear momentum is actually Newton’s Second Law. For a system
mass, the equation can be stated as,
X d(mv)
F = (2.3)
dt
where F is the force applied on the system, m is the mass, v is the velocity and t is the time.
Applying the conservation principle to a differential element volume of a fluid, a general
conservation equation for momentum can be derived:
Dv ∂v
ρ = ρ( + v · ∇v) (2.4)
Dt ∂t
where ρ is the density of the fluid. The applied forces to the system could include body
forces and surface force. Body forces include gravity, electromagnetic force (Lorentz force)
and others. Pressure is the most common surface force in fluid mechanics. For incompress-
ible fluid motion with uniform viscosity, the conservation of momentum can be written as,
Dv
ρ = f − ∇p + µ∇2 v (2.5)
Dt
where f is the body forces, p is the pressure and µ is the dynamic viscosity.
2.3.3 Brownian Motion
The discovery to Brownian motion was credited to Robert Brown, a botanist who ob-
served the random movement of pollen grains under a microscope. Brownian motion was
thus a continuous stochastic process which was normally distributed, and whose variance
14
increases the farther it gets from the origin. It was independent of its remote history, and a
change in its location at any increment in time was independent of a change anywhere else
over the same time increment. Albert Einstein gave the first mathematical demonstration
of the underlying normal distribution of Brownian motion (Einstein 1905). Though other
mathematicians described many of detailed mathematical properties of Brownian motion,
the present research limited its investigation to Einstein’s approach to understand the diffu-
sion phenomena in microfluidic devices.
Einstein applied molecular-kinetic theory and well-known diffusion theory - Fick’s

laws to prove that the probability distribution of the resulting displacements during an ar-
bitrary time t was thus the same as the distribution of random errors,
x2
n e− 4Dt
Fp (x, t) = √ √ (2.6)
4πD t
where Fp is the probability displacement function in x direction, n is number of samples

and D is the diffusion coefficient. With this equation, Einstein related the diffusion coeffi-
cient to displacement by giving the root-mean-square displacement as,
l2
tD = , (2.7)
2D
where l is the diffusion distance, and tD is the time required to diffuse across distance
l. This equation had been widely used in interpreting microfluidic mixing. Reducing the
length in the equation could reduce the diffusion time by a second order, and the relation
revealed one important method to improve mixing. This was that reducing the diffusion
distance between fluid streams improves mixing.
15
2.3.4 Taylor-Aris dispersion
Taylor (Taylor 1953, Taylor 1954) described the dispersion of solute in a circular pipe
along the longitudes (axial) direction rather than transverse direction. In a circular pipe,
the flow follows Poiseuille’s law - the flow pattern was a parabolic shape for a viscous
fluid. The fluid close to the wall had more time to diffuse than the fluid in the center of the
channel. When the flow was slow enough for a complete radial diffusion, but much faster
than axial diffusion, diffusion in axial direction was negligible. Therefore, the dispersion
of solute in such a flow acted like a slug dispersed axially without the obvious parabolic
flow pattern. The virtual dispersion could be written into mathematical form as,
a2 ū2 4L
D∗ = , À Pe À 7 (2.8)
48D a
where D∗ is the virtual dispersivity or Taylor’s dispersivity, a is the radius of the circular
pipe, u is the mean velocity, D is the molecular diffusion coefficient, L is the longitudinal
length and P e is Péclet number, which is a measurement of convective mass transfer to
molecular diffusion, and P e = 2ūa/D (Catchpole & Fulford 1966). Aris extended Taylor’s
research and gave a more general expression (Aris 1956),
a2 ū2 4L
D∗ = D + , Pe < (2.9)
48D a
Equation 2.9 includes convective, axial and radial diffusion, and was particularly use-
ful in microfluidics, in which the Péclet number was very often larger than 100. And when
P e À L/a, mixing was by pure convective dispersion.
There has also been extensive discussions about the Taylor-Aris dispersion in a rectan-
gular channel (Doshi, Daiya & William 1978, Chatwin & Sullivan 1982, Dutta & Leighton
2001). From these discussions, it was generally accepted that the apparent diffusion coef-
ficient D∗ (Taylor’s dispersivity) in a rectangular channel could be described in equation
16
(2.10),
2 ū2 a2 b
D∗ = D + fg ( ) (2.10)
105 D a
where D is the molecular diffusion coefficient, ū is the mean velocity in the channel, a
is half the width of the channel, b is half of the height of the channel and fg is the geometric
function of a rectangular channel. fg is determined by the depth to width aspect ratio b/a.
The geometric function fg is not unity, which is the case when ignoring the side walls
of the channel. However, the side walls of a microchannel cannot be omitted from the
above equation, therefore, instead of being unity, fg = 7.95 while b/a ¿ 0. For a Taylor
dispersivity D∗ , the length between the dispersed leading and trailing edges of slug, lm ,
was estimated by Equation 2.11, which is rewritten from Equation 2.7.
√
lm = 2tD∗ (2.11)
Taylor dispersivity D∗ was much larger than the molecular diffusivity D, therefore,
for a fixed width lm of a channel, the mixing time t was reduced.
Taylor-Aris dispersion mechanism was found to be an important factor in understand-

ing mixing, and is studied in more detail in the analysis of section 4.3.
2.3.5 Chaotic Advection
In section 2.3.3, Einstein’s interpretation of Brownian motion stated that the goal of
mixing was reduction of length scales to achieve the uniformity of concentration. It was
well known that the combined action of stretching and folding produces exponential area
growth and is a desirable goal in achieving efficient mixing (Figure 2.2). The Taylor-Aris
dispersions can be interpreted as the combined effects of convection and diffusion. The
17
actual dispersion was several orders faster than pure molecular diffusion due to convection.
In this section, chaotic advection is introduced. Technically, how chaotic advection could
be created and the theoretical background are presented.
t0 t1 t2
t3 t4
t5 t6
Figure 2.2: Stirring of a ‘blob ’of marked fluid by 2D, time-dependent (time series from t0
∼ t6), laminar flow
(Ottino 1989)
Chaotic was an apt description of stretching in truly turbulent flows. Chaotic advection
in laminar flow was a term introduced by Aref (Aref 1984). He stated in his paper: “Essen-
tially what was being proposed was the existence of a new advective region, intermediate
between turbulent and laminar advection, which one might call ‘chaotic advection’.”
18
In Aref’s study, the passive advection to emphasize that the particle was so light and
inert that it could do nothing but follow the fluid, and at any instant, its own velocity was
equal to that of the ambient flow, up = uf , where up was the velocity of particle and uf was
the velocity of the fluid. Therefore, the particle advection could be calculated by solving the
Navier-Stokes equations for Newtonian flows. The condition that particle velocity (u, v, w)
equals fluid velocity then leads to the advection equations:
dx
= u(x, y, z, t)
dt
dy
= v(x, y, z, t) (2.12)
dt
dz
= w(x, y, z, t)
dt
From the vantage point of dynamic system theory, Equation 2.12 was more than
enough for producing non-integrable or chaotic dynamics. From Equation 2.12, to produce
chaotic particle motion, 3D steady flows or 2D with time-dependent flows (Figure 2.2))
were the necessary conditions for chaotic advection. Steady 2D advection was integrable,
which means non-chaotic.
For an incompressible 2D flow, Hamiltons canonical equations for a one-degree-of-

freedom system can be stated as (Pan 1988):
dx ∂ϕ
=
dt ∂y
(2.13)
dy ∂ϕ
=−
dt ∂x
where ϕ is the stream function.
The coordinates of the particle, x and y, are the conjugate variables. Use of either
one as generalized coordinate, and the other Cartesian coordinate is then the conjugate
generalized momentum. Phase space in this problem was configuration space. Briefly,
2D kinematics of advection for an incompressible flow is equivalent to the Hamiltonian
dynamics of a one-degree-of-freedom system. This applies regardless of whether the dy-
19
namics of the fluid itself is viscous or inviscid. A time-dependent Hamiltonian system with
one degree of freedom can be chaotic (Aref 1984, Metcalfe & Ottino 1994, Rudman 1998).
A necessary condition for chaos was the “crossing ”of streamlines. In a 2D system,
this could be achieved by time modulation of the flow field, for example by motions of
boundaries or time periodic changes in geometry. In a 3D system, a periodic velocity in
space was necessary to create chaotic advection.
2.4 Generation of Fluid Motion in a Microfluidic Device
Pumping methods have a close relation with microfluidic mixing. For example, con-
tinuous microfluidic mixing could be induced by pulsatile micropumps. In this case, the
mixing was achieved by the combination of micropumps and micro check-valves (Deshmukh
2001). For many electrokinetically (EK) driven flows, the mixing channels could be nar-
rower than those of pressure driven flows so that mixing was sufficient by virtue of molec-
ular diffusion. This implied that it was necessary to study mixing in pressure driven flows
rather than electrokinetically driven flows. On the other hand, the flow rates and pres-
sure heads that could be achieved by current micropumps are two important factors to be
considered in the study of microfluidic mixing.
There were two major methods to drive flow in microfluidic devices: pressure driven
and electrokinetically driven pumping techniques. There were other mechanisms, such
as magnetohydrodynamic (MHD) (Hosokawa, Shimoyama & Miura 1993, Lemoff & Lee
2000, Jang & Lee 2000), which was driven by Lorentz force, electrohydrodynamic (EHD)
pumping (Bart, Tavrow, Mehregany & Lang 1990, Darabi, Rada, Ohadi & Lawler 2002),
which was driven by Coulomb force, and pumping by acoustic streaming, ie. ultrasonic
pumping (Nguyen, Meng, Black & White 2000). MHD, EHD and acoustic streaming
pumping had no moving parts and could deliver continuous flows. These pumping tech-
20
niques had close relation with microfluidic mixing and the mechanisms are reviewed in
section 2.5. A brief review of pressure driven pumping and electrokinetic pumping would
be conducted in this section.
The most critical parameters for micropumps were found to be the flow rates and
pressure heads that micropumps could deliver. Some information about flow rates and
pressure drops of micropumps could be found in early review papers (Shoji & Esashi 1994).
In Figure 2.3, the typical flow rates and pressure heads that could be delivered by current
micropumps are illustrated.
Figure 2.3: Flow rates and pressure heads of micropumps could deliver - a review
(Chen & Santiago 2002)
2.4.1 Mechanically Actuated Micropumps
There have been a few reviews on micropumps (Gravesen et al. 1993, Shoji & Esashi
1994, Nguyen, Huang & Chuan 2002). In this literature review, only a few are covered to
21
demonstrate the working principles of mechanical micropumps.
The fluids were driven by a pressure difference between the inlet and outlet in the
micropumps controlled by mechanical micro check-valves. According to the actuation,
they could be categorized into piezoelectric micropumps, pneumatic micropumps, thermo-
pneumatic micropumps, electrostatic micropumps and others. The maximum output pres-
sure of a micropump depended directly on the available force an actuator could deliver.
Figure 2.4 illustrates a typical micropump. When the membrane was actuated by piezo-
electric force or other actuation forces, the membrane disk was deformed and enlarged the
volume of the cavity to create a vacuum. The inlet check valve was opened passively by
the vacuum and fluid flowed in. When the membrane was released from actuation, the
membrane recovered its shape and reduced the volume of the cavity to create a positive
pressure. The outlet check valve was opened passively and fluid flowed out. Figure 2.5
shows a typical passive check valve used in micropumps.
Figure 2.4: Disk piezoelectric reciprocating micropump

(Shoji & Esashi 1994)
Mechanical micropumps without moving parts featured two complicated series of

loops for the intake and outlet valves designed to favour flow in one direction over the
other. Although this approach was not very efficient in terms of overall flow rate, it had the
advantage of being simple in construction and highly reliable because of no moving parts
in the pump (Figure 2.6). The pump was machined by reactive ion etching (RIE) and pack-
22
Figure 2.5: A typical passive micro check valve

(Shoji & Esashi 1994)
aged with electrical field assisted bonding (anodic bonding) with Pyrex. The driving force
was from a piezoelectric driver glued to the Pyrex. Application of voltage potential across
the piezoelectric driver causes it to change its diameter resulting in bowing of the Pyrex
cover. The maximum flow rate that had been achieved was 1600 µl/minute and pressure
heads of over 7 m water had been achieved.
Figure 2.6: Tesla micropump

(Forster, Bardell, Afromowitz, Sharma & Blanchard 1995)
23
2.4.2 Electrokinetic Micropumps
Electrokinetics was the general term describing phenomena that involved the interac-
tion between solid surfaces, ionic solutions, and applied electric fields. Those phenomena
were collectively defined as electrokinetic phenomena and included two important classes:
electroosmosis and electrophoresis. Electroosmosis described the movement of a liquid
relative to a stationary charged interface under the influence of an electric field. The fixed
surface would typically be a capillary tube or porous plug. Electrophoresis was the induced
movement of a charged interface (usually colloidal particles or macromolecules) relative to
stationary liquids under the influence of an applied electrical field.
It had became clear that electrokinetics was one of the important methods to drive
fluids through microchannels with dimensions in the order of 10 to 100 µm and flow
rates in the nanoliter per second range for rapid mixing and detection (Bousse 1999).
Microfabricated capillary electrophoresis integrated systems had a clear benefit for rapid
separations (Harrison, Fluri, Seiler, Fan, Effenhauser & Manz 1993, Manz, Effenhauser,
Burggraf, Harrison, Seiler & Fluri 1994). Compared to conventional capillary electrophore-
sis, the separation times were reduced from an order of 10 minutes down to an order of
seconds and even less than one millisecond.
Electro-Osmotic (EO) flows in microchannels were critical to the design and process
control of various on-chip microfluidic devices such as modern instruments used in chem-
ical analysis and biomedical diagnostics.
To understand electrokinetically driven flows, it is necessary to describe the electrical

double layer model.
24
Electrical double layer
Most solid surfaces acquire a surface electrical charge when brought into contact with
an electrolyte liquid. When a polar liquid was placed in contact with another phase (like the
wall of the surrounding channel), a potential difference develops at the interface. Molecules
of the polar liquid (dipoles) would re-orient themselves near the phase boundary. Ions
or excess electrons present would cause a charge re-distribution at the interface, causing
formation of the so-called electrical double layer (EDL).
The thickness of electrical double layer was described by several models. An impor-
tant one was the Debye length, which could be stated as,
ψ = ψo exp(−redl z) (2.14)
where ψ was the electrical potential, z is the distance from the wall and redl was identified
as the reciprocal of the thickness of the electrical double layer, also commonly referred to
as the “Debye length”.
Electroosmosis
In the electrical double layer, the ionic concentration near the solid surface was higher
than that in the bulk liquid. In the compact layer, which was about 0.5 nm thick, the ions
were strongly attracted to the wall surface and were immobile. In the diffuse layer the ions
were affected less by the wall and were mobile. The thickness of the diffuse layer ranges
from a few nanometers up to several hundreds of nanometers, depending on the electric
potential of the solid surface, the bulk ionic concentration and other properties of the liquid.
When a liquid was forced through a microchannel under hydrostatic pressure, the ions in
the mobile part of the EDL were carried towards one end. This causes an electrical current,
called streaming current, to flow in the direction of the liquid flow. The accumulation of
ions downstream sets up an electrical field with an electrical potential. This field causes
25
a current, called conduction current, to flow back in the opposite direction. When the
conduction current was equal to the streaming current a steady state was reached. It was
easy to understand that, when the ions were moved in the diffuse double layer, the induced
viscous shear stress would pull the liquid along with them. This was the so-called Electro-
Osmotic (EO) flow (Rice & Whitehead 1965, Burgreen & Nakache 1964, Molho, Herr,
Kenny, Mungal, Deshpande, Gilbert, Garguilo, Paul, St. John, Woudenberg & Connell
1998).
When the Debye length is small, the equation of motion for steady flow with low
Reynolds number, and incompressible flow is:
µ∇2 u = −ρE E + ∇p (2.15)
where µ is the dynamic viscosity, u is velocity, E is electric field, p is pressure and ρE is

the electric charge density.
In the limit of a small Debye length, solving the above the equation yields the Helmholtz-
Smoluchowski equation for the electroosmotic velocity (Equation 2.16).
−εζEO Ex
u= (2.16)
µ
where ε is the permittivity, ζEO is the electroosmotic zeta potential and Ex is the electrical
field.
The Helmholtz-Smoluchowski equation predicts a ”plug-like” velocity profile when

the Debye length was much less than the capillary diameter. Electrophoresis describes the
motion of a charged surface submerged in a fluid under the action of an applied electric
field. The electrophoretic velocity of the dye was again described by the HelmholtzSmolu-
chowski equation (Equation 2.17).
−εζEP Ex
u= (2.17)
µ
26
where ζEP is the electrophoretic zeta potential.
The advantages of using the EO flow as a pumping method to transport liquid in mi-
crochannels in on-chip microfluidic devices include no moving part, no noise and easy to
control the electrical field. A key in the applications of EO pump was to be able to deliver
the precise amount of liquid. Two critical disadvantages for EO flow are the very high
electric voltage required and bubble formations inside microchannels due to electrolysis.
Electroviscous Effect
The EDL field at the solid surface exerts electrical forces on the ions in the liquid,
and hence, restricts the motion of these ions. Consequently, the presence of the EDL field
would reduce the liquid flow in comparison with the cases of no EDL effects. Further-
more, the apparent viscosity could be several times higher than the bulk viscosity of the
liquid (Kulinsky, Wang & Ferrari 1999, Li 2001).
2.5 Active Micro-Mixers
The above sections introduced the flow driven mechanisms and flow behavior which
have a close relation with the mixing mechanisms in a microchannel. In this section, the
utilization of active flow driven mechanisms to create secondary flow to reduce diffusion
path is introduced.
Due to the laminar nature of micro flows, applications of active micro-mixers also
try to create secondary flows, and in some cases, chaotic effects could be obtained. For
most active micro-mixers, this could be achieved by periodic perturbation of the flow fields
temporally, or spatially. This section reviews the magneto-hydrodynamic (MHD) (Bau,
Zhong & Yi 2001, Yi, Qian & Bau 2002) and electro-hydrodynamic (ED) (Choi & Ahn
27
2000) micro-mixers. Travelling waves and acoustic or peristaltic waves, were used to drive
flows and improve mixing in several applications (Yang, Goto, Matsumoto & Maeda 2000,
Selverov & Stone 2001), and could be categorized into active mixers. Mechanical agitation
was rarely used in micro-scaled devices, however, it was still possible to fabricate stirring
rods inside a micro-device (Lu, Ryu & Liu 2002).
In this section, micro-mixers with periodic perturbation of fluids spatially or tempo-

rally are reviewed first in sections 2.5.1 to 2.5.3. Then, mixing mechanisms by the contin-
uously generation of secondary flows are reviewed in sections 2.5.4 and 2.5.5. Finally, a
design of mechanical stirrer is described in section 2.5.6.
2.5.1 Planar Laminar Bubble Mixer
A microfluidic planar bubble mixer was fabricated with a five mask process (Fig-
ure 2.7). It was composed of bubble pumps and bubble check valves. Thermo-capillary
bubble valves operate by creating a vapor bubble within chamber, and then, using these
bubbles prevent flow though the chamber. Bubble pumps operate by turning on and off
micro-heaters to create or eliminate bubbles within the chamber, thus the volume thermo-
expansion or shrinkage created piston-type pumping. The pumping process created chaotic
advection to mix fluid streams in the planar, laminar chamber. A chaotic flow field was one
in which the path and final position of a particles placed within the field were extremely
sensitive to their initial position. In the chaotic flow field, particles initially close together
might become widely separated, and the flow as a whole becomes well mixed.
2.5.2 Periodic Perturbation Applied to a Main Fluid Stream
Periodic pressure or electrokinetic perturbation was sometimes declared as chaotic ad-

vection (Figure 2.8). Without further investigation, it was not clear that it was in the chaotic
28
Figure 2.7: 2D chaotic advection microfluidic mixer

(Evans, Liepmann & Pisano 1997)
region. Similar designs had been studied with a different interpretation (Deshmukh 2001).
The idea was to introduce slugs of solutes into the main stream. In accordance with the
concept of Taylor-Aris dispersion in a microchannel, the actual dispersion was much faster
than pure molecular diffusion (section 2.3.4). Under the term of continuous microfluidic
mixing, pulsatile micropumps were successfully microfabricated. The periodic perturba-
tion by the pulsatile pumping produced alternated fluids slugs that created flow patterns
similar to Taylor-Aris dispersion. It was a versatile method to solve microfluidic mixing
problem. However, it should be noted that the claim about high flow rate mixing was not
clearly verified to fulfill the inequality of Taylor-Aris dispersion. For flow rates that were
too high, there was not sufficient residual time for the fluids to diffuse into each other.
High flow rates generate high pressure drops, which was a disadvantage for microfluidics
and should be avoided. Though, this kind of mixer contained moving parts, it was com-
plicated for fabrication and difficult for flow control. Nevertheless, mixing by applying
Taylor-Aris dispersion was found important for further study and is addressed in Chapter
4.
29
Figure 2.8: Microfluidic mixer by periodical perturbation

(Lee, Deval, Tabeling & Ho 2001)
2.5.3 Peristaltically Driven Micromixing
In the previous section (2.5.2), the periodic perturbations to the main flows were un-
der low frequencies. When the frequencies were high and the magnitudes of the waves
were low, this formed another mixing category - peristaltically travelling waves. Most of
the peristaltically type driven micro-mixing was induced by acoustic waves, and some by
mechanically driven peristaltic waves.
Although these mixing methods showed promising results, the mixing mechanisms
were described differently. Turbulence could possibly exist in such a system, as declared
in a previous paper (Yang et al. 2000). However, most researchers agreed that the fluid
motion induced by the acoustic waves caused certain recirculations that favour mixing (Zhu
& Kim 1998, Monnier, Wilhelm & Delmas 1999, Rife, Bell, Horwitz, Kabler, Auyeung &
Kim 2000, Vivek, Zeng & Kim 2000, Selverov & Stone 2001, Liu, Yang, Pindera, Athavale
& Grodzinski 2002, Liu, Lenigk & Grodzinski 2003).
In acoustic microfluidic mixers, piezoelectric acoustic transducers were attached to

the fluidic channels. The transducers converted radio-frequency (RF) electrical energy into
30
an ultrasonic acoustic waves. The acoustic waves were absorbed in the fluids, but generated
mechanical travelling waves ( flexible plate waves or surface acoustic waves (SAW)) at the
attached solid structures (Figure 2.9). The mechanical travelling waves induce fluid mo-
tion. By controlling the direction of the streaming, pumping or mixing could be achieved.
For microfluidic mixing, mixing time across channels with characteristic length l would
be O(l2 /D), where D = kB T /6πµr was the Stokes-Einstein translational diffusion coeffi-
cient, kB was the Boltzmann constant, T was the absolute temperature and r was the radius
of the diffusing particle. Typical times for Peristaltic convective transport were O(1/²2 ω),
where ² is the dimensionless amplitude of the travelling wave and ω is the frequency of the
oscillation of the wall (see Figure 2.9). This implies that frequencies should be O(D/²2 l2 )
or larger. For typical microfluidic mixers driven by piezoelectric waves, convective trans-
port in aqueous solutions requires frequencies of ω À 103 Hz.
Figure 2.9: Mechanical travelling waves induced by acoustic streaming

(Nguyen et al. 2002)
2.5.4 Magnetohydrodynamic (MHD) Micro-mixers
The potential to apply electromagnetic force for fluids pumping had attracted atten-
tion (Hosokawa et al. 1993, Heng, Huang, Wang & Murphy 1999, Lemoff & Lee 2000, Jang
31
& Lee 2000), as mentioned in section 2.4. MHD could be used not only for the purpose
of pumping fluids, but also to control fluid flow in microchannels and to induce secondary
flows that maybe beneficial for mixing. By depositing an array of transverse electrodes on
the bottom of a microfluidic conduit (Bau et al. 2001), the electrodes were connected alter-
nately to two terminals of a DC power supply. As a result, electric currents were induced
in opposite directions between adjacent pairs of electrodes. The interaction between these
currents and uniform magnetic field perpendicular to the bottom of the channel led to the
formation of Lorentz forces in opposing directions between adjacent pairs of electrodes.
This, in turn, led to the formation of convection. This convection could be used to fold and
stretch material lines and enhance mixing (Figure 2.10). Yi et al. extended the study exper-
imentally and theoretically, and demonstrated chaotic stirring effects by MHD in a cylinder
cavity (Yi et al. 2002). However, steam bubbles induced in the MHD micro-mixers was a
major problem.
Figure 2.10: A minute MHD micro-mixer

(Yi et al. 2002)
32
2.5.5 Electrohydrodynamic (EHD) Micro-mixers
In electrically dominated flows, there is always some residual free charge present,
typically residing at the interface between two fluids with different electric conductivities
flowing side by side in a microchannel. The external electrical fields that are tangential to
the interface result in a tangential stress. Dragged by this shear stress, fluid motion can be
formed with the moving charges. This model which deals with fluid motion induced by
electric fields was firstly introduced by Taylor, and Melcher used it extensively to develop
electrohydrodynamics (Brenner & Stone 2000).
In the case of microfluidic mixing, the flow motion induced by an electrical field
transverse to the main stream, forms the secondary flow (Figure 2.11). This secondary
flow stretches and folds material lines, and therefore, enhances microfluidic mixing. The
operation voltage was of a range from 7 to 27 volts.
In Figure 2.11, the electrical field was perpendicular to the surface charges, so it could
form secondary flow. While the external electric field was parallel to the surface charges,
transverse secondary flow was not induced. However, when the working voltage was larger
than 1kV, it creates electrokinetic instability. Although the mechanism for the cause of the
electrokinetic instability was still not clear, its appearance was very much like turbulent
flow, even though the Reynolds number was low. By applying this instability, a micro-
mixer was fabricated in Stanford University (Figure 2.12).
2.5.6 Magnetic Micro-Stirrer
Inspired by the macro-scaled magnetic bar mixers, a micro-scaled mechanical stirring

mixer was fabricated by surface micro-machining (Lu et al. 2002). In the process, a mag-
netic actuated rotary permalloy bar, or an array of these bars, on Pyrex 7740 glass wafer
was encapsulated in the microchannel formed on polydimethylsiloxane (PDMS) substrate
33
b
Figure 2.11: (a). Schematic illustration of the active micro-mixer (b). Cross sectional view
of the active micro mixer
(Choi & Ahn 2000)
34
0 sec 0.3 sec 1 sec 2.3 sec
Figure 2.12: Mixing processes in the micro-mixer by applying electrokinetic instability

(Oddy, Santiago & Mikkelsen 2001)
(Figure 2.13). When the stir bar was rotated, it stirred the fluids and increased the contact
area between the two fluids, hence increased the rate of diffusive mixing between them.
The mixing efficiency would be further improved, if the rotating speed was increased (Fig-
ure 2.14). However, viscosity in a microfluidic device plays a dominant role (Purcell 1977),
and the magnetic force on the stirrer had to overcome the viscosity of the fluids. Hence, it
is difficult to achieve very high rotating speeds. Nevertheless, it was predicted by Purcell
that the mechanical agitation for mixing was not effective at very low Reynolds numbers.
However, this micro-stirrer mixer might have an unique advantage on mixing fluids in a
chamber.
Figure 2.13: Magnetic micro stirrer

(Lu et al. 2002)
35
Figure 2.14: Mixing in microchannel by magnetic micro stirrer, (a). 0 rpm, (b). 150 rpm,
and (c). 300 rpm
(Lu et al. 2002)
36
2.6 Passive Microfluidic Mixers
By definition, passive micro-mixers have no moving parts and other external mecha-
nisms to create secondary flows. Therefore, they are relative easy for fabrication. Because
no electrical and magnetic parts are involved, there is no special requirements for electric
and/or magnetic material properties, and passive micro-mixers have much broader materi-
als to choose from. Less material requirements also reduce the complexity of microfabri-
cation.
A passive micro-mixer can improve mixing by reducing the diffusion path sufficiently,
and by continuously moving fluids away from the interfacial areas, which is actually bulk
mass convection. In a microscopic scale, convection is always coupled with diffusion, due
to the low Reynolds number of microfluidics.
In this section, diffusion dominated passive mixing is first reviewed. This includes T-
type mixers, splitting and recombining mixing. Then, passive mixers coupled with convec-
tion and diffusion type mixers are reviewed. This category includes focusing fluid streams,
sub-injection, collision and passive chaotic mixers.
2.6.1 T-type mixers
In many microfluidic devices, to introduce reagents and samples into the devices, the
T-type shape of channels were often used. The “Y ”shape, cross, and other similar shapes
follow the same design principles and could be categorized into T-type channels. When
the size of microchannel was a few tens of microns, the contacting of two substreams in a
T-type channel was sufficient enough to achieve adequate mixing. Many of electroosmotic
driven microfluidic devices use channels that had the sizes less than 100µm, and could
contact reagent and sample directly in T-type mixer without any other mixing assistance.
37
Mixing Tee in Electrokinetically Driven Microchips
In Figure 2.15, it was shown that in a compact chip architecture, and with a single
voltage source, parallel mixing with different mixing ratios of sample to buffer could be
accomplished by using a series of T-mixers (Figure 2.15 a), and the serial mixing device
was based on an array of cross-intersection and sample shunting (Figure 2.15 b). The
channels had similar cross-sectional areas, which made channel resistances proportional
to length. Applying the same potential to both the sample and buffer reservoirs allowed
accurate and fast mixing.
a b
Figure 2.15: Mixing T structures, (a). Parallel configuration; (b). Serial configuration
(Jacobson, McKnight & Ramsey 1999)
38
Collision of Fluid Streams in Mixing-Tees
When collision of two fluid streams in a T-junction occurs, a flow pattern termed
engulfment flow could be created (Figure 2.16). The swirling flow increased the contact
area between fluid streams and so did the mixing performance.
Figure 2.16: Streamlines of engulfment flow

(Kockmann, Foll & Woias 2003)
Similar in principle, a double mixing tee micro device was applied to improve mixing
(Figure 2.17). The outlet of the first tee formed one inlet of the second mixing tee. This
geometry provided a so-called residence time channel of defined length separated by two
mixing tees. In the idealized case of short mixing length, this design would enable imme-
diate start of a reaction and control of the reaction time by variation of flow rate or length
of flow passage before stopping the reaction by a second mixing process. This structure
needed a high volumetric flow rate to ensure the good mixing by collision of two streams.
T-sensors
A T-type of micro-mixer, sometime termed T-sensor, is one of the micro Total Anal-
ysis System (µTAS) components. T-sensors had been investigated extensively by many
researchers (Weigl & Yager 1999, Kamholz, Weigl, Finlayson & Yager 1999, Kamholz &
Yager 2000, Hatch, Kamholz, Hawkins, Munson, Schilling, Weigl & Yager 2001, Kamholz,
Schilling & Yager 2001, Kamholz & Yager 2001, Kamholz & Yager 2002, Ismagilov, Ros-
marin, Kenis, Chiu, Zhang, Stone & Whitesides 2001, Ismagilov, Stroock, Kenis, White-
sides & Stone 2000). Two or more fluid layers enter the channel side by side, and mixing
39
Figure 2.17: Double-T micro-reactor

(Bokenkamp, Desai, Yang, Tai, Marzluff & Mayo 1998)
takes place between the fluids streams only by diffusion. The mixing zone could be quan-
titatively measured by optical methods (Figure 2.18). For a pressure-driven viscous flow in
a microchannel, the experimental results showed that the flow patterns were parabolic-like
profiles in one direction or even both directions that were perpendicular to the flow axis
and the velocity was zero at the wall and maximum velocity in the center of the channel.
Due to the parabolic flow profiles, the fluids that were close to the wall moved slower,
hence, had more residual time for diffusion and were mixed first. To reduce the incon-
venience of the Taylor dispersion effect for sensing, the depth to width aspect ratio of
T-mixer was usually kept small for the solute to diffuse rapidly across the depth direction
to achieve uniformity before the detection zone. Scaling rules were also studied (Ismagilov
et al. 2000, Kamholz & Yager 2002). Near the wall of microchannel, the width of region
mixed by diffusion scales as the one-third power of the ratio of the axial distance, whereas
in the center of the channel, the diffusion scales by the one-half law. For a small depth to
width aspect ratio, reasonably accurate estimations using Einstein equation (Equation 2.7)
could be made (Kamholz & Yager 2002).
40
Figure 2.18: A schematic working principle of T-Sensor

(Kamholz & Yager 2001)
41
The study of research on T-mixers provided a general understanding of how diffusion

acts in microfluidic devices and many novel applications were found, such as lab-on-a-chip
and medical diagnostics.
2.6.2 Splitting and Recombining
T-type micro-mixers provided a way to introduce two fluid streams into microchan-
nels, and the mixing was by virtue of diffusion. Mixing by diffusion requires a long mi-
crochannel to achieve complete mixing. To reduce the channel length and the space oc-
cupied by the mixing units, the diffusion path between fluid streams must be reduced to a
sufficiently small scale.
It was intuitive to split fluid streams into tens or hundreds of small substreams to
increase diffusion. In this way, large contact surfaces and small diffusion paths were gener-
ated (Koch, Chatelain, Evans & Brunnschweiler 1998, Koch, Witt, Evans & Brunnschweiler
1999, Koch, Schabmueller, Evans & Brunnschweiler 1999, Ehrfeld, Hessel & Lowe 2000)
(Figure 2.19)). These methods were simple, but also require a large amount of space,
which makes the device larger. Another concern was the uneven pressure distribution of
such mixers, which might cause partial mixing at some regions of the mixing unit.
In a microfluidic analytical system with flow rates of the order of nanoliters and/or
even picoliters, the splitting and recombining of fluid streams could be fabricated in-situ.
He et al. described the in-situ micro static mixer with a flow rate as low as 100pL, and
liquid A and liquid B enter a Y-channel along one side (He, Burke, Zhang, Zhang &
Regnier 2001, He, Tait & Regnier 1998). In this design, the flow was redirected later-
ally with shrinkage of channel width to introduce convective effect, and consequently, the
diffusion distance was reduced. The fluid streams were further divided partially by side
bypass channels and recombined at the other end (Figure 2.20). For the electrokinetically
driven microfluidic device, the channel size was small, and the mixing efficiency was high.
42
Figure 2.19: Parallel dividing into substreams.

(Koch et al. 1999)
However, it may be difficult to apply this design in a pressure driven micro device because
of possible high pressure drop. Therefore, to determine the suitability of the micro-mixer
design for pressure driven flow, it was decided that an investigation be carried out and the
results reported in this thesis.
By 3D geometrically design of the micro-mixer (Schwesinger, Frank & Wurmus 1996,

Hessel, Hardt, Lowe & Schonfeld 2003, Nielsen, Branebjerg, Gravesen, Dyhr-Mikkelsen
& Gade 2001), the diffusion process could be arranged to form a multi-layer sandwich-like
structure to reduce the diffusion path (in Figure 2.21). The fluids were split, stretched into
substreams and recombined repeatedly by the fork-shaped divided channels. The lamella
thickness for the diffusion could be reduced to a sufficiently small scale. Moreover, the
structures could also be used for immiscible liquids to produce emulsified mixtures.
Caterpillar micro-mixers were developed by IMM-Mainz. The splitting and recom-

bining was repeated periodically by using 3D structures (Figure 2.22 a). The complex
geometric structures divided the streamlines repeatedly (Figures 2.22 b and c). However,
the complexity of the structures made it very difficult to fabricate and optimise.
The Danfoss design of a micro-mixer used very much a similar mixing method to the
Caterpillar micro-mixers (Nielsen et al. 2001), however, one stream passed underneath an-
43
Figure 2.20: In-situ static micro-mixer by subdividing streams

(He et al. 2001, He et al. 1998)
Figure 2.21: Divided flow into substreams (a) vertical arrangement; (b) horizontal
arrangement.
(Schwesinger et al. 1996)
44
a b
c
Figure 2.22: (a). Channel structures; (b). Schematic drawing of its split-recombine-
principle; (c). Streamlines in a caterpillar mixer.
(Hessel et al. 2003)
45
other stream and joined in one mixing chamber (Figure 2.23). The recombining of streams
was in a vertical manner. The diffusion distance was reduced dramatically, due to the shal-
low mixing chamber and the broadening of channel width, which also reduced the thickness
of the fluid streams.
Figure 2.23: Danfoss mixer

(Nielsen et al. 2001)
Figure 2.24 illustrates the vertical splitting and recombining of fluids in a static micro-
mixer (Klaus Schubert, Wilhelm Bier, Erhard Herrmann, Thomas Menzel & Gerd Linder
2000). The guiding structures divided fluid into multiple smaller filaments that overlap
vertically. The patent claims to achieve uniformed mixture over the whole outlet area of
this mixer.
The above review showed that the diffusion distance could be reduced by fabricating
mixing unit in 3D structures, and the use of device space was decreased. The disadvantage
of this mixer was its complexity for fabrication and the high pressure drops that prohib-
ited their practical applications. For the structures that were made on two separate Silicon
wafers, the packaging process also required high accuracy to align two or more parts to-
gether. To reduce pressure drop in such devices, a modified structure was proposed (Mahe,
Tranchant, Tromeur & Schwesinger 2003). The modified design demonstrated comparable
mixing efficiency and much lower pressure drop.
46
Figure 2.24: Static Micro-mixer

(Klaus Schubert et al. 2000)
2.6.3 Focusing Fluid Streams
Reducing the diffusion distance between fluid streams could be also achieved by so
called “focusing ”. Focusing the fluid streams into a very thin layer could be done by
geometrically shrinking the width of the channel (Veenstra, Lammerink, Elwenspoek &
van den Berg 1999), or by hydrodynamic focusing (Knight, Vishwanath, Brody & Austin
1998). A super-focusing micro-mixer had been developed by combining the geometrically
splitting and recombining technique together with the focusing technique. First the main
fluid streams were divided into many substreams with reduced channel width, then the
substreams were converged into a channel with much smaller width than their total width
of the sub-channels (Figure 2.25). Mixing in such a micro-mixer could be completed within
only a few milliseconds.
In a hydrodynamic focusing system for ultra-fast micromixing (Knight et al. 1998),

47
Figure 2.25: IMM super-focusing micro-mixer

(Hessel et al. 2003)
one stream was squeezed into a main channel by two side flows hydrodynamically (Figure
2.26). For incompressible flow, to keep the conservation of mass (continuity), the flow
rates must be constant. Therefore, by controlling the velocity of the middle fluid stream, in
accordance with v = (flow rate) / (Cross-section area), the cross-section area of the middle
stream was stretched into a thin layer to keep the flow rate constant. While the thickness
of this thin filament was reduced to a few microns, mixing could be achieved rapidly by
diffusing through this thin layer.
2.6.4 Inject Small Fluid Substreams into a Main Stream
In Figure 2.27, the two fluids flowed towards each other and subdivided by the sinu-
soidal shaped slits. The substreams were recombined in a mixing chamber above these
slits, and reduced diffusion paths could be achieved. The main concern of this structure
was the control of the uniformity of the flow. If the mixing zone was significantly larger
than the microchannel width, a non-uniform flow would occur, resulting only in a partial
overlap of the sheets. The injection mixing schemes had another arrangement illustrated in
48
Figure 2.26: Hydrodynamic focusing micro-mixer

(Knight et al. 1998)
Figure 2.28. However, instead of combining fluid streams above the channel, the sub-jets
from the bottom of the structure were guided so they could be combined with the main fluid
stream and flowed side-by-side into the mixing chamber.
Figure 2.29 illustrates injection of sub-jets into a mixing chamber through a micro-
fabricated square-shaped array of nozzles (Elwenspoek, Lammerink, Miyake & Fluitman
1994). The injection of multiple sub-jets formed micro-plumes, which increased the con-
tact surfaces between two fluids.
2.6.5 Passive Microfluidic Chaotic Mixers
If one defines the axial flow as the primary flow, any modification of the speed or
direction of primary flow was termed secondary flow. In some cases, strong secondary
flows can create chaotic advection (covered in section 2.3.5).
49
Figure 2.27: IMM-Mainz Mixing unit

(Ehrfeld, Golbig, Hessel, Lowe & Richter 1999)
Figure 2.28: Injection fluid into a mixing chamber

(Larsen 2001)
50
a b
Figure 2.29: (a). Mixing Chamber with Micro-plumes (b). Simulation of sub-jets of
reagent into main stream
(Elwenspoek et al. 1994)
Chaotic advection could stretch and fold the fluids dramatically. As a result, the con-
tacting interfacial area was increased and the thickness of the fluids layer was decreased
exponentially (Ottino 1989). In other words, the diffusion distance was reduced exponen-
tially. As the process continues, rapid mixing could be achieved.
One important strategy to solve microfluidic mixing problem was to create chaotic
advection. Chaotic advection, sometime called Lagrange chaos, had been studied inten-
sively in the macro-world (Aref 1984, Jones, Thomas & Aref 1989, Ottino 1989). The
knowledge could be borrowed for the design of microfluidic mixers. For instance, using a
twisted pipe, secondary flow occurs at the bends of the pipe (Jones et al. 1989). Chaotic
advection could be a result of these secondary flows. Inspired by the twisted pipe, 3D ser-
pentine channels were fabricated (Beebe, Adrian, Olsen, Stremler, Aref & Jo 2001, Liu,
Stremler, Sharp, Olsen, Santiago, Adrian, Aref & Beebe 2000) (Figure 2.30). The flow in
3D serpentine channels demonstrated chaotic advection at high flow rates (Re ∼ 70) and
good mixing performance. One of the advantages of the serpentine design was to keep the
51
height-to-width ratio near one so that it minimizes the chances of clogging, fouling and
loss of sample in the biological analytical applications. However, the challenges of ser-
pentine mixers were the microfabrication of complex 3D structures and the need for a high
Reynolds number to stir the fluids to generate chaotic advection. Another drawback of this
design is the dead volume created in the perpendicular bends, as the fluid at the corner of
these bends keeps still, the only mass transport is diffusion.
Figure 2.30: Serpentine-Shaped Micro-mixers

(Liu et al. 2000)
Chaotic microfluidic mixers made by patterning grooves in microchannels were com-

patible with microfabrication processes and could be realized in general micromachining
52
centers (Johnson, Ross & Locascio 2002, Stroock, Dertinger, Ajdari, Mezic, Stone &
Whitesides 2002, Stroock, Dertinger, Whitesides & Ajdari 2002). The anisotropic pat-
terned grooves in the micro-channel could create a pressure component laterally to the
primary flow, which induced spiral circulation around the flow axis, and could stretch and
fold streams of liquids to complete mixing in a shorter channel length (Figure 2.31). The
mixing effects were investigated over a range of Reynolds numbers. Both the serpentine
and patterned grooves chaotic mixing strategies showed that the design for chaos was one
of the effective solutions to enhance mixing in a micro-scaled device.
Figure 2.31: Microfluidic chaotic mixer with grooved surfaces

(Stroock et al. 2002)
Further study revealed that deeper grooves create stronger chaotic effects (Johnson
& Locascio 2002), and that the flow patterns for this type of micro-mixer were indepen-
dent of Reynolds number (Stroock et al. 2002). Similar to the slant grooved surface de-
sign, transverse recirculation could also be created by the so called transverse electrokinet-
ics (Ajdari 2002)
The obvious drawback of patterning grooves in microchannels for mixing was the
dead volume created by the grooves. While reducing dead volume is desirable, the deep
grooves were critical for generating helicity and research to date has not yet found any
53
direct solution to reducing the dead volume in such mixers. From an engineering design
point of view, the depth of the grooves is one of the design parameters of this type of mi-
cromixer and needs to be further investigated. This work is carried out in this research.
Nevertheless, the design was compatible with microfabrication processes, and works effi-
ciently, especially at low Reynolds number. Hence, this mixing strategy could be used in
many applications. For instance, disposable microfluidic devices for medical diagnosis and
other applications that did not need the grooves to be cleaned after using them. Under these
conditions, a lower flow rate normally means lower pressure drop, which makes it easier to
drive the flow through a microfluidic device.
2.7 Summary
In a typical microfluidic device, viscosity dominates flow, and as a result, the Reynolds
number is low and the flow is laminar. Therefore, the mixing of two or more fluid streams
in microfluidic devices is by virtue of diffusion, which is a slow process. On the other hand,
for microfluidic applications towards biomedical or chemical diagnostics, such as lab-on-a-
chip, it was very important to mix two or more reagents and test samples thoroughly before
detection can take place. This requires, a) the mixing needs to be done in a confined area;
b) The volume of fluids for detection should be adequate and sufficiently mixed.
Therefore, microfluidic mixing is challenging and device design issues need to be

solved, and many mixing methods have been developed. Based on the mixing mechanism,
micro-mixers have pure molecular diffusion at one end and chaotic mixing on the other.
Based on structures, micro-mixers can be categorized in either the active or the passive
mixing methods. Most of active micro-mixers enhance mixing by stirring flow. The stirring
can be done mechanically or by the means of magneto-hydrodynamic (MHD), electro-
hydrodynamic (EHD) or by acoustic streaming to create secondary flows. The secondary
flows can stretch and fold material lines to reduce the diffusion path between fluid streams,
54
and hence, enhance mixing. Most active micro-mixers are complex to fabricate and require
an external power source.
Some passive mixers reduce diffusion path between fluid streams by splitting and re-
combining. Recently, passive chaotic micro-mixers were investigated by using 3D serpentine-
type channels. The flow in 3D serpentine channels demonstrated chaotic advection at high
flow rates (Re ∼ 70) and good mixing performance. However, the challenges of serpentine
mixers are the microfabrication of complex 3D structures and the need of high Reynolds
number to stir the fluids to generate chaotic advection.
Passive micro-mixers with planar structures are relatively easy to fabricate and can be
made of most engineering materials. The mixers were found to be efficient over a range
of Reynolds numbers and compatible with microfabrication processes. However, most of
the research to date focused on one or a few passive structures, and a full understanding
of geometric structures, and their influence on mixing and flow patterns had not yet been
achieved.
In this research, in-situ passive microfluidic mixers with various planar or semi-planar
geometric structures were systematically studied to understand the influence of these struc-
tures on mixing performance and pressure drop.
It would be impractical to apply solely experimental method to undertake the system-

atic investigation of a series geometric variations. Furthermore, the nonlinear behaviour of
some optics used in visualisation of flow in the existing experimental methods also added
uncertainty to the accuracy of measurement. Therefore, to achieve the objective of this
research, advanced computational fluid dynamics (CFD) techniques were chosen as the
major tools for the investigation. In addition, a simplified experimental method was also
setup to validate selected numerical models used in the CFD simulations. By the experi-
mental validation, the numerical models could be verified, and confidence to use CFD for
further investigation could be achieved. Where applicable, published experimental and an-
55
alytical results were used for the validation of simulation results. The methodology for the
investigation is described in Chapter 3.
For the mixing mechanisms reviewed in this chapter, there was not a clear interpreta-
tion of mixing in relation to geometric variations. To initiate the present research, a brief
description of mixing in a T-channel is introduced in Chapter 4. The study of two-adjacent
fluid streams diffusing into each other was investigated previously. In Chapter 4, further
discussion of Taylor-Aris dispersion is presented. Because pressure drop was also a criti-
cal parameter to the performance of a micro-mixer, the relation between pressure drop and
channel width to depth ratio was also investigated for the better understanding of microflu-
idic mixing.
In regards to periodic geometric structures, a serpentine shaped channel was investi-

gated previously. Serpentine shaped channel structure increases the total length of channel
on an area constrained microfluidic chip, and hence, increases residual time for diffusion.
However, for further improvement of mixing efficiency, it requires high Reynolds num-
ber and 3D serpentine shaped micro-mixers are complicated for fabrication. Chapter 5
presents three methods to modify serpentine and similar structures to improve the mixing
performance of this type of mixers. The first method was to partition flow in a single-side
serpentine to create mixing zones with different concentrations. The second approach was
to short-circuit the serpentine by a straight branch to generate partial splitting and recom-
bining at the interface between two fluid streams. The last approach in this chapter was
inspired by the structures machined by Excimer laser, which had periodic ramping struc-
tures at the bottom of the channel. This structure was similar to 3D serpentine structures.
While the periodic ramping structures inspired further research, it was difficult to de-
fine the geometry of a 3D ramping in a numerical model. Therefore, the ramping structures
are simplified to cylindrical obstacles placed inside the microchannel. Obstacles could
be positioned in the microchannel arbitrarily so that mixing performance versus different
layout of obstacles could be studied. This work is reported in Chapter 6.
56
Although cylindrical obstacles could provide the information about how the positions
act on mixing, other aspects of obstacles are still of the interest to this research. Therefore,
rectangular shaped obstacles were added and the results presented in Chapter 7. Rect-
angular shaped obstacles provided other parameters for the study, such as size, degree of
blockage, offset, angle and height of the obstacles versus the mixer performance. When the
height of the obstacles became negative, the obstacles decayed into grooves at the bottom
of the channel, which resulted in the further study of grooved surfaces in Chapter 8.
Past research on mixers with grooved channels covered experimental and analytical
prediction. Due to the relevance to the study of obstacle structures, numerical investigation
of micro-mixers with slantly patterned grooves was added in this research to determine
guidelines for the engineering design of this type of micro-mixer.
The research presented in this thesis links the existing knowledge on passive mixing
mechanisms with geometric configurations with simplified obstacle structures. The results
of the present research provide useful information for the design of passive micro-mixers
and identify further areas of research in this field. These are summarized in Chapter 9.
57
Chapter 3
Numerical Modelling and Experimental

Methodology
3.1 Introduction
In chapter 2, the literature review found that it was important and necessary to study
passive micromixers with in-situ geometric variations to enhance mixing. These geomet-
ric variations included channel shape, structures inside channel and geometric topology of
channel walls. For the large number of parameters involved, it was difficult to conduct
the research by experiments. Though some researchers delicately utilized micro particle
imaging velocimetry (µPIV) to measure the whole flow field, these techniques are limited
to particular applications, and typically, to slow shallow flow field due to the depth of field
of optical lens (Santiago, Wereley, Meinhart, Beebe & Adrian 1998). Moreover, Computa-
tional Fluid Dynamics (CFD) can neglect un-important factors in the real environment, and
focuses on the fundamental principles of the problem, therefore, CFD can be a versatile
tool for studying laminar incompressible flows.
58
The existing CFD softwares can be divided into three categories. One is the general
CFD packages for macroscopic scale analysis, such as Fluent, CFX, which are based on
the Navier-Stokes equations. One way to apply this kind of CFD softwares to microfluidics
is to understand the principles and the structures of the packages, and then modify some
parameters to model a microfluidic device by utilizing the user program interface embedded
in the software. This category of CFD package, in general, is well tested, and usually
implemented with a versatile CAD system for pre-processing and visualization tools for
post-processing. Therefore, the users can treat the packages as black boxes and concentrate
on the analysis. However, if the structure of the package needs to be changed for simulating
the special flow behavior in a microscopic scale, they may face problems. Because the users
normally cannot access the source code of the packages for modification, and even if they
can, it is a time consuming and a difficult task for the individual researcher.
The second group of CFD package is specially designed for microfluidics, such as
MemCFD from Coventor Inc. and CFDACE+ from CFDRC Inc. These softwares enable
users to do some quick simulations and designs, and hence, save time. However, com-
pared to the general CFD package, they have a much smaller user group, and a shorter
history. Furthermore, the theory of microfluidics was still under development. Therefore,
the users should be aware of the restrictions of these softwares. In addition, because the
surface/volume ratio is high in microfluidics devices, a CFD package should be able to
simulate the interfacial effects, such as capillary flow and electrokinetic flow.
The third group, eg. DiffPack and Matlab, contains some general tools, for exam-
ple, partial differential equation (PDE) solvers, which can complement CFD packages for
developing microfluidic models. This group has the greatest flexibility. Nevertheless, the
development and validation of numerical coding were expensive, time consuming, and
the performance of the simulation depends very much on the experience of individual re-
searchers.
TM
For the above reasons, a microfluidic CFD package, MemCFD , was chosen as
59
the primary tool to study the microfluidic mixing. MemCFD embeds Fluent and Fidap
solvers to solve Navier-Stokes equations and has functions to simulate multi-species mix-
ing, pressure-driven flows as well as electro-osmotic flows. Fluent and Fidap are the state-
of-art CFD solvers to date, and MemCFD also has integrated advanced modelling and
visualization tools, therefore, MemCFD can provide reliable and high quality simulations.
To complement and validate the CFD simulation results and the algorithms developed
as part of this study, experiments were conducted to validate typical models as well as the
use of published experimental data and analytical predictions. In the research, customized
codes, to improve simulation quality, and particle tracing algorithms were also developed
for an extensive investigation of 3D complex advection in micro-mixers with patterned
grooves (section 8.3).
In this chapter, section 3.2 introduces the CFD fundamentals. Section 3.3 describes
the special requirements of microfluidic modelling. Complementary experiments were nec-
essary to understand the physical principles behind the numerical simulations, therefore, a
simple experimental method is introduced in Section 3.4.
3.2 Fundamentals of Computational Fluid Dynamics (CFD)
3.2.1 Overview
The technological value of computational fluid dynamics has become undisputed. In

microfluidics, it might be very difficult to experimentally diagnose the flow behavior at
micron or sub-micron scales. However, it is still feasible for CFD to be used to simulate
micro-flow as long as the validation of continuity is maintained (in Section 2.3.1). For the
understanding of the significance and limitation of CFD in microfluidics, a basic knowledge
is introduced in this section.
60
CFD simulation can be divided into three stages, pre-processing (modelling), analysis
(solve Navier-Stokes equations) and post-processing (results output). Among these CFD
modules, the solvers are the core of the problem.
3.2.2 CFD Solvers
The Navier-Stokes equations for flows of practical interest are usually so complicated
that an exact solution is unavailable, and it is necessary to seek a computational solution.
Computational techniques replace the governing partial differential equations with systems
of algebraic equations, so that a computer can be used to obtain the solution. Since the
governing equations for most cases of fluid dynamics are nonlinear, the computational so-
lution usually proceeds iteratively. Therefore, if the numerical algorithm that performs the
iteration is stable, and the discrete equations are faithful representations of governing equa-
tions, then the computational solution can be made arbitrarily close to the true solution of
governing equations by refining the grid. However, the refinement is restricted by com-
puter resources and accumulation of error during the iteration. Therefore, optimized grid
size should be applied to obtain sufficient accuracy.
In general, there are many methods for discretizing the physical domain, including fi-
nite difference method (FDM), finite volume method (FVM), finite element method (FEM)
and spectral method. FDM is very well developed and very straightforward to understand
the discretization of PDEs. However, FDM requires regular shape of mesh elements and
cannot be used to solve problems with complex geometric shape. In general, FDM is rarely
used in most of the contemporary CFD codes. FVM and FEM are found commonly used,
and are based on the theory of weighted residual methods (WRM).
A weighted residual method assumes that the solution can be represented analytically.
The starting point for a weight residual point is to assume an approximate solution, and
then the partial differential equations can be replaced by a system of algebraic equations or
61
a system of ordinary differential equations for time-dependent problems. By substituting

the exact solution, an equation residual is generated. The integral of the weighted residual
over the computational domain should be zero to determine the necessary coefficients of the
system of algebraic equations or the system of ordinary differential equations in general.
Different choices for the weight (trial) function give rise to different methods in the class
of methods of weight residuals. This forms the theoretical foundation of the finite element
and finite volume methods.
Finite Element Analysis
In the Finite Element Method, the solution domain can be discretized into a number
of uniform or non-uniform finite elements that are connected via nodes. The change of
the dependent variable with regard to location is approximated within each element by an
interpolation function. The interpolation function is defined relative to the values of the
variable at the nodes associated with each element. The original boundary value problem
is then replaced with an equivalent integral formulation. The interpolation functions are
then substituted into the integral equation, integrated, and combined with the results from
all other elements in the solution domain. The results of this procedure can be reformulated
into a matrix equation of the form, which is subsequently solved for the unknown variable.
The finite element method was developed initially as an ad hoc engineering proce-
dure for constructing matrix solutions for stress and displacement calculations in structural
analysis. The method was placed on a sound mathematical foundation by considering the
potential energy of the system and giving the finite element method a variational interpre-
tation. However, very few fluid dynamic (or heat transfer) problems can be expressed in a
variational form. For many situations, the Galerkin method is used to define the trial func-
tion or shape functions in the engineering literature for solving variational problems (Li,
Wang & Yi 1999). Consequently, most of the finite element applications in fluid dynamics
62
have used the Galerkin finite element formulation. The Galerkin finite element method has
two different features from traditional weighted residual methods. Firstly, the approximate
solution is written directly in terms of the nodal unknowns. The second important fea-
ture is that the approximating functions are chosen exclusively from low-order piecewise
polynomials restricted to contiguous elements.
The finite element method discretizes flow domain in two stages. First a piecewise
interpolation is introduced over discrete or finite elements to connect the local solution
to the nodal values. The second stage uses a weighted residual method to obtain alge-
braic equations connecting the solution nodal values. Both stages introduce errors. Small
element size can achieve more accurate interpolation, however, this induces more interpo-
lation steps and the error of each interpolation is accumulated, therefore, a restriction on
the smallest element size should be placed in the computation.
Finite Volume Analysis
Finite volume analysis divides the computational domain into many sub-domains, and
the weighted function is assigned unity if inside a sub-domain or zero if outside. The fi-
nite volume method calculates the values of the conserved variables averaged across the
sub-domain and provides an appropriate framework for enforcing conservation at the dis-
cretized equation level. In this way the conservation properties inherent in the governing
equations are preserved. This is a particular advantage in obtaining accurate solutions for
internal flows.
Finite volume method, which has been used for both incompressible and compress-
ible flow, has two major advantages. First it has good conservation properties, as mentioned
above. Second it allows complicated computational domains to be discretized in a struc-
tured and/or unstructured manner.
63
In summary, the finite volume method is preferable to other methods as a result of

the fact that boundary conditions can be applied non-invasively. This is true because the
values of the conserved variables are located within the volume element, and not at nodes
or surfaces. Finite volume methods are especially powerful on coarse non-uniform grids
and in calculations where the mesh moves to track interfaces. For the above reasons, the
finite volume method was used extensively in this research.
3.3 CFD Modelling in Microfluidics
For the dimensions of the microfluidic devices and the fluids investigated in this re-
search, the Knudsen number (Section 2.3.1) is approximately zero and continuum holds
true. Therefore, the computational fluid dynamic tools can be applied in the investigation
of passive microfluidic mixing. However, viscosity is critical and cannot be ignored for a
successful simulation.
In this study, a CFD package, known as MemCFDTM from CoventorWare, was ap-
plied. MemCFD provides both FVM and FEM solvers to simulate microfluidic problems.
For solving Navier-Stokes equations coupled with convection-diffusion equations, finite
volume method was used extensively for its advantage of mass conservation.
Although there is little significant change of governing equations in CFD solvers for
micro-flows, there are practical issues needed to be considered for a simulation in a micro-
channel. First, it is convenient to use a standard system of units and a modified system of
units for microfluidics was introduced in MemCFD. In this section, an analytical bench-
mark for CFD modelling and mesh sensitivity study are also introduced to calibrate an
accurate CFD simulation. To improve the performance of MemCFD, a mesh converting
program was developed to use an advanced meshing package to create a uniform struc-
tured mesh. The analysis of the simulation results are done by the post-processing tools
64
within MemCFD. Fortran programs were developed for further analysis of the results ob-
tained from the CFD simulations.
All simulations were performed on Win NT4.0 with Pentium III 800MHz CPU and
512MB memory. The maximum number of mesh elements (cells) in the models was up to
300,000 for the limitation of this class of workstation.
3.3.1 System of Units in Microfluidics
As the characteristic lengths of microfluidic devices are of the order a few microns to
a few millimeters, it is practical to use micron as the unit of length. The units of length (L),
time interval (T ), temperature (Θ) and electrical current (A) are regarded as the reference
(M-L-T-Θ-A) system, and other units are derived from them. Therefore, the dimensional
system for microfluidic are mass (kg), length (µm), time interval (second), temperature
(Kelvin) and current (pico Ampere). When applying the reference magnitudes to derive
other ones, the dimensional homogeneity principles must be followed. Then, we can derive
the following units used in microfluidics (Table 3.1). Other units, which are not listed in
this table, can be derived by the same rules.
3.3.2 Boundary Conditions
Boundary conditions are necessary for the CFD solvers to initialize the simulation
and achieve rational results. In microfluidic applications with liquids inside the channel,
the Reynolds number is low and the flow is generally laminar as stated in section 2.3. As
viscosity dominates the flow domain, the velocity close to the wall is zero. For most of
the simulations run in this research, the boundary conditions were set as steady, laminar,
Newtonian, with one fluid (or two fluids to solve diffusion-convection problems). The
inlets were assigned volumetric flow rate/velocity, or pressure boundary conditions.
65
Parameter Symbol Microfluidic Units

Mass m kg
Current C pA
Length l µm
Time interval t s
Temperature T K
Pressure p MP a
Density ρ kg µm−3
Diffusivity D µm2 s−1
Force F µN
Moment of Inertia I kg µm2
Power P pW
Stress σij MP a
Thermal Coefficient of Expansion TCE K −1
Velocity v µm s−1
Viscosity µ MP a s
Table 3.1: Systems of units in microfluidics
3.3.3 Building Numerical Models for Micro-Flows
MemCFD provides an integrated environment to include all of the three basic CFD
modules (Section 3.2.1). However, it does not give users sufficient control of the software,
especially the pre and post-processes. For some complex structures, MemCFD cannot cre-
ate a uniform structured mesh. However, a high quality CFD simulation depended very
much on the pre-processing stage, which includes meshing and boundary condition set-
tings. Therefore, mesh quality is very critical to the accurate simulations. More uniform
mesh and aspect ratio of mesh elements can improve the accuracy of numerical simulations.
To calibrate the CFD models, mesh sensitivity and uniformity studies were performed.
66
Mesh Sensitivity Study
For a general CFD solution of the governing equations over a range of significantly
different grid resolutions should be presented to demonstrate grid-independent or grid-
convergent result. Though it is not possible to use this method to address the truncation
error of the MemCFD solver, it is critical for an accurate CFD simulation. The empiri-
cal mesh density was estimated at least 20 element per edge for flows in a microchannel.
Denser mesh elements can give better accuracy, however, the solving and post-processes
would be slowed down and sometimes halt the system’s operation. To further determine
the size of mesh elements, an analytical benchmark of viscous flow in a rectangular micro-
channel was used to calibrate the mesh quality. Viscous flow in a rectangular micro-channel
also has practical value in studying microfluidics, and would be studied in Chapter 4.
To give a sample of benchmarking process, a viscous flow in a rectangular micro-

channel with width to depth aspect ratio (w/H) of 4 was used to calibrate the CFD model
(Figure 3.1). Figure 3.2 shows the structured mesh of the numerical model. The mesh was
constructed uniformly and had denser elements at walls.
For a given flow rate, the velocity field for a fully developed flow in a rectangular
channel could be predicted analytically by Equation 4.4 described in Section 4.2 or it can be
simulated by CFD solvers. The comparison of analytically prediction and CFD simulation
result gives good indication of the quality of numerical model. The comparisons could
be found in the line plots in Figures 3.3 and 3.4. The analytical prediction agreed with
the simulation results very well. The maximum error existed close to the centre of the
channel, and was less than 2.5%. As both the analytical predictions and CFD simulation
could contribute to errors, it was acceptable to determine the mesh density to be applied
when the difference was less than 5% in this research.
67
Fluid in
2b
2a
Fluid out
Figure 3.1: Viscous flow in a rectangular micro-channel
Mesh Uniformity Study
Because of the limitation for user to control the pre-processes in MemCFD, a state-of-
art mesh generation package, Gambitr from Fluent Inc., was used to improve the meshing
quality. Gambit could build a structured mesh by mapping or sub-mapping techniques.
However, the mesh topology in Gambit is different from that used in MemCFD (Figure
3.5). To use numerical models generated by Gambit, the mesh topologies needed to be
converted to the format that MemCFD could import. A Fortran program was coded to
convert the Gambit neutral file to universal neutral file compatible with MemCFD. The
algorithm in words is described below,
68
Figure 3.2: A typical meshed fluid volume in a rectangular micro-channel
1
0.8
u(y, zmid)
0.6
0.4
CFD simulation
0.2 Analytical prediction
0
-0.5 0 0.5
y/2a
Figure 3.3: Normalized velocity profile u(y, zmid ) in xOy plane

69
1.0
0.8
u(ymid, z)
0.6
0.4
CFD simulation (z plot)

0.2 analytical prediction
0.0
-0.5 0 0.5
z/2b
Figure 3.4: Normalized velocity profile u(ymid , z) in xOz plane

70
8 6 7 6
5 8 5
7
4 3 2
2
3 4
1 1
a b
TM
Figure 3.5: Mesh topology, (a). Gambitr 8-Hexahedral mesh; (b). MemCFD 8-
Hexahedral mesh
a b
Figure 3.6: (a). Structured 8-node Hexahedral mesh by sub-mapping technique; (b). Irreg-
TM
ular 8-node structured mesh by MemCFD
start
read Gambit neutral file
while ( < total number of elements N )
exchange node 3 and node 4
exchange node 7 and node 8
endwhile
write to universal neutral file
stop
The mesh illustrated in Figure 3.6a is more uniform than the mesh in Figure 3.6b
for a complex structure. The meshing uniformity of mesh is critical for an accurate CFD
simulation. Uniform mesh elements provide more accurate simulations and use less CPU
time.
71
Although the CFD prediction of velocity field with an irregular or unstructured mesh
was accurate, solving of convection-diffusion problem could be a problem. For instance,
to simulate a convection-diffusion problem in a microchannel of size H × w = 200µm ×
100µm, two fluids were assigned to two inlets respectively, and each with a volumetric flow
rate of 2 ×108 µm3 s−1 (Re = 5). The diffusion coefficient was 4 × 10−6 cm2 s−1 , which
was calculated according to the chemical groups of a food colouring dye (Reid, Prausnitz
& Poling 1987). The simulation results with uniform and irregular mesh are shown in Fig-
ure 3.7. Even with the same boundary conditions, Figure 3.7a shows inaccurate simulation
results using irregular mesh elements. Two cross sections ((1) and (2) in Figure 3.7b) of
mass distribution are shown in Figure 3.8, and are compared with the published experimen-
tal measurement by Stroock et al. (Stroock et al. 2002). The flow patterns of experimental
measurements and simulation results have the same trend, which indicate the simulation is
accurate and reliable.
3.3.4 Simulation Results Analysis
CFD simulations create a large amount of data to be analyzed. MemCFD provided a

visualization tool to post-process the results directly and a feature to export the results for
further analysis. The relevant functions included scalar plots of flow velocity, mass fraction
distribution, vector plots of velocity field, and flow trajectory (streaklines). Generally, the
visualizer displays the flow pattern inside a micro-channel, and the performance of micro-
mixer is evaluated by the comparison of mixing efficiency integrated from the mass fraction
distributions.
For some applications, the existing post-processing functions in MemCFD were not
adequate. Nevertheless, MemCFD provided a feature to export the entire simulation results
into a text file. This text file can be analyzed in other visualization tools or analyzed by
customized program. In Chapter 8, the velocity fields in a periodic grooved micro-channel
72
(a)
(b)
Figure 3.7: Mesh uniformity study, (a). irregular mesh elements; (b). uniform mesh
elements
73
Section 1 Section 2
Figure 3.7b Figure 3.7b
Figure 3.8: Mixing of two fluids in a T-channel with grooved surfaces, top: simulation
result, two fluids are represented by blue and red, green colour is the mixed region, bottom:
cited experimental results by Stroock, et al. 2002a, two fluids are represented by gold and
black colours.
were exported into text files and analyzed by a particle tracing program that was developed
to calculate the flow trajectories to investigate the folding and stretching of fluid streams.
3.4 Experimental Method
No matter how powerful the CFD solver is, if the numerical model or the bound-
ary condition is not correctly defined, the simulation result will not reflect the reality. In
this aspect, experiments are necessary to verify the numerical model. Because the various
numerical models used for simulating convection-diffusion problems in this thesis were
consistent with each other, only selected numerical models were tested by experiments.
Moreover, the simulation results of flow patterns in microchannels with patterned grooves
were compared with published experimental data, which were used to verify the CFD sim-
ulations of this type of micro-mixers.
In this section, an overview of micro-fabrication of micro-mixers is provided, followed

by description of an experimental setup used for the visualization of micro-mixing.
74
3.4.1 Excimer Laser Micro-machining
Excimer laser micro-machining was used to fabricate the micro-mixers in this re-
search, however, other micro-machining techniques would also be able to create the same
structures in a compatible manner but with silicon as the material. The micro-mixers were
fabricated in the microtechnology laboratory in the Industrial Research Institute Swinburne.
An Excimer Laser (ExiTech) was used to ablate channels onto polymer substrates. The
standard operation procedures could be found in early papers published by the labora-
tory (Harvey, Rumsby, Gower & Remnant 1995). In brief, the pulsed UV laser beam, with
a wavelength of 248nm, is projected through the designed features on a chrome-on-quartz
mask, and then the laser beam is focused to pattern the feature on the substrate (Figure 3.9
). The Excimer laser machined structures had a depth resolution of the order of 0.1µm and
spatial resolution of the order of 1µm.
Input Beam
248 nm (KrF) or 193nm (ArF)
Lens Array Homogenizer
Plane of uniformity
& Mask
Projection Lens
Workpiece step motion
Figure 3.9: A schematic illustration of Excimer laser micromachining

(Harvey et al. 1995)
75
Figure 3.10: Section of Y-channel machined in Polycarbonate by Excimer Laser (SEM

image)
3.4.2 Microfluidic Packaging
Figure 3.10 shows the Excimer laser ablated Y-channel on a polycarbonate substrate.
Following the ablation process, the channel was packaged in a lamination process. A thin
PET plastic foil coated with a melting adhesive film layer was pressed by a heated lami-
nation roller (at 150o C) onto the structure as the lid. Melinex°
R Polyester film type 301
(Dupont Teijin films) having a thickness of 30µm was used in our fabrication, and the
lamination equipment was a MEGA dry film laminator model 305 (Mega Electronics).
3.4.3 Microfluidic Mixing Visualization
Most of the conventional fluid field measurement methods were not suitable for mi-
crofluidic environments. Even though the µPIV technique has been researched since 1998
(Santiago et al. 1998), its complexity forbids it being widely used in general microfluidic
research centers. For studying microfluidic mixing, fluorescent dyes, food colouring dyes,
and ultraviolet absorption techniques have been used for measurement and evaluation of
76
CCD camera
Optical microscope Image capture computer
Observation
window
Figure 3.11: Laboratory setup for micro-mixing visualization

77
mixing performance. In this research, food colouring dyes are used to visualize the flow
and mixing process and is reported in Section 6.3.2.
Two colouring aqueous solutions have been applied in the study. One was mixture of
water and 2.4% food dye E102 and E122 (Yellow). Another one was mixture of water and
2.1% food dye E133 (blue) (Queen Fine Foods Pty Ltd, Australia). The two aqueous solu-
tions were respectively introduced to the channel from two inlets by capillary effect. The
viscosity of the food dye/water solution was estimated as being approximately the same as
water, and the diffusion coefficient was calculated according to the chemical groups (Reid
et al. 1987). The wastes from the outlet were collected by porus paper, and by measuring
the net weight before and after wetting in a known interval of time, hence, flow rates could
be calculated.
Figure 3.11 shows the laboratory testing devices for visualizing mixing performance.
The micro-mixer was filled with two different food colouring dye aqueous solutions, and
the broadening of mixed region was observed under an optical microscope. The image
was captured by a charge-coupled device (CCD) camera, and then transferred to the image
capture card in the computer.
3.5 Summary
To study passive mixing with various geometric structures, there were many param-
eters. Because the design of experiments became too complicated for all the experiments
to be carried out, it was impractical to investigate all these parameters by an experimental
approach.
On the other hand, most of micro flows are laminar with non-slip boundary conditions.
For laminar incompressible viscous flow, advanced microfluidic CFD techniques can be
applied to provide accurate solutions to microfluidic problems. To achieve a reliable CFD
78
prediction of fluid behaviour, the model needed to be built carefully with the assistance
of a sensitivity study. To calibrate the CFD models, Hele-Shaw flow in a finite aspect
ratio rectangular micro-channel was used. As the Hele-Shaw model was well-developed,
it could be used to validate the velocity field predicted by CFD simulation. To achieve a
reliable CFD simulation, the mesh density should no less than 20 elements per edge. Even
the velocity field was accurately predicted, and the simulation of mass transport was very
sensitive to the mesh quality. A uniform, evenly distributed aspect ratio mesh needs to be
applied with sufficient mesh density. To achieve this, an advanced mesh generator was used
and a mesh could be translated and imported into the MemCFD solvers using a program
developed in Fortran. Typical simulations were selected to be validated by experiments and
existing published experimental data. Simple experimental testing equipment was setup
for carrying out tests. The equipment included an optical microscope, a CCD camera and
an image capture computer. The simple laboratory configuration could give qualitative
comparisons between experiment visualizations and numerical simulations. A fabrication
technique to produce simple micromixers in microchannels was also described.
With the numerical and experimental methods developed, the subsequent research on
microfluidic mixing is presented in the following chapters. In Chapter 4, typical microflu-
idic flows in a rectangular microchannel were investigated to provide a necessary under-
standing of flow behaviours and mixing performance in a microscopic scale.
79
Chapter 4
Understanding Microfluidic Mixing
4.1 Introduction
To commence the research on microfluidic mixing, it is necessary to understand the

fundamentals of the scaling effects on the behaviour of microfabricated fluidic devices. The
scaling down of the characteristic length in microfluidic device resulted in low Reynolds
numbers, which initiated research to study fluid dynamics under low Reynolds numbers in
microfluidic devices.
The scaling down of a micro-device gave rise to the diffusion-controlled systems (Manz
et al. 1990). In such a device, the down-scale to 1/10 of size reduced the diffusion time
to 1/100. If the Reynolds number remained constant, the pressure drop increased 100
times. Therefore, the characteristic length of a micro-device was critical to the diffusion-
controlled system.
Taylor-Aris dispersion is an effect for dispersion of slugs of solute in a small tube with
a diameter of 500µm. It is an unwanted effect for most analytical processes, however, it
had unique ability to enhance mixing.
80
The aim of this chapter is to identify the mechanisms behind microfluidic flow be-
haviour, and therefore determine the necessary parameters, such as characteristic length
scale and flow rate, for the design of a micro-mixer. Section 4.2 introduces flow at low
Reynolds number, which is a characteristic feature of microfluidics. The understanding of
Taylor-Aris dispersion is essential to the design of miniaturized mixers, and is covered in
Section 4.3. Critical design parameters of microfluidic mixers were discussed in Section
4.4. Finally, Section 4.5 summarizes the outcome of this chapter.
4.2 Flow at Low Reynolds Number
The Reynolds number is used to relate the inertial forces to the viscous forces, and is
given by equation 4.1,
ρūl ūl
Re = = , (4.1)
µ ν
where ρ is the density, ū is the average velocity, µ is the dynamic viscosity, l is the charac-
teristic length of the microchannel and ν is the kinematical viscosity.
In microfluidic devices, due to the low velocity of the flow and the micron size chan-
nels, the Reynolds number is low. This means that the viscosity plays a dominant role
in microchannels rather than inertia, and turbulence cannot exist in such a viscous flow.
One classic paper described fluid dynamic differences at the low Reynolds number (Purcell
1977). Purcell stated that when Re ¿ 1, the inertia dies off from the Navier-Stokes equa-
tion, and gives the governing equation at such low Reynolds number as,
µ∇2 u = ∇p (4.2)
where p is the pressure and u is the velocity.
Equation (4.2) indicates that the velocity is determined only by the pressure distribu-
tion and all motion is symmetric in time at low Reynolds number, as described in a paper by
81
Brody (Brody, Yager, Goldstein & Austin 1996). This means that the flow can be reversed
back to its original status.
Fluid in
y
Fluid out
Figure 4.1: Rectangular microchannel
Microchannels with a rectangular section were widely used in the microfluidic appli-
cations. The slow viscous flow in a rectangular channel is very similar to a Hele-Shaw flow.
The literature review found vast literature on the Hele-Shaw flow, however, this was out of
the scope of this research. Nevertheless, the results from studying Hele-Shaw flow could
be used to study the flow behaviour in rectangular microchannels, which were the most
commonly used shape in microfluidic applications. Considering the velocity distribution
for the flow in a rectangular channel shown in Figure 4.1, based on conservation of mass,
82
or continuity in other words, the velocity remains the same at the x direction for a fully
developed incompressible flow. The 2D velocity distribution, u(y, z), satisfies Poisson’s
equation 4.3, which is rewritten from equation 4.2,
1 dp
∇2 u(y, z) = − (4.3)
µ dx
where µ is the dynamic viscosity and dp/dx is pressure gradient and assumed an absolute
constant (Lamb 1993).
With necessary wall boundary conditions, uwall = 0, a Fourier series solution of ve-
locity field u(y, z) in a rectangular channel size of −w/2 ≤ w/2 and −H/2 ≤ H/2 can be
written as:
µ ¶X∞ · ¸
4w2 dp (n−1) cosh((2n − 1)πz/w) cos((2n − 1)πy/w)
u(y, z) = − (−1) 1− ,
µπ 3 dx n=1
cosh((2n − 1)πH/2w) (2n − 1)3
(4.4)
where the velocity varies in the channel width to depth aspect ratio w/H. As hyperbolic
function cosh approaches infinity with n, it is non-trivial to achieve exact solution of above
equation 4.4. However, as the sum converge very rapidly, a finite but sufficient large n
value can be used to estimate the solution by numerical computation (see Appendix B.2).
By integrating equation 4.4 over the area of the channel section, the volumetric flow
rate can be obtained and stated as,
µ ¶ ∞ · ¸
8Hw3 dp X 1 2w
q= − − tanh((2n − 1)πH/2w) (4.5)
µπ 4 dx n=1 (2n − 1)4 (2n − 1)5 πH
where q is the volumetric flow rate, dp/dx is the pressure gradient along x. For conve-
nience, the pressure gradient can be related to the mean velocity ū and stated as,
dp 4kµū
=− 2 (4.6)
dx H
where k is a constant related to the aspect ratio of a rectangular channel, and the value of
k, which can be calculated by equations 4.5 and 4.6 with a finite n, is given in Table 4.1.
83
Using equations 4.5 and 4.6, for a known volumetric flow rate, the velocity field for
a fully developed flow in a rectangular channel was computed. Figure 4.2 shows the flow
field in a rectangular channel with width to depth aspect ratio (w/H = 4).
0.8
0.6
u*
0.4
0.2
0
0.5
0.5
z* 0
0
y*
−0.5 −0.5
Figure 4.2: Normalized velocity profile in a rectangular microchannel, w/H = 4; u∗ is

normalized velocity, y ∗ = y/w and z ∗ = z/H
By cutting the normalized velocity profile from Figure 4.2 along the mid of y and
project it on the XZ plane, Figure 4.3 shows that this velocity profile is very close to a
parabolic shape in shallow channel. Using the same method but for various w/H ratios,
Figure 4.4 shows the velocity profiles projected on the XY plane, and that the velocity
profile moves to a plug-like flow when the w/H ratio becomes larger. From Figure 4.4, it
can be concluded that the side-wall effects on velocity profile are negligible when the ratio
of w/H → ∞, which is the Hele-Shaw flow.
It is important to notice that the velocity profiles are not restrained by Reynolds num-
ber when flow is in the laminar region, therefore, the results can be applied to low Reynolds
number microfluidic flows.
84
w/H 1 2 3 4 5 10 20 100 ∞
k 28.454 17.492 15.191 14.244 13.731 12.807 12.3905 12.076 12.002
Table 4.1: Values of k in a rectangular-cross-section channel for various width to depth

aspect ratios
1.0
0.8
u(ymid , z)
0.6
0.4
0.2
0.0
-0.5 0.0 0.5
z/w
Figure 4.3: Normalized velocity profile by the maximum flow velocity plotted in XZ plane
for y = ymid
85
y/H
u*(ymid,z)
Figure 4.4: Normalized velocity profile plotted in XY plane for w/H = 2 (O), w/H = 5
(×) and w/H = 10 (+).
4.3 Taylor-Aris Dispersion
4.3.1 Overview
The nature of microfluidic mixing is similar to a process outlined by Taylor (Taylor

1953). Therefore, it is important to understand the background knowledge of Taylor Dis-
persion. Section 4.3.2 reviews dispersion of slowly moving solute in a rectangular cross-
section microchannel, which can be termed Taylor-Aris dispersion. Section 4.3.3 investi-
gates the application of Taylor-Aris dispersion to enhance mixing in a microchannel with
rectangular cross-section.
86
4.3.2 Dispersion of a Solute in Solvent Flowing Slowly Through a Rect-

angular Cross-Section Microchannel
In Taylor’s work, a section of one fluid was introduced into the flow of another fluid
in a small circular cross-section tube (φ500 µm). The flow was laminar, with a parabolic
velocity profile, where the velocity in the center of channel was maximum and the velocity
at walls was zero. Therefore, the introduced fluid streams in the tube are stretched. Be-
cause the velocity closest to wall is slower, the fluid streams have more time to diffuse into
each other. When the mean velocity is slow and tube is small enough, the diffusion will
sample all the streamlines in the tube. That is the reason that the fluid spread around a
cross-sectional plane moving at the mean fluid velocity is similar to fluid diffusing from
such a plane in still fluid. This spreading of solute could be characterized by a dispersion
coefficient instead of a pure molecular diffusion coefficient. While Taylor’s work dealt
with a round tube, a similar analysis based on his ideas was applied to a 2-D rectangular
channel (Doshi et al. 1978, Chatwin & Sullivan 1982), which can be useful for mixing in a
microchannel.
The concentration distribution in a rectangular channel can be described by the con-

vective diffusion equations,
· 2 ¸
∂C ∂C ∂ C ∂ 2C ∂ 2C
+ u(y, z) =D + + (4.7)
∂t ∂x ∂x2 ∂y 2 ∂z 2
87
along with the following boundary conditions:



0 0 ≤ y ≤ w/2,
C(0, y, z) =

1 −w/2 ≤ y < 0,
∂C
=0
∂t ¯
∂C ¯¯
=0
∂x ¯x→∞ (4.8)
¯
∂C ¯¯
=0
∂y ¯y=0,±w/2
¯
∂C ¯¯
=0
∂z ¯z=0,±H/2
C(∞, y, z) = constant
where C is concentration distribution.
The analytical solution to the above problem is non-trivial. Although there was attempt
to use Taylor’s approach to give the analytical prediction, the result created controversy
because the Taylor-Aris dispersion in rectangular microchannel was different from that
used in the original study of dispersion in a small circular pipe (Beard 2001, Dorfman &
Brenner 2001). Most of the microfluidic mixers researched introduced confluent streams
into the microchannel, and flows were steady (∂C/∂t = 0). The concentration gradient at
any position of the channel does not change with time, and this is a fundamental difference
to the Taylor-Aris dispersion.
The Taylor-Aris dispersion could and did inspire the study of microfluidic mixing,
due to the virtual diffusion coefficient D∗ (Taylor’s dispersivity) being several orders larger
than the pure molecular diffusion coefficient. Micro-mixers had been introduced to apply
Taylor-Aris dispersion using an active pulsed source (Niu & Lee 2003), and was claimed
as chaotic mixer. Though it is true the Taylor dispersion is a chaotic system with time
variable in a 2D Hamilton system, it is much easier to comprehend with well-developed
Taylor-Aris dispersion theories. Nevertheless, the Taylor-Aris dispersion is a special case
88
of the convection-diffusion problem and points out that introducing lateral convection is a
way to enhance microfluidic mixing.
Taylor’s dispersivity in a rectangular channel can be described in equation 4.9,

1 u2 w 2 H
D∗ = D + f( ) (4.9)
210 D w
where D is the molecular diffusion coefficient, w is the width of the channel, ū is the mean
velocity in the channel, H is the depth of the channel and f is the geometric function of a
rectangular channel, decided by the height to width aspect ratio H/w. While H/w ¿ 0,
assuming no side-walls, it gives that the geometric function f a unity value, however, with
the consideration of side-walls, f = 7.95. Therefore, the side-walls of a microchannel
have to be considered while investigating mixing problems. For a Taylor dispersivity D∗ ,
the length between the dispersed leading and trailing edges of slug, lm , is estimated by the
√
Einstein-Smoluchowski relation (lm = 2D∗ t).
10000
1000
D*/D
100
10
1
0 200 400 600 800
u (µm/s)
Figure 4.5: Ratio of Taylor dispersivity to pure molecular diffusion coefficient for different
Reynolds numbers
Taylor dispersivity D∗ is much larger than the molecular diffusivity D (Figure 4.5),
89
therefore, for a fixed width of channel w, the mixing time t is reduced. Analytical results
can evaluate Taylor’s dispersivity for rectangular channels with different depth to width
aspect ratio (Doshi et al. 1978). The transient behavior of most practical applications are
difficult to interpret by analytical approaches, however, it is much easier to interpret through
numerical simulations.
4.3.3 Applying Taylor-Aris Dispersion to Enhance Mixing
The simulation results in Figure 4.6 demonstrate a time series of mass fraction of
a slug of fluid 1 dispersed in fluid 2 for a time period of 1 to 20 seconds. Boundaries
conditions were kept the same as for Taylor-Aris’s analytical interpretation, the flow was
very slow (Re = 0.005) and all wall boundaries were assigned zero velocity. Therefore,
the simulation could be verified by comparing with the analytical prediction.
Figure 4.6 clearly shows that the length of the dispersed fluid is longer than a pure
diffusion length, because of the convection. As Taylor’s dispersivity coefficient is a func-
tion of mean velocity, D∗ = f (ū), the inference is to have a large flow rate for fast mixing.
However, as velocity increases, there will be less than sufficient time for the diffusion to
sample all the streamlines.
In practice, the velocity in a microchannel is mostly too large to satisfy Taylor-Aris

requirement for an analytical solution. However, the effects of Taylor-Aris dispersion still
exist in a faster moving stream with a more complex flow pattern.
Simple experiments were carried out to demonstrate mixing of two fluids in a Y-type
microchannel, illustrated in Figure 3.10. Initially, the two fluids flowed parallel to each
other without significant mixing. When one fluid supply was stopped, the remaining of
this fluid was embedded into the main stream and its dispersion into the second fluids
90
1sec
2sec
5sec
8sec
10sec
12sec
15sec
18sec
20sec
Figure 4.6: Taylor dispersion along a straight channel: numerical simulation for various
time intervals
was enhanced dramatically. Later in the research, it was found that the phenomena was
termed Taylor dispersion in the literature. Nevertheless, the experimental results are shown
in Figure 4.7, where the green region was the mixed fluids between fluid 1 (yellow) and
2 (blue). The dispersed solute was initially a narrow parabolic shape at time t1, and then
developed into a wider parabolic shape at time t2 in a short distance. Although there is lack
of quantitative measurement for evaluating the mixed volume of the two fluids, the enlarged
area of mixed green region and lengthened parabolic shape indicated the enhancement of
mixing. In addition, this phenomena was repeated several times during the experiments
and demonstrated similar results.
91
500µm
500µm
(a) (b)
Figure 4.7: Experimental results: the solute dispersed at, (a) time t and (b) time t + δt
The simulation results shows the same as experimental results that the broadening of
the slug of fluid1 in the axial direction, and the result demonstrated the actual dispersivity
was much faster than pure molecular diffusion. The matching of simulation results and
experimental results also shows the accuracy of the simulation and prove CFD can be a
reliable tool in studying microfluidic mixing.
4.4 Two Important Parameters for Designing a Microflu-

idic Mixer
Unlike chamber mixing, passive microfluidic mixing is an in-situ process in a mi-

crochannel. The mixing performance is measured by the length to achieve complete mix-
ing, which can be measured by the Péclet number. From the literature review, it was found
that the Péclet number is a critical parameter for the design of a micro-mixer, and hence,
discussed first. In addition, the aspect ratio of a rectangular microchannel is another impor-
tant parameter regarding the diffusion path and pressure drop, and is also discussed in this
section.
92
4.4.1 Péclet number
For mixing at low Reynolds number, diffusion is dominant and mechanical agitation
is ineffective. Molecular diffusion is a random process, and the diffusion time tD to cross
a distance l was given by Einstein in his paper on Brownian motion (Einstein 1905). For
typical small molecular liquids with a diffusion coefficient D of 103 µm2 /s, ie. H2 O, the
diffusion time to cross a 100 µm distance is about 5 seconds. However, for 1mm distance,
the mixing time is more than 8 minutes. For most bio-molecules with diffusion coefficients
one order or even two orders lower than that of water, diffusion time is of the order of hours,
which is not considered realistic in a microscaled device. For mixing to be completed in
milliseconds, the diffusion path should be of the order of 1 µm.
Considering microfluidic devices with microchannels of the order of a few hundred

microns, diffusion is not very efficient due to the large diffusion path. However, for low
Reynolds number flow in a microchannel (Re < 10), there is no turbulence to assist mixing,
therefore convection is the only method for mass transport in laminar flows other than pure
molecular diffusion (Purcell 1977). The convection time through the channel can be noted
as,
Lm
tc = (4.10)
ū
where ū is the mean velocity, and Lm is the length of the channel for mixing.
The molecular diffusion time tD through a distance of a can be stated as,
a2
tD = (4.11)
D
where a is diffusion path (half of the width of channel in this research) and D is the diffu-
sion coefficient. Therefore, if tc ≥ tD , the fluids can be diffused into each other completely
by the end of the channel. Hence, the least channel length for complete mixing can be
stated as, Lm ,
a2
Lm ≥ ū (4.12)
D
93
Lm can be refined as,

Lm ≥ ū · tD = P e · a. (4.13)
where Péclet number, P e, is the ratio of bulk mass transfer (Convection or Stirring) to
molecular diffusivity,
aū
Pe = (4.14)
D
According to Equation 4.13, it would be desirable to have a small Péclet number in

order to have a better mixing. When the flow velocity is very slow to have a small Péclet
number, therefore, fluids have sufficient residual time for diffusion to take place. A com-
plete mixing can be achieved without further assistant. However, it can be found that flow
velocity is a constrained design parameter to supply sufficient mixture for certain applica-
tion. In such cases, to reduce a in Equation 4.14 is one way to reduce Péclet number. From
previous review, convection induced virtual diffusion coefficient D∗ (Taylor’s dispersiv-
ity) is several folds higher than pure molecular diffusion coefficient, to introduce Taylor’s
dispersivity D∗ is another way to reduce Péclet number in Equation 4.14.
Using the Péclet number, microfluidic mixing can be divided into three zones (Figure
4.8). To clarify Figure 4.8, the diffusion zone means diffusion dominates and convection
is neglected, and the convection zone means the diffusion can be neglected. Figure 4.8
was used to define the scope of this research. If the maximum Péclet number for effective
mixing was defined as P et . When P e > P et , complete mixing at end of the channel is not
achievable. In this case, a method has to be introduced to reduce the Péclet number. If a
(in equation 4.14) is redefined as the distance for fluids to diffuse through and D is defined
as the virtual diffusion coefficient D∗ , the upper limit of Péclet number is restricted by
the ability of a micromixer to reduce a and enlarge D∗ to bring the Péclet number back to
the convection-diffusion region in Figure 4.8. Therefore, for a specific microchannel with
known L/w ratio and velocity, the design of a micromixer should be able to reduce a and
enlarge D∗ to make P e ≤ P et , where P et is defined as the threshold Péclet number.
94
L/w
104
Convection-
Diffusion
103
Diffusion
102
Convection
101
Pet Pe
10-6 10-4 10-2 100 102 104 106
Péclet Number
Figure 4.8: Three microfluidic mixing zones defined by using the Péclet number
When the flow is sufficiently slow, the Péclet number is less than unity, the flow is
in the diffusion zone, mixing can be achieved without further assistance. And when the a
Péclet number that is larger than 106 , because the flow velocity is so high that it is difficult to
transport in micro-scale and the dispersion of solute is by pure convection for Pe À a/D
according to the Taylor-Aris description (section 2.3.4). For such a high flow velocity,
microfluidics loses its advantage for mixing. The relation between Reynolds number and
Péclet number can be calculated by,
Pe × D
Re = . (4.15)
ν
From Equation 4.15, range of Reynolds numbers can be calculated, and Reynolds numbers
0.01 to 100 were chosen as the scope of this research for the above reasons.
95
4.4.2 Aspect Ratio of a Rectangular Microchannel
The width to height aspect ratio is one of the critical technical specifications for fabri-
cation of microchannels. For a given flow rate, the mixing performance acts differently for
different aspect ratio. For two fluids flowing parallel in a channel, the diffusion occurs at
the interface of the two fluids (Figure 4.1). The fabrication of low aspect ratio microstruc-
tures is generally easier, however, for mixing, low aspect ratio microfluidic devices give
small contact area between fluids.
Gobby, et al. investigated the aspect ratio with constant width (Gobby, Angeli &
Gavriilidis 2001), and we extend this investigation to a constrained area. For simplicity, we
assume the bi-fluids travel at same velocity ū or flow rate, q, through a channel with cross
section area A. Thus,
q = Aū, (4.16)
where
A = wH, (4.17)
and the channel depth to width ratio β is
H
β= , (4.18)
w
then,
q
u= (4.19)
βw2
Einstein’s expression on dispersion of particles in a fluid stream can be stated as equa-

96
tion 4.20,
(w/2)2 w2
t= = (4.20)
2D 8D
Using equation 4.20 and 4.19, we can obtain the travel distance Lm for mixing in time t,
q
Lm = ūt = . (4.21)
8βD
Therefore, if the flow rate is constant, the minimum length of the channel for a com-
plete mixing, Lm , is identified as the reciprocal of the aspect ratio of the channel.
Equation 4.21 is an estimation of mixing length, and satisfies equation 4.13. Because
chemical analysis requires low noise to signal ratio, therefore, a small height to width as-
pect ratio for sensing applications is widely used to achieve low optical signal distortion,
as a result, a long mixing length (Lm ) is required (Kamholz & Yager 2001). Conversely,
large channel height to depth aspect ratio increases the contact area of confluent fluids to
mix fluids rapidly to have high productivity without concerning about any signal distor-
tion. Therefore, large channel height to depth ratio is preferable for chemical synthesis
productivity.
Using equation 4.6 and the values in Table 4.1 computed by the program described
in Appendix B.2, the relation between pressure drop and channel aspect ratio can be cal-
culated. Figure 4.9 shows how the aspect ratio of a rectangular microchannel affected the
pressure drop. When designing a microchannel for mixing applications, the aspect ratio
should be optimized to favour the mixing, however, the pressure drop needs also to be con-
sidered. Under any condition, the least pressure drop happens at unity aspect ratio. Pressure
drops for channel aspect ratios β < 1 have a symmetric relation with β > 1 against β = 1,
and were omitted in Figure 4.9.
When the mixing is due to pure molecular diffusion, equation 4.21 holds true when the
flow is fully developed. Figure 4.10 shows shorter channel required for large channel width
to depth (β) ratio, because the diffusion works efficiently for flow in narrower channels.
97
1.0E-05
Pressure gradient (MPa / µm) 1.0E-06
1.0E-07
1.0E-08
Re = 0.2
1.0E-09 Re = 2
1.0E-10 Re = 20
1.0E-11
1 10 100
Channel depth to width ratio, β
Figure 4.9: Pressure drop along a rectangular microchannel with the same cross-section
area, for different height to width aspect ratio, at various flow rates q
1.0E+06
Minimum channel length for
complete mixing ( µm)
1.0E+05
1.0E+04
1.0E+03
Re = 0.2
1.0E+02
Re = 2
Re = 20
1.0E+01
1 10 100
Channel depth to width ratio, β
Figure 4.10: Aspect Ratio and flow rate affect the total length to achieve complete mixing
98
For other arbitrary cross-section shaped microchannels, the pressure gradient, dp/dx,
for fully developed laminar flow is approximately (Pan 1988),
dp 128µq
= , (4.22)
dx πd4h
where the hydraulic diameter dh is given by the equation 4.23
4 × Area
dh = (4.23)
wetted perimeter
4.5 Summary
In this chapter, mass transport at low Reynolds number was introduced. At low
Reynolds number, the inertial component in the Navier-Stokes equation is relaxed and vis-
cosity is predominant in the flow. Because rectangular Microchannels were used regularly
in Microfluidic devices, the behaviour of a viscous incompressible flow in a rectangular
was studied.
The Taylor-Aris dispersion in a small channel is well-known in the analytical chem-

istry for its unwanted features, however, its effect may be a way to disperse solutes con-
vectively in the channel. The Taylor dispersivity is normally a few orders greater than
pure molecular diffusion coefficient, therefore, it is effective to create Taylor dispersion
in a micro-mixer. Section 4.3.3 described the dispersion of a solute slug in a straight mi-
crochannel by simulation and experiment. The result demonstrated that it was a promising
method to broaden the dispersion of solute to a fluid in a microchannel. To apply Taylor
dispersion in a micromixer, solute slugs need to be introduced into the microchannel con-
tinuously to form a sandwiched structure, which may be achievable by an active mechanical
shutter inside a microchannel. However, any active mechanism inside a micro-device will
unavoidably increase the complexity of the system, and the complicated control system
99
to cope with diffusion coefficient with the pulse frequencies and external energy source
requirement make it a non-preferred mixing method. In addition, for most of the microflu-
idic mixers which introduce confluent fluid streams into the channels, it is difficult to apply
Taylor dispersion using an active mechanism.
Taylor-Aris dispersion is a special case of convection enhanced mixing, and the intro-
duction of convection effects forms an important microfluidic mixing method. To evaluate
the convection effect, the Péclet number, which is the bulk mass transport to diffusion mass
transport ratio, was introduced. For a microfluidic device, it was found that only a range
of Péclet numbers need to be investigated. At very low Péclet numbers, mixing is in the
pure molecular diffusion zone, and mixing methods other than diffusion are not required.
Also, at very high Péclet numbers, mixing is in the pure convection zone. Because all the
mixing processes require diffusion to complete the mixing, therefore, a very high Péclet
number is also not preferable for mixing in microfluidic devices. Moreover, a high Péclet
number induces high pressure drop and can make it difficult to drive the fluids through the
microchannel. Péclet number is an important parameter for the design of a micromixer,
and was also used to form the scope of this research.
Another important parameter is the geometry of a micromixer. Because many mi-

cromixers are structured with rectangular shaped cross-sections, the microchannel depth
to width aspect ratio, (β), was investigated. It was shown in Section 4.4 that mixing with
a constrained flow rate and channel transversal area was affected by β, and the increas-
ing value of β benefited mixing. Therefore, one can expect good mixing performance in
high depth to width aspect ratio microchannels. On the other hand, it was found that mi-
crofluidic analytical devices, ie. µTAS, required a small β value to reduce the distortion of
flow pattern. A small β value increases the width for the fluids to travel, and hence, it has
poor mixing performance. Therefore, the attempts in the following chapters were made to
find methods to enhance mixing performance in a microfluidic diagnostic device with low
aspect ratio microchannels.
100
In Chapter 5, three variations in design of serpentine structures or serpentine-like

structures are investigated to try to overcome the shortcoming of using high Reynolds
numbers in the 3D serpentine micro-mixers. In these designs, transverse convection was
brought into the main stream in the longitudinal direction. Hence, the convection could
disturb the main flow to induce bulk mass dispersion.
On the other hand, Taylor-Aris dispersion is a special case of convection enhanced

mixing. Inspired by the Taylor-Aris dispersion, methods need to be found to introduce
convection into enhance mixing.
101
Chapter 5
Microfluidic Mixing by Partitioning

Flows
5.1 Introduction
It may be intuitive to use a narrow channel to mix two confluent fluid streams rapidly,
because of the reduction of diffusion distance, however, for pressure driven flow, the pres-
sure drop in a narrow channel is high. Therefore, wide channels are usually used and
sometimes divided into multiple narrow sub-channels. Then, fluid streams could be intro-
duced into these multiple narrow sub-channels alternately like sandwich structures (Koch
et al. 1998). This method of mixing was termed the splitting and recombining technique,
and it occupied a large area of a microchip and such designs were too complicated to fab-
ricate. In some of these structures, the pressure drop was so high that it was very difficult
to drive the fluids through the channels (Schwesinger et al. 1996).
In this chapter, three types of serpentine-like designs were studied to improve the split-
ting and recombining techniques to reduce diffusion path. In addition, a guideline for these
102
designs was to introduce the transverse convection. The studies were carried out by sim-
ulating mixing in all three types of mixers using MemCFD. In section 5.2, a single-side
serpentine shaped channel was divided into a number of mixing zones, where concentra-
tions were different. In section 5.3, a dollar-shaped serpentine channel is introduced. This
serpentine channel was short-circuited by a straight channel. Inspired by the complex 3D
structures that were machined by an Excimer laser, a 3D periodic ramping structure is
investigated in section 5.4. Section 5.5 summarizes the findings of Chapter 5.
5.2 Dividing Fluids into Discrete Concentration Zones
In many chemical or biological diagnostic or analysis applications, it was not always

necessary to achieve a homogenous solution. In such applications, a mixed zone that satis-
fies the concentration requirement might be sufficient for conducting an analysis. Inspired
by the splitting and recombining techniques in microfluidic mixing (section 2.6), a single-
sided serpentine channel with splitting zones was designed. However, instead of dividing
streams parallel to the flow direction, the main streams were divided into multiple sub-
streams laterally to the main flow. For the continuity of incompressible flow, the velocity
in the divided branches had to slow down to keep a constant volumetric flow rate. There-
fore, the design provided more residual time for diffusion to take place in certain divided
zones. The divided zones also reduced the diffusion path, and homogenous mixing can be
achieved in some of the zones.
The results in Figure 5.1 show that two fluid streams coloured red and blue were
pumped into a T-type mixer via two inlets. The mass fraction of fluid 1 was determined
as a relative percentage to fluid 2. A mass fraction of 0.0 indicates 100% concentration of
fluid 1, and a mass fraction of 0.5 denotes a complete mixture with equal quantities of fluid
1 and 2. A mass fraction of 1.0 denotes a 100% concentration of fluid 2. The main streams
were divided into eight zones. As the concentration of two streams equaled each other, the
103
100.0%
90.0%
80.0%
4 5
Mass Fraction
70.0%
60.0% 3
50.0%
40.0% 2 6
30.0% 7
20.0%
10.0% 1 8
0.0%
0 500 1000 1500 2000 2500 3000 3500
position ( µ m)
(a)
1 2 3 4 5 6 7 8
300µm
(b)
Figure 5.1: Different concentration gradient separated in the branches, (a). relative mass
fraction distribution in percentage; (b). corresponding mass fraction in contour plot
104
mass fraction of homogenous solution would be 50%. Zones 2, 6 & 7 were approximately
homogenous, and other zones had different mass fraction of fluid 1.
5.3 Dollar Sign Shaped Micromixers
In the fork shaped micromixers reviewed in section 2.6.2, there were two major con-
cerns: pressure drop and pressure distribution. The high pressure drop limited its practical
application and it also had difficulty in distributing pressure evenly to the fork modules. In
addition, the serpentine shape microchannels were commonly used to increase the residual
time of fluids for sufficient diffusion. The serpentine and the fork-shaped mixing designs
were combined and a by-pass channel was introduced to short circuit the fluid streams. This
by-pass channel disrupted the main flow to create lateral convection, and also, partially split
the main fluid streams into the by-passing shortcut channel. This design is referred to as
the dollar sign shaped micromixer ($-mixer) and is shown in Figure 5.2, and was developed
as part of the present research.
The width of the bypassing channel was smaller than the width of serpentine channel,
therefore, only some of the fluids flowed into the by-passing channel. The remaining fluid
stream formed a thinner layer close to the wall. Figure 5.3 shows the splitting and recom-
bining of fluids streaklines inside the channel, mapped with mass fraction distribution. As
explained in Section 4.2, the flow pattern in a rectangular channel was nearly parabolic.
The velocity closest to the wall was slower than the flow in the center of the channel. This
gave fluid streams that were closest to the walls more residual time for diffusion to take
place. At the next turn of serpentine, the other fluid stream was partially split and the pro-
cedure repeated itself. In the simulation, because of the limitation of computer memory
and CPU speed, the size of the numerical model was restricted and only one section of
serpentine was illustrated in this research. Practically, many turns of serpentine could be
applied to achieve further mixing.
105
200µm
Figure 5.2: Mass fraction in a planar $-mixer

106
Figure 5.3: Streaklines mapped with mass fraction in the planar $-mixer
107
To further improve the performance of $-mixer, a multi-layer structure could be ap-

plied. Figure 5.4 shows the mass fraction distribution with a 3D variation of planar dollar-
mixer. The by-passing fluid stream could pass underneath the main stream. By doing this,
the by-passing fluid could join the other fluid at other side, which is impossible in the pla-
nar structures. Illustration of the flow using streakines, where the splitting and recombining
process was very much like plaiting fluid filaments into a pigtail is shown in Figure 5.5.
The plaiting blended two fluid streams into each other, and reduced the diffusion path as
well as increasing residual time, which was similar to planar structures.
200µm
Figure 5.4: Mass fraction in 3D variation of the planar $-mixer

108
200µm
Figure 5.5: Streaklines mapped with mass fraction in the 3D $-mixer

109
5.4 3D Ramping Micro-mixers
5.4.1 Overview
Using the Excimer Laser micromachining process, very complex channels structures
could be fabricated. Figure 5.6, shows two ramping structures made by Excimer laser mi-
cromachining in the microlithography laboratory in the Industrial Research Institute Swin-
burne (Hayes, Harvey, Ghantasala & Dempster 2001). These ramping structures were
created by periodically moving the machine platform and mask rotation. Therefore, the
ramping structures also had the periodic configurations. To some extent, the ramping struc-
tures could be recognized as serpentine-like structures. While the design intention was
to improve mixing, it had not been demonstrated whether these novel ramping structures
could actually benefit microfluidic mixing or not. Therefore, a study was carried out in this
section.
a b
Figure 5.6: Two samples of 3D ramping structures made by Excimer Laser
micromachining.
(courtesy of Dr. J.P. Hayes)
110
5.4.2 Simulations of 3D Ramping Micro-mixers
The ramping structures described above were very complex for computer modelling.
Therefore, the geometric model was simplified but the critical features were kept. Then,
the geometric model was discretized for CFD simulations to examine its effect on mixing
(Figure 5.7). To achieve rational interpretation of simulation results, mesh elements must
be structured and have high uniformity. The mesh elements were created by mapping and
sub-mapping techniques (described in Section 3.3).
100µm
Figure 5.7: Computer model of a ramping structure meshed with 8-node Hexahedrons
5.4.3 Discussion of Results
The ramping structures created 3D convective effects. The interfacial region between
two fluid streams was moved by convection created by the ramping structures (Figure 5.8).
When the interfacial area was in the slow velocity region closest to the wall, the fluid had
more residual time to diffuse across the streamlines. Because diffusion is a irreversible
111
Figure 5.8: A 3D ramping mixer, at Re = 5

112
process against time, hence, mixed fluids remain mixed even when the region is moved
back to the high velocity region closer to the center of the channel due to the convection. A
periodic structure was necessary to repeat this process to achieve a required concentration.
When the Reynolds number was high (Re ∼ 50), this created strong secondary convection
(secondary flows). If the periodic structure of 3D ramping could generate periodic velocity
field, the strong secondary flows may lead to recirculations around the flow axis (Figure
5.9). Figure 5.10 shows the irregularity of flow pattern at a high Reynolds number. The
flow pattern demonstrated that the mixing not only relied on the diffusion across stream-
lines, but also and more importantly relied on the cross-streamline advection. In some
literature reviewed in Chapter 2, it was mentioned that the high flow velocity may give rise
to chaotic advection. However, as stated in the previous Section 4.4, this research concen-
trated on investigating mixing at low Reynolds number. Moreover, chaotic advection is not
a necessary condition for ideal passive microfluidic mixing, because while it reduces the
diffusion residual time, it requires longer channels for complete mixing.
The results suggest that chaotic mixing in a channel with periodic 3D ramping struc-
ture provide a basis for further research.
5.5 Summary
The work in this chapter showed that it was straight forward to split and recombine
fluid streams directly by using geometric topology. Three design concepts were presented
that could generate different concentration gradients.
A modification of the serpentine microchannel design concept showed that by-passing

channels could be added to create partial flows. The partial flow reduced the amount of
fluid entering the main channel, the fluid stream became thinner and closer to the wall.
Under this condition, the diffusion path was reduced and residual time was increased.
113
Figure 5.9: Mass fraction distribution in a 3D ramping mixer, at Re = 50

114
Figure 5.10: Streaklines in 3D ramping mixer mapped with mass fraction, at Re = 50

115
The last approach was to use 3D ramping structure to enhance mixing. It was difficult
to fabricate 3D ramping structures using techniques other than Excimer Laser microma-
chining. The ramping structures could create convective effects to enhance mixing, and at
high Reynolds numbers (Re = 50), irregular flow pattern was created. While mixing using
a high Reynolds number was not the focus of this research, for the sake of completeness,
the results were included and considered useful for future research.
The 3D ramping structure was novel and could be used to enhance mixing in a mi-
crochannel. However, it was difficult to define parameters for the study. Therefore, the 3D
ramping structures needed to be simplified.
In Chapter 6, these ramping surface structures are further simplified to cylindrical

obstacles. It was envisaged that the cylindrical obstacles could provide disturbances to
the flow that was similar to the ramping structures, and resistance and disruption could be
represented by number and position of obstacles. By studying distribution of mass fraction
and velocity field affected by the obstacles in the channels, the effect of disturbance could
be revealed. Ideally, the disruption to flow velocity field alters the flow direction to create
lateral convection, and increases the contact area between the confluent fluid streams, and
hence, enhance mixing.
116
Chapter 6
Optimizing Layout of Cylindrical

Obstacles in Microchannels for
Enhanced Mixing
6.1 Introduction
In the previous chapter, it was found that the geometric topology of a microchannel
could create useful disturbance for microfluidic mixing. This chapter focuses on several
aspects of passive micro-mixer design and the optimization of layout of cylindrical obsta-
TM
cles. MemCFD was used as the principal tool to study the mixing of two liquids in
a Y-channel (or T-channel) with complementary experimental verification of mixing per-
formance. Obstacle structures used in this study were totally compatible with the current
micromachining techniques and were realized by Excimer laser.
It is emphasised here that obstacles in microchannels do not generate turbulence when

the Reynolds numbers are much lower than the threshold number for a turbulent flow,
117
however, the obstacles still disrupt the fluid to create lateral bulk mass transport to improve
the mixing.
In this chapter, the simulations of cylindrical obstacles inside a microchannel are first
described, followed by the experimental validation. In section 6.2, the detailed numerical
models and the method to calculate mixing efficiency from mass fraction distribution are
introduced. In section 6.3, the simulation results and experimental validation are discussed.
Section 6.4 summarises the outcomes of Chapter 6.
6.2 Numerical Modelling of Mixing in a Microchannel with

Cylindrical Obstacles
The obstacle structure in the T-type channel is illustrated in Figure 6.1. The T-type
channel has two inlet and one outlet ports. The width of the inlets is 200 µm, the width of
the outlet channel is 300 µm, and the height of the channel is 100 µm. In order to reduce
computing time, we used a channel length of only 1.2mm, except for the two obstacle
arrays (configuration no.8) with a length of 2 mm. The diameter of the obstacles is 60 µm.
Table 6.1 gives the number of the obstacles in eight different designs. Figure 6.2 illustrates
the layout of design no.6, no.7 and no.8, and the spacing for other designs are the same.
Configuration number No.1 No. 2 No. 3 No. 4 No.5 No.6 No.7 No.8
Number of obstacles 0 1 1 2 3 9 9 18
Table 6.1: Configuration of microchannel with Cylindrical obstacles
Using the numerical modelling techniques described in section 3.3.3, the geometric
models could be discretized into finite element meshes (Figure 6.3). Boundary conditions
for flow at the walls of the microchannel and walls of obstacles were set to zero velocity,
due to the nature of viscous flow in a microchannel. Volumetric flow rate boundary con-
118
Fluid 2
Fluid 1 Cylindrical
obstacles Mixture
Figure 6.1: 3D model of a T-channel with two arrays of cylindrical obstacles
φ φ
(a). (b).
(c).
Figure 6.2: Layout of square and triangular configuration of obstacle array. (a). configura-
tion No.6. (b). configuration No.7 (c). configuration No.8
119
ditions, q = 2 × 107 µm3 /sec, were assigned to both two inlets, which is equivalent to
Reynolds number, Re = 0.4.
100 µm
Figure 6.3: Meshed fluid volume in a T-type channel with two arrays of cylindrical
obstacles
TM
Simulations were performed using MemCFD on Win NT4.0 with Pentium III 800
MHz CPU and 128MB memory. It was non-trivial to apply uniform structured mesh with
multiple cylindrical obstacles, and to improve the simulation quality, parabolic elements
were used instead of linear elements. The finite element solver required considerable com-
puter memory, and even for modest size 3D CFD models, it consumed the memory rapidly
120
and used the hard disk as virtual memory. This made the simulation very slow. For the
series of simulations to be run, 3D finite element analysis was very time consuming. How-
ever, the structures of the microchannels were planar with low aspect ratio, and diffusion is
rapid in the depth direction when it is a slow flow. In addition, the velocity and the deriva-
tion of velocity in the depth direction was zero. Hence, the consideration of dispersion in
the depth direction could be relaxed for these reasons. Therefore, it was possible to ap-
ply 2D finite analysis to achieve accurate simulation without losing the focus on the major
problem. The accuracy of the simulation will depend on the height to width aspect ratio.
For a slow flow in a low aspect ratio channel, the diffusion is rapid across the streamlines
in the depth direction (z-axis in Figure 4.1) of the fluids due to the parabolic shape (see
Taylor dispersion in section 4.3). The mixing problem would be on the interfacial region
between two fluid streams (Figure 4.1). The width to depth aspect ratio of the channel was
1/3 for the purpose of demonstrating the concept. The sample fluids used in the simulation
were water and ethanol (Table 6.2) . Steady flow with a Peclet number of 200 was assumed
in the simulation.
Fluids Viscosity (kg µm−1 s−1 ) Diffusivity (µm2 s−1 ) Density (kg µm−3 )
Water 9.0e-10 1.2e3 9.998e-16
Ethanol 1.2e-9 1.2e3 7.89e-16
Table 6.2: Properties of water and ethanol at 20 ˚ C
6.3 Results and Discussions
Using the numerical models described above, CFD simulations were carried out. Se-
lected experiments provided qualitative verification and also presented in this section.
121
6.3.1 Simulation Results
The simulation results demonstrated the evolution of the design layouts (Figures 6.2)
using cylindrical obstacles in the microchannels. In Figure 6.4, the red colour represented
one fluid (Ethanol) and blue colour represented another fluid (water) and the green colour
was for the mixed solution. The wider the green area at the outlet, the better the mixing.
Figure 6.4 shows the simulation result of relative mass fraction distribution with cylin-
drical obstacles inside microchannels. Figure 6.5 shows the line plots of mass fraction
distribution at different position along the microchannel width at location Y-Y.
Figure 6.4: Relative mass fraction distribution in solution in a microchannel with two arrays
of cylindrical obstacles.
The mass fraction distribution illustrated in Figure 6.5 indicates the degree of mixing.
However, to quantitatively evaluate the performance of mixing, the mixing efficiency can
122
100.0%
90.0%
mass fraction (%)

80.0%
70.0%
60.0%
50.0%
40.0%
30.0%
20.0%
10.0%
0.0%
0 50 100 150 200 250 300
channel y - axis (µm)
Figure 6.5: Relative mass fraction distribution in solution in a microchannel with two arrays
of cylindrical obstacles.
be integrated from Figure 6.5 and calculated by equation 6.1 (Jeon, Dertinger, Chiu, Choi,
Stroock & Whitesides 2000),
 Rw 
|c − c∞ |dx
 
mef f = 0
1 − Rw
 × 100%,
 (6.1)
|c0 − c∞ |dx
0
where c is the mass concentration distribution across the transverse direction at the out-
let, c∞ is the concentration of a complete mixing, and c0 is the initial distribution of the
concentration before any mixing.
Because of the disruption of obstacles in a microchannel, the pressure drop is affected

by the blockage of obstacles. The pressure drop is a crucial factor to the design of a mi-
crofluidic device, therefore, the overall performance of microfluidic mixer should include
the evaluation of pressure drop. This overall mixing performance is termed mixing index,
midx , and used in this research to evaluate the overall performance of a passive microfluidic
mixer. The mixing index is defined in Equation 6.2.
∆pmax
midx = mef f × (6.2)
∆p
123
where mef f is mixing efficiency, ∆p is pressure drop between outlet and inlet of a mi-
cromixer and ∆pmax is the maximum pressure drop in the same set of simulation.
Figure 6.6 shows the evolution of the design with cylindrical obstacles for eight con-
figurations.
1 2 3 4
5 6 7 8
Figure 6.6: Evolution of design with obstacles in the Y-channel for configurations 1 to 8
The simulation results of mixing efficiency and pressure drop in the Y-channel were
calculated and plotted in Figure 6.7a, and the mixing index is shown in Figure 6.7b.
The simulation results (Figures 6.6 and 6.7) show that mixing efficiency increases as
the number of obstacles increases. The increased mixing can be interpreted as the obstacles
causing the parabolic velocity distribution to be evened out and give more time for diffusion
at the interface of the two liquids. However, the improvement of mixing performance is not
proportional to the increased pressure drop or the number of obstacles. Configurations no.2
and 3 have approximately the same pressure drop, but no.3 has much more mixing. Con-
figurations no.6 and 7 have the same number of obstacles, but no.7 has lower pressure drop
124
60.0% 1.6E-03
8
Pressure Drop (MPa)

7 8 1.4E-03
50.0%
Mixing Efficiency (%)
4 5 1.2E-03
40.0% 3 6
1.0E-03
6
30.0% 8.0E-04
7
2 6.0E-04
20.0%
1 4.0E-04
10.0%
4 5 2.0E-04
1 2 3
0.0% 0.0E+00
Number of obstacles
mixing efficiency Pressure Drop

a
5.0
2
4.0 4
3.0
midx
5
3
2.0 1
1.0 7
6 8
0.0
0 2 4 6 8 10 12 14 16 18
Number of obstacles
b
Figure 6.7: Finite element simulation results of (a) mixing efficiency, pressure drop be-
tween inlet and outlet versus number of obstacles, and (b) Mixing index versus number of
obstacles(configuration no.1 ∼ no.8).
125
and higher mixing efficiency. These results show that the asymmetric layout of cylindrical
obstacles were favouring mixing, but symmetric ones were not.
Figure 6.7b shows the mixing index (defined in equation 6.2), and the optimized con-
figuration would be one cylindrical obstacle with asymmetric layout (configuration no.3).
The criterion for designing a passive micro-mixer with cylindrical obstacles should com-
prise the capability of a micro-pump to overcome pressure drop with an adequate mixing
efficiency. Therefore, the mixing efficiency has to compromise the pressure drop to opti-
mize the design of a micro-mixer.
For a given flow rate, the total time for the flow through the channel remains the same,
because the average velocity at inlets and outlet remains the same for incompressible fluids
(conservation of mass). The local velocity, where the obstacles are, is increased to keep the
constant flow rate. The asymmetric layout of obstacles gives different resistance to the flow
in the lateral direction, so the fluids find their path through the lower resistance area, similar
to an electric circuit. This means that part of the fluid flow is distorted with the redirection
of the flow and this convective effect is shown in Figure 6.8. The lateral convection created
bulk mass transport other than diffusion, and this was the reason that mixing relied largely
on the lateral convection generated by obstacles with asymmetric configurations.
Figure 6.8: Section of the channel illustrated the convective effect using velocity vector
field
126
6.3.2 Experimental Results
To validate the simulation results by experiment, Excimer laser machined microchan-

nels with cylindrical obstacles were filled with miscible liquids to visualize the mixing
phenomena (Figure 6.9). Figure 6.10 shows that mixing performance was poor when there
were no obstacles in the channel. Figure 6.11 shows that the asymmetric arrangement im-
proved mixing performance. In Figure 6.11, the fluids were Ferric Ammonium Sulfate
and Ammonium Thiocyanate. The apparent diffusion was faster than it should have been
due to chemical reaction. However, there was still an obvious broadening of mixing re-
gion. Food colouring aqueous solutions had one-third diffusion coefficient compared to
the diffusion coefficient of water and ethanol. Nevertheless, Figure 6.12 shows that the
experimental visualization of the mixed area (green region) using food colouring dyes is
broadened compared to that of the two configurations in Figures 6.10 and 6.11.
200 µm
15kvolts, x180, 04/03/2002
Figure 6.9: Channel with cylindrical obstacles (SEM image)
Although quantitative measurement was not possible in these experiments, by com-

paring the mass fraction distribution between simulation results and experimental results,
the visualization provided qualitative agreement with the simulation results. Therefore, the
127
Figure 6.10: Configuration no.1, channel without cylindrical obstacles (food colouring),
Re = 0.4
Figure 6.11: Configuration no.5, channel with 3 cylindrical obstacles (Ferric Ammonium
Sulfate and Ammonium Thiocyanate solutions), Re= 0.4
Figure 6.12: Configuration no.8, channel with 18 cylindrical obstacles (food colouring),
Re = 0.4
128
MemCFD simulation software was demonstrated to be an appropriate method to evaluate

the performance of passive microfluidic mixer designs.
6.4 Summary
With viscosity dominating flow in microchannels, mixing of two fluid streams mainly
depends on diffusion. As a simplification of the ramping structure, obstacles were placed
in the microchannels, and mixing performance was investigated by simulation and experi-
ment. Obstacles in a microchannel with a low Reynolds number cannot generate eddies or
turbulence. However, the results demonstrated that obstacles could improve mixing perfor-
mance by affecting the flow pattern, and that the asymmetric arrangement of obstacles alter
the flow directions and force one fluid to flow into another to create lateral mass transport.
In this chapter, this phenomenon was termed convective effect in order to differentiate it
from pure molecular diffusion.
A very important outcome from the studies was that asymmetric layout of obstacles
could enhance microfluidic mixing, while symmetric layouts had little effect on mixing.
Placing obstacles in the microchannels was a novel method for mixing in microfluidic de-
vices, and properly designed geometric parameters, such as layout and number of obstacles,
could improve the mixing performance without significantly increasing the pressure drop.
The use of cylindrical obstacles provided only a limited number of parameters for
the study. For a more comprehensive study on obstacles, more design parameters for this
type micro-mixer needed to be investigated. Therefore, cylindrical obstacles would be
replaced by rectangular obstacles to study their effects on mixing in Chapter 7. Rectangular
shaped obstacles provide more design parameters than cylindrical shaped obstacles do, and
these parameters are useful for optimizing the design of passive microfluidic mixers using
geometric variations.
129
Chapter 7
Enhancing Mixing in Microchannels

Using Square and Rectangular Obstacles
7.1 Introduction
The previous chapter revealed that although cylindrical obstacles in microchannel

could broaden the mixed region, more detailed parameters needed to be examined. To
continue the research on how structures in microchannels affect mixing, passive mixing by
rectangular obstacles was studied and reported in this chapter. One reason to use a rectan-
gular shape was the number of parameters that could be investigated, such as width, height,
gap between obstacles, position in the channel and obstacle orientation. In addition, rect-
angular structures were much easier to mesh into uniform structured elements to achieve
high simulation quality with reduced computing time.
In respect to the nature of laminar flow in a microchannel, the geometric variations

were designed to try to improve the lateral convection. By doing this, the dispersion of
one fluid into another was not only assisted by diffusion, but also, and more importantly,
130
by the convection in the lateral direction. Geometric parameters versus the mixing perfor-
mance were investigated systematically in T-type channels by applying MemCFD solvers
for microfluidics. Various obstacle shapes, sizes and layouts were studied, and the mixing
performance of microchannels with obstacles were evaluated by mass fraction.
In practice, the gap, size and other parameters associated with obstacles were not sep-
arated. In this chapter, the author provides linkages between these parameters and mixing
performance.
7.2 Numerical Modelling of Square and Rectangular Ob-

stacles
The channels were meshed into 8-nodes hexahedral elements. The simulation was
run as 3D, steady, laminar, Newtonian, with two fluids to evaluate the dispersion of one
fluid into another. The inlets were assigned flow rate (ranging from 2× 107 µm3 s−1 to 2×
109 µm3 s−1 in this research) boundary condition, and all the channel walls were assigned
wall boundary condition (velocity components in Cartesian coordinates: vx = vy = vz =
0). Two fluids entered the two inlets respectively (red dye aqueous solution and blue dye
aqueous solution), with a diffusion coefficient, D = 4×10−6 cm2 s−1 , which was calculated
according to their chemical groups (Reid et al. 1987).
For all the studies in this chapter, the T-type channels were used, with two inlets and
one outlet. In Figure 7.1, the T-channel with four rectangular obstacles is illustrated. The
length, width and depth of the main channels were 5mm, 200 µm and 100µm respectively.
The length of the inlet, L1 , is not critical but should be sufficient for full developed flow,
which requires L1 À 0.01wRe (Bird et al. 1960). The main channel length is divided
into three parts, where L0 is the entry channel length for the full development of flow, Lr
is length of the section to create lateral convection and Lm is the channel length for the
131
broadening of solution by diffusion.

Fluid 2
a
Mixture
Fluid 1
b
L0 Lr Lm
g g g
L1
w
δl
Figure 7.1: T-channel with 4 rectangular obstacles, (a) 3D view; (b) layout.
The simulation results in Figure 7.2 show how the obstacles affected mixing. To eval-
uate the performance of a micromixer, several important factors need to be considered:
mixing efficiency (Equation 6.1), pressure drop and overall mixing performance, which
was termed mixing index and defined in Equation 6.2.
In this section, the effect of number, size, angle and other important parameters of
132
1 2 3 4
5 6 7 8
a
9 10 11 12
13 14 15
b
Figure 7.2: Mass fraction distributions, (a) symmetric layout of obstacles, (b) asymmetric
layout of obstacles
133
obstacles to mixing and pressure drop were evaluated using mixing efficiency, pressure
drop and mixing index respectively.
7.3.1 Number of Obstacles
In Chapter 6, it was found that symmetric placement of obstacles does not improve
mixing, whereas asymmetric layout of obstacles can enhance mixing, and the mixing effi-
ciency increases with the number of obstacles. By repeating the simulations using square
obstacles, Figure 7.2 shows similar results. The mass fractions were extracted close to the
end of the channel and the pressure drops were calculated between the inlet and outlet.
Mass fractions and pressure drops were used to compute the mixing efficiency and mixing
index as defined in Section 6.3.1. Figure 7.3 shows the mixing efficiency, and it clearly
indicates that symmetric layouts have little effect on enhancing mixing and asymmetric
layouts can improve mixing efficiency.
50.00%
45.00%
mixing efficiency (%)
40.00%
35.00%
30.00%
asymmetric
25.00%
symmetric
20.00%
15.00%
10.00%
5.00%
0.00%
0 2 4 6 8 10 12 14 16 18 20
number of obstacles
Figure 7.3: Mixing efficiency versus number of square obstacles
As explained in the previous chapter, pressure drop is a critical parameter for a mi-
crofluidic device. Figure 7.4 shows the pressure drop versus the number of obstacles, and
134
Figure 7.5 shows the mixing efficiency measured by mixing index, which was defined as
the mixing efficiency divided by the normalized pressure drop (Section 6.3.1). Although
the mixing efficiency increased with the increase of number of obstacles, the improvement
was not significant when obstacles reached a certain number, because the increase in num-
ber of obstacles also blocks the lateral convection. The pressure drop in Figures 7.4 and
7.5 shows that the asymmetric configuration with six obstacles gave optimized overall per-
formance. In Figure 7.3, the mixing efficiency with two obstacles is better than three
2.5E-05
Pressure drop (MPa)
2.0E-05
1.5E-05
1.0E-05
5.0E-06
0.0E+00
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
number of obstacles
asymmetric symmetric
Figure 7.4: Pressure drop versus number of square obstacles
obstacles. It does not generally agree with the view that the mixing efficiency improves
with the number of obstacles. Therefore, a further study of these two configurations was
necessary. Figure 7.6 uses the angle between the flow axis and obstacle axis.
The further simulation results showed that two obstacles inside a microchannel had a
better mixing efficiency than that of three obstacles (Figure 7.7). This may be interpreted as
the mid obstacle actually reducing the convective effect and partially blocking the diffusion
process at the interfacial area between two fluid streams. Figure 7.8 shows that there is no
significant difference on pressure drop. Therefore, the mixing index shown in Figure 7.9
135
1.4
1.2
1.0
midx
0.8
0.6
0.4
0.2
0.0
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
number of obstacles
asymmetric symmetric
Figure 7.5: Mixing index, midx , versus number of square obstacles
had the same trend as in Figure 7.7, and two obstacles configuration demonstrated better
overall performance than three obstacles.
7.3.2 Size of Obstacles
The width of the obstacles was found to be a more important dimension than length,
because width was the most effective length scale to influence the flow patterns (Figure
7.10). The configuration and the simulation results are illustrated in Figure 7.11. Two
streams of fluids flowed from left to right along the x-axis. All the parameters including
flow rate were fixed, except the width of the obstacles. The height of the obstacles equaled
the depth of the channel. By changing the width of the obstacles, this changed the degree
of blockage to the flow. The velocity component lateral to the flow increased, and so did
mixing efficiency, when the size of obstacles became bigger. However, at a certain width
(80% of blockage in Figure 7.11) when the gaps between obstacles and wall equaled each
other, the flow pattern became symmetric. As shown previously, the mixing efficiency
136
Figure 7.6: Angle between flow axis and obstacle axis, two and three obstacles
137
30.00%
Mixing efficiency (%)
20.00%
2dots
3dots
10.00%
0.00%
0 20 40 60 80
Angle (degree)
Figure 7.7: Mixing efficiency with two and three obstacles inside a microchannel versus
angle of layout
8.0E-06
Pressure drop (MPa)
6.0E-06
4.0E-06 2dots
3dots
2.0E-06
0.0E+00
0 20 40 60 80
Angle (degree)
Figure 7.8: Pressure drop versus two and three square obstacles
138
0.5
0.4
0.3
midx
2dots
0.2
3dots
0.1
0.0
0 20 40 60 80
Angle (degree)
Figure 7.9: Mixing index, midx , versus two and three square obstacles
Figure 7.10: Mass fraction distribution for two rectangular obstacle layouts with different
size of obstacles
139
60.0%

50.0%
40.0%
30.0%
20.0%
10.0% 1
0.0%
0.0% 20.0% 40.0% 60.0% 80.0% 100.0%
degree of blockage (b/w) (%)
Figure 7.11: Mixing efficiency versus size of obstacles
was poor for any symmetric arrangement. When the width of the obstacles continued to
increase until it reached the wall, the configuration was similar to the serpentine shaped
microchannels (Liu et al. 2000).
Larger size obstacles block the flow more (Figure 7.12), and Figure 7.13 shows the
overall performance was maximum around the 40% blockage, which meant the size of
obstacle was 40% of the channel width.
7.3.3 Angle of Obstacles
From the configuration described in section 7.3.2, we rotate the obstacles against their
center of mass simultaneously (Figure 7.14). The orientation of the obstacles causes a dif-
ferent type of disturbance to the flow (Figure 7.15). The results showed that the maximum
mixing efficiency was located in the range from 60o to 75o and 105o to 120o .
The lateral projection area of obstacle decided the resistance to the flow, therefore
the obstacles perpendicular to the channel wall presented the highest pressure drop (Figure
140
1.0E-02
Pressure drop (MPa) 1.0E-03
1.0E-04
1.0E-05
1.0E-06
0.0% 20.0% 40.0% 60.0% 80.0% 100.0%
degree of blockage (%)
Figure 7.12: Pressure drop versus size of obstacles
140.0
120.0
100.0
80.0
midx
60.0
40.0
20.0
0.0
0.0% 20.0% 40.0% 60.0% 80.0% 100.0%
degree of blockage (%)
Figure 7.13: Mixing index, midx , versus size of obstacles
7.16). Figure 7.17 shows the optimized overall performance is around 60o and 120o to 150o .
141
Figure 7.14: Mass fraction distribution, various angles of obstacle layouts
60.0%
50.0%
40.0%
30.0%
20.0%
10.0%
0.0%
0 50 100 150 200
Angle (degree)
Figure 7.15: Mixing efficiency versus angle of rectangular obstacles

142
1.5E-04
1.0E-04
5.0E-05
0.0E+00
0 50 100 150 200
Angle (degree)
Figure 7.16: Pressure drop versus angle of obstacles
0.8
0.7
0.6
midx
0.5
0.4
0.3
0.2
0 50 100 150 200
Angle (degree)
Figure 7.17: Mixing index, midx , versus angle of obstacles

143
7.3.4 Obstacle Offsets from Channel Walls
When two or more rectangular obstacles are placed in a microchannel, the relative
position of the obstacles to each other and the channel walls needs to be studied. In this
section, the parameter investigated was measured as a to b ratio and ten obstacles were
placed in the channel (Figure 7.18). In Figure 7.19, when a/b = 1, it indicated that the
obstacles lay in the center of the channel symmetrically, which did not improve the mixing
performance. When a/b = 0, the obstacles were connected to the wall, which actually
formed a serpentine-shaped channel. However, it was preferable to have a gap between the
obstacles and the wall, as it could eliminate the dead volume created at corners between
obstacles and channel walls.
a b
Figure 7.18: Mass fraction distribution, for different offsets of obstacles from channel
walls, (a). no offset from the center; (b). with offset.
When there was no space between the obstacles and channel walls, this gave the best
mixing efficiency, however, this raised the potential for dead volume at the corners between
obstacles and channel walls. An offset of obstacles from the channel walls removes dead
volumes. Figure 7.20 shows there is no significant difference in pressure drop, and Figure
7.21 demonstrates that obstacles closer to the wall give better overall performance of the
144
a
b
120.00% …
100.00% b

a
80.00%
60.00%
40.00%
20.00%
0.00%
0 0.2 0.4 0.6 0.8 1 1.2
a/b ratio
Figure 7.19: Mixing efficiency versus offset from the wall

Pressure drop (MPa)
2.0E-04
1.5E-04
1.0E-04
5.0E-05
0.0E+00
0 0.2 0.4 0.6 0.8 1
a/b ratio
Figure 7.20: Pressure drop versus offset of obstacles from wall

145
1.2
1.0
0.8
midx
0.6
0.4
0.2
0.0
0 0.2 0.4 0.6 0.8 1
a/b ratio
Figure 7.21: Mixing index, midx , versus offset of obstacles from wall
device. A space between obstacle and the channel wall removes the dead volume, and the
a to b ratio between 0.1 to 0.6 gave the optimized results.
7.3.5 Gap between Obstacles
Figures 7.22 and 7.23 show the results of mixing of two fluids for two obstacle gap
arrangements. The green colour is the complete mixture. The smaller gap shows less
mixing for the same overall length. The mixing efficiency is expected to increase with the
reduction of gap between the obstacles, due to the diffusion distance being reduced by the
gap. However, a smaller gap tends to offer greater resistance to flow, and hence, reduces the
flow in the gap and increases the flow through the upper and lower spaces. Therefore, this
may mean a better mixed fluid in the smaller gap, but less of it, and the overall contribution
to mixing at the outlet by the fluid flowing through the smaller gap is then offset by reduced
flow.
Because the lateral projected area of obstacles was the same, the pressure drop re-
mained the same for the different configurations (Figure 7.24). Figure 7.25 demonstrates
146
a b
Figure 7.22: Mass fraction distribution for different obstacle gaps, (a). 160 µm gap; (b).
400µm gap;
70.0%
60.0%
50.0%
40.0%
30.0%
20.0%
10.0%
0.0%
0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
Gap to width ratio (g/w)
Figure 7.23: Mixing efficiency versus gap between obstacles

147
the overall performance was improved with the increase of gap. However, the improvement
rate slowed as the gap approached the same scale as the width of the channel.
1.0E-04
1.0E-06
1.0E-07
100 200 300 400
gap (µ m)
Figure 7.24: Pressure drop versus gap between two obstacles
0.7
0.6
0.5
0.4
midx
0.3
0.2
0.1
0.0
100 200 300 400
gap (µ m)
Figure 7.25: Mixing index, midx ,versus gap, g, between two obstacles
148
7.3.6 Flow Rate
When a certain flow rate (Re ∼ 1) was reached, the mixing efficiency stopped de-
creasing (Figure 7.26). For the optimized obstacle structures, which have the maximum
mixing index, i.e. the structures described in Figure 7.13 with 40% blockage, the mixing
efficiency increased with the flow rate: 2 to 1000nl/sec (Re = 0.005 ∼ 2.5). As there was
no chaotic effects observed from the flow pattern, the interpretation of this would be that
there was an increase of convective components in the lateral direction to the flow.
The pressure drop would be expected to increase with Reynolds number (Figure 7.27).
However, the slower flow improved the overall mixing performance more, which in turn
proved that slow fluid motion in a microfluidic device was preferred (Figure 7.28).
80%
70%
60%
50%
40%
0 0.5 1 1.5 2 2.5 3
Reynolds number
Figure 7.26: Mixing efficiency versus flow rates

149
4.0E-04
2.0E-04
1.0E-04
0.0E+00
0 0.5 1 1.5 2 2.5 3
Reynolds number
Figure 7.27: Pressure drop versus Reynolds number
1000.0
100.0
midx
10.0
1.0
0.1
0 0.5 1 1.5 2 2.5 3
Reynolds number
Figure 7.28: Mixing index, midx , versus Reynolds number, Re

150
7.3.7 Ratio of Height of Obstacles to Depth of Channel
Figure 7.29 illustrated that obstacles with different heights in microchannels had an
influence on mixing performance. Higher obstacles created stronger 3D convection (Figure
7.30). The height, h, of obstacles was one variable and width, b = 280µm, was another.
Other parameters were the same as in Figure 7.11. The most useful disturbance to improve
a b
Figure 7.29: Mass fraction versus obstacle height to channel depth aspect ratio, h/H, (a).
h/H = 0.2, (b). h/H = 0.9
mixing is that caused by obstacles perpendicular to the interface of two fluid streams. The
full height obstacle proved to be most efficient.
When the ratio of obstacle height to channel depth was zero, the channel was with-
out any obstacles inside, and the mixing efficiency was obviously the lowest due to pure
diffusion based mixing mechanism. However, when this ratio became ’negative’, which
indicated that the obstacles became grooves, the mixing was improved again. Mixing in
microchannels having slantly patterned grooves are reported in Chapter 8. Figure 7.31
shows that the pressure drop increases quickly with height of the obstacles, and Figure
7.32 shows that the negative height, which illustrates grooved channel, gives the best over-
151
100.0%
80.0%
60.0%
40.0%
20.0%
0.0%
-0.2 0 0.2 0.4 0.6 0.8 1
Aspect ratio (h/H)
Figure 7.30: Mixing efficiency versus height to depth aspect ratio of obstacles
2.0E-04
Pressure drop (MPa)
1.5E-04
1.0E-04
5.0E-05
0.0E+00
-0.2 0 0.2 0.4 0.6 0.8 1
Aspect ratio (h/H)
Figure 7.31: Pressure drop versus obstacle height to channel depth ratio
152
500.0
400.0
300.0
midx
200.0
100.0
0.0
-0.2 0 0.2 0.4 0.6 0.8 1
Aspect ratio (h/H)
Figure 7.32: Mixing index, midx , versus obstacle height to channel depth ratio, h/H
all performance. This means that when pressure drop is a restriction to the design of a
micromixer, using grooved structures in the microchannel can provide an alternative, in-
stead of using obstacles. This provided a very important design guideline for micromixers.
7.4 Summary
In this chapter, square and rectangular obstacle structures in a microchannel were used
to improve mixing performance. Geometric parameters, such as number, size, gap between
obstacles, position, orientation in channel, and height ratio of obstacles, were studied for
their effect on mixing. The investigation of how flow rate affected mixing performance
illustrated promising results, and mixing efficiency remained constant or even increased
with flow rate when Re > 1. This suggested that convective effects dominated mixing over
diffusion, and diffusion over a short distance is negligible. The results were used to identify
feasible strategies for improving mixing by placing obstacles in microchannels.
When the height of obstacles to depth of channel ratio became ‘negative’, the obstacles
153
became grooves patterned on the channel wall. These grooves provided slip-like bound-
ary conditions, because of the lower resistance to flow due to fluid to fluid friction other
than fluid to wall friction. Therefore, a lateral pressure gradient was built up and created
secondary flow. A detailed investigation of slantly patterned grooves inside a rectangular
microchannel are carried out in Chapter 8 using a different numerical approach, termed
particle tracing algorithm, to trace the flow trajectories in the microchannel.
154
Chapter 8
Numerical Investigation of Mixing in

Microchannels with Patterned Grooves
8.1 Introduction
In Chapter 7, square and rectangular obstacles were used to disrupt main flow to cre-
ate lateral convection, which can be also referred to as secondary flow. Stronger secondary
flow was found to improve mixing performance more. As with many other microfabri-
cated structures, the obstacles were planar, and the structure is sometimes referred to as
2.5-Dimension,which can be defined by height, width and length. The flow patterns were
mostly planar, with only two non-zero velocity components. When the height of rectan-
gular obstacles to channel depth ratio was reduced, the convective secondary flows were
not just lateral, but showed signs of of 3D convection, which indicated that all the three
velocity components were non-zero. When the height ratio became negative, which meant
the obstacles were actually decayed into grooves, 3D convection created by the grooves
could also improve mixing performance.
155
Several recently published papers investigated microfluidic mixers with slanted grooves
in microchannel walls by theoretical predictions and experimental measurement (Johnson
et al. 2002, Stroock et al. 2002). Unlike obstacles, the grooves created secondary flows
by providing local slip boundary conditions, and this type of mixer was declared a chaotic
microfluidic mixer (Stroock et al. 2002). However, the actual flow pattern in this type of
mixer was not clearly verified, and the grooves also provided splitting and recombining
effects to the main fluid streams. Therefore, grooves in microchannels needed to be further
studied to provide design guidelines for this type of micromixer.
In this chapter, the numerical model was critical for accurate simulations and is ad-
dressed first in section 8.2. Then, a set of particle tracing algorithms is developed and
presented in section 8.3. In section 8.4, the flow patterns are analyzed using the developed
particle tracing algorithms. Streaklines and Poincaré maps are computed to demonstrate
3D helical flow patterns. The splitting and folding fluid material lines by direct simulations
of two fluids diffusion-convection problems is also illustrated in this section. Section 8.5
summarizes the results and concludes that the results obtained in the chapter could be used
as guidelines for designing this type of micromixers.
8.2 Numerical Setup
The performance of microfluidic mixers with patterned grooves was evaluated using
several numerical approaches. A Computational Fluid Dynamics (CFD) package for Mi-
crofluidic applications was used to simulate the 3D velocity field as well as two fluids
mixing, and algorithms were developed to compute rate of shear and particle trajectories to
study the stretching and folding of fluids.
The 3D geometric models (Figure 8.1) were kept consistent with a paper by Stroock et
al. (2002a). In order to obtain high quality structured mesh elements, a CFD preprocessing
156
z inlet -1
y
x Outlet
H
L0 Lm W
Lr
inlet -2 Ltotal
Note: not to scale
Figure 8.1: Fluid volume in microchannel with grooved bottom surface. L0 is the entry
channel length, Lr is the grooved channel length and Lm is the channel length to complete
the mixing by diffusion
package, Gambitr v2.04, was used to build 3D models and 8-node Hexahedral elements
(Figure 8.2). It should be emphasized that high quality uniform mesh elements are crit-
ical to achieving accurate simulation results, especially for simulating two fluids mixing.
The uniformity of the mesh can be achieved by applying advanced mesh mapping or sub-
mapping techniques (covered in section 3.3). The density of the mesh elements was first
calibrated by comparing the simulation results with analytical solutions of viscous flow in
a rectangular channel. By this approach, the size of the element could be decided for mod-
elling channels with patterned grooves. Then, by mesh sensitivity study, the mesh density
can be increased or decreased until a stable solution is reached, which either can be de-
termined by comparing with an analytical solution or when the simulation reaches a point
where accuracy is insensitive to mesh density. Because the size of the numerical model
is directly related to mesh density, sometimes it is necessary to reduce the unnecessary
density of mesh. MemCFD v2001.3d from Coventor was used to solve the Navier-Stokes
equations and the diffusion-convection problem.
The simulations were run as steady, laminar, Newtonian, and with one fluid (or two
fluids in the case of solving diffusion-convection problems). In Figure 8.1, both inlet1
and inlet2 were assigned the same volumetric flow rate boundary conditions. In the case
of simulations using two fluids, 100% fluid1 was assigned to inlet1 and 100% fluid2 was
157
z
y
x
Figure 8.2: Meshed fluid volume
assign to inlet2 with the same volumetric flow rate value. In this section, the mean ve-
locities used ranged from 100 µm s−1 to 50000 µm s−1 , which could be worked out
from the flow rates. All the channel walls were assigned “WALL” boundary condition
(vi = 0, i = 1, 2, 3, represented three orthogonal velocity components). The diffusion co-
efficient, D, was 4×10−6 cm2 s−1 for simulating two fluids diffusion-convection flow. A set
of particle tracing algorithms was developed and coded in Fortran77 to compute the particle
trajectories and Poincaré maps using the 3D velocity field exported from the CFD simu-
lations (Appendix B.2). The length, width and depth of the channels were 5mm, 200µm
and 100µm respectively. Simulations were performed on Win NT4.0 with Pentium III
800MHz CPU and 512MB memory. The number of mesh elements (cells) in the models
ranged from 250,000 to 300,000.
158
8.3 A Particle Tracing Algorithm
To evaluate the performance of a mixer, one can place non-diffusive particles in the
velocity fields. To observe the advection of these particles, the convection of mass transport
can be evaluated. In experiments, measurement of velocity can be carried out by injecting
very fine coloured ink streams in the flow fields or by particle image velocimetry (PIV).
A numerical approach was developed that places virtual particles, which have no physical
properties, in the velocity fields. The trajectories of these particles can be computed by the
Lagrangian method, which records the spatial positions of the particles at each time step.
In a steady flow, particle trajectories can be integrated with the system of equations (8.1):
dp(t)
= v(p(t), t) (8.1)
dt
where p(t) is the particle position at time t, and v is the velocity field. Integrating equation
(8.1) yields: Z t+δt
p(t + δt) = p(t) + v(p(t), t)dt (8.2)
t
where δt is the time step.
The integral term on the right hand side can be evaluated numerically using a multi-
stage method or a multi-step method. Regardless of how it is solved, the end result is
a displacement, which when added to the current position, p(t), gives the new particle
location at time t + δt.
In a discrete velocity field, the velocity of a particle at certain position is interpolated

by the element containing this position. Hence, this element needs to be located first.
After the element is located, the particle velocity can be interpolated. Then, a time step is
assigned to integrate its new position. The procedure is repeated until the particle leaves
the flow domain. The algorithm is described in words below:
159
find cell containing initial position (point location)

while ( particle in grid )
determine velocity at current position (interpolation)
calculate new position (integration)
find cell containing new position (point location)
endwhile
8.3.1 Point location
Point location is required to find the cell that contains a specified point. For the com-
plex mesh of microchannels with grooved surface, it is difficult to locate particle position
in the 8-node hexahedral elements. So it is necessary to split hexahedral elements into 5 or
6 tetrahedral elements, and 5 tetrahedral elements are used commonly (in Figure 8.3). For
each tetrahedron, assume that we have a point p with coordinates p(x, y, z). The simplest
way to find the element into which point p falls is to perform a loop over all the tetrahe-
drons, evaluating their shape-functions, Ni , with respect to p(x, y, z) using Equation 8.3.
X X X X
x= Ni x i ; y = Ni y i ; z = Ni z i ; 1 = Ni (8.3)
i i i i
Tetrahedrization is only performed in the cells along the path of the line and the results
do not have to be stored. Each tetrahedral element has 4 shape-functions (i = 4). The above
equations can be written into a linear system (Equation (8.4)), and the shape functions can
be evaluated by the standard Gaussian elimination method.
x = X · N, −→ N = X−1 · xp (8.4)
min(Ni , 1 − Ni ) > 0, ∀i (8.5)
The criterion set forth in Equation (8.5) is used to determine whether a point lies within
160
7
6
8 5
3
2
4
1
Figure 8.3: Hexahedral cell decomposed into 5-Tetrahedrons
the confines of an element. All the shape functions must have positive values but less than
unity to determine a point within an element. Negative values or values larger than unity
mean that the point is outside the element. This can be interpreted geometrically as the
shape functions are the ratios of four volumes divided by point p(x, y, z) to the volume
of the tetrahedron (in Figure 8.4) (Li et al. 1999). The sub-volumes cannot be negative
or larger than the volume of the element, therefore, any negative value of Ni or 1 − Ni
means that the point p(x, y, z) is outside this tetrahedral element. Then, the point location
proceeds by advancing to, and crossing the respective face into the adjoining tetrahedron.
The worst violator of the four conditions is used to predict which tetrahedron to try next.
Even if the bounding tetrahedron is not found in the immediate neighbor, then by always
moving in the direction of the worst violator the search technique will converge upon the
correct cell.
8.3.2 Interpolate velocity
One of three techniques may be used for the spatial interpolation of velocity: phys-
ical space linear interpolation, volume weighted interpolation, and linear shape function
interpolation. All three are mathematically equivalent and produce identical interpolation
161
p(x,y,z)
N4
N1 N3
N2
Figure 8.4: Tetrahedron geometry. Ni , i = 1, 4. are the shape functions
functions (Kenwright & Lane 1996). The linear shape function was the most efficient
technique for this application because it reused the shape function Ni , i = 1, 4 computed
during point location. The linear shape function for spatial velocity interpolation is:
4
X
v= Ni × vi (8.6)
i=1
where, vi are the velocity vectors at the 4 vertices of the tetrahedron.
Integration
Many integration methods are shown in the literature, ranging from the simple first-
order Euler scheme to the fourth-order Runge-Kutta scheme or even higher-order methods,
162
applied with fixed or variable time steps. The fourth-order Runge-Kutta scheme becomes:
p(t + δt) = p(t) + (k1 + 2k2 + 2k3 + k4 )/6
k1 = δt · v(p(t), t)
k2 = δt · v(p(t) + k1 /2, t + δt/2) (8.7)
k3 = δt · v(p(t) + k2 /2, t + δt/2)
k4 = δt · v(p(t) + k3 /2, t + δt/2)
where ki , i = 1 ∼ 4 are the parameters decided by particle position p, velocity v, time t

and time step δt. Thus, using equation (8.7), a particle can be numerically integrated and
traced through the flow field. In this section, an adaptive variable time step fourth-order
Runge-Kutta scheme was used.
With the algorithms developed in the section (8.3), one arbitrary periodical flow sec-
tion (Figure 8.5) of the exported velocity field for a single fluid flow was extracted to be
used to compute the Poincaré map. Mixing of two fluid streams was also presented to
demonstrate the twisting and stretching of the interfacial area between them.
8.4.1 Flow patterns
One periodical flow was illustrated for vector planes in Figure 8.6, which corre-
sponded to surface of sections ( π2 i, i = 0 to 4) in Figure 8.5. The void areas in the vector
planes indicated the cross section of the solid parts of the grooves. The aspect ratio of
grooves (α) can be defined as half of the groove depth to channel depth (H). The flow pat-
terns for high aspect ratio (α = 0.30) grooves were more complex than that of low aspect
163
200
αH
H
0 π/2 π 3π /2 2π
sec tion pla nes
Figure 8.5: One periodic section of patterned grooves, with section planes labelled from
0 ∼ 2π, and particles advected from plane 0 to plane 2π
ratio (α = 0.05). The rotating like flow patterns can be called secondary flow. Therefore, a
higher aspect ratio created stronger secondary flow.
8.4.2 Streaklines
A streakline is defined as a line formed by the particles, which pass through a given
location in the flow field. A streakline can be made visible by injecting a dye into the fluid
at the given location. In a steady flow, the streaklines coincide with the particle traces, and
the streamlines are lines to which the velocity vectors are tangent at all points.
To visualize the effects of the mixing numerically, a particle may be traced along the
streakline. The algorithms to determine streaklines are given in Appendix B.3. In Figure
8.7, the set of streaklines are twisting like a helical shape. This indicates folding and
stretching of fluids, which favours mixing. For the patterned grooves with the periodic
configuration, the velocity field can be used repeatedly. CFD simulations could simulate
only a section of the channel length due to the limitation of computer speed and memory.
However, by reusing the periodic velocity field repeatedly, streaklines may, in principle, be
164
130
105
0 0 200
0
200
130
105
π/2 0 200
0
200
130
105
π 0 200
0
200
130
105
3π/2 0 200
0
200
130
105
2π 0 200
0
200
α = 0.05 α = 0.30
Figure 8.6: Cross section velocity vector planes, 0 ∼ 2π, Re = 5

165
120
Z (µm)
Fluid 2
60
Fluid 1
0
200 Direction of fluid
150 Start points
100 5000
4000
y (µm) 50 3000
2000
0 1000
0 x (µm)
Figure 8.7: Streaklines of two neighboring fluid particles from two fluid streams respec-
tively, α = 0.20, Re = 5
computed for a longer channel length.
The distances between neighboring fluid elements’ traces (Figure 8.7) were calculated
along the flow direction. A plot of the normal distance between the two traces along the
length of the channel showed that the distance between the traces varied randomly (Figure
8.8). The plot indicates that stretching and folding of the fluid material lines occurs.
8.4.3 Poincaré maps
The same algorithms to compute streaklines can be used to record Poincaré maps,
which are useful for evaluating mixing performance. To generate a Poincaré map in a
spatially periodic system, one or a series of passive particles are advected by the periodic
velocity field and pass through a series of periodic planes in the system. For the mixer
considered here the plane is located at the exit of each mixing segment, and each position of
the particle which hits this plane is recorded. The plane is called the Poincaré map. Regular
patterns in the Poincaré map indicate integrable (non-chaotic) behavior, while jumbled
166
80.0
Gap between two streaklines (µm)

60.0
40.0
20.0
0.0
3000 3500 4000 4500 5000 5500 6000
Distance in flow direction (µm)
Figure 8.8: Gap between the two Streaklines in Figure 8.7
patterns may indicate the existence of chaotic behavior. In Figure 8.5, a spatially periodic
system between section plane 0 and plane 2π was clipped from the exported velocity field,
and the Poincaré maps were actually recorded on plane 2π.
In Figure 8.9a, a small aspect ratio (α = 0.05) had a regular recirculation pattern,
which implied that it was not chaotic. With high aspect ratio (α = 0.30), the flow pattern
illustrated in Figures 8.9b and 8.9c became more irregular. However, there was still not
sufficient evidence to verify it was chaotic advection. Nevertheless, the particles trajectories
were circling around the flow axis. Hence, by simply counting the dots per circle (illustrated
with the dark lines in Figure 8.9) in the Poincaré map, the necessary length of the grooved
channel to complete one recirculation can be calculated. The folding and stretching of
material lines is closely related to the recirculations. In this aspect, the necessary length
for one circulation may be used as a criterion to evaluate the performance of micromixers
with patterned grooves. In Figure 8.10, these lengths had an exponential relation with the
groove aspect ratio. When α = 0.3, the necessary lengths to complete one circulation were
also calculated for different Reynolds numbers, Re = 5 and 0.01 (Figures 8.9b and 8.9c).
The results were the same and this indicated that the recirculation had little to do with flow
velocity, but mostly was a function of geometric parameters, especially the aspect ratio of
167
Figure 8.9: Poincaré map, 45o patterned grooves: (a) α = 0.05, Re = 5; (b) α = 0.30, Re =
5; (c) α = 0.30, Re = 0.01.
168
grooves.
Length of grooved 102

channel (mm)
101
100
10-1
0.1 0.2 0.3
Groove Aspect Ratio (D)
Figure 8.10: Length of grooved channel to complete one circulation, for H/w = 0.5
Moreover, the recirculations can reorient fluids layers, hence, the reduction of the
diffusion path can be more effective in the depth direction, because most of microchannels
have low depth to width ratios due to the planar nature of micro-devices. In this aspect,
shallow channels with patterned grooves can achieve rapid mixing. In addition, the length
of channel required to achieve complete mixing is a criterion to evaluate the performance
of a passive microfluidic mixer. In Figure 8.1, Lr is a function of groove aspect ratio, and
Lm is proportional to the mean flow velocity. Therefore, when designing a micromixer,
instead of increasing the velocity to create chaotic advection to improve mixing, the use
of wide shallow channels with patterned grooves can offer an alternative solution to mix
fluids in a constrained channel length.
169
8.4.4 Simulation of Mixing in Microchannels with Patterned Grooves
In the previous sections, the particle tracing technique was used to optimize the design
of micromixers with patterned grooves. In this section, MemCFD was used as an alternative
approach to evaluate mixing performance.
To initiate the simulations, the geometric model and boundary conditions were set
in accordance with the description in section 8.2. Two fluids were assigned to two inlets
respectively to simulate diffusion-convection problems, with ū = 50000µm/s (Re = 5),
α = 0.1 and D = 4 × 10−6 cm2 /s. Due to anisotropic patterned grooves, the interfacial
area between two fluid streams was shifted from the flow axis and stretched. The shift
of interfacial area and relative mass fraction (percentage of fluid1, or fluid2, distributed
spatially in the channel) are shown in Figure 8.11. In the figure, both patterns were derived
from the same simulation result. In Figure 8.11 (a), only a thin interfacial layer between
two fluids is illustrated in a 3D view, which shows that the area of this layer is enlarged by
stretching and twisting. Figure 8.11b shows the mass fraction pattern. The two fluids were
represented by two colours, the green area (or bright area in a grey print) indicated where
the fluids became mixed. Figure 8.12 shows a similar simulation run for groove aspect
α = 0.30 at Re = 5. The twisting effect was much stronger and the interfacial regions
between two fluid streams were enlarged by the stretching, and even torn apart.
The shifting of fluid streams could be analyzed by calculating the rate of shear just
below the flat plate (Stroock et al. 2002). A measure of the helicity of the flow can be
made by comparing the transverse and longitudinal rate of shear. In this section, the rate of
shear was computed from the 3D velocity field exported from the CFD simulations. The
measurement of helicity can be provided by the angle Ω between the channel axis x and
interfacial line of two fluid streams, which is shifted by the helical flow pattern (Figure
8.11).
In Figure 8.13, the mean helicity computed numerically in this research is compared
170
100% fluid 1 Mixed 100% fluid 2
Figure 8.11: Two fluid streams flow into a T-type channel with patterned grooves, α = 0.10
and Re = 5. (a). The interfacial line is indicated by the bright mixed region. The helicity
was measured by the angle Ω between this interfacial line and the channel axis, (b). The
visualization of mixing
171
Figure 8.12: Interfacial regions when α = 0.30 and Re = 5
16.0
Rate of Shear (tan( Ω )x100)
12.0
8.0
4.0
0.0
0 0.1 0.2 0.3
Aspect ratio, α
analytical prediction numerical prediction measured
Figure 8.13: Mean helicity measured by the angle Ω for H/w = 0.5
172
with published analytical predictions and three experimental measured points from Stroock’s
paper (Stroock et al. 2002). When α > 0.2, the numerical prediction does not agree with
the analytical one due to the finite channel depth to width ratio (H/w = 0.5) used in this
study.
Stroock’s analytical prediction gave the mean helicity, and the periodic terms were not
counted in the calculation. In Figure 8.14, the helicity was computed numerically along
the channel longitudinal direction for different Reynolds numbers. The amplitudes of the
waves were proportional to the Reynolds number, however, the waves vibrated about the
same value, which was the mean used in Figure 8.13. The number of peaks in wave form
coincide with the number of periodic grooves in the channel. We can conclude that the
mean value of helicity only relies on groove aspect ratio α, width of the grooves and depth
of channel, and is independent of Reynolds number. Aspect ratio, α, contributes the most
to the helicity, and deeper grooves improve the mixing performance.
4 200 µm
Rate of Shear (tan(Ω)x100)
3.5
3
2.5
2
1.5
1
0.5
0
1200 1700 2200 2700 3200
Distance in the x direction (µm)
Re = 0.25 Re = 0.5 Re = 2.5 Re = 5
Figure 8.14: Helicity at different Reynolds numbers, numerical estimations for α = 0.10
173
8.5 Summary
A set of particle tracing algorithms was developed and used together with CFD sim-
ulations to evaluate the performance of micromixers with patterned grooves. The particle
tracing technique was used to construct Poincaré maps, which showed at low groove aspect
ratio (α = 0.05), the flow pattern was regular, and became less regular when the groove
aspect ratio became larger (α = 0.3). There was insufficient evidence to indicate chaotic
advection existed in such systems. However, the ability to create recirculation was analyzed
by calculating the necessary channel length to generate one recirculation. It was shown that
this length had an exponential relationship with the groove aspect ratio and was not signif-
icantly affected by the flow velocity. The diffusion path in the micromixer can be reduced
by the folding and stretching, and is directly related to the number of recirculations.
Simulation of diffusion-convection problems using a CFD package provides an under-

standing of mass transport in the micromixers with patterned grooves. Although the length
of channel that can be simulated is limited by the computer memory and speed, it is possi-
ble to demonstrate the stretching and folding of processes of material lines. Also, the rate of
shear computed from the simulations can be used to evaluate the performance of this type
of mixer. The mean value of helicity calculated numerically is independent of Reynolds
number (for Re ≤ 5), which agrees with the published results (Stroock et al. 2002).
Passive mixers with patterned grooves in microchannels have the drawback of creating
dead volumes, however, deeper grooves improve the mixing efficiency and reduce the chan-
nel length required for complete mixing. Also, these mixers work at relatively lower flow
velocity (Re ≤ 5), which reduces pressure drop, and are compatible with microfabrication
processes.
The investigation of micromixers with patterned grooves concluded the systematic

research of geometric variations on microfluidic mixing. The next chapter will summarize
174
the results and contributions to the knowledge, and make recommendations for further
research.
175
Chapter 9
Conclusions and Further Research
9.1 Overview
This study has investigated the ability of various geometric configurations to improve
the performance of microfluidic mixing, by using computational fluidic dynamic (CFD)
simulations, experimental investigations and developed particle tracing algorithms for so-
phisticated geometric configurations. This chapter presents the findings, contributions and
limitations of this research on passive mixing in microchannels. The last section makes
recommendations to extend this study on microfluidic mixing.
9.2 Contributions of the Research
9.2.1 Clarification of Related Theories to Microfluidic Mixing
The clarification of theories that could be applied to microfluidics is necessary and

helpful to this research. Therefore, the first contribution of this research was the identifica-
176
tion of applicable theories to micro-flow and microfluidic mixing.
This research started from the review of existing mixing techniques in a micro-device.
In microfluidic devices, it is very important to mix two or more reagents together. However,
in a typical microfluidic device, viscosity dominates flow, and as a result, mixing of fluids in
microfluidic devices is by virtue of diffusion, which is a slow process. This attracted many
researchers attention to study methods to improve microfluidic mixing. These methods
can be categorized in either the active or the passive mixing methods. Most active micro-
mixers are complex to fabricate and require an external power source. Passive mechanisms
have also been studied extensively, however, previous researchers have concentrated on
typically one of geometric factors, and an overall view of geometric configurations and their
effect on mixing was not clear. Therefore, it was the intention of this research to provide a
comprehensive research on passive mixing with various geometric configurations.
The influence of various geometric configurations of microchannels to mixing per-

formance could be studied experimentally or numerically. Because of the various com-
binations of design parameters of a microfluidic mixer, it is impractical to apply solely
experimental approaches. In this study, advanced CFD packages were used as the pri-
mary research method. Experimental approaches were also used as the complementary
method to prove the numerical approaches faithfully represented the physical principles of
microfluidic mixing. Work was undertaken to correctly and accurately carry out numeri-
cal simulations. In chapter 3, several methods were introduced to improve the quality of
numerical simulation, such as mesh density and sensitivity study. Then, CFD simulations
results were calibrated by comparing with the analytical results of Hele-Shaw flow in a
finite aspect ratio rectangular micro-channel. It was found that velocity field predicted by
CFD simulation could be validated by this method. However, in the case of simulating a
convection-diffusion problem, even though the velocity field was accurately predicted by
using sufficient mesh density, the simulation of mass transport was very sensitive to the
mesh quality and could not accurately predict fluids concentration distributions. It was
177
found that a uniform high quality mesh needs to be applied with sufficient mesh density.
In the case of modelling fluids in a complex configured microchannel, it is difficult to
achieve a uniform mesh. With the help of an advanced mesh generator, and mesh conver-
sion program developed in this research (Appendix B.1), high quality uniform mesh could
be translated into the MemCFD solvers. By comparing with the published experimental re-
sults, the numerical result using improved mesh could simulate a real convection-diffusion
problem. A simple experimental arrangement was also setup for carrying out tests. The
equipment included an optical microscope, a CCD camera and an image capture computer.
The simple laboratory configuration could give qualitative comparisons between experi-
ment visualizations and numerical simulations. A fabrication technique to produce simple
micro-mixers in microchannels was also described.
After the methodology was decided, several characteristics of microfluidic flow and
mixing were clarified. Three aspects were investigated, the influence of viscosity, the
convection-diffusion in a micro-device and the geometric characteristics of a common mi-
crochannel. It is widely accepted that the flow in a microchannel is slow viscous incom-
pressible flow. The investigation of micro-flow is to study the flow at very low Reynolds
number. Because the majority of microchannels are structured in a rectangular or similar
cross-section shape, viscous flow in a rectangular channel was studied. Because there was
a lack of theory to study micro-flow, it was the intention of this research to link microflu-
idics to the study of viscous flow in a duct, such as Hele-Shaw flow. Because the theory for
viscous flow in a duct has been well-developed, this linkage is important to provide theory
to the study of microfluidics. The linkage between Taylor-Aris dispersion and microflu-
idic mixing was also provided. The theory of Taylor-Aris dispersion in a small channel is
well-developed in analytical chemistry. By theoretical analysis and experimental proof, the
conditions of microfluidic mixing are very much similar to that of the Taylor-Aris disper-
sion. The Taylor dispersivity is normally a few orders greater than pure molecular diffusion
coefficient, therefore, this provided a guide to this research to introduce a convective effect
178
to improve microfluidic mixing.
To evaluate the convective effect, the Péclet number, which is the bulk mass transport
to diffusion mass transport ratio, was introduced. For a microfluidic device, it was found
that only a range of Péclet numbers need to be investigated. At very low Péclet numbers,
mixing is in the pure molecular diffusion zone, and mixing method other than diffusion is
not required. Also at very high Péclet numbers, a high Péclet number induces high pressure
drop and can make it difficult to drive the fluids through the microchannel. Therefore, the
Péclet number is an important parameter for the design of a micromixer, and was also used
to form the scope of this research.
Another important factor is the geometry of a micromixer. Because many micromix-

ers are structured with rectangular shaped cross-sections, the depth to width aspect ratio
of a microchannel was also investigated. The influences of channel depth to width aspect
ratio on pressure drop in the channel and mixing performance were investigated. It was
found that mixing benefited from increasing of the channel aspect ratio, and pressure drop
increased if the aspect ratio is either larger or less than unity. However, because of the dif-
ficulties to fabricate high aspect ratio channels, and also high channel depth to width aspect
ratio creates optical distortion for analyzing mixing results, therefore, shallower channels
are normally used in microfluidic applications. A small channel aspect ratio increases the
distance for the fluids to travel, and hence, it has poor mixing performance. Therefore,
the intention of this research was to find methods to enhance mixing performance in a
microfluidic analytical device with low aspect ratio microchannels.
9.2.2 Modifications to Serpentine-Shaped Microchannels
Three modifications of serpentine structures or serpentine-like structures were inves-

tigated in Chapter 5. In these designs, transverse convection was brought into the main
stream in the longitudinal direction. Hence, the convection could disturb the main flow to
179
induce bulk mass dispersion.
Some microfluidic applications do not require a complete mixing to achieve homoge-

nous solution. A concept of micromixer could generate different concentration gradients
was proposed in this chapter. The mixer was divided into several mixing zones by parti-
tion walls. The mixing simulation results showed that different concentration zones were
created in these divided zones.
The second modification to the serpentine micromixer design added by-passing chan-
nels. Partial flows were induced in the by-passing channels. The partial flow reduced the
amount of fluid entering the main channel from one side of the channel, therefore, one
fluid stream became thinner and closer to the wall. The slower flow motion at the wall and
thinner diffusion distance improves the mixing performance.
The last approach was to use 3D ramping structure to enhance mixing. It was inspired
by the complex 3D micro-structures created by Excimer Laser micromachining. At the
beginning of this research, the effects of these 3D structures on microfluidic mixing were
not clear. Using CFD modelling technique, the diffusion-convection simulation results
showed the positive influence of these structures on microfluidic mixing. At high Reynolds
numbers (Re = 50), irregular flow pattern in a 3D ramped channel was also shown. The 3D
ramping structure was novel and could be used to enhance mixing in a microchannel.
9.2.3 Design Guidelines for Obstacle Structures
3D ramping structures were difficult to parameterize for further investigation using

CFD. Therefore, the 3D structures were simplified to obstacles. The study of obstacle
structures were carried out in Chapters 6 and 7. In Chapter 6, it was found that cylindrical
obstacles could provide disturbances to the flow that was similar to ramping structures,
and resistance and disruption could be controlled by number and position of obstacles. By
180
studying distribution of mass fraction and velocity field affected by the obstacles in the
channels, the effect on mixing could be revealed. Obstacles in a microchannel with a low
Reynolds number cannot generate eddies or turbulence. However, the results demonstrated
that obstacles could improve mixing performance by affecting the flow pattern, and that the
asymmetric arrangement of obstacles alter the flow directions and force one fluid to flow
into another to create lateral mass transport. In Chapter 6, this phenomenon was termed
convective effect in order to differentiate it from the pure molecular diffusion.
A very important outcome from the studies was that asymmetric layout of obstacles
could enhance microfluidic mixing, while symmetric layouts had little effect on mixing.
Placing obstacles in the microchannels was a novel method for mixing in microfluidic de-
vices, and properly designed geometric parameters, such as layout and number of obstacles,
could improve the mixing performance without significantly increasing the pressure drop.
In Chapter 7, rectangular shaped obstacles provided more design parameters to be

studied, and these parameters are useful for optimizing the design of passive microfluidic
mixers using geometric variations. These parameters, such as number, size, gap between
obstacles, position, orientation in channel, and height ratio of obstacles, were studied for
their effect on mixing. The investigation of how flow rate affected mixing performance
illustrated promising results, and mixing efficiency remained constant or even increased
with flow rate when Re > 1. This suggested that convective effects dominated mixing over
diffusion, and diffusion over a short distance is negligible. It was found the greater the
lateral convection created by obstacles, the better the mixing performance. However, pres-
sure drop is also a critical performance criteria. A mixing index was defined to combine
the mixing performance and pressure drop for the design of a microfluidic mixer. Using the
mixing index, it was found that grooved patterns gave better overall mixing performance.
The results were used to identify feasible strategies for improving mixing by placing ob-
stacles in microchannels.
181
9.2.4 Design of Microfluidic Mixers with Patterned Grooves
When the height of obstacles to depth of channel ratio became ‘negative’, the ob-
stacles became grooves patterned on the channel wall. These grooves provided slip-like
boundary conditions, because of the lower resistance to flow due to fluid to fluid friction
rather than fluid to wall friction. Therefore, a lateral pressure gradient was built up to create
secondary flow. A detailed investigation of slantly patterned grooves inside a rectangular
microchannel was carried out in Chapter 8 using a numerical approach, termed particle
tracing algorithm, to trace the flow trajectories in the micro-channel.
A set of particle tracing algorithms was developed in this research and used together
with CFD simulations to evaluate the performance of micromixers with patterned grooves.
The particle tracing technique was used to construct Poincaré maps, which showed at low
groove aspect ratio (α = 0.05), that the flow pattern was regular, and became less regular
when the groove aspect ratio became larger (α ∼ 0.3). There was insufficient evidence
to indicate chaotic advection existed in such systems. However, the ability to create re-
circulation was analyzed by calculating the necessary channel length to generate one re-
circulation. It was shown that this length had an exponential relationship with the groove
aspect ratio and was not significantly affected by the flow velocity. The diffusion path in
the micromixer can be reduced by the folding and stretching, and is directly related to the
number of recirculations.
Simulation of diffusion-convection problems using a CFD package provides an under-

standing of mass transport in the micromixers with patterned grooves. Although the length
of channel that can be simulated is limited by the computer memory and speed, it is possi-
ble to demonstrate the stretching and folding processes of material lines. Also, the rate of
shear computed from the simulations can be used to evaluate the performance of this type
of mixer. The mean value of helicity calculated numerically is independent of Reynolds
number (for Re ≤ 5), which agrees with the published results (Stroock et al. 2002).
182
Passive mixers with patterned grooves in microchannels have the drawback of creating
dead volumes, however, deeper grooves result in better mixing efficiency and reduce the
channel length required for complete mixing. Also, these mixers work at relatively lower
flow velocity (Re ≤ 5), which reduces pressure drop, and are compatible with microfabri-
cation processes.
9.3 Limitations of the Research
The scope of this research was restricted within a limited range of Reynolds num-
bers. Mixing two or more fluid streams with very slow motion does not need assistance
other than by pure molecular diffusion, and fluids with very high velocity can be mixed
by conventional methods. However, microfluidic mixing provides a uniform mixture that
other technologies cannot achieve, and the study of high output micromixers for chemical
synthesis is another important aspect of research on microfluidic mixing.
One limitation of the research was the computational resources. It was impractical
to simulate the complete length of the microchannel in many cases. Instead, most of the
numerical models were built to simulate a section of the device. Although this did not
affect the evaluation of mixing performance of individual mixers, it can limit its application
to complex designs. In addition, because the size of the numerical model is limited by the
computer memory and main frequency of the computer central processing unit (CPU), the
quality of computer simulations is also affected. This is another reason why large models
could not be simulated.
Experimental measurements of mixing efficiency can provide overall information and

is useful for evaluating prototype micromixers. This research designed a simple testing
setup, which could provide qualitative measurement, in comparison to sophisticated mi-
crofluidic diagnostic technology used by some researchers (Santiago et al. 1998), including
183
µPIV and other advanced high power laser optical diagnostic systems.
9.4 Recommendations for Further Research
9.4.1 Experimental Methods
Numerical simulation provides a versatile method to study flow and mixing in mi-
crochannels. Numerical simulations can represent many practical problems when the mod-
els are correctly setup. Numerical approaches omit many less important conditions as in a
real environment, such as temperature, humidity and other environmental noise, and focus
on principal parameters. Numerical modelling and simulation may reveal more informa-
tion about the problem, especially when the incompressible viscous flows can be predicted
by numerical simulations. However, numerical modelling and simulation can represent the
reality only when the theories are well-developed and validated, such as in the case of in-
compressible viscous flow studied in this research. However, because the governing equa-
tions for µCFD involving electrohydrodynamics (EHD), magnetohydrodynamics (MHD)
and driving forces other than pressure, have not been fully developed, CFD was restricted
to be widely applied to microfluidic application. For example, the instability in an EHD
flow under high voltage has not built its theoretical foundation (Oddy et al. 2001), hence,
its governing equations have not been developed. In such cases, direct experimental mea-
surements of flow field are desirable. Nevertheless, even for the flow that can be predicted,
experimental measurements are required to validate prediction and build up confidence in
the simulations.
A personal communication in January 2003, with Professor Juan Santiago of Stanford

University suggested that the micro particle image velocimetry (µPIV) was now a mature
technique and could provide rich information in micro flows. Hence, 2D or even 3D veloc-
ity field of micro-flow can be measured and recorded. µPIV can be a complementary tool
184
to study microfluidic mixing.
There are still no µPIV algorithms for unsteady micro-flows, because it is difficult, if
not impossible, to direct a laser sheet beam into the micro-flow field. Instead, the depth-
of-field of the microscope has been used in current µPIV systems and this also limits its
application to measure unsteady or transient flows. Intensified CCD cameras also put re-
striction onto the development µPIV systems for its high cost. However, relatively low
cost CCD cameras still can measure flow velocity up to several hundred micrometers per
second, and this is sufficient for many microfluidic devices.
9.4.2 Complex 3D Geometry
3D Ramping
Although it is difficult for other microfabrication techniques to build very complex

3D structures, it is relatively simple for Excimer Laser micromachining. Applying mask
dragging and rotation techniques, periodic or non-periodic 3D structures in microchannels
can be rapidly fabricated. In Chapter 5, a 3D ramping structure was briefly introduced
and it was shown that it was possible to create 3D convection for high Reynolds number.
A microfluidic mixer with grooved surfaces has dead volume, which restricts its applica-
tion, and 3D ramping can provide a solution to the dead volume problem. However, more
research needs to be done to reveal the parameters of 3D structures that could create 3D
helical-shape secondary flows. This research has shown that periodic geometric structures
were the necessary conditions to create periodic secondary flows, and further research of
3D structures may provide another solution to microfluidic mixing.
185
Mixer with Patterned Grooves
The results in Chapter 8 have shown that the analysis of Poincaré maps to be a useful
approach for the study of mixing in a system that involves periodic recirculations. However,
for Re ≤ 5, the use of the Poincaré mapping technique did not provide a sound justification
for labelling micromixers reported in this research as chaotic or not. Chaotic mixing can
be verified by demonstrating that at least one positive Lyapunov exponent exists in the
system (Wolf, Swift, Swinney & Vastano 1985), and it is recommended as a topic of further
work. Chaotic mixing, which defines a flow region intermediate between turbulent and
laminar advection, has been recently introduced into microfluidic mixing literature.
Moreover, Figure 8.10 in Section 8.4.3 provides the relation between length of channel
to complete one recirculation and groove aspect ratio, for the channel depth to width ratio
H/w = 0.5. A chart of this relation for a range of H/w ratios, ie. 0.1 to 1, would be a
useful guideline for the design of the T-type micromixer reported here, in a field where at
present only few design guidelines have been established.
The periodic structures investigated in Chapter 8 were compared with the results in
one of Stroock’s paper (Stroock et al. 2002). In another Stroock’s similar paper (Stroock
et al. 2002), the patterned grooves changed orientation several times down the channel.
However, whether these variants are better for mixing than pure periodical structures or
worse is not clear, and is an interesting topic for further research.
9.5 Summary
This research developed feasible methods to enhance microfluidic mixing by using

geometric configurations in microchannels of microfluidic devices.
In a typical microchannel, viscosity dominates flow, and as a result, the Reynolds

186
number is low and the flow is laminar. Therefore, mixing of two or more fluid streams
in microfluidic devices is by virtue of diffusion, which is a slow process. On the other
hand, mixing two or more reagents together is a critical part of sample preparation in a
micro Total Analysis System (µTAS), and mixing elements need to be integrated on the
microchip. Therefore, microfluidic mixing is challenging and attracts the attention of many
researchers. Based on the mixing mechanism, micromixers have pure molecular diffusion
at one end (Kamholz & Yager 2001) and chaotic mixing on the other (Liu et al. 2000,
Stroock et al. 2002). Based on structures, micromixers can be categorized in either the
active or the passive mixing methods. Most active micromixers enhance mixing by stirring
flow. The stirring can be done mechanically (Lu et al. 2002) or by the means of magneto-
hydrodynamic (MHD) (Lemoff & Lee 2000), electro-hydrodynamic (EHD) (Choi & Ahn
2000) or by acoustic streaming (Zhu & Kim 1998, Yang et al. 2000) to create secondary
flows. The secondary flows can stretch and fold material lines to reduce the diffusion path
between fluid streams, and hence, enhance mixing. Active micromixers are particularly
suitable for chamber mixing, however, most active micromixers are complex to fabricate
and require an external power source.
Some passive mixers reduce diffusion path between fluid streams by splitting and re-
combining (Schwesinger et al. 1996, Koch et al. 1998, Ehrfeld et al. 1999). Recently, pas-
sive chaotic micromixers were investigated by using 3D serpentine-type channels (Beebe
et al. 2001, Liu et al. 2000). The flow in 3D serpentine channels demonstrated chaotic
advection at high flow rates (Re ∼ 70) and good mixing performance. However, the chal-
lenges of serpentine mixers are the micro-fabrication of complex 3D structures and the
need of high Reynolds number to stir the fluids to generate chaotic advection.
In this thesis, the author started the study by understanding passive microfluidic mix-
ing in a rectangular microchannel to provide the necessary background theory. Because
most of the effort in this research was on the study of various geometric structures, the
geometric variations were then defined as the channel geometric topologies and structures
187
inside a microchannel. Inspired by the serpentine shaped micromixers, modified serpentine

or serpentine-like structures were applied to improve their mixing performance. The 3D
ramping structures indicated promising mixing performance. However, 3D ramping struc-
tures had limited and arbitrary parameters, and were difficult to model, therefore, the 3D
ramping structures were simplified into cylindrical obstacles inside a microchannel. The
study on obstacles was extended to square and rectangular obstacles, which provided a rich
source of parameters for the study. During the research, overlapping studies by other re-
searchers were found that investigated grooved microchannels. Because of the relevance to
the present study of geometric structures of this research, extended research on mixers with
slantly grooved surfaces was carried out with particle tracing algorithms developed by the
author.
The vast parameters to be studied meant that there was a lack of experimental tools,
and therefore, the study was carried out principally by computational methods. To comple-
ment the numerical approach, selected samples were verified by experiments qualitatively.
The experiments conducted in this research were visualization of mixing in channels with
cylindrical obstacles and the visualization of Taylor-Aris dispersion in a straight rectangu-
lar channel. The numerical computations of mixers with grooved surfaces were compared
with the published analytical and measured data.
In conclusion, the principal objective of this research was to systematically investi-

gate the relation between various geometric configurations and their stirring mechanisms
towards microfluidic mixing by using computational fluid dynamic (CFD) simulations and
experimental validations. The results of this research provide a link between individual
geometric configuration of a passive mixer, and useful guidelines for the design of mi-
crofluidic passive mixers.
188
References
Ajdari, A. (2002). Transverse electrokinetic and microfluidic effects in micropatterned

channels: Lubrication analysis for slab geometries, Physical Review E 6501(1): art.
no.–016301.
Aref, H. (1984). Stirring by chaotic advection, Journal of Fluid Mechanics 143: 1–21.
Aris, R. (1956). On the dispersion of a solute in a fluid flowing through a tube, Proceed-
ing of the Royal Society of London, series A (Mathematical and Physical Science)
235: 67–77.
Bart, S. F., Tavrow, L. S., Mehregany, M. & Lang, J. H. (1990). Microfabricated electro-
hydrodynamic pumps, Sensors and Actuators, A: Physical 21(1-3): 193–197.
Bau, H. H., Zhong, J. H. & Yi, M. Q. (2001). A minute magneto hydro dynamic (MHD)
mixer, Sensors and Actuators B-Chemical 79(2-3): 207–215.
Beard, D. A. (2001). Taylor Dispersion of a Solute in a Microfluidic Channel, Journal of

Applied Physics 89(8): 4667–4669.
Beebe, D. J., Adrian, R. J., Olsen, M. G., Stremler, M. A., Aref, H. & Jo, B.-H. (2001).
Passive mixing in microchannels: Fabrication and flow experiments, Mecanique &
Industries 2(4): 343–348.
189
Bird, B. R., Stewart, W. E. & Lightfoot, E. N. (1960). Transport Phenomena, John Wiley
and Sons, Inc., USA.
Bokenkamp, D., Desai, A., Yang, X., Tai, Y., Marzluff, E. & Mayo, S. (1998). Microfab-
ricated Silicon Mixers For Submillisecond Quench-Flow Analysis, Analytical Chem-
istry 70(2): 232–236.
Bousse, L. (1999). Electrokinetic microfluidic systems, Proceedings of SPIE - The Inter-

national Society for Optical Engineering 3879: 2–8.
Brenner, M. P. & Stone, H. A. (2000). Modern Classical Physics Through the Work of G.i.
Taylor, Physics Today 53(5): 30–35.
Brody, J. P., Yager, P., Goldstein, R. E. & Austin, R. H. (1996). Biotechnology at low
Reynolds numbers, Biophysical Journal 71(6): 3430–3441.
Burgreen, D. & Nakache, F. R. (1964). Electrokinetic Flow in Ultrafine Capillary Slits, The
Journal of Physical Chemistry 68(5): 1084.
Catchpole, J. P. & Fulford, G. (1966). Dimensionless Groups, Industrial and engineering

chemistry 58(3): 46–60.
Chatwin, P. & Sullivan, P. (1982). The effect of aspect ratio on longitudinal diffusivity in
rectangular channels, Journal of Fluid Mechanics 120: 347–358.
Chen, C. H. & Santiago, J. G. (2002). A planar electroosmotic micropump, Journal of

Microelectromechanical Systems 11(6): 672–683.
Choi, J.-W. & Ahn, C. H. (2000). An active microfluidic mixer using electrohydrodynamic
(EHD) convection, Solid-State Sensor and Actuatoer Workshop, Hilton Head Island,
SC, pp. 52–55.
190
Darabi, J., Rada, M., Ohadi, M. & Lawler, J. (2002). Design, fabrication, and testing
of an electrohydrodynamic ion-drag micropump, Journal of Microelectromechanical
Systems 11(6): 684–690.
Deshmukh, A. A. (2001). Continuous Microfluidic Mixing Using Pulsatile Micropumps,

PhD thesis, University of California, Berkeley.
Dorfman, K. D. & Brenner, H. (2001). Comment on ”Taylor Dispersion of a Solute in a

Microfluidic Channel”, Journal of Applied Physics 90(12): 6553–6554.
Doshi, M. R., Daiya, P. M. & William, N. G. (1978). Three dimensional laminar dispersion
in open and closed rectangular conduits, Chemical Engineering Science 33: 795–804.
Dutta, D. & Leighton, D. T. (2001). Dispersion Reduction in Pressure Driven Flow Through
Microetched Channels, Analytical Chemistry 73(3): 504–513.
Ehrfeld, W., Golbig, K., Hessel, V., Lowe, H. & Richter, T. (1999). Characterization of
mixing in micromixers by a test reaction: Single mixing units and mixer arrays, In-
dustrial & Engineering Chemistry Research 38(3): 1075–1082.
Ehrfeld, W., Hessel, V. & Lowe, H. (2000). Microreactors: New Technology for Modern
Chemistry, John Wiley & Sons.
Einstein, A. (1905). On the Motion of Small Particles Suspended in Liquids, Annalen der
Physik 17: 549–560.
Elwenspoek, M., Lammerink, T., Miyake, R. & Fluitman, J. (1994). Towards integrated
microliquid handling systems, Journal of Micromechanics and Microengineering
4: 227–245.
Evans, J., Liepmann, D. & Pisano, A. P. (1997). Planar laminar mixer, Proceedings of the
IEEE Micro Electro Mechanical Systems (MEMS) pp. 96–101.
191
Forster, F. K., Bardell, R. L., Afromowitz, M. A., Sharma, N. R. & Blanchard, A. (1995).
Design, fabrication and testing of fixed-valve micro-pumps, American Society of Me-
chanical Engineers, Fluids Engineering Division (Publication) FED 234: 39–44.
Gad-el Hak, M. (1999). Fluid mechanics of microdevices-the Freeman scholar lecture,

Journal of Fluids Engineering, Transactions of the ASME 121(1): 5–33.
Gobby, D., Angeli, P. & Gavriilidis, A. (2001). Mixing Characteristics of T-Type Microflu-
idic Mixers, Journal of Micromechanics and Microengineering 11(2): 126–132.
Gravesen, P., Branebjerg, J. & Jensen, O. S. (1993). Microfluidics - a review, Journal of

Micromechanics and Microengineering 3(4): 168–182.
Harrison, D. J., Fluri, K., Seiler, K., Fan, Z., Effenhauser, C. S. & Manz, A. (1993). Micro-
machining a Miniaturized Capillary Electrophoresis-Based Chemical Analysis Sys-
tem on a Chip, Science 261(5123): 895–897.
Harvey, E. C., Rumsby, P. T., Gower, M. C. & Remnant, J. L. (1995). Microstructuring

by Excimer Laser, SPIE Proceedings, Micro- and Nanofabricated Electro-Optical
Mechanical Systems for Biomedical and Environmental Applications 2639: 266–277.
Hatch, A., Kamholz, A. E., Hawkins, K. R., Munson, M. S., Schilling, E. A., Weigl, B. H.
& Yager, P. (2001). A Rapid Diffusion Immunoassay in a T-Sensor, Nature Biotech-
nology 19(5): 461–465.
Hayes, J., Harvey, E., Ghantasala, M. & Dempster, B. (2001). A laser-ablation and lift-
off technique for fabricating simple microfluidic devices, Proceedings of SPIE - The
International Society for Optical Engineering 4236: 115–125.
He, B., Burke, B. J., Zhang, X., Zhang, R. & Regnier, F. E. (2001). A picoliter-volume
mixer for microfluidic analytical systems, Analytical Chemistry 73(9): 1942–7.
He, B., Tait, N. & Regnier, F. (1998). Fabrication of nanocolumns for liquid chromatogra-
phy, Analytical Chemistry 70(18): 3790–7.
192
Heng, K.-H., Huang, L., Wang, W. & Murphy, M. C. (1999). Development of a dif-
fuser/nozzle type micropump based on magnetohydrodynamic (MHD) principle, Pro-
ceedings of SPIE - The International Society for Optical Engineering 3877: 66–73.
Hessel, V., Hardt, S., Lowe, H. & Schonfeld, F. (2003). Laminar mixing in different in-
terdigital micromixers: I. Experimental characterization, AIChE Journal 49(3): 566–
577.
Hosokawa, K., Shimoyama, I. & Miura, H. (1993). Study of MHD (magnetohydrody-

namic) micropump, Nippon Kikai Gakkai Ronbunshu, C Hen/Transactions of the
Japan Society of Mechanical Engineers, Part C 59(557): 205–210.
Ismagilov, R. F., Rosmarin, D., Kenis, P. J. A., Chiu, D. T., Zhang, W., Stone, H. A. &
Whitesides, G. M. (2001). Pressure-Driven Laminar Flow in Tangential Microchan-
nels: an Elastomeric Microfluidic Switch, Analytical Chemistry 73(19): 4682–4687.
Ismagilov, R. F., Stroock, A. D., Kenis, P. J. A., Whitesides, G. & Stone, H. A. (2000).
Experimental and Theoretical Scaling Laws for Transverse Diffusive Broadening in
Two-Phase Laminar Flows in Microchannels, Applied Physics Letters 76(17): 2376–
2378.
Jacobson, S. C., McKnight, T. E. & Ramsey, J. M. (1999). Microfluidic Devices for Elec-
trokinetically Driven Parallel and Serial Mixing, Analytical Chemistry 71(20): 4455–
4459.
Jang, J. S. & Lee, S. S. (2000). Theoretical and Experimental Study of MHD (Magnetohy-
drodynamic) Micropump, Sensors and Actuators A-Physical 80(1): 84–89.
Jeon, N. L., Dertinger, S. K. W., Chiu, D. T., Choi, I. S., Stroock, A. D. & Whitesides, G. M.
(2000). Generation of Solution and Surface Gradients Using Microfluidic Systems,
Langmuir 16(22): 8311–8316.
193
Johnson, T. J. & Locascio, L. E. (2002). Characterization and optimization of slanted well

designs for microfluidic mixing under electroosmotic flow, Lab on a Chip 2(3): 135–
140.
Johnson, T. J., Ross, D. & Locascio, L. E. (2002). Rapid Microfluidic Mixing, Analytical
Chemistry 74(1): 45–51.
Jones, S. W., Thomas, O. M. & Aref, H. (1989). Chaotic advection by laminar flow in a
twisted pipe, Journal of Fluid Mechanics 209: 335–357.
Kamholz, A. E., Schilling, E. A. & Yager, P. (2001). Optical Measurement of Transverse

Molecular Diffusion in a Microchannel, Biophysical Journal 80(4): 1967–1972.
Kamholz, A. E., Weigl, B. H., Finlayson, B. A. & Yager, P. (1999). Quantitative Analysis of
Molecular Interaction in a Microfluidic Channel: the T-Sensor, Analytical Chemistry
71(23): 5340–5347.
Kamholz, A. E. & Yager, P. (2001). Theoretical analysis of molecular diffusion in pressure

driven laminar flow in microfluidic channels., Biophysical Journal 80(1): 155–160.
Kamholz, A. E. & Yager, P. (2002). Molecular Diffusive Scaling Laws in Pressure-Driven

Microfluidic Channels: Deviation From One-Dimensional Einstein Approximations,
Sensors and Actuators B-Chemical 82(1): 117–121.
Kamholz, A. & Yager, P. (2000). Performing quantitative diffusion-based assays in a mi-

crofluidic sensor, Annals of Biomedical Engineering 28(SUPPL.1): S97.
Kenwright, D. N. & Lane, D. A. (1996). Interactive Time-Dependent Particle Tracing

Using Tetrahedral Decomposition, IEEE Transactions on Visualization and Computer
Graphics 2(2): 120–129.
Klaus Schubert, K., Wilhelm Bier, E.-L., Erhard Herrmann, L., Thomas Menzel, H. &
Gerd Linder, K. (2000). Static Micromixer. Patent Number: US6082891.
194
Knight, J. B., Vishwanath, A., Brody, J. P. & Austin, R. H. (1998). Hydrodynamic fo-
cusing on a silicon chip: Mixing nanoliters in microseconds, Physical Review Letters
80(17): 3863–3866.
Koch, M., Chatelain, D., Evans, A. & Brunnschweiler, A. (1998). Two simple micromixers
based on silicon, Journal of Micromechanics and Microengineering 8(2): 123–126.
Koch, M., Schabmueller, C. G. J., Evans, A. G. R. & Brunnschweiler, A. (1999). Microma-

chined Chemical Reaction System, Sensors and Actuators A-Physical 74(1-3): 207–
210.
Koch, M., Witt, H., Evans, A. & Brunnschweiler, A. (1999). Improved characteriza-
tion technique for micromixers, Journal of Micromechanics and Microengineering
9(2): 156–158.
Kockmann, N., Foll, C. & Woias, P. (2003). Flow Regimes and Mass Transfer Characteris-
tics in Static Micro Mixers, in H. Becker & P. Woias (eds), Proceedings of SPIE, Mi-
crofluidics, BioMEMS, and Medical Microsystems, Vol. 4982, SPIE, San Jose, USA,
pp. 319–329.
Kulinsky, L., Wang, Y. & Ferrari, M. (1999). Electroviscous effects in microchannels,

Proceedings of SPIE - The International Society for Optical Engineering 3606: 158–
168.
Lamb, Horace, S. (1993). Hydrodynamics, 6th edn, Cambridge University Press.
Larsen, U. D. (2001). Micromixer Having a Mixing Chamber for Mixing Two Liquids
Through The Use of Laminar Flow. Patent Number: US6241379.
Lee, Y.-K., Deval, J., Tabeling, P. & Ho, C.-M. (2001). Chaotic Mixing in Electrokinetically
and Pressure Driven Micro Flows, The 14th IEEE Workshop on MEMS, Interlaken,
Switzerland, pp. 483–486.
195
Lemoff, A. V. & Lee, A. P. (2000). An AC magnetohydrodynamic micropump, Sensors

and Actuators B-Chemical 63(3): 178–185.
Li, D. Q. (2001). Electro-viscous effects on pressure-driven liquid flow in microchannels,

Colloids and Surfaces a-Physicochemical and Engineering Aspects 195(1-3): 35–57.
Li, Q., Wang, N. & Yi, D. (eds) (1999). Numerical Analysis, 3rd edn, Huazhong University
Press.
Liu, R. H., Lenigk, R. & Grodzinski, P. (2003). Acoustic micromixer for enhancement of
DNA biochip systems, Journal of Microlithography, Microfabrication, and Microsys-
tems 2(3): 178–184.
Liu, R. H., Stremler, M. A., Sharp, K. V., Olsen, M. G., Santiago, J. G., Adrian, R. J.,
Aref, H. & Beebe, D. J. (2000). Passive Mixing in a Three-Dimensional Serpentine
Microchannel, Journal of Microelectromechanical Systems 9(2): 190–197.
Liu, R. H., Yang, J. N., Pindera, M. Z., Athavale, M. & Grodzinski, P. (2002). Bubble-
induced acoustic micromixing, Lab on a Chip 2(3): 151–157.
Lu, L. H., Ryu, K. S. & Liu, C. (2002). A magnetic microstirrer and array for microfluidic
mixing, Journal of Microelectromechanical Systems 11(5): 462–469.
Mahe, C., Tranchant, J. F., Tromeur, M. & Schwesinger, N. (2003). Emulsions on demand
using microstructured devices, in H. Becker & P. Woias (eds), SPIE: Microfluidics,
BioMEMS, and Medical Microsystems, Vol. 4982, San Jose, CA, USA, pp. 290–296.
Manz, A., Effenhauser, C. S., Burggraf, N., Harrison, D. J., Seiler, K. & Fluri, K. (1994).
Electroosmotic pumping and electrophoretic separations for miniaturized chemical
analysis systems, Journal of Micromechanics and Microengineering 4: 257–265.
Manz, A., Graber, N. & Widmer, H. (1990). Miniaturized total chemical analysis systems:
A novel concept for chemical sensing, Sensors and Actuators, B: Chemical B1(1-
6): 244–248.
196
Metcalfe, G. & Ottino, J. M. (1994). Autocatalytic Processes in Mixing Flows, Physical

Review Letters 72(18): 2875–2878.
Molho, J. I., Herr, A. E., Kenny, T. W., Mungal, M. G., Deshpande, M., Gilbert, J. R.,
Garguilo, M. G., Paul, P. H., St. John, P. M., Woudenberg, T. M. & Connell, C. (1998).
Fluid transport mechanisms in microfluidic devices, American Society of Mechanical
Engineers, Dynamic Systems and Control Division (Publication) DSC 66: 69–76.
Monnier, H., Wilhelm, A. M. & Delmas, H. (1999). Influence of ultrasound on mixing

on the molecular scale for water and viscous liquids, Ultrasonics Sonochemistry 6(1-
2): 67–74.
Nguyen, N. T., Huang, X. Y. & Chuan, T. (2002). MEMS-Micropumps: A Review, Trans-

actions of ASME 124: 384–392.
Nguyen, N. T., Meng, A. H., Black, J. & White, R. M. (2000). Integrated Flow Sensor
for in Situ Measurement and Control of Acoustic Streaming in Flexural Plate Wave
Micropumps, Sensors and Actuators A-Physical 79(2): 115–121.
Nielsen, C. R., Branebjerg, J. A., Gravesen, P., Dyhr-Mikkelsen, P. & Gade, N. (2001).
Micro-Mixer and Mixing Method, US Patent US6,190,034 B1.
Niu, X. & Lee, Y.-K. (2003). Efficient spatial-temporal chaotic mixing in microchannels,
Journal of Micromechanics and Microengineering 13(3): 454–462.
Oddy, M. H., Santiago, J. G. & Mikkelsen, J. C. (2001). Electrokinetic instability mi-

cromixing, Analytical Chemistry 73(24): 5822–5832.
Ottino, J. (1989). The Kinematics of mixing: stretching, chaos, and transport, Cambridge
University Press.
Pan, W. (1988). Engineering Fluid Dynamics, Tsing Hua University Press.
Petersen, K. E. (1982). Silicon as a Mechanical Material, Proceeding IEEE 70(5): 420–457.

197
Purcell, E. M. (1977). Life at low Reynolds number., American Journal of Physics 45(1): 3–
11.
Reid, R. C., Prausnitz, J. M. & Poling, B. E. (1987). The Properties of Gases and Liquids,
4th edn, McGraw-Hill, Inc.
Rice, C. L. & Whitehead, R. (1965). Electrokinetic Flow in a Narrow Cylindrical Capillary,

The Journal of Physical Chemistry 69(11): 4017.
Rife, J. C., Bell, M. I., Horwitz, J. S., Kabler, M. N., Auyeung, R. C. Y. & Kim, W. J.
(2000). Miniature Valveless Ultrasonic Pumps and Mixers, Sensors and Actuators
A-Physical 86(1-2): 135–140.
Rudman, M. (1998). Mixing and particle dispersion in the wavy vortex regime of Taylor-
Couette flow, Aiche Journal 44(5): 1015–1026.
Santiago, J. G. (2002). Fluid Flow in Micro Devices, Lecture notes on microfluidics, De-
partment of Mechanical Engineering, Stanford University .
Santiago, J. G., Wereley, S. T., Meinhart, C. D., Beebe, D. J. & Adrian, R. J. (1998). A Par-
ticle Image Velocimetry System for Microfluidics, Experiments in Fluids 25(4): 316–
319.
Schwesinger, N., Frank, T. & Wurmus, H. (1996). A Modular Microfluidic System With an
Integrated Micromixer, Journal of Micromechanics and Microengineering 6(1): 99–
102.
Selverov, K. P. & Stone, H. A. (2001). Peristaltically driven channel flows with applications
toward micromixing, Physics of Fluids 13: 1837–859.
Shoji, S. & Esashi, M. (1994). Microflow devices and systems, Journal of Micromechanics
& Microengineering 4: 157–171.
198
Stroock, A. D., Dertinger, S. K., Ajdari, A., Mezic, I., Stone, H. A. & Whitesides, G. M.
(2002). Chaotic Mixer for Microchannels, Science 295(5555): 647–651.
Stroock, A. D., Dertinger, S. K., Whitesides, G. M. & Ajdari, A. (2002). Patterning flows
using grooved surfaces, Analytical Chemistry 74(20): 5306–5312.
Taylor, G. (1953). Dispersion of soluble matter in solvent flowing slowly through a tube,
Proceeding of the Royal Society of London, series A (Mathematical and Physical
Science) 219: 186–203.
Taylor, G. (1954). Conditions under which dispersion of a solute in a stream of solvent can
be used to measure molecular diffusion, Proceeding of the Royal Society of London,
series A (Mathematical and Physical Science) 225: 473–477.
Veenstra, T., Lammerink, T., Elwenspoek, M. & van den Berg, A. (1999). Characterization
method for a new diffusion mixer applicable in micro flow injection analysis systems,
Journal of Micromechanics and Microengineering 9(2): 199–202.
Vivek, V., Zeng, Y. & Kim, E. S. (2000). Novel acoustic-wave micromixer, Proceedings of
the IEEE Micro Electro Mechanical Systems (MEMS) pp. 668–673.
Weigl, B. H. & Yager, P. (1999). Microfluidics - Microfluidic Diffusion-Based Separation

and Detection, Science 283(5400): 346–347.
Wolf, A., Swift, J. B., Swinney, H. L. & Vastano, J. A. (1985). Determining Lyapunov
exponents from a time series, Physica D: Nonlinear Phenomena 16(3): 285–317.
Yang, Z., Goto, H., Matsumoto, M. & Maeda, R. (2000). Active Micromixer for Mi-
crofluidic Systems Using Lead- Zirconate-Titanate(Pzt)-Generated Ultrasonic Vibra-
tion, Electrophoresis 21(1): 116–119.
Yi, M. Q., Qian, S. Z. & Bau, H. H. (2002). A magnetohydrodynamic chaotic stirrer,

Journal of Fluid Mechanics 468: 153–177.
199
Zhu, X. & Kim, E. S. (1998). Microfluidic motion generation with acoustic waves, Sensors
and Actuators, A: Physical 66(1-3): 355–360.
200
Appendix A
List of Publications
A.1 Journal Publications
Wang, H., Iovenitti, P., Harvey, E. C. & Masood, S. (2003). Numerical investigation of
mixing in microchannels with grooved surfaces, Journal of Micromechanics and Micro-
engineering 13(6):801-8.
Wang, H., Iovenitti, P., Harvey, E. & Masood, S. (2002). Optimizing layout of obstacles
for enhanced mixing in microchannels, Smart Materials and Structures 11(5): 662-7.
201
A.2 International Conference Publications
Wang, H., Iovenitti, P., Harvey, E. & Masood, S. (2003). Passive mixing in microchannels
by applying geometric variations, SPIE, Micromachining and Microfabrication, San Jose,
USA, 4982: 282-9.
Wang, H., Iovenitti, P., Harvey, E. & Masood, S. (2002). Mixing of two fluid streams in a
microchannel using Taylor-Aris dispersion effect, Proceedings of SPIE, Smart Electronics
and MEMS, Melbourne, Australia, 4937: 158-163.
Wang, H., Iovenitti, P., Harvey, E. & Masood, S. (2001). Mixing of liquids using obsta-
cles in microchannels, Proceedings SPIE, BioMEMS and Smart Nanostructures,Adelaide,
Australia, 4590: 204-212.
Wang, H., Masood, S., Iovenitti, P. & Harvey, E. (2001). Application of fused deposition
modelling rapid prototyping system to the development of microchannels, Proceedings of
SPIE, Smart Electronics and MEMS, Adelaide, Australia, 4590: 213-220.
Wang, H., Iovenitti, P., Harvey, E. & Masood, S. (2000). A simple approach for modelling
flow in a microchannel, SPIE, Smart Electronics and MEMS, Melbourne, Australia, 4236:
99-106.
202
A.3 Local Publications and Presentations
Wang, H., Iovenitti, P., Harvey, E. C. & Masood, S. (2002). Passive mixing in a microchan-
nel, in D. Toncich (ed.), Profiles in Industrial Research Knowledge and Innovation 2002,
Industrial Research Institute Swinburne, Hawthorn, pp. 261-268, ISBN 1876 567 04X
Wang, H., Iovenitti, P., Harvey, E. C. & Masood, S. (2001). Mixing of fluids in microchan-
nels, in D. Toncich (ed.), Profiles in Industrial Research Knowledge and Innovation 2001,
Industrial Research Institute Swinburne, Hawthorn, pp. 225-236, ISBN 1876 567 03 1
Wang, H., Iovenitti, P., Harvey, E. C. & Masood, S. (2000). Flow in micro-domain, in D.
Toncich (ed.), Profiles in Industrial Research Knowledge and Innovation 2000, Industrial
Research Institute Swinburne, Hawthorn, ISBN 1876 567 023
203
Appendix B
Selected Fortran Codes Developed to

Study Microfluidic Mixing
In this appendix, programs were developed for the purpose of numerical simulations.
These programs included Fortran code for converting neutral mesh to an I-Deas universal
mesh format (Appendix B.1), that was used to improve mesh quality through simulations
carried out in this research, a program developed in Section 4.2 to calculate viscous flow
velocity distribution in a rectangular duct (Appendix B.2) and a main program for particle
tracing developed in Section 8.3 (Appendix B.3).
B.1 Fortran Codes for Mesh Conversion
TM
For complex geometric structures, MemCFD could not generate uniform structured
mesh. Another finite element mesh generation package, Gambitr 2.04, was used. The
TM
following Fortran code was developed to import Gambitr neutral mesh into MemCFD .
This mesh converting program had been used throughout the whole research project pre-
204
sented in this thesis. The flow chart is shown in Figure B.1
C*******************************************************************
C** *
C** convert Gambit neutral file to ideas universal neutral file *
C** for import mesh file to CoventorWare *
C** *
C** by Hengzi Wang *
C** *
C** 10-July-2002 *
C** *
C** hengzi_wang@yahoo.com.au *
C** *
C** IRIS, SWINBURNE UNIVERSITY OF TECHNOLOGY *
C*******************************************************************
C** *
C** PROGRAM GAMBIT2UNV *
C** *
C*******************************************************************
C**
C**
PROGRAM GAMBIT2UNV
USE DFPORT
IMPLICIT NONE
CHARACTER*80 NEUFILE,UNVFILE,CONTRINFO,GFILE,HED,PROG,REC5,
& ENDOFSECTION,MTHYEAR,GROUP,ELEMENT,MATERIAL,NSFLAGS,
& ELMAT,TODAY,MATNAM
INTEGER NUMNP,NELEM,NGRPS,NBSETS,NDFCD,NDFVL,ND,NE,I,
& NTYPE,NDP,NODE(30),NGP,NELGP,MTYP,NFLAGS,ISOLVE(10),
& NELT(1000000),REDATE(8),MATNUM,NODET
REAL(8) REVL,X,Y,Z
C**
C** Give the name for both input and output mesh files
C**
WRITE(*, ’(A,$)’) ’ Enter the Input mesh file, .neu: ==>’
READ(*, ’ (A) ’) NEUFILE
WRITE(*, ’(A,$)’) ’ Enter the OUTPUT mesh file: .unv==>’
READ(*, ’ (A) ’) UNVFILE
C**
C** Open Gambit neutral file, and ideas universal neutral file
C**
OPEN(UNIT=10,FILE=NEUFILE,FORM=’FORMATTED’,STATUS=’OLD’)
OPEN(UNIT=20,FILE=UNVFILE)
C**
205
Read gambit neutral mesh
Write universal datasets

1710, 164, 2400, 2420, 2411,
and 2412.
no
End of writing?
yes
Covert mesh topology to
universal dataset 2412
no
End of converting?
yes
Export mesh to a universal file.
Figure B.1: Flow chart to convert gambit neutral mesh to universal mesh
206
C** Write IDEAS universal dataset 1710

C**
WRITE(20,’(I6)’) -1
WRITE(20,’(I6)’) 1710
WRITE(20,*)’=================================================
&=====================’
WRITE(20,*) ’MATERIAL’
WRITE(20,*)’=================================================
&=====================’
MATNUM=2
WRITE(*,’(A,$)’) ’INPUT MATERIAL NAME: ’
READ(*,’(A) ’) MATNAM
MATNAM=’WATER’
WRITE(20,10) MATNUM,MATNAM
10 FORMAT(I10,2X,A40)
WRITE(20,’(I10,A)’) 0,’ LINE(S) OF TEXT’
WRITE(20,’(I10,A)’) 0,’ MATERIAL CLASS(ES)’
WRITE(20,’(I10,A)’) 0,’ MATERIAL ATTRIBUTE(S)’
WRITE(20,’(I10,A)’) 0,’ MATERIAL COMPONENT(S)’
WRITE(20,’(I10,A)’) 0,’ MATERIAL SPECIFICATION(S)’
WRITE(20,*) ’-------------------------------------------------
&---------------------’
WRITE(20,’(I10,A)’) 0,’ MATERIAL VARIABLE(S)’
WRITE(20,*) ’-------------------------------------------------
&---------------------’
WRITE(20,’(I10,A)’) 0,’ MATERIAL PROPERT(IES)’
WRITE(20,*) ’-------------------------------------------------
&---------------------’
WRITE(20,’(I10,A)’) 1,’ REFERENCE ENTITIES’
WRITE(20,’(I10,A)’) 1,’ MATERIAL TYPES’
WRITE(20,*) ’FEM ISOTROPIC
& MATERIALS’
WRITE(20,*)’=================================================
&=====================’
WRITE(*,*) ’end of writing dataset 1710...’
WRITE(*,*) fdate()
WRITE(*,*) ’start writing dataset 164...’
WRITE(20,’(I6)’) -1
C**
C**
WRITE(20,’(I6)’) -1
WRITE(20,’(I6)’) 164
WRITE(20,’(I10,A20,I10)’) 5,’mm (milli-newton) ’,2
207
WRITE(20,’(3D25.17)’) 1000.0,1000.0,1.0,273.149999999999977
WRITE(20,’(I6)’) -1
C**
C**
WRITE(20,’(I6)’) -1
WRITE(20,’(I6)’) 2400
WRITE(20,’(I12,2I6,I12)’) 1001,7,1,0
WRITE(20,’(A)’) ’FEM_1’
WRITE(20,*)
WRITE(20,’(32I2)’) 0,1,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,
&0,0,0,0,0,0,0,0,0,0
TODAY=FDATE()
TODAY=TODAY(5:)
PRINT *, TODAY
WRITE(20,’(A20,3I12)’) TODAY,30,30,17
WRITE(20,’(I12)’) 0
WRITE(20,’(I6)’) -1
C**
C**
WRITE(20,’(I6)’) -1
WRITE(20,’(I6)’) 2420
WRITE(20,’(I10)’) 1002
WRITE(20,’(A)’) ’PART_1’
WRITE(20,’(3I10)’) 1,0,8
WRITE(20,’(A)’) ’CS1’
WRITE(20,’(1P3D25.16)’) 1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0,
& 0.0,0.0,0.0
WRITE(20,’(I6)’) -1
C**
C** Read Gambit neutral file, and Write Ideas dataset 2411
C**
C**
C** Read the header
READ(10,100) CONTRINFO
READ(10,100) GFILE
WRITE(*,100) GFILE
READ(10,110) HED
WRITE(*,110) HED
READ(10,110) PROG
WRITE(*,110) PROG
READ(10,110) MTHYEAR
208
READ(10,110) REC5
C** Write the header of 2411
WRITE(20,’(I6)’) -1
WRITE(20,’(I6)’) 2411
READ(10,130) NUMNP,NELEM,NGRPS,NBSETS,NDFCD,NDFVL
WRITE(*,130) NUMNP,NELEM,NGRPS,NBSETS,NDFCD,NDFVL
READ(10,100) ENDOFSECTION
WRITE(*,100) ENDOFSECTION
100 FORMAT(A40)
110 FORMAT(A80)
120 FORMAT(A44,F5.2)
130 FORMAT(6I10)
140 FORMAT(5X,’NUMNP’,5X,’NELEM’,5X,’NGRPS’,4X,’NBSETS’,5X,’NDFCD’,
& 5X,’NDFVL’)
C**
C** Read and write nodal coordinates
C**
READ(10,’(A30)’) CONTRINFO
WRITE(*,’(A30)’) CONTRINFO
ND=1
DO WHILE(ND.LT.NUMNP)
READ(10,150) ND,X,Y,Z
WRITE(20,’(4I10)’) ND,1,1,11
WRITE(20,’(1P3D25.16)’) X,Y,Z
END DO
150 FORMAT(I10,1X,3E20.12)
READ(10,’(A40)’) ENDOFSECTION
IF (ENDOFSECTION.EQ.’ENDOFSECTION’) THEN
WRITE(20,’(I6)’) -1
ELSE
WRITE(*,*) ’Reading error, exit 2411’
ENDIF
C**
C** Read and write elements (Ideas 2412)
C**
WRITE(20,’(I6)’) -1
WRITE(20,’(I6)’) 2412
NE=1
DO WHILE(NE.LT.NELEM)
READ(10,160) NE,NTYPE,NDP,(NODE(I),I=1,NDP)
C** 115 for linear elements, and 116 for parabolic elements
C** always use 115 transfer mesh from Gambit to CoventorWare
C** unless, FEMTOOL is called!
209
WRITE(20,’(6I10)’) NE,115,1,2,7,NDP
C**
C** Convert gambit topology to ideas universal topology
C**
C** Exchange no.3 & 4 elements
NODET=NODE(3)
NODE(3)=NODE(4)
NODE(4)=NODET
C** Exchange no.7 & 8 elements
NODET=NODE(7)
NODE(7)=NODE(8)
NODE(8)=NODET
WRITE(20,’(8I10)’) (NODE(I),I=1,NDP)
END DO
160 Format(I8,1X,I2,1X,I2,1X,7I8:/(15X,7I8:))
C**
WRITE(20,’(I6)’) -1
ELSE
WRITE(*,*) ’Reading error, exit 2412’
ENDIF
C**
C** Read and write element group control information record
C**
WRITE(*,’(I6)’) -1
READ(10,170) GROUP,NGP,ELEMENT,NELGP,MATERIAL,MTYP,NSFLAGS,NFLAGS
WRITE(*,180) GROUP,NGP,ELEMENT,NELGP,MATERIAL,MTYP,NSFLAGS,NFLAGS
170 FORMAT(A8,I10,A11,I10,A10,I10,A8,I10)
180 FORMAT(A8,I10,A11,I10,A10,I10,A8,I10)
READ(10,’(A40)’) ELMAT
WRITE(*,’(A40)’) ELMAT
READ(10,190) (ISOLVE(I),I=1,NFLAGS)
WRITE(*,190) (ISOLVE(I),I=1,NFLAGS)
190 FORMAT(10I8)
READ(10,190) (NELT(I),I=1,NELGP)
WRITE(*,190) (NELT(I),I=1,NELGP)
C**
WRITE(*,’(I6)’) -1
ENDIF
CLOSE(10)
210
CLOSE(20)
STOP
END
211
B.2 Fortran Codes for Velocity Field in a Rectangular Duct
The analytical prediction of viscous flow in a rectangular channel comprises Fourier

series. Because these Fourier series converge very quickly, a finite number of Fourier
series should be calculated (Figure B.2). This program gives the analytical solution of flow
field in a rectangular microchannel. This program was developed and used in Section 4.2
to calculate the velocity field in a rectangular microchannel, and was used throughout this
research to find out proper mesh uniformity and density to achieve high quality simulations.
program heleshaw
USE DFPORT
c
c This function was used to calculate the velocity distribution in
c a microchannel with finite aspect ratio - you know I am talking
c about the velocity field a microchannel.
c
real*8 muu,ymin,ymax,zmin,zmax,dy,dz,z,W,b,L,G,u(100,100),qfr,PI
1,u0,u1(100,100),u2(100,100),za(100),ya(100),ku
character*40 uyz,TODAY
integer*4 npts,n,ny,nz
C**
C** Read flow (qfr) and the output file name.
C**
WRITE(*,’(A,$)’)’Enter the flow rate in (mu3/sec): ==>’
READ(*, *) qfr
WRITE(*, ’(A,$)’) ’ Enter the OUTPUT velocity field: .txt==>’
READ(*, ’(A) ’) uyz
OPEN(UNIT=20,FILE=uyz)
c qfr = 2e6
muu = 1.002e-9
ymin = -200
ymax = 200
c width of the channel
zmin = -50
zmax = 50
c thickness of the channel
npts = 40
dy = (ymax-ymin)/npts
dz = (zmax-zmin)/npts
212
Input channel geometry
Calculate pressure gradient

at certain channel width to
depth aspect ratio
Calculate velocity field, with

the above pressure gradient
Export velocity field
End
Figure B.2: Flow chart for viscous flow in a rectangular channel

213
a = (ymax-ymin)/2
b = (zmax-zmin)/2
L = 5000
PI = 3.14159265359
G1 = 0.0
c aspect ratio = 4
ku = 14.2
c input of flow rate to calculate pressure gradient G

G = ku*muu*qfr/(ymax-ymin)/(zmax-zmin)/(zmax-zmin)**2
write(20,’(1p1e16.6)’) G
n = 1000
ny = int((ymax-ymin)/dy)
nz = int((zmax-zmin)/dz)
do i = 1,ny+1
do j = 1,nz+1
u (i,j) = 0.0
enddo
enddo
do i=1,nz+1
za(i) = zmin+(i-1)*dz
enddo
do i=1,ny+1
ya(i) = ymin+(i-1)*dy
enddo
do i = 1,ny+1
do j=1,nz+1
do k = 1,800
u(i,j) = u(i,j)+(16*G*a**2)/(muu*PI**3)*((-1)**(k-1)*
1 (1-cosh((2*k-1)*PI*za(j)/2/a)/cosh((2*k-1)*
1 PI*b/2/a))*cos((2*k-1)*PI*ya(i)/(2*a))/(2*k-1)**3)
enddo
enddo
enddo
write(*,*) nz+1,ny+1
write(20,’(1P1E20.6)’ ) (za(i),i=1,nz+1)
write(20,*)
write(20,*)
write(20,’(1P1E20.6)’ ) (ya(i),i=1,ny+1)
write(20,*)
214
write(20,*)
write(20,’(1P41E20.6)’ ) ((u(i,j),i=1,ny+1),j=1,nz+1)
write(20,*)
C**
TODAY=FDATE()
PRINT *, TODAY
WRITE(20,*) TODAY
C**
CLOSE(20)
stop
end
215
B.3 Fortran Codes for Particle Tracing - Main Program
This program was developed to compute the particle trajectories in a micromixer. One
periodic flow field was clipped from the exported simulation results and used repeatedly to
record streaklines or Poincaré map. It was a lengthy program, and only the main program
is included here to describe the logical flow of the algorithms. The detailed subroutines are
not included in this thesis but can be worked out from Section 8.3, where the algorithms
were described. This program was developed to calculate Poincaré maps and streaklines
so that mixing performance of micromixers with grooved surface could be evaluated. The
flow chart is shown in Figure B.3
C*******************************************************************
C** *
C** extract velocity vector field from CoventorWare exported *
C** result text file(s). The extracted velocity field is used to *
C** calculate the POINCARE MAP TO EVALUATE PARTICLE CONVECTION. *
C** *
C** by Hengzi Wang *
C** *
C** 02-Sept-2002 *
C** *
C** hengzi_wang@yahoo.com.au *
C** *
C** IRIS, SWINBURNE UNIVERSITY OF TECHNOLOGY *
C*******************************************************************
C** *
C** PROGRAM POINCARE version 3.0 BETA *
C** *
C*******************************************************************
C** version 3, change the interpolation strategy. Use eight nodes
C** to interpolate the velocity.
C**
PROGRAM POINCARE
USE DFPORT
USE IMSLF77
IMPLICIT NONE
CHARACTER*80 TECFILE,XSECFILE,CONTRINFO,NODES,ELEMS,RESTS,TODAY
CHARACTER*150 GFILE
INTEGER NNODE,NELEM,ND,FNODE(100000),NODEX,NANG,NSTR,I,J,K,L,M,N,
216
Read simulation result
Determine and clip one

periodic grooved section
no
End of writing?
yes
Assigning initial points locations
Integrating new position
no
In the flow domain?
yes
Locate the mesh element
containing this particle
Interpolating velocity of this

particle
Record the exit point to Poincare map
Initial point <= exit point
Stop when cycles are completed
Figure B.3: Flow chart to compute particle trajectories

217
& YI,POIN,L1,L2,L3,NDD,NEIDX,NIIDX,NJIDX,NKIDX,NE,NL,
3 NIMAX,NJMAX,NKMAX,NS,IDX(8),NMAX,NP,N5TET(5,4),NEWNS,
4 NI,NJ,NK,NEIN(8),CELLIDX(2,2,2),IELE(500000,8),NNODEUP
LOGICAL FLAG
DOUBLE PRECISION CORD(500000,3),V(500000,3),PRESURE(500000),
& MF1(500000),MF2(500000),VISO(500000),XMIN,XMAX,XRES,
& YMAX,YMIN,ZMAX,ZMIN,XMID,YMID,ZMID,XINC,XM,XS,
& DELTZ,DELTA,FLORAT,GRIDX,AREAO,STRP(20,3),XN1(3),DX,
& XSTART,XSTOP,VN(3),VP(3),STRLIN(2000,3),POICAR(10000,3),
& K1(3),K2(3),K3(3),K4(3),PV(100,100,3),
& P1(5,100,100,3),
& P1V(100,100,3),P2V(100,100,3),P3V(100,100,3),P4V(100,100,3),
& T,DT,P,XP(0:2),YP(0:2),ZP(0:2),EPSILON,
& P5(100,100,3),P5V(100,100,3),PT(3),XI(5),ETA(5),ZETA(5)
COMMON /CO/CORD,V/EL/IELE
1 /MAXS/XMAX,XMIN,YMAX,YMIN,ZMAX,ZMIN,XSTART,XINC
C**
C** Give the name for both input and output results files
C**
WRITE(*,’(A,$)’)’Enter the ASCII tecplot result file, .TXT: ==>’
READ(*, ’ (A) ’) TECFILE
WRITE(*, ’(A,$)’) ’ Enter the OUTPUT velocity field: .txt==>’
READ(*, ’ (A) ’) XSECFILE
XINC=200.0
C**
C** Open TECPLOT file, and export velocity vector file
C**
OPEN(UNIT=10,FILE=TECFILE,FORM=’FORMATTED’,STATUS=’OLD’)
OPEN(UNIT=20,FILE=XSECFILE)
OPEN(UNIT=30,FILE=’POIMAP.TXT’)
C**
C** Start to write data into vecfile
C**
TODAY=FDATE()
PRINT *, TODAY
WRITE(30,*) TODAY
C**
C** Read the header
C**
READ(10,’(A)’) CONTRINFO
READ(10,’(A)’) GFILE
WRITE(*,’(A)’) GFILE
READ(10,’(A)’)
READ(10,110) NODES,NNODE,ELEMS,NELEM,RESTS
218
WRITE(*,110) NODES,NNODE,ELEMS,NELEM,RESTS
110 FORMAT(A7,I6,A4,I6,A30)
C**
C** Read and write nodal coordinates
C** Read and write velocity field
C**
WRITE(*,*) ’Now reading the file, please wait...’
DO ND=1,NNODE
READ(10,’(10E16.6)’) CORD(ND,1),CORD(ND,2),CORD(ND,3),V(ND,1),
& V(ND,2),V(ND,3),PRESURE(ND),MF1(ND),MF2(ND),VISO(ND)
C & V(ND,2),V(ND,3),PRESURE(ND)
IF (ND.EQ.NNODE) THEN
WRITE(*,’(A)’) ’END OF READING NODES’
READ(10,*)
ENDIF
END DO
C
DO N=1,NELEM
READ(10,’(8I8)’) (IELE(N,J),J=1,8)
IF (N.EQ.NELEM) THEN
WRITE(*,’(A)’) ’END OF READING ELEMENTS’
ENDIF
END DO
C
WRITE(*,’(A)’) ’CLOSING THE INPUT FILE’
CLOSE(10)
DO I=1,10
WRITE(*,’(8I8)’) (IELE(I,J),J=1,8)
ENDDO
WRITE(*,’(8I8)’) (IELE(NELEM,J),J=1,8)
C**
C** FIND THE BOUNDRARY (XMIN,XMAX,YMIN,ZMIN,YMAX ZMAX)
C**
NNODEUP=1
DO I=1,NELEM
DO J=1,8
IF (IELE(I,J).GE.NNODEUP) NNODEUP=IELE(I,J)
ENDDO
ENDDO
WRITE(*,’(I8)’) NNODEUP
XMIN=CORD(1,1)
XMAX=CORD(1,1)
DO I=1,NNODEUP
IF (CORD(I,1).GT.XMAX) THEN
219
XMAX=CORD(I,1)
ELSE IF (CORD(I,1).LT.XMIN) THEN
XMIN=CORD(I,1)
ENDIF
END DO
WRITE(*,*) ’XMIN = ’,XMIN,’XMAX = ’,XMAX
C
YMIN=10000.0
YMAX=0.0
DO I=1,NNODEUP
IF (ABS(CORD(I,1)-XMAX).LE.0.01) THEN
IF (CORD(I,2).GT.YMAX) THEN
YMAX=CORD(I,2)
ELSE IF (CORD(I,2).LT.YMIN) THEN
YMIN=CORD(I,2)
ENDIF
ENDIF
END DO
WRITE(*,*) ’YMIN = ’,YMIN,’YMAX = ’,YMAX
C
ZMIN=CORD(1,3)
ZMAX=CORD(1,3)
DO I=1,NNODEUP
IF (CORD(I,3).GT.ZMAX) THEN
ZMAX=CORD(I,3)
ELSE IF (CORD(I,3).LT.ZMIN) THEN
ZMIN=CORD(I,3)
ENDIF
END DO
WRITE(*,*) ’ZMIN = ’,ZMIN,’ZMAX = ’,ZMAX
XMID=(XMIN+XMAX)/2
YMID=(YMAX+YMIN)/2
ZMID=(ZMAX+ZMIN)/2
C**
NANG=INT(XMAX/XINC)
C**
XSTART=0.0
XSTOP=0.0
DO L=1,NANG
IF(XSTART.LT.XMID) THEN
XSTART=XSTART+XINC
ENDIF
END DO
XSTOP=XSTART+XINC
220
C
XS=XMIN
DO I=1,NNODEUP
IF(ABS(CORD(I,2)-YMIN).LE.0.1.AND.ABS(CORD(I,3)-
& ZMIN).LE.0.1) THEN
IF ((XSTART-CORD(I,1)).GE.0.0.AND.CORD(I,1).GE.XS) THEN
XS=CORD(I,1)
NS=I
ENDIF
ENDIF
END DO
CALL CS(NELEM,NS,NMAX,IDX,NEIN)
DO I=1,NMAX
IF (IDX(I).EQ.1) THEN
NS=IELE(NEIN(I),4)
ENDIF
ENDDO
XSTART=CORD(NS,1)
XSTOP=XSTART+XINC
c write(20,*) ’main program’
c write(20,*) XMAX,XMIN,YMAX,YMIN,ZMAX,ZMIN,XSTART,XINC
C
C*** SEARCH THE DOMAIN FOR ALL NODES GREAT THEN XSTART AND LESS
C*** THAN XSTOP
C
C SERACH THE ELEMENTS CONTAINING NODE (NS)
C
NODEX=0
DO I=1,NNODEUP
IF((CORD(I,1)-XSTART).GE.-20.0.AND.(CORD(I,1)-XSTOP)
1 .LE.20.0) THEN
NODEX=NODEX+1
FNODE(NODEX)=I
ENDIF
END DO
C
FLORAT=1E9
GRIDX=8.0
AREAO=ABS((ZMAX-ZMIN)*(YMAX-YMIN))
WRITE(30,’(A)’) ’Write the Poincare map to a file:’
WRITE(30,’(2A20)’) ’CORD Y’,’CORD Z’
DO I=1,10
STRP(I,1)=XSTART
221
STRP(I,2)=YMIN+(YMAX-YMIN)/11*I
c STRP(I,2)=YMID+ABS(YMAX-YMID)/10*(I-1)
c STRP(I,2)=6.75E+02
c STRP(I,3)=ZMID
STRP(I,3)=50.0
END DO
WRITE(20,*) ’ START POINT’
WRITE(20,’(3F16.6)’) (STRP(I,1),STRP(I,2),STRP(I,3),I=1,10)
C
C STREAKLINE dye injected from fixed position for period of time
C - tracer of dye show the fluid flow.
C initial position (x0,y0,z0),t=0, find the path (x(t),y(t),z(t))
C motion of a particle is give by dx/dt=vx; dy/dt=vy;dz/dt=vz
C
C==================================================================
C*** START TO GET PARTICLE TRACKS ***
C==================================================================
POIN=0
EPSILON = 2.0
C GIVET DELTA T (TIME STEP)
C
DO 200 I=1,10
C CALL INTERPOLATION ROUTINE TO CALCULATE THE VLOCITY
C TO COMPUTE THE VELOCITY FOR NEXT POINT
C
C Define the time step for x=vx*dt
C
C INITIALIZING START POINT OF STREAKLINE
DO L=1,3
STRLIN(1,L)=STRP(I,L)
END DO
C
DO 300 J=1,200
NSTR=1
DT=GRIDX/(FLORAT/AREAO)
C
C GET THE VELOCTIY OF THE START POINT
DO K=1,3
PT(K)=STRLIN(NSTR,K)
ENDDO
CALL PS(NODEX,FNODE,NS,PT)
CALL CELIDX(NELEM,NMAX,IDX,NEIN,CELLIDX)
CALL ITP(NMAX,CELLIDX,VN,PT,FLAG)
222
IF (FLAG) GOTO 333

IF (NMAX.EQ.8) THEN
IF (CORD(NS,3).GT.PT(3)) THEN
NEWNS = IELE(CELLIDX(1,1,1),2)
CALL CS(NELEM,NEWNS,NMAX,IDX,NEIN)
ELSE
ENDIF
ENDIF
333 if(abs(vn(1))+abs(vn(2))+abs(vn(3)).eq.0.0) GOTO 330
C
C*** GET THE PARTICLE TRACK
C RESTORE DELTA T
DO 350 WHILE (STRLIN(NSTR,1).LE.XSTOP)
C
C*** USE OF 4TH-ORDER RUNGE-KUTTA TO INTEGRATE VELOCITY
C
NSTR = NSTR+1
C
210 CALL RK4V(NODEX,FNODE,NSTR,STRLIN,VN,DT)
DO L=1,3
XN1(L)=STRLIN(NSTR,L)
ENDDO
DT = DT/2
C
C CALL RK4V(NDD,NSTR,STRLIN,VN,DT)
CALL RK4V(NODEX,FNODE,NSTR,STRLIN,VN,DT)
C
DX = 0.0
DO L=1,3
DX=DX+(XN1(L)-STRLIN(NSTR,L))**2
ENDDO
DX=SQRT(DX)
IF (DX.GT.0.0.AND.DX.LT.EPSILON/2.0) THEN
DT=DT*4
GOTO 210
ENDIF
C
DO 400 WHILE (DX.GE.EPSILON)
223
C
C
DO L=1,3
XN1(L)=STRLIN(NSTR,L)
ENDDO
C
DT = DT/2
DX = 0.0
DO L=1,3
DX = DX + (XN1(L)-STRLIN(NSTR,L))**2
ENDDO
DX=SQRT(DX)
400 CONTINUE
C
C CALL VVEC(NODEX,FNODE,NSTR,STRLIN,VN)
DO K=1,3
PT(K)=STRLIN(NSTR,K)
ENDDO
CALL PS(NODEX,FNODE,NS,PT)
IF (FLAG) GOTO 334
IF (NMAX.EQ.8) THEN
IF (CORD(NS,3).GT.PT(3)) THEN
ELSE
ENDIF
ENDIF
C
334 if(abs(vn(1))+abs(vn(2))+abs(vn(3)).eq.0.0) GOTO 330
C
350 CONTINUE
C
IF (ABS(STRLIN(NSTR,1)-XSTOP).LE.0.01) THEN
224
ELSE
DO L=0,2
XP(L)=STRLIN(NSTR-3+L,1)
YP(L)=STRLIN(NSTR-3+L,2)
ZP(L)=STRLIN(NSTR-3+L,3)
END DO
STRLIN(NSTR,1)=XSTOP
STRLIN(NSTR,2)=P(XP,YP,2,XSTOP)
STRLIN(NSTR,3)=P(XP,ZP,2,XSTOP)
END IF
C
WRITE(*,’(A,$)’) ’ CYCLING...’
WRITE(*,’(4X, I, 4X, I)’) J,NSTR
C
C WRTIE TO POINCARE MAP
C
POIN=POIN+1
POICAR(POIN,1)=STRLIN(NSTR,2)
POICAR(POIN,2)=STRLIN(NSTR,3)
WRITE(20,’(A,$)’) ’Write the streakline for point (x,y,z):’
WRITE(20,’(1P3E16.6)’) (STRP(I,L),L=1,3)
WRITE(20,’(3A20)’) ’CORD X’,’CORD Y’,’CORD Z’
C+++ WRITE THE STREAKLINE TO FILE
C DO L1=1,NSTR
C WRITE(20,’(1P3E20.6)’) (STRLIN(L1,L2),L2=1,3)
C END DO
C RESAMPLE THE START POINT
STRLIN(1,1)=XSTART
STRLIN(1,2)=STRLIN(NSTR,2)
STRLIN(1,3)=STRLIN(NSTR,3)
C+++
C
300 CONTINUE
C+++ WRITE THE STREAKLINE TO FILE
330 WRITE (30,’(A,$)’) ’POINCARE MAP, POINT: ’
WRITE (30,’(2X,I8)’) I
DO L1=1,POIN
WRITE(30,’(1P2E20.6)’) (POICAR(L1,L2),L2=1,2)
END DO
C+++
POIN=0
200 CONTINUE
C**
TODAY=FDATE()
225
PRINT *, TODAY
WRITE(30,*) TODAY
C**
CLOSE(20)
CLOSE(30)
C**
STOP
END

Passive Mixing in Microchannels - Thesis - 2004

Uploaded by

Copyright:

Available Formats

Passive Mixing in Microchannels - Thesis - 2004

Uploaded by

Document Information

Original Description:

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Passive Mixing in Microchannels - Thesis - 2004

Uploaded by

Copyright:

Available Formats

Passive Mixing in Microchannels with Geometric

B.S. (Shanghai Jiao Tong University) 1990

A thesis submitted to Industrial Research Institute Swinburne (IRIS) in

SWINBURNE UNIVERSITY OF TECHNOLOGY,

Hengzi Wang September 19, 2004

Generally, before performing chemical or biological analysis, samples and reagents

Three types of geometric variations were researched. First, partitioning in channels

I visited microfluidic diagnostics laboratory in Stanford University, and I wish to thank

List of Tables xviii

List of Variables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xix

List of abbreviations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxii

1.1 Background and Significance of the Research . . . . . . . . . . . . . . . . 1

1.2 Statement of the Problem . . . . . . . . . . . . . . . . . . . . . . . . . . . 2

1.3 Objectives and Scope of the Research . . . . . . . . . . . . . . . . . . . . 3

1.4 Outline of this thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

2.3 Fundamental Theories for Microfluidics . . . . . . . . . . . . . . . . . . . 11

2.3.2 Conservation of Momentum . . . . . . . . . . . . . . . . . . . . . 13

2.3.3 Brownian Motion . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

2.3.4 Taylor-Aris dispersion . . . . . . . . . . . . . . . . . . . . . . . . 15

2.3.5 Chaotic Advection . . . . . . . . . . . . . . . . . . . . . . . . . . 16

2.4 Generation of Fluid Motion in a Microfluidic Device . . . . . . . . . . . . 19

2.4.1 Mechanically Actuated Micropumps . . . . . . . . . . . . . . . . . 20

2.4.2 Electrokinetic Micropumps . . . . . . . . . . . . . . . . . . . . . . 23

2.5 Active Micro-Mixers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26

2.5.1 Planar Laminar Bubble Mixer . . . . . . . . . . . . . . . . . . . . 27

2.5.2 Periodic Perturbation Applied to a Main Fluid Stream . . . . . . . 27

2.5.3 Peristaltically Driven Micromixing . . . . . . . . . . . . . . . . . 29

2.5.4 Magnetohydrodynamic (MHD) Micro-mixers . . . . . . . . . . . . 30

2.5.5 Electrohydrodynamic (EHD) Micro-mixers . . . . . . . . . . . . . 32

2.5.6 Magnetic Micro-Stirrer . . . . . . . . . . . . . . . . . . . . . . . . 32

2.6 Passive Microfluidic Mixers . . . . . . . . . . . . . . . . . . . . . . . . . 36

2.6.1 T-type mixers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36

2.6.2 Splitting and Recombining . . . . . . . . . . . . . . . . . . . . . . 41

2.6.3 Focusing Fluid Streams . . . . . . . . . . . . . . . . . . . . . . . . 46

2.6.4 Inject Small Fluid Substreams into a Main Stream . . . . . . . . . 47

2.6.5 Passive Microfluidic Chaotic Mixers . . . . . . . . . . . . . . . . . 48

3 Numerical Modelling and Experimental Methodology 57

3.2 Fundamentals of Computational Fluid Dynamics (CFD) . . . . . . . . . . . 59

3.2.2 CFD Solvers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60

3.3 CFD Modelling in Microfluidics . . . . . . . . . . . . . . . . . . . . . . . 63

3.3.1 System of Units in Microfluidics . . . . . . . . . . . . . . . . . . . 64

3.3.2 Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . . . . 64

3.3.3 Building Numerical Models for Micro-Flows . . . . . . . . . . . . 65

3.3.4 Simulation Results Analysis . . . . . . . . . . . . . . . . . . . . . 71

3.4 Experimental Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73

3.4.1 Excimer Laser Micro-machining . . . . . . . . . . . . . . . . . . . 74

3.4.2 Microfluidic Packaging . . . . . . . . . . . . . . . . . . . . . . . . 75

3.4.3 Microfluidic Mixing Visualization . . . . . . . . . . . . . . . . . . 75

4 Understanding Microfluidic Mixing 79

4.2 Flow at Low Reynolds Number . . . . . . . . . . . . . . . . . . . . . . . . 80

4.3 Taylor-Aris Dispersion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85

4.3.2 Dispersion of a Solute in Solvent Flowing Slowly Through a Rect-

4.3.3 Applying Taylor-Aris Dispersion to Enhance Mixing . . . . . . . . 89

4.4 Two Important Parameters for Designing a Microfluidic Mixer . . . . . . . 91

4.4.1 Péclet number . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92

4.4.2 Aspect Ratio of a Rectangular Microchannel . . . . . . . . . . . . 95

5 Microfluidic Mixing by Partitioning Flows 101

5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101