- Email: [email protected]
Electroviscous flow through nanofluidic junctions
Accepted Manuscript Electroviscous flow through nanofluidic junctions J.D. Berry, M.R. Davidson, D.J.E. Harvie PII: DOI: Reference:
S0307-904X(14)00075-4 http://dx.doi.org/10.1016/j.apm.2014.02.018 APM 9876
To appear in:
Appl. Math. Modelling
Received Date: Revised Date: Accepted Date:
22 July 2013 28 January 2014 12 February 2014
Please cite this article as: J.D. Berry, M.R. Davidson, D.J.E. Harvie, Electroviscous flow through nanofluidic junctions, Appl. Math. Modelling (2014), doi: http://dx.doi.org/10.1016/j.apm.2014.02.018
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Electroviscous flow through nanofluidic junctions
J. D. Berry a,1 , M. R. Davidson a , D. J. E. Harvie a,∗
a Department
of Chemical and Biomolecular Engineering, University of Melbourne, Melbourne, Victoria 3010, Australia
Abstract The steady-state electroviscous (pressure-driven) flow of an aqueous solution of potassium chloride through a nanofluidic borosilicate glass junction is modelled using a finite volume approach. Electric-field propagation through the glass section of the junction is included in the analysis. Channel half-widths of 20 − 70 nm and salt concentrations of 10−6 − 10−1 M are considered. Two types of common nano/microfluidic flow through the junction are analysed: a separating flow and a mixing flow. The hydrodynamic (pressure) and electrokinetic (electric potential) resistances of the junction are quantified in terms of equivalent lengths. For both flows, the hydrodynamic equivalent length of a 90◦ junction is found to be between ∼ 0.85 and ∼ 1.05 channel half-widths, while the electrokinetic equivalent length varies from ∼ 0.3 to ∼ 0.6 channel half-widths. Decreasing the junction angle decreases the equivalent lengths of the junction.
Preprint submitted to Elsevier Science
6 March 2014
1
Introduction
Nanofluidic and microfluidic lab-on-a-chip (LOC) devices can be utilised to perform many applications, including the analysis and processing of biological fluids [1, 2, 3]. Many LOC devices use pumps to drive flow through networks of channels. Nano/microfluidic pressure-driven flow can be subject to significant electrokinetic effects; depending on the channel dimensions, the conductivity of the fluids present and the materials used to manufacture the device. An effect known as electroviscosity describes the increased flow resistance due to the presence of mobile charge carriers within a fluid flowing through a nano/microfluidic device with charged walls [4]. In order to optimise the design of LOC devices, quantitative understanding of the pressure and potential changes occurring across network elements due to electroviscosity is required. Our group has previously characterised these changes in microfluidic contractions and nanofluidic bends [5, 6, 7, 8]. Another common component in LOC devices is a junction, used to passively mix multiple species [9], and to sort drops and cells based on physical properties [10, 11]. The main difficulty of mixing in LOC devices is that diffusion, the dominant mechanism for mixing at the nano and micro scales, is an inherently slow process and requires long channel lengths. As a result, large pressure drops can develop, potentially leading to problems integrating these ∗ Corresponding author. Tel.: +61 3 8344 6421 Email addresses: [email protected] (J. D. Berry), [email protected] (M. R. Davidson), [email protected] (D. J. E. Harvie). 1 Present address: CSIRO Computational Informatics, Clayton, Victoria 3169, Australia
2
components onto a LOC device [12]. To determine if the use of a junction is a feasible approach to mixing or sorting on a particular LOC device, simple methods of calculating changes in pressure and potential over junctions are required. One such concept, commonly used for high Reynolds number pipe design, is equivalent lengths based on changes in pressure and potential due to the junction. In this study, a finite volume numerical method is used to simulate two types of flow through a junction. The effects of the channel dimensions and the salt concentration on the performance of the junction are quantified using equivalent lengths.
Nomenclature
A
measured quantity (either U or P )
c+
Positive ion concentration
c−
Negative ion concentration
c0
Mean molar concentration
D
ion diffusivity (D = D+ = D− )
∆Ajunction
change in quantity A due to junction
dA ds
gradient of quantity A
dA ds FD
fully-developed gradient of quantity A
3
e
elementary charge
FD
fully-developed flow
K
dimensionless inverse Debye length
k
Boltzmann constant
Lin
length of inlet channel
Lout
length of outlet channel
Leq.,P
equivalent length based on pressure drop
Leq.,U
equivalent length based on potential increase
n+
positive ion number density
n−
negative ion number density
n0
geometric mean ion number density
nwall
surface normal to wall, pointing into solid region
p
pressure
Pe
Peclet number
s
coordinate pointing in flow direction (parallel to wall)
Re
Reynolds number
T
temperature
4
Uf
fluid potential
Uw
solid potential
v
velocity
V¯
primary channel mean flow velocity
W
channel half-width
x
position
z
ion valence
ε0
permittivity of free space
εf
fluid permittivity
εw
solid permittivity
µ
dynamic viscosity
ρ
density
ρE
charge density
σ
wall surface charge
θ
junction angle
Ωf
fluid domain
Ωw
solid domain
5
a) Separating flow
b) Mixing flow
θ Ωw
Glass
Secondary channel
Ωf
Ls
2W
Glass
Primary channel
3
Lp
W
c)
W
Fig. 1. a) - b) Schematic of flow scenarios considered. c) Schematic of computational domain. The dashed-dotted line denotes the line of symmetry used in the simulations.
2
Model Formulation
The geometry considered in this study is a two-dimensional nanofluidic junction. Two types of flow are considered: a separating flow, whereby the flow 6
through one channel is split into two equal flows, and a mixing flow, where equal flow through two channels is combined into one flow (Figure 1a),b)). Both flows are assumed to be steady-state. The channel contains a symmetric 1:1 electrolytic fluid which is driven through the junction by a pressure gradient. The liquid is assumed to be an aqueous solution of potassium chloride with constant density, viscosity and permittivity; and anions and cations of equal and opposite valence and equal diffusivity.
In general, the surface charge density of silica varies in response to local ion concentrations, in particular the proton concentration (pH) which is itself influenced by other chemical species present in the liquid (e.g. buffers) [13, 14]. As out focus is on the electroviscous performance of a simple binary electrolyte, we neglect these more complex chemical effects and assume that the surface charge density along the channel walls is constant and uniform, and approximately equal to the surface charge density of glass at a local pH= 7 (-10 mC/m2 ) [13, 15]. Despite the simplification, previous studies have shown that such a constant wall charge condition can accurately model flow through glass nanochannels, albeit under electro-osmotic flow conditions [16, 17].
The physical properties used in the simulation are representative of typical microfluidic flows past borosilicate glass (Table 1) [8]. Also included in the model is the solid glass section of the junction material, because electricfield propagation through the walls of nano/microfluidic components has been shown to have a significant effect on the flow dynamics of these systems [6]. 7
Table 1 Physical properties used in the model. (Based on flow of an aqueous solution of potassium chloride past borosilicate glass.) Notes: a [18], Property
[19],c [15],d [20] Value
Temperature (K)
298.0
Viscosity (10−4 Pa s)
8.91a
Density (kg m−3 )
1000.0a
Relative fluid permittivity (-)
80.1a
Relative solid permittivity (-)
4.6a
Channel half-width (nm)
20.0 – 70.0
Concentration (M)
10−6 – 0.1
Ion diffusivity (10−9 m2 s−1 ) Wall charge density (mC m−2 ) Primary channel flow velocity (mms−1 )
2.1
b
2.0b −10.0c 5.0d
Governing Equations
For nanochannels of size greater than 5 nm, the ion transport and fluid flow can be analysed using continuum dynamics [21]. The continuum equations governing the steady flow in the fluid region Ωf are given by the continuity equation:
∇ · v = 0,
(1) 8
the momentum equation:
h
i
∇ · (ρvv) = −∇p + ∇ · µ ∇v + (∇v)T − ρE ∇Uf ,
(2)
the ion-transport equation for the anions and cations present in the fluid:
z± eD± n± ∇ · (n± v) = ∇ · D± ∇n± + ∇Uf , kT
(3)
and the Poisson equation relating the electric potential to the local charge density:
∇ · εf ∇Uf = −
X ezi ni ρE =− . ε0 i=+,− ε0
(4)
In the solid region, the electric potential satisfies ∇2 Uw = 0.
2.2
(5)
Boundary Conditions
The computational domain is shown in Figure 1c). Across the fluid/solid interface between regions Ωf and Ωw the electric field normal to the wall undergoes a jump, given by
εf ε0 ∇Uf
wall
· nwall − εw ε0 ∇Uw
wall
· nwall = σ,
(6)
where nwall is the unit normal to the wall pointing into the solid region. Across the fluid/solid interface the electric potential is continuous, meaning
Uf
wall
= Uw 9
wall
.
(7)
Along the external walls of the device, the gradient of electric potential is constrained in the direction normal to the wall by
∇Uw
wall
· nwall = 0,
(8)
representing a boundary with zero surface charge and no field propagation. The fluid/solid interface was considered to be clean, smooth and hydrophilic (appropriate for water on clean glass), and as such was modelled as a no-slip boundary [22, 23, 24]. A zero ion flux normal to the interface was also imposed. Symmetry conditions were imposed on all quantities along the centreline of the junction. At the inlet and outlet of the channel, the velocity gradient normal to the flow and the tangential velocity component were set to zero. At the √ inlet, a flow-rate and the geometric mean ion number density n0 = n+ n− (calculated from the inlet molar concentration c0 ) were specified, along with a zero normal gradient condition on the counter-ion concentration. At the outlet the normal gradient of both ion concentrations was specified as zero. The total current through the channel was set to zero, the necessary condition for electroviscous flow.
2.3
Numerical Method
The governing equations Equations 1-5 were solved in conjunction with the requisite boundary conditions on a composite structured/unstructured mesh generated with the program gmsh [25]. To ensure fully-developed flow at the inlet and outlet of the junction, Lin = Lout = 22.5W . The fluid region Ωf was meshed using rectangular elements, with 50 elements per channel half-width W . When the junction angle was other than 90◦ , a small region of triangular 10
a)
b)
Fig. 2. Coarse representation of the mesh for a) a 90◦ junction and b) a 45◦ junction. The dashed white lines delineate the boundary between the fluid region Ωf , and the solid region Ωw . Note that only a portion of the mesh used in the simulations is shown.
elements was used (Figure 2b)). A typical mesh used in the simulation had ∼ 3.1 × 105 elements. The distribution of elements across the channel is nonuniform (Figure 2), with fewer elements near the centreline and more elements near the wall, in order to capture the high gradients of electric potential and salt concentrations in the electric double layer adjacent to the wall. The open-source multi-physics solver arb[26] was used to solve the system of equations and boundary conditions. In the context of electrokinetic and nano/microfluidic flows, arb has been used to calculate the electrokinetic flow of a binary electrolyte through parallel nanochannels[27] as well as the interfacial mass transfer and reaction in Y-Y microfluidic channels[28, 29, 30]. arb discretises transport equations using an implicit finite-volume method, and solves the resulting numerical system via a backstepped, Newton Raphson method. Interpolation and differentiation of cell-centred values to element faces is accomplished via a variation of the Moving Least Squares method, as 11
detailed in [31]. For this study 2nd order polynomials were used as the basis for this method. Velocity interpolation to element faces references the local pressure and potential fields evaluated at the face, in the spirit of Rhie-Chow interpolation [32]. Ion concentrations are advected using a high-order, limited upwind method. Input files which fully detail the discretisation technique for the flow of a binary electrolyte through a charged channel, and upon which the present simulations are based, are available for download 2 . The matrix at each Newton-Raphson iteration is inverted using the sparse linear solver PARDISO that is packaged within the Intel Math Kernel Library. The non-dimensional convergence criterion employed was = 10−10 [26].
Testing of the method has shown that the pressure and potential gradients for fully-developed electrokinetic channel flow on a structured mesh are second order accurate (that is, the error decreases proportionally to the average element dimension squared). For an unstructured mesh, used only within the junction region of the non 90◦ junctions, the order of accuracy order is more difficult to quantify, but significantly, the pressure and potential gradient errors for fully-developed channel flow do decrease consistently with mesh refinement. For the mesh resolution used in this study, mesh sensitivity testing indicates that the reported potential and pressure equivalent lengths have errors of less than 1% over the entire parameter range studied.
2
http://people.eng.unimelb.edu.au/daltonh/downloads/arb/code/
latest/examples/electrokinetic channel flow/
12
3
Results
Figure 3a) and b) depict the co-ion distribution for both the separating and mixing flows through a junction of channel half-width W = 50 nm and geometric ion ion number density n0 corresponding to molar concentration c0 = 10−4 M. Electric double layers are present, with a depletion of co-ion concentration near the junction walls. For the physical properties used in the simulations the Peclet number is low (Pe = V¯ W/D = 0.125). Hence diffusion and conduction are the dominant means of ion transport, meaning the size of the double layer is largely unaffected by the presence of the junction and the direction of the flow, and that symmetry of the double layer is present about the leading and trailing edge of the junction corner.
Figure 3c) and d) show the contours of velocity magnitude and streamlines. Because of the apparently symmetric nature of the ion distribution about the junction and the fact that the Reynolds number is extremely low (Re = ρV¯ W/µ = 2.8 × 10−3 ), the velocity distribution at the junction could be expected to be insensitive to the flow direction. However, there are subtle differences between the mixing and the separating flow, clearly highlighted in the contours of geometric mean ion concentration (Figure 3e) and f)). The low but finite Pe number effect means that the distribution of ions is biased very slightly towards the outlet of the junction. Thus, due to the non-linear electric body-force term present in Equation 2, the flow through the teepiece is irreversible. The extent of irreversibility is quantified below.
Figure 3g) and h) compare the potential contours for separating and mixing flow. The junction has an effect on the potential distribution in the region of 13
c− (mM) 0.000 0.008 0.016 0.024 0.032 0.040
24
24
22
22
-4
-2
0
2
4
a)
-4
-2
0
2
4
-2
0
2
4
-2
0
2
4
-2
0
2
4
b) |u| (mm/s) 0
1
2
3
4
5
6
7
8
24
24
22
22
-4
-2
0
2
4
c)
-4
d) 0.096 0.097 0.098 0.099 0.100 0.101
c0 (mM)
24
24
22
22
-4
-2
0
2
4
e)
-4
f) U (mV) -130 -114 -97 -81 -65 -49 -32 -16
0
24
24
22
22
-4
g)
-2
0
2
4
-4
h)
Fig. 3. Distribution of co-ions for a) separating flow, and b) mixing flow. Contours of velocity magnitude and streamlines in the junction for c) separating flow, and d) mixing flow. Distribution of geometric mean ion number density for e) separating flow, and f) mixing flow. Contours of electric potential for g) separating flow, and h) mixing flow. The geometric ion number density n0 corresponds to molar concentration c0 = 10−4 M, the channel half-width W = 50 nm, and the axes have been normalised with length-scale W . All physical properties are listed in Table 1.
the geometry change for both mixing and separating flows. The gradient of potential along the channel is opposite for the mixing flow in comparison to the separating flow, due to the change in flow direction. The changes in pressure and potential caused by the presence of the junction can be quantified using [8] ∆Ajunction
dA = ∆A − Li , ds FD,i i=in,out X
(9)
where A is either the pressure or the potential. This definition is valid only if ∆A is measured over a length beginning and ending in fully-developed flow regions in the inlet and outlet respectively. The pressure drop and potential increase caused by the presence of a junction with separating flow are given in Figure 4a) for a range of channel half-widths and salt concentrations. Over the range of concentrations and channel half-widths modelled, the pressure drop due to the presence of the junction varies by almost an order of magnitude. The magnitude of the pressure drop increases with decreasing channel half-width W . Figure 4 compares the pressure drop of a junction with W = 50 nm for both mixing and separating flows. Due to the distribution dependence of the ions on the flow direction the flow is irreversible, and the pressure drop caused by the junction also depends on the flow direction. However, the flow is only weakly irreversible, and the pressure drop for mixing flow is higher than the separating flow value by only ∼ 2 Pa (< 1%) over the range of concentrations studied. There is also significant variation in the potential increase due to the junction (Figure 5). For low concentrations, the increases in potential due to the presence of the junction are constant, but dependent on the size of the channel half-width W . For very high concentrations the increase in potential due 15
-100 -200
∆P (Pa)
-300 -400
W = 70 nm W = 60 nm W = 50 nm W = 40 nm
-500 W = 30 nm -600 -700 -800 -6 10
W = 20 nm -5
10
a)
-4
-3
10 10 Concentration (M)
-2
10
-1
10
-280
∆P (Pa)
-285
-290
-295
-300
-305 -6 10
b)
Separating flow Mixing flow No electrokinetics -5
10
-4
-3
10 10 Concentration (M)
-2
10
-1
10
Fig. 4. The pressure drop due to the junction as a function of average salt concentration for: a) separating flow over a range of channel half-widths, and b) separating and mixing flow for channel half-width W = 50 nm. All physical properties are listed in Table 1. Each symbol represents a simulation result.
W= W= W= W= W= W=
0.6 0.5
20 nm 30 nm 40 nm 50 nm 60 nm 70 nm
∆U (mV)
0.4 0.3 0.2 0.1 0 -6 10
-5
10
a)
-4
-3
10 10 Concentration (M)
-2
10
-1
10
Separating flow Mixing flow
0.6
∆U (mV)
0.5 0.4 0.3 0.2 0.1 0 -6 10
b)
-5
10
-4
-3
10 10 Concentration (M)
-2
10
-1
10
Fig. 5. a) The potential increase due to the junction as a function of average salt concentration for: a) separating flow over a range of channel half-widths, and b) separating and mixing flow for channel half-width W = 50 nm. All physical properties are listed in Table 1. Each symbol represents a simulation result.
to the junction drops to zero. The potential increase due to the junction for mixing flow is higher than that for separating flow (Figure 5b)). Two useful quantities can be introduced to explain the behaviour of the pressure drop, and potential increase, caused by the presence of the junction. The first quantity is the equivalent length of the junction based upon the pressure drop or the potential increase. The equivalent length is defined as the length of uniform inlet channel (with half-width W) needed to match the pressure drop (potential increase) caused by the effect of the junction [7, 8]. Formally, these equivalent lengths are defined as Leq.,A =
∆Ajunction , dA ds FD,primary
(10)
where A is either the pressure or the potential. The second quantity is the dimensionless inverse Debye length, defined as s
K=
2z 2 e2 n0 W 2 . εf ε0 kT
(11)
This parameter can also be interpreted as the square root of a dimensionless concentration. When K is small, the thickness of the electric double layer is greater than the channel dimension, and significant double-layer overlap occurs. At this limit, there is a nearly constant amount of charge (opposite and equal to the fixed surface charge) across the channel due to the exclusion of co-ions, and no region of electrically neutral fluid. When K is large, the double layer thickness is small in comparison with the channel dimension, and all the charge present in the flow is located in a thin region adjacent to the wall. As a result, the bulk flow is electrically neutral. When the pressure drops caused by the junction given in Figures 4 and 5 are 18
reformulated in terms of equivalent lengths and dimensionless inverse Debye lengths K, there is significant collapse of the data (Figures 6 and 7). Over the parameter space simulated, the equivalent length based on the pressure drop varies from ∼ 0.85W to ∼ 1.05W , while the equivalent length based on the potential increase varies from ∼ 0.3W to ∼ 0.6W . For values of K < 1 both equivalent lengths are weakly dependent upon the channel half-width W ; but independent of K due to the exclusion of co-ions in the channel, and the consequent fixed charge therein. As the double layer thickness becomes comparable to the channel dimension however, the resistance of the junction begins to vary with K. As K increases, so does Leq.,P . For values of K > 10 the fluid bulk is electrically neutral, and as a consequence the equivalent length based on the pressure is the same as the equivalent length calculated for purely hydrodynamic flow (found by an additional simulation). Interestingly, for values of K > 3 and W > 50 nm, the equivalent length based on the pressure is independent of W . The equivalent length of the junction based on the potential displays similar behaviour for K < 1. For increasing values of K above 1, Leq.,U decreases until a minimum resistance is reached at K ∼ 5. As K increases further, the double layer decreases in thickness, and electrokinetic effects decrease. Consequently, the potential increase caused by the junction (numerator of Equation 10), and the potential increase per unit length of uniform channel (denominator of Equation 10) go to zero as K → ∞. However, upon inspection of the data, the potential increase caused by the junction goes to zero more slowly than the potential increase per unit length of uniform channel, resulting in Leq.,U increasing as the double layer becomes thinner. The effect of flow direction on the resistance of the junction to the flow is 19
1.1
Leq.,P
1
W = 20 nm W = 30 nm W = 40 nm W = 50 nm W = 60 nm W = 70 nm No electrokinetics
0.9
0.8 0.1
1
a)
10 K
1.1
Separating flow Mixing flow No electrokinetics Leq.,P/W
1
0.9
0.8 0.1
b)
1
10 K
Fig. 6. Equivalent length of the junction based on the pressure drop over the junction as a function of dimensionless inverse Debye length K for a) separating flow for a range of channel half-widths, and b) separating flow and mixing flow for channel half-width W = 50 nm.
0.7
0.6
Leq.,U/W
0.5
0.4
0.3
0.2
W= W= W= W= W= W=
20 nm 30 nm 40 nm 50 nm 60 nm 70 nm
0.1
1
a)
10 K
0.7
Separating flow Mixing flow
Leq.,U/W
0.6
0.5
0.4
0.3 0.1
b)
1
10 K
Fig. 7. Equivalent length of the junction based on the potential increase over the junction as a function of dimensionless inverse Debye length K for a) separating flow for a range of channel half-widths, and b) separating flow and mixing flow for channel half-width W = 50 nm.
1
Separating flow Mixing flow
Leq.,P/W
0.75
0.5
0.25
0 0
10
20
30
40 50 60 Junction angle θ
70
80
90
Fig. 8. Equivalent length of the junction based on the pressure drop over the junction as a function of junction angle θ for separating flow and mixing flow. The channel half-width W = 50 nm, and the geometric ion number density n0 corresponds to molar concentration c0 = 3.8 × 10−5 M (K = 1).
minor, with the mixing flow exhibiting slightly higher equivalent lengths due to the pressure drop (Figure 6b)), and the potential increase (Figure 7b)), in comparison to the separating flow. The effect of junction angle θ is shown in Figure 8 for mean ion number density n0 corresponding to molar concentration c0 = 3 .8 × 10−5 M (K = 1). The resistance of the junction decreases quasilinearly with decreasing junction angle. The effect of flow direction is again very small, with the resistance of separating flow and mixing flow differing by less than ∼ 4%. 22
4
Conclusions
The steady-state pressure-driven flow of an aqueous solution of potassium chloride through a junction has been studied for a range of channel dimensions and salt concentrations. Two types of flow through the junction are compared: a separating flow, commonly used in nano/microfluidic devices to sort drops and cells, and a mixing flow, used to mix multiple species and reagents. The resistance of the junction to the each flow has been quantified in terms of equivalent lengths based upon the pressure drops, and potential increases, caused by the junction. Over the parameter space simulated, the equivalent length based on the pressure drop varies from ∼ 0.85W to ∼ 1.05W , while the equivalent length based on the potential increase varies from ∼ 0.3W to ∼ 0.6W . For the parameters examined, the flow through the teepiece is weakly irreversible because of a slight readjustment of the double-layer downstream to the junction. As a consequence the effect of flow direction (mixing or splitting) on both equivalent lengths is minor, with differences of less than 1% observed. Decreasing the junction angle acts to decrease the resistance of the junction to the flow quasi-linearly.
Thus, for lab-on-a-chip design purposes, the excess pressure drop of a 90◦ junction can be estimated by multiplying the fully-developed one-dimensional electroviscous pressure gradient by 0.95 ± 0.1, and the excess potential increase can be approximated by multiplying the fully-developed one-dimensional electroviscous potential gradient by 0.45 ± 0.15 . 23
5
Acknowledgments
The authors thank Mr Christian Biscombe for enlightening discussions. This research was supported by the Australian Research Council Grants Scheme and by computational resources on the National Computational Infrastructure Facility through the National Computational Merit Allocation Scheme.
References [1] S. B. Cheng, C. D. Skinner, J. Taylor, S. Attiya, W. E. Lee, G. Picelli, D. J. Harrison, Development of a multichannel microfluidic analysis system employing affinity capillary electrophoresis for immunoassay, Analytical Chemistry 73 (7) (2001) 1472–1479. [2] A. Hatch, A. E. Kamholz, K. R. Hawkins, M. S. Munson, E. A. Schilling, B. H. Weigl, P. Yager, A rapid diffusion immunoassay in a T-sensor, Nature Biotechnology 19 (5) (2001) 461–465. [3] T. M. Squires, S. R. Quake, Microfluidics: fluid physics at the nanoliter scale, Reviews of Modern Physics 77 (3) (2005) 977–1026. [4] R. J. Hunter, Zeta potential in colloid science, Academic Press, London, 1981. [5] M. Davidson, D. Harvie, Electroviscous effects in low Reynolds number liquid flow through a slit-like microfluidic contraction, Chemical Engineering Science 62 (16) (2007) 4229–4240. [6] J. D. Berry, M. R. Davidson, D. J. E. Harvie, Electrokinetic development length of electroviscous flow through a contraction, ANZIAM Journal 52 (2011) C837–C852. [7] J. D. Berry, M. R. Davidson, R. P. Bharti, D. J. E. Harvie, Effect of wall 24
permittivity on electroviscous flow through a contraction, Biomicrofluidics 5 (2011) 044102. [8] J. D. Berry, A. E. Foong, C. Lade, C. J. C. Biscombe, M. R. Davidson, D. J. E. Harvie, Electroviscous resistance of nanofluidic bends, submitted to Chemical Engineering Science. [9] S. Wong, M. Ward, C. Wharton, Micro T-mixer as a rapid mixing micromixer, Sensors and Actuators B: Chemical 100 (3) (2004) 359–379. [10] K. Ahn, C. Kerbage, T. P. Hunt, R. M. Westervelt, D. R. Link, D. A. Weitz, Dielectrophoretic manipulation of drops for high-speed microfluidic sorting devices, Applied Physics Letters 88 (2006) 024104. [11] B. Ahn, K. Lee, R. Panchapakesan, K. W. Oh, On-demand electrostatic droplet charging and sorting, Biomicrofluidics 5 (2) (2011) 024113. [12] A. A. S. Bhagat, E. T. K. Peterson, I. Papautsky, A passive planar micromixer with obstructions for mixing at low Reynolds numbers, Journal of Micromechanics and Microengineering 17 (5) (2007) 1017–1024. [13] L. H. Allen, E. Matijevic, Stability of colloidal Silica. II. Ion exchange, Journal of Colloid and Interface Science 33 (3) (1970) 420–429. [14] P. J. Scales, F. Grieser, T. W. Healy, L. R. White, D. Y. C. Chan, Electrokinetics of the silica-solution interface: a flat plate streaming potential study, Langmuir 8 (3) (1992) 965–974. [15] R. B. Schoch, H. van Lintel, P. Renaud, Effect of the surface charge on ion transport through nanoslits, Physics of Fluids 17 (10). [16] D. Stein, M. Kruithof, C. Dekker, Surface-charge-governed ion transport in nanofluidic channels, Physical Review Letters 93 (3) (2004) 1–4. [17] D. J. E. Harvie, C. J. C. Biscombe, M. R. Davidson, Microfluidic circuit analysis I: Ion current relationships for thin slits and pipes, Journal of Colloid and Interface Science 365 (2012) 1–15. 25
[18] H. Ellis (Ed.), Book of Data, 4th Edition, Longman Group Limited, Harlow, 1986. [19] W. M. Haynes (Ed.), CRC Handbook of Chemistry and Physics (Internet version), 91st Edition, CRC Press/Taylor and Francis, Florida, 2011. [20] O. J¨annig, N.-T. Nguyen, A polymeric high-throughput pressure-driven micromixer using a nanoporous membrane, Microfluidics and Nanofluidics 10 (3) (2010) 513–519. [21] H. Daiguji, Ion transport in nanofluidic channels, Chemical Society Reviews 39 (2010) 901–911. [22] C. Honig, W. Ducker, Thin film lubrication for large colloidal particles: experimental test of the no-slip boundary condition, Journal of Physical Chemistry C 111 (44) (2007) 16300–16312. [23] C. Honig, W. Ducker, No-slip hydrodynamic boundary condition for hydrophilic particles, Physical Review Letters 98 (2) (2007) 028305. [24] O. I. Vinogradova, A. V. Belyaev, Wetting, roughness and flow boundary conditions., Journal of Physics. Condensed Matter 23 (18) (2011) 184104. [25] C. Geuzaine, J.-F. Remacle, Gmsh: A 3-d finite element mesh generator with built-in pre- and post-processing facilities, International Journal for Numerical Methods in Engineering 79 (11) (2009) 1309–1331. [26] D. J. E. Harvie, An implicit finite volume method for arbitrary transport equations, ANZIAM Journal 52 (2012) C1126–C1145. [27] C. Biscombe, M. Davidson, D. Harvie, Electrokinetic flow in parallel channels: circuit modelling for microfluidics and membranes, Colloids and Surfaces A: Physicochemical and Engineering Aspects, 440 (2014) 63–73. [28] L. R. Mason, D. Ciceri, D. J. E. Harvie, J. M. Perera, G. W. Stevens, Modelling of interfacial mass transfer in microfluidic solvent extraction. Part I. Heterogeneous transport, Microfluidics and Nanofluidics, 14 (2013) 26
197–212. [29] D. Ciceri, L. R. Mason, D. J. E. Harvie, J. M. Perera, G. W. Stevens, Modelling of interfacial mass transfer in microfluidic solvent extraction. Part II. Heterogeneous transport with reaction, Microfluidics and Nanofluidics, 14 (1-2) (2013) 213–224. [30] D. Ciceri, L. R. Mason, D. J. E. Harvie, J. M. Perera, G. W. Stevens, Extraction kinetics of Fe (III) by di-(2-ethylhexyl) phosphoric acid using a Y–Y shaped microfluidic device, Chemical Engineering Research and Design, In press. [31] L. Cueto-Felgueroso, I. Colominas, X. Nogueira, F. Navarrina, M. Casteleiro, Finite volume solvers and moving least-squares approximations for the compressible Navier–Stokes equations on unstructured grids, Computer Methods in Applied Mechanics and Engineering 196 (4548) (2007) 4712–4736. [32] C. Rhie, W. Chow, Numerical study of the turbulent flow past an airfoil with trailing edge separation, AIAA Journal 21 (11) (1983) 1525–1532.
27
S0307-904X(14)00075-4 http://dx.doi.org/10.1016/j.apm.2014.02.018 APM 9876
To appear in:
Appl. Math. Modelling
Received Date: Revised Date: Accepted Date:
22 July 2013 28 January 2014 12 February 2014
Please cite this article as: J.D. Berry, M.R. Davidson, D.J.E. Harvie, Electroviscous flow through nanofluidic junctions, Appl. Math. Modelling (2014), doi: http://dx.doi.org/10.1016/j.apm.2014.02.018
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Electroviscous flow through nanofluidic junctions
J. D. Berry a,1 , M. R. Davidson a , D. J. E. Harvie a,∗
a Department
of Chemical and Biomolecular Engineering, University of Melbourne, Melbourne, Victoria 3010, Australia
Abstract The steady-state electroviscous (pressure-driven) flow of an aqueous solution of potassium chloride through a nanofluidic borosilicate glass junction is modelled using a finite volume approach. Electric-field propagation through the glass section of the junction is included in the analysis. Channel half-widths of 20 − 70 nm and salt concentrations of 10−6 − 10−1 M are considered. Two types of common nano/microfluidic flow through the junction are analysed: a separating flow and a mixing flow. The hydrodynamic (pressure) and electrokinetic (electric potential) resistances of the junction are quantified in terms of equivalent lengths. For both flows, the hydrodynamic equivalent length of a 90◦ junction is found to be between ∼ 0.85 and ∼ 1.05 channel half-widths, while the electrokinetic equivalent length varies from ∼ 0.3 to ∼ 0.6 channel half-widths. Decreasing the junction angle decreases the equivalent lengths of the junction.
Preprint submitted to Elsevier Science
6 March 2014
1
Introduction
Nanofluidic and microfluidic lab-on-a-chip (LOC) devices can be utilised to perform many applications, including the analysis and processing of biological fluids [1, 2, 3]. Many LOC devices use pumps to drive flow through networks of channels. Nano/microfluidic pressure-driven flow can be subject to significant electrokinetic effects; depending on the channel dimensions, the conductivity of the fluids present and the materials used to manufacture the device. An effect known as electroviscosity describes the increased flow resistance due to the presence of mobile charge carriers within a fluid flowing through a nano/microfluidic device with charged walls [4]. In order to optimise the design of LOC devices, quantitative understanding of the pressure and potential changes occurring across network elements due to electroviscosity is required. Our group has previously characterised these changes in microfluidic contractions and nanofluidic bends [5, 6, 7, 8]. Another common component in LOC devices is a junction, used to passively mix multiple species [9], and to sort drops and cells based on physical properties [10, 11]. The main difficulty of mixing in LOC devices is that diffusion, the dominant mechanism for mixing at the nano and micro scales, is an inherently slow process and requires long channel lengths. As a result, large pressure drops can develop, potentially leading to problems integrating these ∗ Corresponding author. Tel.: +61 3 8344 6421 Email addresses: [email protected] (J. D. Berry), [email protected] (M. R. Davidson), [email protected] (D. J. E. Harvie). 1 Present address: CSIRO Computational Informatics, Clayton, Victoria 3169, Australia
2
components onto a LOC device [12]. To determine if the use of a junction is a feasible approach to mixing or sorting on a particular LOC device, simple methods of calculating changes in pressure and potential over junctions are required. One such concept, commonly used for high Reynolds number pipe design, is equivalent lengths based on changes in pressure and potential due to the junction. In this study, a finite volume numerical method is used to simulate two types of flow through a junction. The effects of the channel dimensions and the salt concentration on the performance of the junction are quantified using equivalent lengths.
Nomenclature
A
measured quantity (either U or P )
c+
Positive ion concentration
c−
Negative ion concentration
c0
Mean molar concentration
D
ion diffusivity (D = D+ = D− )
∆Ajunction
change in quantity A due to junction
dA ds
gradient of quantity A
dA ds FD
fully-developed gradient of quantity A
3
e
elementary charge
FD
fully-developed flow
K
dimensionless inverse Debye length
k
Boltzmann constant
Lin
length of inlet channel
Lout
length of outlet channel
Leq.,P
equivalent length based on pressure drop
Leq.,U
equivalent length based on potential increase
n+
positive ion number density
n−
negative ion number density
n0
geometric mean ion number density
nwall
surface normal to wall, pointing into solid region
p
pressure
Pe
Peclet number
s
coordinate pointing in flow direction (parallel to wall)
Re
Reynolds number
T
temperature
4
Uf
fluid potential
Uw
solid potential
v
velocity
V¯
primary channel mean flow velocity
W
channel half-width
x
position
z
ion valence
ε0
permittivity of free space
εf
fluid permittivity
εw
solid permittivity
µ
dynamic viscosity
ρ
density
ρE
charge density
σ
wall surface charge
θ
junction angle
Ωf
fluid domain
Ωw
solid domain
5
a) Separating flow
b) Mixing flow
θ Ωw
Glass
Secondary channel
Ωf
Ls
2W
Glass
Primary channel
3
Lp
W
c)
W
Fig. 1. a) - b) Schematic of flow scenarios considered. c) Schematic of computational domain. The dashed-dotted line denotes the line of symmetry used in the simulations.
2
Model Formulation
The geometry considered in this study is a two-dimensional nanofluidic junction. Two types of flow are considered: a separating flow, whereby the flow 6
through one channel is split into two equal flows, and a mixing flow, where equal flow through two channels is combined into one flow (Figure 1a),b)). Both flows are assumed to be steady-state. The channel contains a symmetric 1:1 electrolytic fluid which is driven through the junction by a pressure gradient. The liquid is assumed to be an aqueous solution of potassium chloride with constant density, viscosity and permittivity; and anions and cations of equal and opposite valence and equal diffusivity.
In general, the surface charge density of silica varies in response to local ion concentrations, in particular the proton concentration (pH) which is itself influenced by other chemical species present in the liquid (e.g. buffers) [13, 14]. As out focus is on the electroviscous performance of a simple binary electrolyte, we neglect these more complex chemical effects and assume that the surface charge density along the channel walls is constant and uniform, and approximately equal to the surface charge density of glass at a local pH= 7 (-10 mC/m2 ) [13, 15]. Despite the simplification, previous studies have shown that such a constant wall charge condition can accurately model flow through glass nanochannels, albeit under electro-osmotic flow conditions [16, 17].
The physical properties used in the simulation are representative of typical microfluidic flows past borosilicate glass (Table 1) [8]. Also included in the model is the solid glass section of the junction material, because electricfield propagation through the walls of nano/microfluidic components has been shown to have a significant effect on the flow dynamics of these systems [6]. 7
Table 1 Physical properties used in the model. (Based on flow of an aqueous solution of potassium chloride past borosilicate glass.) Notes: a [18], Property
[19],c [15],d [20] Value
Temperature (K)
298.0
Viscosity (10−4 Pa s)
8.91a
Density (kg m−3 )
1000.0a
Relative fluid permittivity (-)
80.1a
Relative solid permittivity (-)
4.6a
Channel half-width (nm)
20.0 – 70.0
Concentration (M)
10−6 – 0.1
Ion diffusivity (10−9 m2 s−1 ) Wall charge density (mC m−2 ) Primary channel flow velocity (mms−1 )
2.1
b
2.0b −10.0c 5.0d
Governing Equations
For nanochannels of size greater than 5 nm, the ion transport and fluid flow can be analysed using continuum dynamics [21]. The continuum equations governing the steady flow in the fluid region Ωf are given by the continuity equation:
∇ · v = 0,
(1) 8
the momentum equation:
h
i
∇ · (ρvv) = −∇p + ∇ · µ ∇v + (∇v)T − ρE ∇Uf ,
(2)
the ion-transport equation for the anions and cations present in the fluid:
z± eD± n± ∇ · (n± v) = ∇ · D± ∇n± + ∇Uf , kT
(3)
and the Poisson equation relating the electric potential to the local charge density:
∇ · εf ∇Uf = −
X ezi ni ρE =− . ε0 i=+,− ε0
(4)
In the solid region, the electric potential satisfies ∇2 Uw = 0.
2.2
(5)
Boundary Conditions
The computational domain is shown in Figure 1c). Across the fluid/solid interface between regions Ωf and Ωw the electric field normal to the wall undergoes a jump, given by
εf ε0 ∇Uf
wall
· nwall − εw ε0 ∇Uw
wall
· nwall = σ,
(6)
where nwall is the unit normal to the wall pointing into the solid region. Across the fluid/solid interface the electric potential is continuous, meaning
Uf
wall
= Uw 9
wall
.
(7)
Along the external walls of the device, the gradient of electric potential is constrained in the direction normal to the wall by
∇Uw
wall
· nwall = 0,
(8)
representing a boundary with zero surface charge and no field propagation. The fluid/solid interface was considered to be clean, smooth and hydrophilic (appropriate for water on clean glass), and as such was modelled as a no-slip boundary [22, 23, 24]. A zero ion flux normal to the interface was also imposed. Symmetry conditions were imposed on all quantities along the centreline of the junction. At the inlet and outlet of the channel, the velocity gradient normal to the flow and the tangential velocity component were set to zero. At the √ inlet, a flow-rate and the geometric mean ion number density n0 = n+ n− (calculated from the inlet molar concentration c0 ) were specified, along with a zero normal gradient condition on the counter-ion concentration. At the outlet the normal gradient of both ion concentrations was specified as zero. The total current through the channel was set to zero, the necessary condition for electroviscous flow.
2.3
Numerical Method
The governing equations Equations 1-5 were solved in conjunction with the requisite boundary conditions on a composite structured/unstructured mesh generated with the program gmsh [25]. To ensure fully-developed flow at the inlet and outlet of the junction, Lin = Lout = 22.5W . The fluid region Ωf was meshed using rectangular elements, with 50 elements per channel half-width W . When the junction angle was other than 90◦ , a small region of triangular 10
a)
b)
Fig. 2. Coarse representation of the mesh for a) a 90◦ junction and b) a 45◦ junction. The dashed white lines delineate the boundary between the fluid region Ωf , and the solid region Ωw . Note that only a portion of the mesh used in the simulations is shown.
elements was used (Figure 2b)). A typical mesh used in the simulation had ∼ 3.1 × 105 elements. The distribution of elements across the channel is nonuniform (Figure 2), with fewer elements near the centreline and more elements near the wall, in order to capture the high gradients of electric potential and salt concentrations in the electric double layer adjacent to the wall. The open-source multi-physics solver arb[26] was used to solve the system of equations and boundary conditions. In the context of electrokinetic and nano/microfluidic flows, arb has been used to calculate the electrokinetic flow of a binary electrolyte through parallel nanochannels[27] as well as the interfacial mass transfer and reaction in Y-Y microfluidic channels[28, 29, 30]. arb discretises transport equations using an implicit finite-volume method, and solves the resulting numerical system via a backstepped, Newton Raphson method. Interpolation and differentiation of cell-centred values to element faces is accomplished via a variation of the Moving Least Squares method, as 11
detailed in [31]. For this study 2nd order polynomials were used as the basis for this method. Velocity interpolation to element faces references the local pressure and potential fields evaluated at the face, in the spirit of Rhie-Chow interpolation [32]. Ion concentrations are advected using a high-order, limited upwind method. Input files which fully detail the discretisation technique for the flow of a binary electrolyte through a charged channel, and upon which the present simulations are based, are available for download 2 . The matrix at each Newton-Raphson iteration is inverted using the sparse linear solver PARDISO that is packaged within the Intel Math Kernel Library. The non-dimensional convergence criterion employed was = 10−10 [26].
Testing of the method has shown that the pressure and potential gradients for fully-developed electrokinetic channel flow on a structured mesh are second order accurate (that is, the error decreases proportionally to the average element dimension squared). For an unstructured mesh, used only within the junction region of the non 90◦ junctions, the order of accuracy order is more difficult to quantify, but significantly, the pressure and potential gradient errors for fully-developed channel flow do decrease consistently with mesh refinement. For the mesh resolution used in this study, mesh sensitivity testing indicates that the reported potential and pressure equivalent lengths have errors of less than 1% over the entire parameter range studied.
2
http://people.eng.unimelb.edu.au/daltonh/downloads/arb/code/
latest/examples/electrokinetic channel flow/
12
3
Results
Figure 3a) and b) depict the co-ion distribution for both the separating and mixing flows through a junction of channel half-width W = 50 nm and geometric ion ion number density n0 corresponding to molar concentration c0 = 10−4 M. Electric double layers are present, with a depletion of co-ion concentration near the junction walls. For the physical properties used in the simulations the Peclet number is low (Pe = V¯ W/D = 0.125). Hence diffusion and conduction are the dominant means of ion transport, meaning the size of the double layer is largely unaffected by the presence of the junction and the direction of the flow, and that symmetry of the double layer is present about the leading and trailing edge of the junction corner.
Figure 3c) and d) show the contours of velocity magnitude and streamlines. Because of the apparently symmetric nature of the ion distribution about the junction and the fact that the Reynolds number is extremely low (Re = ρV¯ W/µ = 2.8 × 10−3 ), the velocity distribution at the junction could be expected to be insensitive to the flow direction. However, there are subtle differences between the mixing and the separating flow, clearly highlighted in the contours of geometric mean ion concentration (Figure 3e) and f)). The low but finite Pe number effect means that the distribution of ions is biased very slightly towards the outlet of the junction. Thus, due to the non-linear electric body-force term present in Equation 2, the flow through the teepiece is irreversible. The extent of irreversibility is quantified below.
Figure 3g) and h) compare the potential contours for separating and mixing flow. The junction has an effect on the potential distribution in the region of 13
c− (mM) 0.000 0.008 0.016 0.024 0.032 0.040
24
24
22
22
-4
-2
0
2
4
a)
-4
-2
0
2
4
-2
0
2
4
-2
0
2
4
-2
0
2
4
b) |u| (mm/s) 0
1
2
3
4
5
6
7
8
24
24
22
22
-4
-2
0
2
4
c)
-4
d) 0.096 0.097 0.098 0.099 0.100 0.101
c0 (mM)
24
24
22
22
-4
-2
0
2
4
e)
-4
f) U (mV) -130 -114 -97 -81 -65 -49 -32 -16
0
24
24
22
22
-4
g)
-2
0
2
4
-4
h)
Fig. 3. Distribution of co-ions for a) separating flow, and b) mixing flow. Contours of velocity magnitude and streamlines in the junction for c) separating flow, and d) mixing flow. Distribution of geometric mean ion number density for e) separating flow, and f) mixing flow. Contours of electric potential for g) separating flow, and h) mixing flow. The geometric ion number density n0 corresponds to molar concentration c0 = 10−4 M, the channel half-width W = 50 nm, and the axes have been normalised with length-scale W . All physical properties are listed in Table 1.
the geometry change for both mixing and separating flows. The gradient of potential along the channel is opposite for the mixing flow in comparison to the separating flow, due to the change in flow direction. The changes in pressure and potential caused by the presence of the junction can be quantified using [8] ∆Ajunction
dA = ∆A − Li , ds FD,i i=in,out X
(9)
where A is either the pressure or the potential. This definition is valid only if ∆A is measured over a length beginning and ending in fully-developed flow regions in the inlet and outlet respectively. The pressure drop and potential increase caused by the presence of a junction with separating flow are given in Figure 4a) for a range of channel half-widths and salt concentrations. Over the range of concentrations and channel half-widths modelled, the pressure drop due to the presence of the junction varies by almost an order of magnitude. The magnitude of the pressure drop increases with decreasing channel half-width W . Figure 4 compares the pressure drop of a junction with W = 50 nm for both mixing and separating flows. Due to the distribution dependence of the ions on the flow direction the flow is irreversible, and the pressure drop caused by the junction also depends on the flow direction. However, the flow is only weakly irreversible, and the pressure drop for mixing flow is higher than the separating flow value by only ∼ 2 Pa (< 1%) over the range of concentrations studied. There is also significant variation in the potential increase due to the junction (Figure 5). For low concentrations, the increases in potential due to the presence of the junction are constant, but dependent on the size of the channel half-width W . For very high concentrations the increase in potential due 15
-100 -200
∆P (Pa)
-300 -400
W = 70 nm W = 60 nm W = 50 nm W = 40 nm
-500 W = 30 nm -600 -700 -800 -6 10
W = 20 nm -5
10
a)
-4
-3
10 10 Concentration (M)
-2
10
-1
10
-280
∆P (Pa)
-285
-290
-295
-300
-305 -6 10
b)
Separating flow Mixing flow No electrokinetics -5
10
-4
-3
10 10 Concentration (M)
-2
10
-1
10
Fig. 4. The pressure drop due to the junction as a function of average salt concentration for: a) separating flow over a range of channel half-widths, and b) separating and mixing flow for channel half-width W = 50 nm. All physical properties are listed in Table 1. Each symbol represents a simulation result.
W= W= W= W= W= W=
0.6 0.5
20 nm 30 nm 40 nm 50 nm 60 nm 70 nm
∆U (mV)
0.4 0.3 0.2 0.1 0 -6 10
-5
10
a)
-4
-3
10 10 Concentration (M)
-2
10
-1
10
Separating flow Mixing flow
0.6
∆U (mV)
0.5 0.4 0.3 0.2 0.1 0 -6 10
b)
-5
10
-4
-3
10 10 Concentration (M)
-2
10
-1
10
Fig. 5. a) The potential increase due to the junction as a function of average salt concentration for: a) separating flow over a range of channel half-widths, and b) separating and mixing flow for channel half-width W = 50 nm. All physical properties are listed in Table 1. Each symbol represents a simulation result.
to the junction drops to zero. The potential increase due to the junction for mixing flow is higher than that for separating flow (Figure 5b)). Two useful quantities can be introduced to explain the behaviour of the pressure drop, and potential increase, caused by the presence of the junction. The first quantity is the equivalent length of the junction based upon the pressure drop or the potential increase. The equivalent length is defined as the length of uniform inlet channel (with half-width W) needed to match the pressure drop (potential increase) caused by the effect of the junction [7, 8]. Formally, these equivalent lengths are defined as Leq.,A =
∆Ajunction , dA ds FD,primary
(10)
where A is either the pressure or the potential. The second quantity is the dimensionless inverse Debye length, defined as s
K=
2z 2 e2 n0 W 2 . εf ε0 kT
(11)
This parameter can also be interpreted as the square root of a dimensionless concentration. When K is small, the thickness of the electric double layer is greater than the channel dimension, and significant double-layer overlap occurs. At this limit, there is a nearly constant amount of charge (opposite and equal to the fixed surface charge) across the channel due to the exclusion of co-ions, and no region of electrically neutral fluid. When K is large, the double layer thickness is small in comparison with the channel dimension, and all the charge present in the flow is located in a thin region adjacent to the wall. As a result, the bulk flow is electrically neutral. When the pressure drops caused by the junction given in Figures 4 and 5 are 18
reformulated in terms of equivalent lengths and dimensionless inverse Debye lengths K, there is significant collapse of the data (Figures 6 and 7). Over the parameter space simulated, the equivalent length based on the pressure drop varies from ∼ 0.85W to ∼ 1.05W , while the equivalent length based on the potential increase varies from ∼ 0.3W to ∼ 0.6W . For values of K < 1 both equivalent lengths are weakly dependent upon the channel half-width W ; but independent of K due to the exclusion of co-ions in the channel, and the consequent fixed charge therein. As the double layer thickness becomes comparable to the channel dimension however, the resistance of the junction begins to vary with K. As K increases, so does Leq.,P . For values of K > 10 the fluid bulk is electrically neutral, and as a consequence the equivalent length based on the pressure is the same as the equivalent length calculated for purely hydrodynamic flow (found by an additional simulation). Interestingly, for values of K > 3 and W > 50 nm, the equivalent length based on the pressure is independent of W . The equivalent length of the junction based on the potential displays similar behaviour for K < 1. For increasing values of K above 1, Leq.,U decreases until a minimum resistance is reached at K ∼ 5. As K increases further, the double layer decreases in thickness, and electrokinetic effects decrease. Consequently, the potential increase caused by the junction (numerator of Equation 10), and the potential increase per unit length of uniform channel (denominator of Equation 10) go to zero as K → ∞. However, upon inspection of the data, the potential increase caused by the junction goes to zero more slowly than the potential increase per unit length of uniform channel, resulting in Leq.,U increasing as the double layer becomes thinner. The effect of flow direction on the resistance of the junction to the flow is 19
1.1
Leq.,P
1
W = 20 nm W = 30 nm W = 40 nm W = 50 nm W = 60 nm W = 70 nm No electrokinetics
0.9
0.8 0.1
1
a)
10 K
1.1
Separating flow Mixing flow No electrokinetics Leq.,P/W
1
0.9
0.8 0.1
b)
1
10 K
Fig. 6. Equivalent length of the junction based on the pressure drop over the junction as a function of dimensionless inverse Debye length K for a) separating flow for a range of channel half-widths, and b) separating flow and mixing flow for channel half-width W = 50 nm.
0.7
0.6
Leq.,U/W
0.5
0.4
0.3
0.2
W= W= W= W= W= W=
20 nm 30 nm 40 nm 50 nm 60 nm 70 nm
0.1
1
a)
10 K
0.7
Separating flow Mixing flow
Leq.,U/W
0.6
0.5
0.4
0.3 0.1
b)
1
10 K
Fig. 7. Equivalent length of the junction based on the potential increase over the junction as a function of dimensionless inverse Debye length K for a) separating flow for a range of channel half-widths, and b) separating flow and mixing flow for channel half-width W = 50 nm.
1
Separating flow Mixing flow
Leq.,P/W
0.75
0.5
0.25
0 0
10
20
30
40 50 60 Junction angle θ
70
80
90
Fig. 8. Equivalent length of the junction based on the pressure drop over the junction as a function of junction angle θ for separating flow and mixing flow. The channel half-width W = 50 nm, and the geometric ion number density n0 corresponds to molar concentration c0 = 3.8 × 10−5 M (K = 1).
minor, with the mixing flow exhibiting slightly higher equivalent lengths due to the pressure drop (Figure 6b)), and the potential increase (Figure 7b)), in comparison to the separating flow. The effect of junction angle θ is shown in Figure 8 for mean ion number density n0 corresponding to molar concentration c0 = 3 .8 × 10−5 M (K = 1). The resistance of the junction decreases quasilinearly with decreasing junction angle. The effect of flow direction is again very small, with the resistance of separating flow and mixing flow differing by less than ∼ 4%. 22
4
Conclusions
The steady-state pressure-driven flow of an aqueous solution of potassium chloride through a junction has been studied for a range of channel dimensions and salt concentrations. Two types of flow through the junction are compared: a separating flow, commonly used in nano/microfluidic devices to sort drops and cells, and a mixing flow, used to mix multiple species and reagents. The resistance of the junction to the each flow has been quantified in terms of equivalent lengths based upon the pressure drops, and potential increases, caused by the junction. Over the parameter space simulated, the equivalent length based on the pressure drop varies from ∼ 0.85W to ∼ 1.05W , while the equivalent length based on the potential increase varies from ∼ 0.3W to ∼ 0.6W . For the parameters examined, the flow through the teepiece is weakly irreversible because of a slight readjustment of the double-layer downstream to the junction. As a consequence the effect of flow direction (mixing or splitting) on both equivalent lengths is minor, with differences of less than 1% observed. Decreasing the junction angle acts to decrease the resistance of the junction to the flow quasi-linearly.
Thus, for lab-on-a-chip design purposes, the excess pressure drop of a 90◦ junction can be estimated by multiplying the fully-developed one-dimensional electroviscous pressure gradient by 0.95 ± 0.1, and the excess potential increase can be approximated by multiplying the fully-developed one-dimensional electroviscous potential gradient by 0.45 ± 0.15 . 23
5
Acknowledgments
The authors thank Mr Christian Biscombe for enlightening discussions. This research was supported by the Australian Research Council Grants Scheme and by computational resources on the National Computational Infrastructure Facility through the National Computational Merit Allocation Scheme.
References [1] S. B. Cheng, C. D. Skinner, J. Taylor, S. Attiya, W. E. Lee, G. Picelli, D. J. Harrison, Development of a multichannel microfluidic analysis system employing affinity capillary electrophoresis for immunoassay, Analytical Chemistry 73 (7) (2001) 1472–1479. [2] A. Hatch, A. E. Kamholz, K. R. Hawkins, M. S. Munson, E. A. Schilling, B. H. Weigl, P. Yager, A rapid diffusion immunoassay in a T-sensor, Nature Biotechnology 19 (5) (2001) 461–465. [3] T. M. Squires, S. R. Quake, Microfluidics: fluid physics at the nanoliter scale, Reviews of Modern Physics 77 (3) (2005) 977–1026. [4] R. J. Hunter, Zeta potential in colloid science, Academic Press, London, 1981. [5] M. Davidson, D. Harvie, Electroviscous effects in low Reynolds number liquid flow through a slit-like microfluidic contraction, Chemical Engineering Science 62 (16) (2007) 4229–4240. [6] J. D. Berry, M. R. Davidson, D. J. E. Harvie, Electrokinetic development length of electroviscous flow through a contraction, ANZIAM Journal 52 (2011) C837–C852. [7] J. D. Berry, M. R. Davidson, R. P. Bharti, D. J. E. Harvie, Effect of wall 24
permittivity on electroviscous flow through a contraction, Biomicrofluidics 5 (2011) 044102. [8] J. D. Berry, A. E. Foong, C. Lade, C. J. C. Biscombe, M. R. Davidson, D. J. E. Harvie, Electroviscous resistance of nanofluidic bends, submitted to Chemical Engineering Science. [9] S. Wong, M. Ward, C. Wharton, Micro T-mixer as a rapid mixing micromixer, Sensors and Actuators B: Chemical 100 (3) (2004) 359–379. [10] K. Ahn, C. Kerbage, T. P. Hunt, R. M. Westervelt, D. R. Link, D. A. Weitz, Dielectrophoretic manipulation of drops for high-speed microfluidic sorting devices, Applied Physics Letters 88 (2006) 024104. [11] B. Ahn, K. Lee, R. Panchapakesan, K. W. Oh, On-demand electrostatic droplet charging and sorting, Biomicrofluidics 5 (2) (2011) 024113. [12] A. A. S. Bhagat, E. T. K. Peterson, I. Papautsky, A passive planar micromixer with obstructions for mixing at low Reynolds numbers, Journal of Micromechanics and Microengineering 17 (5) (2007) 1017–1024. [13] L. H. Allen, E. Matijevic, Stability of colloidal Silica. II. Ion exchange, Journal of Colloid and Interface Science 33 (3) (1970) 420–429. [14] P. J. Scales, F. Grieser, T. W. Healy, L. R. White, D. Y. C. Chan, Electrokinetics of the silica-solution interface: a flat plate streaming potential study, Langmuir 8 (3) (1992) 965–974. [15] R. B. Schoch, H. van Lintel, P. Renaud, Effect of the surface charge on ion transport through nanoslits, Physics of Fluids 17 (10). [16] D. Stein, M. Kruithof, C. Dekker, Surface-charge-governed ion transport in nanofluidic channels, Physical Review Letters 93 (3) (2004) 1–4. [17] D. J. E. Harvie, C. J. C. Biscombe, M. R. Davidson, Microfluidic circuit analysis I: Ion current relationships for thin slits and pipes, Journal of Colloid and Interface Science 365 (2012) 1–15. 25
[18] H. Ellis (Ed.), Book of Data, 4th Edition, Longman Group Limited, Harlow, 1986. [19] W. M. Haynes (Ed.), CRC Handbook of Chemistry and Physics (Internet version), 91st Edition, CRC Press/Taylor and Francis, Florida, 2011. [20] O. J¨annig, N.-T. Nguyen, A polymeric high-throughput pressure-driven micromixer using a nanoporous membrane, Microfluidics and Nanofluidics 10 (3) (2010) 513–519. [21] H. Daiguji, Ion transport in nanofluidic channels, Chemical Society Reviews 39 (2010) 901–911. [22] C. Honig, W. Ducker, Thin film lubrication for large colloidal particles: experimental test of the no-slip boundary condition, Journal of Physical Chemistry C 111 (44) (2007) 16300–16312. [23] C. Honig, W. Ducker, No-slip hydrodynamic boundary condition for hydrophilic particles, Physical Review Letters 98 (2) (2007) 028305. [24] O. I. Vinogradova, A. V. Belyaev, Wetting, roughness and flow boundary conditions., Journal of Physics. Condensed Matter 23 (18) (2011) 184104. [25] C. Geuzaine, J.-F. Remacle, Gmsh: A 3-d finite element mesh generator with built-in pre- and post-processing facilities, International Journal for Numerical Methods in Engineering 79 (11) (2009) 1309–1331. [26] D. J. E. Harvie, An implicit finite volume method for arbitrary transport equations, ANZIAM Journal 52 (2012) C1126–C1145. [27] C. Biscombe, M. Davidson, D. Harvie, Electrokinetic flow in parallel channels: circuit modelling for microfluidics and membranes, Colloids and Surfaces A: Physicochemical and Engineering Aspects, 440 (2014) 63–73. [28] L. R. Mason, D. Ciceri, D. J. E. Harvie, J. M. Perera, G. W. Stevens, Modelling of interfacial mass transfer in microfluidic solvent extraction. Part I. Heterogeneous transport, Microfluidics and Nanofluidics, 14 (2013) 26
197–212. [29] D. Ciceri, L. R. Mason, D. J. E. Harvie, J. M. Perera, G. W. Stevens, Modelling of interfacial mass transfer in microfluidic solvent extraction. Part II. Heterogeneous transport with reaction, Microfluidics and Nanofluidics, 14 (1-2) (2013) 213–224. [30] D. Ciceri, L. R. Mason, D. J. E. Harvie, J. M. Perera, G. W. Stevens, Extraction kinetics of Fe (III) by di-(2-ethylhexyl) phosphoric acid using a Y–Y shaped microfluidic device, Chemical Engineering Research and Design, In press. [31] L. Cueto-Felgueroso, I. Colominas, X. Nogueira, F. Navarrina, M. Casteleiro, Finite volume solvers and moving least-squares approximations for the compressible Navier–Stokes equations on unstructured grids, Computer Methods in Applied Mechanics and Engineering 196 (4548) (2007) 4712–4736. [32] C. Rhie, W. Chow, Numerical study of the turbulent flow past an airfoil with trailing edge separation, AIAA Journal 21 (11) (1983) 1525–1532.
27