PROOF COPY 022212PHF
PHYSICS OF FLUIDS VOLUME 14, NUMBER 12 DECEMBER 2002 Nonlinear electrokinetic ejection and entrainment due to polarization at nearly insulated wedges
Sunil Kumar Thamida and Hsueh-Chia Chang
Department of Chemical Engineering, University of Notre Dame, Notre Dame, Indiana 46556 Received 4 June 2002; accepted 16 September 2002 We examine a singular electrokinetic flow around a corner or a wedge in micro-channels constructed from dielectric materials whose permittivity is small but finite compared to that of the electrolyte. When the wedge angle is less than 180�, the applied electric field, which is tangential far from the corner, develops a normal surface component that becomes singular at the corner. This normal field leakage causes opposite polarization at the two sides of the wedge and produces a converging singular tangential electrokinetic flow that ejects liquid from the tip. By expanding in cylindrical harmonics, we estimate this ejecting flow as a function of the permittivity ratio, applied electric field, angle of the wedge and the microscopic corner curvature that suppresses the singularity. The ejecting flow entrains tangential flow on the front side of the wedge and produces a vortex on the downstream side. This entrainment offers a long-range attractive hydrodynamic force that complements short-range electrostatic DLVO Derjaguin�Landau�Verwey�Overbeek and dielectrophoretic forces to enhance corner deposition and aggregation of colloids and proteins during electrophoresis/electro-osmosis. � 2002 American Institute of Physics. DOI: 10.1063/1.1519530 O PR OF CO I. INTRODUCTION Electrokinetic flow is the mechanism of choice for transporting fluids in future generations of biochips, as it allows easy flow control, metering, and maneuvering.1 However, one major problem for electrokinetic flow is that proteins tend to precipitate at channel junctions and colloids aggregate at the junction corners. As an example, a submicron latex colloidal suspension is driven electro-osmotically and electrophoretically in water through the micro-channel junction in Fig. 1 a . After about 20 min, colloids aggregate in large amounts around the inner corner of the channel junction as shown in Fig. 1 b . As shown in the figure, the aggregates form a curious and relatively large spiral-like structure that curls downstream and spans nearly half of the 80 m channel. A small vortex is visible just downstream of the spiral, as shown in Fig. 1 b , when the colloid trajectories are scrutinized with a microscope. The aggregation becomes acute with higher electric fields 30 V/cm and smaller channels. Several physical mechanisms can be responsible for the corner aggregation phenomenon. The spiral aggregate exceeds 40 m in size and is much larger than homogeneous aggregates formed in comparable time by the DLVO Derjaguin�Landau�Verwey�Overbeek mechanism.2 One hence expects a preferential migration of the particles towards the corner that is not present in homogeneous aggregation. As we shall demonstrate in subsequent analysis, both the tangential and normal electric fields are very large at the corner. The latex particles can be polarized by this intense field and migrate by dc dielectrophoresis.3 This dielectrophoretic migration is driven by electric fields that are highly nonuniform spatially, such as those at the junction corner.
1070-6631/2002/14(12)/1/14/$19.00 The nonuniform field will impart a net force on the bipolar particle despite its opposite charges on the two sides perpendicular to the field line. However, the direction of this dielectrophoretic motion is determined by the difference in the medium and particle dielectric permittivities. For our systems, the dielectric constant of the water medium is 81 and that of the latex particle is 2.5. Hence, the water medium is more polar than the particle and the dielectrophoretic migration should be towards regions with low field intensity. This is clearly opposite of our observation in Fig. 1, where particles migrate towards the high-field corner region. Moreover, dielectrophoretic motion cannot explain the observed vortex. The field nonuniformity will also be shown to exist in only a small neighborhood of the corner whose dimension is much smaller than the aggregate size. Particle polarization occurs and is a key mechanism that holds the aggregate together, but it probably does not contribute to its formation by particle migration towards the corner in this transport-limited aggregation process. Instead, the large aggregate dimension, the observed vortex, and the aggregate spiral structure that seems to be shaped by the vortex, all suggest that the particle migration towards the corner is driven by hydrodynamic convection. Such a theory is proposed in this article. The small region of high field intensity near the corner is shown to affect the flow in a much larger region by a nonlinear electrokinetic mechanism. We shall demonstrate through theoretical analysis and numerical simulation that this long-range hydrodynamic entrainment mechanism can drive particles to the corner, produce the observed vortex, and allow the formation of a spiral aggregate. In an electrokinetic flow, fluid is driven by an electric field. Usually, an electric potential drop is applied across the PY
1 02 22 P 12 HF � 2002 American Institute of Physics PROOF COPY 022212PHF PROOF COPY 022212PHF
2 Phys. Fluids, Vol. 14, No. 12, December 2002 S. K. Thamida and H.-C. Chang FIG. 1. Experimental snapshot of the micro-channel junction. The narrow channel is 80 m wide and the latex colloids move towards the narrower channels. The silica channel and the latex colloids are oppositely charged and hence both electrophoretic and electro-osmotic motions are in same direction. a During initial stages of electrokinetic flow of submicron colloidal suspension in the direction indicated by the arrow. b During later stages showing a spiral colloidal aggregation at the inner corner and a small vortex as sketched schematically. ends of a channel using electrodes which establishes an electric field parallel to the walls. The micro-channels are made of dielectric materials like silica, polyester, acrylic, etc. Due to surface functional groups, this material creates a charged Debye layer in the electrolyte neighboring it, such that the counter-ions within screen the surface charges. The extent of charge separation in the electrolyte is usually measured by the electric potential across it, the zeta potential , which would be proportional to the surface charge of the wall. For typical dielectric-electrolyte pairs,2 the potential is in the range of 25�100 mV. Due to the Maxwell body force within the Debye layer, an electric field parallel to the wall drives the bulk fluid at a velocity given by the Smoluchowski slip velocity4 O PR OF
f us Et , 1 where f is the liquid permittivity, is its viscosity, and E t is the tangential field. For dielectric materials with small permittivity compared to the electrolyte, the normal electric field E n in the wall vanishes and the charge distribution within the Debye layer is independent of the applied field. The polarization in the Debye layer is produced by the field generated by the surface charge. As such, both the zeta potential and the Debye layer thickness are uniform in the entire channel of the same material. Since there is no field leakage and no liquid leakage, this linear electrokinetic slip 1 implies that the velocity field everywhere in the channel is simply u
f E, 2 is the electric field. Hence, the electric powhere E tential also becomes the velocity potential and, despite the minuscule Reynolds number in micro-devices ( 10 1 ), linear electrokinetic flow behaves like an inviscid potential flow
PROOF COPY 022212PHF due to the absence of wall shear. The electric field lines coincide with the streamlines in this case for an open channel.4,5 More interestingly, this potential flow is irrotational and the generation of microscopic vortices in electrokinetic devices for mixing purpose is a major challenge.6 Yet, Fig. 1 b indicates that a vortex exists just downstream of the junction. We shall demonstrate that this vortex is generated by a nonuniform zeta potential which is dependent on the applied field. The resulting nonlinear electrokinetic flow produces an ejecting flow from the corner and a back pressure that generates the vortex. This nonlinear electrokinetic phenomenon is a result of field leakage at the corner. If the dielectric wall has a small dielectric permittivity mathematically zero compared to the fluid, then the electric potential around a wedge corner can be solved by decomposing the bulk Ohmic electro-neutral region from the wall Debye charged layer and by the usual harmonic expansion.4 For a wedge angle less than 180�, the tangential field is singular at the corner from classical electrostatic potential theory.7 However, the normal field remains zero everywhere on the boundary and even the tangential field is smooth away from the boundary. Hence, the coincidence between streamlines and electric field lines still holds away from the corner. However, biochips are made with dielectrics with a small but finite permittivity. While field leakage is negligible over most of the chip, it is very important at corners because of convexity and the singular tangential field. Using a perturbation analysis, we show that the singular tangential field now leads to a singular normal field that also blows up at the corner. This gives rise to large but opposite polarization on two sides of the corner see Fig. 2 in addition to the uniform polarization due to surface charges. These normal fields drive two oppositely charged ions to each side. Such polarization produces a destabilizing normal Maxwell stress across a deformable membrane.8 Here, it produces another instability with a converging electro-osmotic flow towards the corner. CO PY 02 22 P 12 HF PROOF COPY 022212PHF
Phys. Fluids, Vol. 14, No. 12, December 2002 Nonlinear electrokinetic ejection and entrainment 3 FIG. 3. Schematic diagram of a 90� bend in a micro-channel with electrodes located at far ends. FIG. 2. Schematic diagram showing polarization and the boundary conditions on the electric field of a wall-electrolyte interface around a wedge. A wedge with a finite radius of curvature is also shown. This instability occurs because the zeta potential near a corner would have opposite contribution from the normal leakage field across the corner, which brings excess ions that overwhelm the original surface charge of the wall. This normal field is singular infinite at the corner, localized, and decays to an insignificant amount over a few microns away from the corner. The resulting nonlinear Smoluchowski slip velocity, which is proportional to the product of zeta potential and tangential field, also becomes nonlinear with respect to the applied field, singular at the corner, and changes sign across the corner over a microscopic length. Thus an electrostatic field leakage, though distributed, causes a momentum dipole that produces an ejecting electrokinetic flow from the corner. Although the field penetration length is typically small 1 m , the ejection distance hydrodynamic penetration or entrainment length is long at tens and hundreds of microns due to the singularity. We estimate the finite strength of this ejecting flow by assigning a finite but microscopic radius of curvature to the corner tip. The ejecting velocity would be directed symmetrically away from the inner corner towards the outer corner and its strength would be a function of the ratio of dielectric permittivities of wall and the fluid, angle of the wedge, the far-away tangential field, zeta potential of the wall, and the radius of curvature of the corner. For atomically sharp corners, the ejecting flow can exceed the bulk electro-kinetic flow. Due to this strong ejecting flow from the corner, the bulk flow is blocked and vortices develop at the corner. The irrotational feature of linear electrokinetic flow is revoked due to the nonuniform and field-dependent potential at the corner. More precisely, if is nonuniform in 2 , the incompressibility condition ""u 0 no longer produces the Laplace equation for potential flow. We demonstrate vortex formation by obtaining a local solution for the stream function from the biharmonic equation. We also develop a lattice-Boltzmannmethod LBM code suitable to simulate the outer microfluidic flow in an actual 90� bent channel of finite width. We O PR have implemented a multi-scale approach in which the electrostatics and nonlinear electrokinetic flow near the singular point is dealt analytically while the electrokinetic microhydrodynamics is simulated numerically with a LBM code. Hence, a very high resolution exists within the thin Debye layer for the field-dependent potential and the slip velocity on the channel walls whereas for flow simulation, the lattice node spacing is only of medium resolution--a few microns. The ejecting flow at the corner stipulates a detailed derivation of the corner boundary conditions and lattice spacing. In fact, the long-range corner ejecting flow entrains most of the in-coming flow. With this large hydrodynamic entrainment length of tens to hundreds of microns, colloids and proteins are convected towards the corner and are exposed to the singular normal field. They would then polarize and aggregate due to the classical DLVO mechanism5 when they are within a submicron neighborhood of the growing aggregate. We propose this long-range hydrodynamic mechanism to be the cause for polarization, alignment, and aggregation of colloids. We obtain the hydrodynamic entrainment fraction through theory and favorably compare it with the multi-scale simulation. Our study then allows us to design channels with minimum hydrodynamic entrainment, protein precipitation and colloid aggregation. The remaining sections of this article are planned as follows. In Sec. II, we present the governing equations for electrokinetic flow along with the electrostatics in electrolyte and the wall. We propose a model to obtain the normal field leakage and effective nonlinear Smoluchowski slip velocity near a wedge corner. In Sec. III, we obtain a local solution to the vortices and overall flow in the vicinity of the wedge corner. In Sec. IV, we confirm the local results with numerical simulation LBM for an electrokinetic flow around a 90� bend.
II. ELECTROSTATICS BY PERTURBATION EXPANSION OF CO PY 02 Consider the schematic diagram in Fig. 3 depicting a micro-channel with a 90� bend. The electrodes are located at the far ends of each branch. When a potential drop is applied across the electrodes, an electric field would be established in the fluid and also the neighboring wall. Since the fluid is an 22 P 12 HF PROOF COPY 022212PHF PROOF COPY 022212PHF
4 Phys. Fluids, Vol. 14, No. 12, December 2002 S. K. Thamida and H.-C. Chang electrolyte, there is a possibility of charge separation and net ionic charge near the channel wall. The electrostatic potential in the electrolyte or the fluid would be governed by the Poisson equation,4
2 f , 3 where is the electrostatic potential, is the net ionic charge density, and f is the dielectric permitivity of the fluid. The potential distribution in the wall, w , would also be governed by a similar equation. There is no charge within the solid dielectric material and its potential obeys the Laplace equation, O PR
w 2 w 0. 4 Note that these equations and the ones to follow are twodimensional. One spatial variable is along the length and another along the width of the channel. The micro-channels are assumed to have uniform and constant slip velocities at the top and bottom. With its small depth-to-width ratio, the variations in the thin direction can be neglected.1 The lateral boundary conditions are OF w f n on side wall,
o w n 5a CO on side wall, 5b due to continuity of potential and field across the two media. Here, o is the uniform surface charge density on the wall since we use the same functional material along the entire channel. To solve the above set of equations, one needs to model the ionic concentration distribution in 3 and estimate the cumulative charge density in the fluid, z ic i , 6 PY 02
FIG. 4. Schematic diagram showing the ionic path lines and zero flux lines around a wedge in a actual geometry and b transformed orthogonal coordinates. 22 i where c i is the concentration of ith ionic species and z i its ionic charge. Ideally, every ionic species should be captured by the full Nernst�Planck convection-diffusion equation with electromigration. However, we will construct specific orthogonal coordinates in a given geometry such that one set of lines is along the direction of ionic flux and the other set of perpendicular lines along which ionic flux is zero. One of these two coordinates correspond to the streamlines. There is hence no convection in the other coordinate orthogonal to the streamlines. The only ion flux mechanisms along them are diffusion and electromigration. Moreover, since the wall is a streamline, the latter orthogonal coordinates must terminate at the wall. Since there is no ion flux into the wall, electromigration must then cancel diffusion exactly everywhere along this set of no-flux lines. Thermodynamic equilibrium must exist along these lines and everywhere the lines are defined on a plane. A simple integration of the convectionless Nernst�Planck equation along these lines then yields the equilibrium Boltzmann distribution,4 c i c i,o exp
PROOF COPY 022212PHF z iF RT , 7 where c i,o is the concentration of ith ionic species in the bulk, is the over-potential relative to the potential in the bulk, F is the Faraday constant, R is the gas constant, and T is the temperature of the fluid. As described earlier regarding the orthogonal coordinates, would be the potential along those set of orthogonal lines that emanate from the impermeable wall and along which there is no ionic flux. These lines would be perpendicular to the wall in straight channels. The thickness of the Debye layer is simply the characteristic decays to the zero length along these lines over which reference value in the bulk. For the case of a wedge, the thickness is measured along the no-flux coordinate. This noflux coordinate is multi-valued at the cusp, as shown in Fig. 4 a . The equilibrium Boltzmann distribution 7 hence becomes undefined exactly at the cusp. However, we shall smooth the cusp with an effective radius of curvature in our theory and the singularity exactly at wedge corner is not significant. Since the leading-order Ohmic field and the wall field are determined by the Laplace equation, we shall make extensive use of the conformal map of the cylindrical (r, ) coordinates of the wedge to rectilinear , coordinates P 12 HF PROOF COPY 022212PHF
Phys. Fluids, Vol. 14, No. 12, December 2002 Nonlinear electrokinetic ejection and entrainment 5 along a flat plane as shown in Fig. 4 b . Away from the corner, the latter coordinates will also be shown to be the desired orthogonal streamlines and no-flux lines for the Boltzmann distribution and the Debye layer thickness. Since the normal field in the Debye layer far exceeds the applied tangential field, classical electrokinetic theory4 decomposes the fluid potential into the Ohmic potential , representing the tangential field for the neutral bulk that sees an insulated wall, and an over potential , representing the normal field in the polarized Debye layer. Moreover, the over obeys the Poisson equation with a wall flux potential specified by the surface charge through 5b . In the classical case, the wall is insulated and w is zero exactly in 5b . Since n t , this decomposition allows an expansion of the quadratic nonlinear Maxwell term in the momentum equation such that it becomes linear in the tangential field t . This expansion leads to the slip velocity 1 directly. Moreover, the decomposition allows a onedimensional resolution of and the momentum equation due to the corresponding length-scale separation between the Debye layer and the channel dimensions that accompanies the separation in the field strengths. Finally, the harmonic, Ohmic potential obeys the Laplace equation with insulated walls and is hence easy to resolve with mathematical techniques like conformal maps. We shall also employ the same convenient decomposition of the fluid potential into an Ohmic component and an over potential component . At the wall polarized layer, as in classical theories. This large field n is n n attributed to the surface charge o in 5b . The new wall leakage term w n w is negligible away from the corner but must be comparable to o near it to reverse the polarization as shown in Fig. 2. We hence associate the wall leakage and o terms on the w n w in 5b to the dominant f n fluid side to produce a new boundary condition. The remaining term f n for the wall Ohmic field on the fluid side is much smaller than these three terms. Hence, to leading order, the fluid Ohmic potential obeys a no-flux condition at the wall. We note that the wall field leakage term w n w has opposite signs across the corner while o has always the same sign. Hence, the wall field on the liquid side f n can have opposite signs across the corner--the polarization and slip velocity can change sign across the corner. The sign change is in fact necessary for back flow and vortex formation. Our theory will produce a new slip velocity that reduces to the classical version 1 when the field leakage is absent. Decomposing the fluid potential as stated, The leading-order Ohmic potential in the fluid is hence governed by the following equations in the cylindrical coordinates centered around the wedge corner as in Fig. 4 a ,
2 0, 0 at . 9a 9b n The large fluid electric field at the wall, as governed by 5b , is assigned to that of the over potential . The solution to the Ohmic potential of 9a in the physical coordinates is an eigenfunction expansion in the cylindrical harmonics of the Laplace equation that satisfy the boundary conditions 9b , Ar sin where 2n 1 2 , n 0, 1, 2,... . 10b , 10a O PR , 8 where is the Ohmic potential governing the tangential transport of neutral electrolyte and is the over potential across the polarized layer near the wall due to field leakage at the corner. As such, the linear boundary condition 5b without the surface charge term is decomposed into a no-flux condition for and a gradient condition for that is determined by the surface charge o and the wall field w n w . We solve for and iteratively by also employing the Boltzmann distribution 7 and using boundary conditions 5 alternatively in the iteration.
PROOF COPY 022212PHF If the wedge were a corner in the bent channel, then the field far from the corner should be a constant in either of the channel branches. Only one eigenfunction corresponding to n 0 would not blow up far away and be zero at the line of 0. Hence for every wedge of certain angle symmetry at there is only one possible eigenfunction with eigenvalue /2 . The constant A would then be obtained from matching with the outer solution. Though the actual value of A is case dependent, it is clear that it scales linearly with the applied field E 0 . Let A A * E o W with A * a dimensionless constant, E o the constant electric field in the channel far from the corner, and W the width of the channel. The dimen^ ^ sionless space variable in 10a is r r /W where r is the actual dimensional radial position. Since the far-way field is constant at E o , it is chosen as the appropriate electric field scaling. An appropriate length scale is the width of the channel W since the effect of a wedge in a bent channel is assumed to persist only over a distance of the width of the channel, which is a reasonable one to assume because beyond a distance equal to width from the corner, both the walls are parallel and the field would be nearly a constant. From numerical solution in bent channels, it will be shown that A (1.0 1.2)E o W, for all angles of the wedge 0 /2. We next determine the wall potential w governed by the Laplace equation 4 . Although the normal field strength of the over potential n is much larger than t or n of the Ohmic field, the actual potential drop across the polarized layer is small compared to the wall Ohmic potential drop across the corner due to leakage. For Debye layers, the zeta potential is the potential drop and is of the order of 25�200 mV. While the over potential with field leakage is higher on one side, we do not expect it to exceed 200 mV, which is far smaller than the drop in the wall Ohmic potential across the corner when leakage exists. We hence neglect in 5a and use the Ohmic potential continuity condition OF CO PY 02
w The solution to the wall potential distribution can then be obtained from 4 with boundary condition 11 as 22
at P 12 HF
, 11 PROOF COPY 022212PHF
6 Phys. Fluids, Vol. 14, No. 12, December 2002 S. K. Thamida and H.-C. Chang w Ar sin sin . 12 2 i z i c i,o exp
f z i F /RT . 14a The normal field leakage part of E n in the adjacent electrolyte at the wall is then
w n f w f w The boundary condition 5b on the normal electric field now produces a condition for the over-potential on the liquid side,
n A W cot r 1 at . 13 En o f w n f w at , 14b It can be seen from the schematic in Fig. 2 that the normal field changes sign in the region near the corner. It is because and rethe field leakage enters from the fluid at enters the fluid at . Without surface charge, this field leakage would establish Boltzmann distributions of oppositely charged ions and create oppositely polarized surfaces on the opposite edges of the wedge, rendering the wedge corner bipolar. We have shown in an earlier work that electrokinetic flow past a large bipolar particle can produce back-flow and vortices.6 Here, back flow and vortices develop because the zeta potential changes sign across the corner and the resulting Smoluchowski slip velocity converges towards the corner. The tangential slip velocity hence decreases and changes sign around a corner and continuity stipulates that a back pressure must appear to eject the flow, as sketched in Fig. 2. This back pressure also generates the downstream vortex. However, the coefficient A will be determined through matched asymptotics with a constant field applied far from the wedge. This far-field behavior depends on the channel geometry and will be resolved numerically in a later section with a finite channel width. The effect of surface charge, ignored in the above argument, will also be included. Although the Ohmic field and wall potentials are solved in the original cylindrical coordinates of the wedge, we shall solve for in the transformed coordinates , of Fig. 4 b . These coordinates are defined by the conformal map that maps the wedge into a plane. They are both harmonic func2 0. Insofar tions that satisfy the Laplace equation 2 as the Ohmic fluid potential is a harmonic function for the wedge geometry, the coordinate simply corresponds to a dimensionless version of . The constant contours are isopotential contours. The orthogonal constant contours correspond to the streamlines of linear electrokinetic flow with an insulated wedge. Hence, the orthogonal coordinates are the desired no-flux lines along which the Boltzmann distribution 7 can be derived and the Debye layer thickness can be defined. The Boltzmann distribution hence holds everywhere the coordinates are defined. Near the corner, wall leakage would mean that coordinates are no longer the no-flux lines but are good approximations of them near the wall. At precisely the corner point, however, no-flux lines cannot be defined at all. Hence, the Boltzmann distribution 7 is valid everywhere except at the corner and a small sector beyond it. We can then use 7 to determine the net charge density i z i c i everywhere except the point on the corner. Since the Ohmic potential is a harmonic function, 3 then yields where E n represents the normal electric field near the wall on the fluid side and includes both surface charge and field leakage contributions. Another boundary condition is that 0 RT/F throughout, then we can furfar from the wall. If ther invoke the Debye�Huckel approximation2 on 14a by expanding in (z i F /RT) and transforming the Poisson� Boltzmann equation to a linear partial differential equation, 1 W2
2 r, 2, O PR 15 OF ( i c i,o F 2 z 2 / f RT) 1/2 is the Debye length and where i use has been made of the electric neutrality condition at large , that is i z i c i 0. We now transform this Poisson equation by the conformal map of the wedge to a plane in Fig. 4 b . The variable is the Ohmic potential function (r, ) in its dimensionless form and corresponds to the field lines. We get from 10 A * r cos A * r sin , , 16a 16b CO PY ^ /W with and being the dimen^ /W and where sionless coordinates. Transforming 15 from (r, ) coordinates to , coordinates gives
2 2 02
W
2 2 where 2 , / 2 2 / 2 is the Laplace operator in rectilinear coordinates. Since is same as the nondimensional 2 E 2 /E 2 potential function, the nonlinear coefficient o where E is the local Ohmic electric field. Since this E t at Ohmic field does not leak into the wall, E t the wall. The field E o is the far field used to scale . In the 2 appropriate variables of , , it can be shown that 2 2 2 (1 )/ ( ) . Since the Debye layer is very close to the wall, that is , it would imply that 0 / which upon substituting in 17 gives an orand / dinary differential equation in as d2 d^2 (1 )/ is the tangential field at the wall where E t which is only a function of and is independent of . Note that ^ is the dimensional form of , that is, ^ W where W is the width of the channel. The boundary conditions in 14 on in the transformed variables are 22
, 2, 17 P 12
E t /E o
2 2, HF
at ^ 0 18 ^ En E t /E o 19a PROOF COPY 022212PHF PROOF COPY 022212PHF
Phys. Fluids, Vol. 14, No. 12, December 2002 Nonlinear electrokinetic ejection and entrainment 7 0 at ^ . 19b Integration of 18 with 19a and 19b then yields the fielddependent , the value of at ^ 0, En . 20 Note that unlike for a plane wall, the Debye layer thickness from 18 near a sharp corner is larger by a factor of E t /E o . But, the normal field in 19a is reduced by the same factor E t /E o . Hence these opposing contributions cancel each other and the formula for zeta potential given in 20 remains the same as for a plane wall which is equal to the product of the normal field at the wall and Debye layer thickness. Since is singular at the wedge, 18 becomes invalid at Et t the corner as shown in Fig. 4 a . However, we expect 20 to be valid away from some neighborhood of the corner. If E n is too large or RT/F, then the above approximation would overshoot the actual value that would be obtained by solving 14a and 14b exactly. Then, we would have to use either Gouy-Chapman theory for symmetrical electrolytes or Grahame's analysis for asymmetrical electrolytes.2 Note from 14b that the normal field at the wall E n and the potential in 20 would consist of two parts, one is due to the field of surface charge ( o / f ) and the second is due to field ( w / f ) n w in the wall neighboring the fluid. With singular field leakage, the second contribution becomes dominant in the vicinity of the corner. We can now estimate the zeta potential from 20 by substituting 13 and 14b to obtain the natural and induced potentials, nates the linear slip. The nonlinear slip velocity has opposite signs on either sides of the wedge which implies that it is directed towards the corner on either sides. Moreover, the nonlinear slip occurs only if 1, that is for wedge angles less than or /2 from 10b . Hence there would be field leakage at the inner corner of a bent channel of Fig. 3. However, for the outer corner, with angle 3 /2 greater than , the eigenvalue will be greater than 1 and there would not be any field leakage, polarization, or vortex formation. This explains why protein and colloid aggregates are never found at the outside corner. Moreover, the induced 1 changes sign across the wedge corner.
III. LOCAL HYDRODYNAMICS O PR
r
1 o When we solve for the hydrodynamic flow, the stream function can also be split into two parts, one for linear and one for nonlinear electrokinetics,
o 1. 23 OF The first one would be due to the linear slip velocity u s,o caused by the uniform surface charge with a constant zeta potential, o . The second part is due to nonlinear slip velocity u s,1 caused by the nonuniform zeta potential 1 due to field leakage. The linear electrokinetic flow biharmonic equation with uniform slip is
4 o CO 0,
t o o 24a 0 u s,o at at , . 24b 24c PY un ut
n o o o f w 1 f 1, 21a , A W 21b With linear slip, it corresponds to potential flow. Hence, of all the cylindrical harmonics of the biharmonic equation, only one Laplace harmonic survives--the same one as the /2 , electro-static harmonic 10a with 02 22 cot at f . 21c 0 Cr cos . 25 The linear and nonlinear parts of Smoluchowski slip velocity can then be estimated by substituting 21 in 1 , u s u s,o u s,1 , u s,o
f w f w f oE t 22a A r W cot cot , r
1 , E2 t
2 1 22b The streamlines of this linear electrokinetic flow are hence exactly identical to the electric field lines, as shown in Fig. 5 a . The magnitude of the coefficient C would be proportional to the bulk flow rate u o W as C A * u o W, where A * is the same constant that would appear in the potential function of 10a . The nonlinear electrokinetic flow stream function 1 must be governed by a biharmonic equation as in 24a and have zero normal flux of fluid at the walls of the wedge, P 12 HF u s,1
f 1E t t 1 0 at , 26 A2 W2 2 and yield opposite tangential velocities on the two sides of the wedge, as shown in Fig. 5 b . Hence, the only permissible harmonics of the biharmonic equation are 22c
1 at Dr sin 27 where u s,o is the linear slip velocity in the absence of field leakage effects and u s,1 is the additional nonlinear slip velocity due to corner field leakage. The linearity and nonlinearity refer to the scalings with respect to the local tangential field E t which is linear for u s,o while quadratic for u s,1 . The nonlinear slip exists only near the corner where it actually domiPROOF COPY 022212PHF n /( ), n 1, 2,... . Since there is an outer with wall, the momentum of the ejected flow should decay due to back pressure and it stipulates that 1 0 as r . Retaining the dominant singular harmonic that decays the slowest, we obtain n 1 and / . The value of this exponent is exactly negative of twice of the exponent for PROOF COPY 022212PHF
8 Phys. Fluids, Vol. 14, No. 12, December 2002 S. K. Thamida and H.-C. Chang FIG. 5. Streamlines of electrokinetic flow around a wedge. a Linear electrokinetic flow o for a perfectly insulating wall. Equation 25 with C 1. b Vortex stream function 1 due to field leakage effect. Equation 27 with D 0.0004. c Overall stream function o 1. the electrokinetic stream function in 10a . However, the exponent is incompatible with that of the nonlinear slip 22c for any positive wedge angle , since 1 2 1 2 1 . 28 In fact, the exponent of the nonlinear slip velocity is incompatible with any of the harmonics in 27 for any wedge /2. Unlike its linear electrokinetic flow counterangle part, the nonlinear creeping electrokinetic flow at a wedge is strictly ill-posed. The difficulty may arise from the absence of a normal and the inability to define a slip direction at the tip. However, real wedges are not infinitely sharp but rather
PROOF COPY 022212PHF O PR OF CO PY 02
have a finite radius of curvature, as seen in Fig. 2. This regularized problem is well-posed and the incompatible harmonics of 27 can still be used to approximate its solution. A finite ejecting flow of strength u c now exists for the smoothed cusp as shown in Fig. 2 and it is taken to be the total flow within the Debye layer on either sides of the wedge converging towards the tip. If the tip is rounded off to a radius r c , then u c can be evaluated from u s,1 of 22c at a r c cot from the tip along the wedge wall: distance r 22 P 12 HF
2 * uc
f w f A2 W2 cot r * W 2 1 . 29a PROOF COPY 022212PHF
Phys. Fluids, Vol. 14, No. 12, December 2002 Nonlinear electrokinetic ejection and entrainment 9 If we invoke the matched asymptotic result, A A * E o W, where A * has to be determined later, a more explicit estimate of u c can be obtained as uc
f w f E 2A *2 o 2 cot r * W 2 1 . 29b The ejecting velocity can be compared to the far-field electrokinetic flow velocity u o given by 1 , uc uo
w f Eo
o A *2 2 cot r * W 2 1 29c 2 /4, For a 90� bend 3 , the ejecting velocity becomes comparable to u o for (r /W) 10 4 for a typical values of * 1 50 nm, o 50 mV, and E o 100 V/cm. For ( w / f ) 20 , a mm-dimension channel, this is reached at r c of roughly 100 nm, a reasonable roughness scale for silica channels fabricated by plasma etching. Due to the cot scaling, this critical ejecting velocity becomes even more severe for sharp edges. As shown in Fig. 2, it produces an ejecting flow of strength u c at a radius roughly equal to that of Debye layer thickness . This gives an estimate of D that is the strength of the circulating vortices and we get 1 in 27 as O PR
2 OF * FIG. 6. The entrainment stream function e scales as a power-law with respect to r c and increases with decreasing wedge angle, . Here is the coefficient in 32 , ( w / f )(E o / 0 )( /W)A * 2 ( 2 /2)cot . CO 1 uc 2 r sin 2 . 30 * e w f Eo
0 2 W A *2 2 cot r c cot W 2 1 . 32 However, we do not expect the flow near the wedge to be dominated by only one Stokes harmonic, especially one with an incompatible exponent from the nonlinear slip. The incompatibility gives rise to an error that grows towards the corner. Rather, 30 should only be valid sufficiently far from the wedge. In fact, we shall show from our subsequent simulation that, with a sufficiently sharp wedge, the flow near the wedge is a radial potential flow. Nevertheless, the existence of vortices and hydrodynamic entrainment can be estimated from 1 . The streamlines of linear electrokinetic flow 0 , nonlinear electrokinetic flow 1 , and the combined flow 0 1 are shown in Figs. 5 a �5 c . It can be seen from Fig. 5 c that the streamlines are drawn towards the corner on the front side of the wedge. All the streamlines approaching the corner are bent inwards and make a half circulation before passing on to other side of the wedge. This hydrodynamic entrainment is due to the ejecting momentum at the corner. The amount of flow that is entrained can be estimated from the value of the overall stream function near the wedge. Near the corner, the circulatory flow due to 1 would be the dominating one. It can be estimated as the maximum value of 1 or total ejecting flow from 30 in the region of singularity, that is, r /2. If e is the fraction of the flow and 0 entrained, then uc . 2 31 PY e Hence the entrainment efficiency to the total flow u o W scales as
PROOF COPY 022212PHF * e that is the ratio of e The channel width is W and a pure electrokinetic flow with a flat velocity profile u o described by the Smoluchowski slip 1 has been assumed far from the corner. We will demonstrate from subsequent matched asymptotic result for the electric field in finite-width channels that A * 1 for all wedge angles less than . The entrainment efficiency, * , e and a power-law has complex scalings with respect to scaling with respect to r c . Nevertheless, it increases with decreasing and r c , as seen in Fig. 6. If the wedge were a corner in a bent channel with a known flow rate wall stream function , then it can be predicted that the entire flow could be entrained if e is equal to the wall stream function. This critical condition depends on the angle of the wedge, radius of curvature or sharpness, dielectric constant of the wall, the mean electric field, Debye layer thickness, width of the channel, and the mean zeta potential. On the front side of the wedge, both main flow stream function o and the ejecting flow stream function 1 are positive whereas on the back side of the wedge, o is positive while 1 is negative. The resultant overall flow would have a back-flow or a vortex in the back side wherever the overall stream function is negative. The size of this vortex is insignificant for smaller strengths of the ejecting velocity, but it increases gradually as the ejecting velocity u c increases. A vortex also appears at the upstream side but tends to be much weaker. The current analysis is a local one near the corner. Only dominant harmonics have been retained and the flow near the wedge can be very different from 30 . The outer flow can also distort the vortices in Fig. 5 c . However, the back vortex and the entrainment condition 31 would be 02 22 P 12 HF PROOF COPY 022212PHF
10 Phys. Fluids, Vol. 14, No. 12, December 2002 S. K. Thamida and H.-C. Chang verified for an actual 90� bend with numerical simulation LBM in the following section, where we remove the nonlinear slip and use the estimated ejecting velocity u c of 29c in an appropriate manner to capture the ejecting flow due to field leakage.
IV. NUMERICAL FLOW SIMULATION The nonlinear ejecting flow due to field leakage is localized to the inner corner and we use 29b or 29c to model its effects. The simulation is done on the remaining domains with a linear slip where field leakage is absent. We hence impose 1 or 22b with 0 of the natural potential. It is necessary to solve the Ohmic electric field without leakage in a bent channel of Fig. 3 with wedge angle and unit width. A numerical Laplace solver9 is used to solve Eq. 9a with no-flux boundary condition 9b on the walls and the far away electric field is set to unity. The magnitude of the tangential electric field along the inner wall and the outer wall of a 90� channel are plotted in Fig. 7 a . As expected, the field near inner corner bears a singularity and decreases monotonically to its far field value of unity over a distance equal to the width of the channel. For the outer wall, field is zero at the corner and increases monotonically to its far field value of unity over a distance equal to the width of the channel. The outer corner does not bear a singularity as the singularity occurs only for wedge angles the angle on the wall side less than 180�. The Ohmic tangential field from 10a is E t A * E o r (1 ) on the inner wall near the wedge where E t must scale linearly with respect to the far-away field E o being set to unity. The magnitude of the coefficient A * is obtained from its asymptotic value as r0 for various wedge angles and it is plotted in Fig. 7 b . It is clear that the coefficient is nearly a constant of order unity A * 1. We have even simulated walls with a finite but small dielectric permittivity. The tangential field away from the corner does not vary much from the solution plotted in Fig. 7 a for perfectly insulating walls. The field profile near the corner shifts downwards as the dielectric constant of the wall increases. We can use A * 1 in 32 to estimate the capture efficiency of an arbitrary wedge. To verify the accuracy of 30 and 31 , we have developed a lattice Boltzmann method LBM code for simulating the flow. It is a mesoscopic simulation technique for fluid flows and other transport phenomena.10 It incorporates the features of density probability distribution, single relaxation time constant, discretized spatial lattice, discrete particle velocities. We refer to the literature for the simulation technique.11 For our purpose, we choose the one that corresponds to an incompressible fluid and that recovers the solution to Navier�Stokes equation. As shown in the schematic Fig. 8, for a 90� bend channel, a square lattice is used with nine velocity directions resolved at each lattice node.11 In total, there are 100 nodes along the width of the channel and 400 nodes along the length. Because the ejecting velocity 29b is obtained via a flow balance within the Debye layer, the slip and ejecting boundary conditions imposed on the wall nodes of the LBM code must be done with care. Linear and nonlinear electro- O PR OF CO PY 02
FIG. 7. a Electric field variation along the length of a 90� channel at inner corner and outer corner. b The magnitude of the harmonic coefficient A * obtained from its asymptotic value at r0 for various wedge angles in radians. kinetic tangential velocities approach constant asymptotic values away from the wall, giving rise to effective universal slip expressions like the linear Smoluchowski slip 1 or 22b and our nonlinear version 22c . This is because the bulk is neutral and the Maxwell stress vanishes away from the polarized layer. The universality means that the slip velocity can be applied anywhere near or on the boundary, as far as the outer flow is concerned. Such asymptotic behavior does not exist for normal ejection velocity from the corner. Due to either radial mass flux considerations or momentum loss, it decays monotonically from the corner. As such a universal effective boundary condition does not exist for the normal ejection velocity, its value depends on where this 22 P 12 HF PROOF COPY 022212PHF PROOF COPY 022212PHF
Phys. Fluids, Vol. 14, No. 12, December 2002 Nonlinear electrokinetic ejection and entrainment 11 FIG. 8. Schematic diagram of slip boundary condition for LBM simulation of electrokinetic flow in a micro-channel with a 90� bend. The channel width W is scaled to unity. boundary condition is imposed for the outer flow and this location cannot be exactly at the wall. It hence becomes a lattice-size-dependent effective normal velocity condition that must be applied more carefully than the universal slip conditions. The wall nodes are located at a half of the lattice spacing from the immediate neighboring nodes as shown on Fig. 8. The first layer of nodes along the wall are assigned a slip velocity given by 22b , after the tangential wall electric field has been computed from the Laplace solver. The only exception to this rule occurs at the corner node. The ejecting velocity u c is assigned to the first node diagonally away from the corner. As seen in Fig. 2, the derived ejecting Debye layer thickness slip velocity u c is at a distance from the tip and tangential on either sides of the wedge. The actual ejecting velocity assigned to the corner node located at a distance / 2 from the corner is not u c but u node such that c the net ejecting flow condition is satisfied. We expect the ejecting flow to be nearly irrotational flow near the corner just outside the Debye layer, as is the case of the classical outer linear electrokinetic flow. However, at a finite distance away, the vortices and back pressure will render it a viscous Stokes flow that is eventually dominated by the harmonic of 30 . Exactly where the transition occurs depends on the ejecting velocity assigned to the corner node, which has yet to be specified. It depends on the local Reynolds number based on the channel width W and the assigned ejecting velocity. Based on these arguments, the radial velocity should decay as r 1 , as is the case for a point source of inviscid flow in two dimensions. Hence, the appropriate ejecting velocity to be assigned to the corner node is u node u c c
PROOF COPY 022212PHF , / 2 33 which results from flow rate balance at a radius of / 2 and assuming that the ejecting flow spans an angle of /2 on either sides of the ejecting line of symmetry as shown in Fig. 2. Its value depends on the grid spacing. In our simulations, is held constant at the value of 0.01W and we then expect from 31 that the entrainment should be approxi)/2 2 u node . mately ( c We impose the ejecting velocity of 29c at the corner node as shown in Fig. 8, directed towards the outer corner. We have performed LBM flow simulation for different values of u node . Hence, each of the solutions and flow profiles c would correspond to different strengths of the corner singularity and field leakage effects. In the LBM simulation, a 1 is used. It corresponds to a relaxation time constant of 1 2 1 /6 3 units as in LBM. The far-field viscosity of velocity is given a small value of u o 0.001. Since the width of the channel consists of 100 nodes, the Reynolds number 1 is the density of of the flow is Re uoW/ 0.3 where the fluid or particles at each node. The creeping flow condition described by the biharmonic equation is hence well approximated in the LBM simulation. Each simulation is run until a steady flow is established or for 104 time steps. Stream function of the flow is constructed from the obtained velocity fields at the end of simulation. The contours of stream lines are plot