Competition between Induced-Charge Electro-Osmosis and Electro-Thermal Effects around a Weakly-Polarizable Microchannel Corner
11 Competition between Induced-Charge Electro-Osmosis and Electro-Thermal Effects around a Weakly-Polarizable Microchannel Corner
Matan Zehavi, Alicia Boymelgreen and Gilad Yossifon*
Faculty of Mechanical Engineering, Micro- and Nanofluidics Laboratory, Technion–Israel Institute of Technology, Technion City 32000, Israel
The microchannel corner is a common inherent component of most planar microfluidic systems and thus its influence on the channel flow is of significant interest. Application of an alternating current electric field enables quantification of the non-linear induced-charge electro-osmosis (ICEO) ejection flow effect by isolating it from linear electro-osmotic background flow which is present under dc forcing. The hydrodynamic flow in the vicinity of a sharp channel corner is analyzed using experimental micro-particle-image-velocimetry ( µ PIV) and numerical simulations for different buffer concentrations, frequencies and applied voltages. Divergence from the purely ICEO flow with increasing buffer conductivity is shown to be a result of increasing electro-thermal effects due to Joule heating. *Corresponding author: [email protected] I. INTRODUCTION
Previous experimental studies of induced-charge-electro-osmosis (ICEO) around weakly polarizable (i.e. dielectric) structures embedded within an electrolyte have focused mainly on the direct current (DC) case [1]–[7]. At the same time, a number of theoretical studies of time dependent ICEO have been performed for perfectly conducting cylinders in both step-wise and alternating current (AC) fields [8], lossless dielectric structures using macroscale time-dependent boundary conditions [9], lossy dielectric cylinders under AC fields [10] and dipolophoresis of Janus particles [11]. In contrast to “alternating-current electroosmosis” (ACEO) which describes hydrodynamic flow generated by polarization of the active electrodes themselves[12], ICEO, whether under DC or AC fields, refers to a floating structure which is polarized by a field applied at external electrodes[8]. In the current work we focus on ICEO processes occurring at a sharp corner of a L-shaped microchannel junction. This work extends previous studies [5]–[7] focused on colloid and fluid flow dynamics in the same microchannel design, but utilizes AC rather than DC applied fields. In contrast to the DC case where the effects of natural and induced zeta potentials are superimposed, in the AC case, only the latter non-linear effect that has a non- vanishing time-averaged effect. In addition it will be shown that Joule heating resulting in electro-thermal effects [13] becomes significant at sufficiently large solution conductivity/applied electric fields. In the following, we describe the experimental methods in sec. II, the numerical model in sec. III, the results and discussion in sec. IV and concluding comments in sec. V.
II.
EXPERIMENTAL METHODS A.
Fluidic device and experimental setup
The fluidic device is made of two reservoirs connected by a L-junction micro-channel, (Figure 1) 140µm in depth (normal to the plane of view). The micro-channels were fabricated from PDMS (Polydimethylsiloxane, Dow-Corning Sylgard 184) using the rapid prototype technique [14]. A high resolution chrome mask with a 3µm resolution was used for the SU8 mold. The nominal L-junction corner has a radius of curvature is 6±2 µm. The channel was bonded to a microscope slide coated with a 30µm layer of PDMS using plasma bonding[15]. An electric function generator (Agilent 33250A) was connected through an amplifier (TREK 2220) to platinum wire electrodes (0.5mm platinum wire, Sigma-Aldrich) which were inserted in the reservoirs. The electric signal was monitored using an oscilloscope (Tektronix TPS-2024).
Figure 1 : (color online) (a) Schematic design of the microfluidic chip and scaled up microscopic image of the L-junction; (b) experimental setup (not to scale). All dimensions are in µ m. B. µ -PIV measurements KCl solutions of 3 different concentrations (solution conductivities of 2.21 ± ± µ PIV system (TSI) installed on a Nikon TI microscope. The PIV images were obtained with a Litron Nano S65-PIV Laser (YAG 532nm wavelength, 400mJ, 4ns at 15Hz) through a Nikon Plan-Fluor 20x/0.50 and a dichroic filter (560nm dichroic/565nm emission). The time interval between images was 10ms, with a repetition rate of 5 image pairs per second, and 50 image pairs per measurement set. The images were then passed through a linear Gaussian filter (Filter size 5x5 pixels, standard deviation parameter 0.5) and background subtraction. The processing interrogation region is 32x32 pixels on a standard Nyquist grid, with the results presented as a sum of correlations of the 50 image pairs (see [16] for more details on µ PIV technique). C. Temperature visualization
Ross et al. [17] developed a Rhodamine B based thermometry technique for direct measurements of the in-channel temperature profiles in microfluidic systems. Rhodamine B is a fluorescent dye whose quantum yield is strongly dependent on temperatures in the range of 0°C to 100°C, making it ideal for liquid based systems. Here, the neutral (zwitterionic) fluorescent dye Rhodamine B (Sigma 83690) at 10 -4 M in concentration was added to KCl solutions of different conductivities 2.1, 50.3, 351µS/cm, corresponding to concentrations of 1.4x10 -5 , 3.3x10 -4 , 2.3x10 -3 M, respectively. The chip was placed on a Nikon TI microscope and illuminated using a fluorescent lamp (Nikon-Intensilight C-HGFL, wavelength 380-600nm, 130Wx50Hz). Images were captured using Andor-NEO camera, through a Nikon Plan-Fluor 10x/0.30 objective lens and an Omega Optical (560nm dichroic/565nm emission) dichroic filter. Ambient temperature of C ° was measured on the glass slide surface using a K-type thermocouple. III.
NUMERICAL MODEL A.
Electrostatics
The AC electrostatic problem for the channel geometry shown in Figure 1(a) is solved numerically using finite-element based commercial software (COMSOL). Herein, we follow the approach taken in Yossifon 2009 [9] wherein macroscale boundary conditions for the time-dependent problem (eqs. 2.15, 2.16 in [9]) were derived in a way that eliminated the need to directly resolve transport within the electric double layers (EDL). A symmetric electrolyte ( z z z + − = − = ) with equal ionic diffusivities (
D D D + − = = ) is considered along with the standard assumptions of a thin EDL ( a δ λ = ≪ ) and weak electric fields ( zF RT ψ φ = ≪ ). Herein, λ denotes the Debye length, a the geometric length scale of the problem chosen to be the corner radius of curvature, F Cmol − = ⋅ is the Faraday number, R JK mol − − = the universal gas constant, T K = ° the absolute temperature and φ the characteristic induced electric potential. The latter was set [9] to ( ) w f w E a a φ ε λ ε ε λ = + , wherein f ε and w ε denote the dielectric constants of the electrolyte solution and solid wall, respectively and E the characteristic magnitude of the external field within the narrow microchannel. For the current geometry, ( ) ( )
21 2 hVL h h E V m ∆ + = ≈ ⋅ , where h m µ = , h m µ = , L mm = so that for V V ∆ = , mV φ ≈ in rough agreement with the weak field assumption. Within this framework, the hydrodynamic and electrostatic problems are decoupled. Furthermore, the electrostatic solution can be separated into three sub-domains: the ‘outer’ electroneutral solid wall and bulk solutions with their harmonic electric potentials w φ and f φ , respectively, and the ‘inner’ domain of the diffuse EDL whose electric potential φ satisfies Poisson’s equation. The macroscale boundary conditions connecting w φ and f φ on the solid-electrolyte interface, constituting the main result of [9], were obtained by focusing on the inner domain. The linearity of the present problem allows for the superposition of the linear equilibrium (exclusively associated with a native surface charge/zeta potential) and nonlinear induced charge solutions. However, when using AC electric forcing, the former linear solution integrates to zero over a given period so that in the sequel we need focus only on the induced-charge distribution. For the special case of harmonic AC forcing, i.e. j t E e ω = within the narrow channel away from the junction and normalized by E , the potentials within the solid wall and bulk solution can be expressed as j tf f e ω φ φ = and j tw w e ω φ φ = , respectively. Herein, f φ and w φ are the non-dimensional complex amplitudes (phasors) normalized by φ , j is the imaginary number, ω is the non-dimensional frequency normalized by τ and t is the time normalized by the EDL diffusion time D τ λ = . Both potentials obey the Laplace equation f φ ∇ = , w φ ∇ = , (1) along with boundary conditions wf w n φαφ φ γ ∂− = ∂ , (2a) w f n n φ φα γδ γ ∂ ∂ = ∂ − ∂ , (2b) which were directly derived from Eqs. 2.15’, 2.16’ [9] for the special case of harmonic forcing in the limit of infinite time (i.e. when the transient phase of the solution vanishes to obtain a periodic solution). Herein, j j j fD D λ λγ ω ω π = + = + = + ɶɶ , f ɶ is the frequency and n is the spatial coordinate normal (in the outward sense) to the surface of the solid, normalized by a . The tilde stands for dimensional quantities. The parameter ( ) w f α ε ε δ = , in the terminology of an equivalent RC-circuit model, represents the ratio of the capacitance of the dielectric solid ~ w a ε to that of the EDL ~ f ε λ . Eqs.(2a,b) also satisfy the extended boundary conditions of Zhao and Yang[10] in the limit of an ideally dielectric wall (i.e. lossless dielectric), which for the case of PDMS used in the current study (material conductivity O(10 -10 µS/cm[18]), constitutes a reasonable assumption. We note that substitution of (2a) into (2b) yields the charging equation for the EDL ( ) fw f n φδζ φ φ γ γ ∂= − = − ∂− , (3) where ζ is the induced part of the zeta potential phasor in terms of the electric potential in the fluid domain. For the DC case ( ) γ = eqs. (2a,b) reduce to f n φ ∂ =∂ , (4a) ww f n φφ α φ ∂+ =∂ , (4b) which is in agreement with previous results [5], when the contribution of the equilibrium zeta potential is neglected. For the sake of comparison with other related studies, the above results can be examined within certain limits. For example, in the case of a conductive cylinder (i.e. w φ = ) subject to an externally applied field ˆ j t e ω z (normalized by E , wherein ˆ z is the normal unit vector of the axisymmetry axis of coordinate ( ) cos z r θ = ), the solution for the potential phasor is of the form [8] ( ) ( ) ( )
20 0 cos 1 / f E a r g r ω φ φ θ = − + , (5) with the coordinate r normalized by the cylinder radius a and ( ) ( ) g j j ωδ ωδ − − = − + . From eqs. (3) and (5) one obtains the induced zeta potential phasor ( ) ( ) cos2 1 E a j θζ φ γ ωδ − = + , (6) which for small enough frequencies (i.e. ω ≪ ), i.e. j γ ω = + ≅ , stands in agreement with eq. 4.9 in Squires & Bazant [8]. B. Thermal analysis
For sufficiently high solution conductivity, [ ] Sm σ − , and/or large electric fields, Joule heating effects [19], [20] may become significant due to generation of a large power density within the fluid domain [ ] q Wm σ − = E ɶɺ , (7) with σ being the solution conductivity, f φ = −∇ E ɶɶ ɶ is the electric field within the fluid domain and the time average. The temperature field is obtained by solving the energy equation with the internal Ohmic heat source q ɺ p p f DT Tc c T k T qDt t ρ ρ ∂= + ⋅ ∇ = ∇ + ∂ u ɶ ɶ ɶ ɶ ɶ ɶɶ ɺɶ ɶ , (8) where kg m ρ = , p c J kgK = ⋅ and f k J mKs = correspond to the fluid mass density, specific heat (at constant pressure) and thermal conductivity of the fluid, respectively, and D Dt ɶ is the material derivative. Here, we have neglected the viscous dissipation term, which is substantially smaller than the Joule heating term [13]. An order of magnitude analysis of (8) gives the diffusion time over which thermal equilibrium is established as diff T t H s α = ≈ ɶ , where T f p k c m s α ρ − = = ⋅ is the thermal diffusivity and H m µ = the channel depth over which the largest temperature gradients exist. Thus it is evident that for sufficiently large frequencies (i.e. diff diff t Hz ω ω π = = ɶɶ ɶ≫ ) the differential temperature change, diff T T ω ω ∆ ≈ ɶ ɶ ɶ ɶ , will be negligible. Hence, for diff t t > ɶ ɶ and applied frequencies diff ω ω ɶ ɶ≫ steady-state conditions can be assumed while solving for the time-averaged temperature field. Additionally, we note that the heat convection term is negligible relative to the diffusion term in Eq.(8) , i.e. p f c T k T ρ ⋅ ∇ ∇ u ɶ ɶ ɶ ɶɶ ≪ . Their ratio can be shown to scale as p f c uH k L ρ and for our system, this value is always much smaller than one - even in the case of a relatively large characteristic velocity of u mm s ≈ and channel length L mm = , p f c uH k L ρ ≈ . Thus, the effect of fluid flow on the temperature field can be neglected and eq.(8) can be simplified to f k T q ∇ + = ɶ ɶ ɺ . (9) Unlike momentum and species transport analysis, which is confined to the fluidic domain, thermal analysis presents a unique challenge as thermal diffusion necessarily extends the simulation domain from the region of interest (i.e. the fluidic domain) to encompass a significant portion, if not the entire chip. The consideration of the thermal diffusion process through the solid however, poses significant computational problems as it introduces larger length and time scales as well as three-dimensionality. In order to circumvent this and avoid the “whole chip” approach [19] in which the thermal analysis is performed on the entire 3D chip structure (solid walls and fluidic channel), we will formulate an equivalent 2D model in terms of ( ) , T T x y = ɶ ɶ (instead of ( ) , , T T x y z = ɶ ɶ ), in which heat transfer in the z-direction and within the channel walls is included via an effective resistance at the top and bottom substrates. Using an infinitesimal cuboid control volume ( ) , , dx dy H the following 2D version of the heat equation (8) can be derived ( ) f US BS T TT Tk qx y H R R ∞ − ∂ ∂+ + − + = ∂ ∂ ɶ ɶɶ ɶ ɺɶ ɶ , (10) where T C ∞ = ° ɶ is the environment temperature, US US PDMS PDMS
R t k h m KW − = + = is the thermal resistance (per unit area) of the upper substrate ( US ) consisting of a conduction ( ) US PDMS t k and convection ( ) h resistances in series, while BS BS PDMS PDMS BS Glass Glass
R t k t k m KW − − = + = ⋅ is the thermal resistance (per unit area) of the bottom substrate ( BS ) consisting of two conduction resistances in series. The thickness of the various layers (see Fig.1) are , US PDMS t mm = , , BS PDMS t m µ = , , BS glass t mm = , with the corresponding thermal conductivities PDMS k Wm K − − = [19] and glass k Wm K − − = [13]. Following [19], a natural convection heat transfer coefficient of ~ 10 h Wm K − − is used based on the classical configuration of a “heated upper plate”. This value is within the range of typical natural heat convection coefficients Wm K − − − [21], all of which would result in US BS
R R ≫ . Since these thermal resistances are connected in parallel (Eq.(10)) the conduction resistance through the bottom substrate BS R is dominant. Note that in contrast to the upper surface where heat transfer by convection was considered, an isothermal condition is implemented at the bottom surface since it is in contact with the microscope stage which can be assumed to be at a constant temperature. Since most of the steady state heat rejection is through the lower substrate [19] instead of the chip side walls, we assume that ( ) SW f
T T R k T n ∞ − = − ∂ ∂ ɶ ɶ ɶ ɶ , (11) where the resistance (per unit area) of the microchannel side walls is approximated SW BS
R R ≈ and n ɶ is the coordinate normal to the wall pointing into the wall domains. C. Hydrodynamics
Governing equations
The time-averaged velocity of the fluid u is governed by the Stokes equation in the low Reynolds number limit ( ~ 0.1 uH ρ µ for typical values for velocity ~ 1 u mm s and dynamic viscosity ~ 10 kgm s µ − − − ) E p µ −∇ + ∇ = u + f ɶɶ ɶɶ ɶ , (12) together with the mass conservation equation for an incompressible fluid ∇ ⋅ = u ɶ ɶ , (13) where p ɶ is the pressure and E f ɶ is the electric force density outside the EDL. The effect of E f ɶ within the thin EDL is captured in the following slip velocity boundary condition. Hydrodynamic boundary conditions
The fluid velocity u satisfies the no-penetration condition and the quasi-steady Helmholtz-Smoluchowski slip condition [22] { } { } || Re Re slip u ζ = ⋅ = − u t E , (14) which is strictly valid for frequencies much smaller than that corresponding to the viscous relaxation within the EDL ( ns λ ν ≈ where for aqueous solutions we take the kinematic viscosity cm s ν − − ≈ and ~ 10 nm λ ). Here the slip velocity is normalized by f E ε ε φ µ , t is the unit vector tangential to the wall, || = ⋅ E E t is the component of the electric field tangential to the wall in the fluid domain and { }
Re ... stands for the real part of { } ... . The induced zeta potential is derived from (3) as ( ) j t j tw f e e ω ω ζ ζ φ φ = = − . (15) Thus the time-averaged slip velocity is ( ) || || cos2 slip u ζ ζ = − − E E ∡ ∡ , (16) with ... being the absolute value and ∡ the phase angle of the complex amplitudes. Electrothermal body force
To first order approximation (obtained by performing a perturbative expansion assuming the relative changes in the permittivity, ε , and conductivity, σ , are small for typical temperature increments), the time averaged electrical force density for harmonic AC forcing, with the electric field expressed as j t e ω = E E ɶɶ ɶ ɶ is ([13], [20]) ( ) ( ) E A B T A Tj σε εσ ωε −= ∇ ⋅ − ∇ + f E E E ɶ ɶ ɶ ɶ ɶ ɶ ɶ ɶɶ , (17) 0 where A T εε ∂ = ∂ , B T σσ ∂ = ∂ and * indicates complex conjugate. For an aqueous KCl solution A K − ≈ − and B K − ≈ [13]. However, within this work, we chose to treat , A B as fitting parameters between the numerical and experimental results (see Results and Discussion section). The force density is frequency dependent and has two distinct limits. At low frequencies ( ω σ ε ɶ ≪ ) the Coulomb force (1 st term of the right hand side of (17)) dominates since, for an aqueous solution, the relative change in the conductivity is greater than that of the permittivity (i.e. B A > ). At high frequencies, ω σ ε ɶ ≫ and the dielectric force dominates. For a typical experimental value of ~ 2 S cm σ µ with f ε ε ε = wherein
12 1 10
CV m ε − − − = ⋅ is the permittivity of vacuum and the relative dielectric constant of water is ~ 80 f ε , one obtains MHz σ ε ∼ . D. Summary of equations
Following the above set of approximations, the electrical, thermal and mechanical problems are decoupled and we solve the governing equations as follows. First we solve the electrical problem consisting of the Laplace equation (1) for the potentials in both the fluid and wall domains, f φ ɶ and w φ ɶ respectively, with the boundary condition (2a,b) to relate them. Second we calculate the temperature field T ɶ in the fluid (eqs.(7), (10) and (11)) using the solution for the electric field f φ = −∇ E ɶɶ ɶ . Finally, we use these two solutions to compute the electrical force density (17) and solve the Navier–Stokes’ equation (12-13) for the fluid velocity, u ɶ , and pressure, p with the ICEO slip velocity boundary condition (16). In order to facilitate the numerical simulations we have chosen to resolve a computational domain in the area of interest that is much smaller than the full chip dimensions (Fig.1), wherein the lengths of the narrow (width h m µ = ) and wide (width h m µ = ) are taken as l mm = and l mm = , respectively, instead of L mm = . The electric potential drop between the opposite microchannel entrances can then be approximated by ( ) ( ) ( ) ' 1 2 2 1 2 1 ~ V V l h l h Lh Lh ∆ ∆ + + , which in the case of an applied potential drop of
V V ∆ = yields ' V V ∆ ≈ . A more accurate value of ' V V ∆ ≈ is obtained through solving the Laplace equation f φ ∇ = in the full chip geometry, which accounts for 1 the curvature of the corner, with an electric potential drop of V V ∆ = between the microchannel ends with an insulation condition at all other solid-electrolyte interfaces. The radius of curvature of the L-junction corner used in the simulations was 6µm.
IV.
RESULTS AND DISCUSSION
Temperature distributions for different solution conductivities
Based on the measured fluorescence intensity (see supporting materials [23]), Figure 2 illustrates that the average temperature within the fluid increases with increasing solution conductivity [24] since the Joule heating effect is linearly dependent on the medium conductivity ( q E σ ∝ ɺ ). More specifically, the time evolution of the averaged normalized temperature within the narrow (Fig.2b) and wide (Fig.2c) channels shows a sharp increase (decrease) upon application (shut off) of the electric field at t=10s (t=30s). It is also clear that the Joule heating effect is more pronounced within the narrow channel compared to the wide one, in accordance with the existence of larger electric fields within the former. The 2D intensity maps at time t=30s (just before shut off) clearly show the existence of sharp axial temperature gradients at the vicinity of the channel corner, whereas an almost uniform temperature distribution exists along the channel’s axial direction away from the corner. This is a clear indication that most of the generated heat dissipates orthogonally to (rather than along) the channel’s axial direction, thus, supporting the main assumption of the above thermal theoretical analysis, i.e., the consideration of heat dissipation through the top and bottom surfaces of the microfluidic channels (3 rd term on the r.h.s. of eq.(10)) as well as through the side walls (boundary condition eq.(11)). An order-of-magnitude estimate of the incremental temperature rise (within the narrow channel) can be made based on eqs. (7) and (9) as ( )
22 1 2 f T Vk H L h h σ ∆ ∆≈ + ɶ or ( )
21 2 f V HT k L h h σ ∆ ⋅∆ ≈ + ɶ , where V ∆ is the potential difference across the electrodes, h m µ = and h m µ = are the narrow and wide channels widths, respectively. For an applied potential amplitude of V V ∆ = at 300Hz and with
2, 50, 350
S cm σ µ = the temperature rise can be approximated as T C ∆ ≈ ° ɶ , respectively. This stands in reasonable agreement with the experimental results in Figure 2. Figure 2: (color online)
Experimental visualization of the temperature field using Rhodamine B. (a) The interrogation windows used for generating the time evolution of the average temperature at the narrow (b) and wide (c) channels at various solution conductivities. The applied voltage amplitude of 2000V and 300Hz was introduced at t=10s and shut off at t=30s. Two-dimensional temperature maps at t=30s (just before shutting off the voltage) for various solution conductivities: (d) 2.1µS/cm; (e) 50.3µS/cm; and (f) 351µS/cm.
The numerically calculated temperature distribution that is shown in Fig.3 qualitatively supports the experimental findings (Fig. 2), such that indeed most of the axial temperature gradients exist at the vicinity of the corner along with elevated temperatures within the narrow channel. The latter are in good quantitative agreement with the experimental observations (Fig.2). The magnitude of these temperature gradients implies that they cannot be neglected and will give rise to a significant electrothermal body force (Fig.3(b)), which in turn, affects the local flow pattern as discussed in more detail in the following sections. T [ ° C ] S [ µ m] Figure 3: (color online)
Numerical simulation results at a fixed frequency of 300Hz and 2000V applied voltage showing the: (a) temperature profiles for various conductivities along the corner wall versus the coordinate S (see inset); (b) electro-thermal body force (vector plot) and temperature (surface plot) distributions for 350µS/cm. Flow patterns for different solution conductivities
The µ PIV generated velocity fields for solutions of varying conductivity, at a fixed voltage (2000V) and frequency (300Hz) are shown in Fig.4. In contrast to previous studies ([6], [7]) on a similar device but performed under DC field conditions, where a superposition of linear electroosmotic and ICEO flows must be accounted for, here only the ICEO contribution is seen as the linear contribution is time-averaged to zero. For a solution of low S (a) (b)
4 conductivity ( S cm µ ), the flow field clearly exhibits the expected non-linear ICEO jetting flow, stemming from the polarized induced zeta-potential around the corner (Fig.5), along with a vortex pair that is formed due to conservation of mass[5]. The maximum value of the induced zeta potential at the largest applied voltage, 2000V in amplitude is ~10mV, hence, smaller than the thermal potential. As the conductivity of the solution is increased, a decrease in the calculated effective zeta potential is observed, which corresponds with previous experimental results [25]. At the same time, a clear divergence of the flow pattern from the expected ICEO jetting is obtained. The onset of this divergence corresponds with the increase in temperature gradients within the electrolyte as illustrated in Figure 2. By incorporating electrothermal effects in the numerical solutions (Fig.4 (d)-(f)) the obtained flow patterns are in good qualitative agreement with the experimental observations. In the sequel, we further assess the influence of Joule heating, by examining the voltage and frequency dependent behavior of electrolytes of varying conductivity. Figure 4 : (color online) (a)-(c) µPIV generated flow fields for a 2000V voltage amplitude, 300Hz signal, and varying solution conductivities exhibit divergence from the expected ICEO corner ejection flow with increasing conductivity. Arrows: velocity vectors. Color: angular momentum parameter - swirl intensity [26]; (d)-(f) The resulting hydrodynamic flow patterns
5 (velocity streamlines) for various electrolyte conductivities indicating a clear divergence from the classical ICEO ejection flow with increasing conductivity. Arrows: velocity vectors.
Figure 5 : (color online) Numercially calculated induced zeta-potential around the corner (inset) for a 2000V voltage amplitude, 300Hz signal, and varying solution conductivities.
Voltage scan
Fig.6a shows that the measured ejection velocity at the lowest conductivity is quadratic in the electric field at a fixed frequency. This stands in agreement with the theoretical predictions of ICEO phenomenon [8] and previous experimental studies [6]. The measured velocity is obtained by averaging the µ PIV velocities in an interrogation box of 30 µ m X 30 µ m in size (see Fig.4a), where we differentiate between its absolute magnitude and its component projected on the corner bisector (with its positive sense pointing into the electrolyte). As seen for σ µ = both experimental and calculated values collapse onto the line of purely ICEO response (i.e. without electro-thermal effects). This result verifies that indeed the ICEO flow consists of strong ejection flow along the corner bisector and that at such low conductivity electro-thermal effect is negligible. It is noted that in order to obtain an agreement between the numerical solution and the experiments, a fitted value for -10-50510-500 -250 0 250 500 ζ [ m V ] S [ µ m] S
6 the effective wall dielectric constant of w ε = (see Fig.S2 in supporting materials [23]) is used rather than the material property for PDMS of w ε = . The need for an artificially increased wall permittivity, implies that the PDMS wall behaves as if it is more polarizable than its material property suggests. The underlying mechanism behind why the numerical simulations underestimate the experimental results is not clear and it is in opposite trend from what is commonly reported in the literature regarding ICEO over conductors (or ideally polarizable), wherein the theoretical and numerical results tend to overestimate the experiments by at least an order of magnitude [25], [27], [28]. Several physical mechanisms have been suggested to address the latter, such as surface conduction [27], non-linear capacitance and steric effects associated with the induced-charge [29], [30]. However, all of these mechanisms add corrections to the numerical solution so as to decrease the resulting ICEO velocities and hence are likely not relevant in the current study wherein we are facing an opposite problem of how to increase the numerically predicted velocities. Extension of the weak field approach ( RT zF φ ≪ ), underlying eqs.(2a,b), to strong fields ( ~ RT zF φ ) about a dielectric surface of zero surface charge [31] also cannot explain this discrepancy, and on the contrary, shows a transition from an E dependence at moderate fields to an essential linear variation with E at strong fields. In line with the decreased induced zeta potential (Figure 5) and the departure of the flow pattern from classical ICEO ejection (Figure 4) observed with increased conductivity, we also see a marked divergence between the absolute value of the measured velocity and its projected component. In fact, for
373 S / cm σ µ = , the projected component even reverses sign and instead of pointing outwards from the corner, it points inward. At the same time, although at low voltages the absolute velocities are smaller than the case of σ µ = , at large enough voltage, they exceed even the velocities of the σ µ = solution. Overall, for
373 S / cm σ µ = , the scaling of the velocity with voltage changes from quadratic to biquadratic (i.e. to the fourth power) as shown in the inset of Fig.6c, in agreement with the theoretical scaling of the electro-thermal effects [32]. The dominance of the electro-thermal effects at high conductivity solutions, which are amplified at high applied voltages are further verified by the numerical simulations. -1.5E-040.0E+001.5E-043.0E-040.0E+00 2.0E+06 4.0E+06 V e l o c i t y [ m / s ] V [V ] Exp: 50uS normExp: 50uS projSim: 50uS no ETSim: 50uS norm ETSim: 50uS proj ET -1.5E-040.0E+001.5E-043.0E-040.0E+00 2.0E+06 4.0E+06 V e l o c i t y [ m / s ] V [V ] Exp: 2uS normExp: 2uS projSim: 2uS no ETSim: 2uS norm ETSim: 2uS proj -1.5E-040.0E+001.5E-043.0E-040.0E+00 2.0E+06 4.0E+06 V e l o c i t y [ m / s ] V [V ] Exp: 373uS normExp: 373uS projSim: 373uS no ETSim: 373uS norm ETSim: 373uS proj ET (a) (b) (c) Figure 6 : (color online) Measured (symbols) and calculated ejection velocity versus increasing electric potential amplitude (squared) for different solution conductivities at 300 Hz. Continuous lines stand for the numerical model without electro-thermal (ET) effects; dashed line stand for the norm of the velocity; and dotted line stand for the projected (along the corner bisector) component of the velocity. As seen there is an increasing divergence between the norm and projected values with increasing conductivity. A fourth power dependence of the velocity on voltage is shown in the inset of part (c) indicating the existence of strong electro-thermal effects. Error bars indicate the measured standard deviation.
Frequency scan
Here we concentrate on studying the frequency dependence of the ICEO ejection velocity. Based on Fig.6 we use the component of the velocity projected onto the corner bisector to differentiate between ICEO and electro-thermal effects (see also supporting materials [23]). Figure 7 depicts the dependency of the ejection velocity on the applied frequency and exhibits for the case of low conductivity solution (2µS/cm) both low and high frequency dispersion with a maxima around 200-400Hz. Applying a numerical model, which considers the wall to be a lossless dielectric with the fitted effective wall dielectric constant of w ε = , we observe similar high frequency dispersion behavior as in the experiments but also record an order of magnitude difference between the experimental and numerical RC times for the low conductivity solution. The greater discrepancy of two orders of magnitude in RC time for the intermediate conductivity case is likely also influenced by the non-negligible electro-thermal effects. 9 Figure 7 : (color online) Ejection velocity versus frequency as measured (symbols with error bars of a single standard deviation) and numerically calculated (continuous lines) for the case of low conductivity (2µS/cm) solution and applied voltage amplitude of 2000V. Both measured and calculated ejection velocities are the projected velocity along the corner angle bisector averaged within an interrogation window of 30X30 µ m (see Fig.4a). V. CONCLUSIONS
As opposed to their DC counterpart [7], the use of AC fields facilitates a direct experimental observation of the non-linear ICEO corner jetting flow by eliminating the linear EOF effects (which have a zero time-average contribution). The ICEO flow pattern was observed using µ PIV, clearly showing the corner ejection flow along with the associated vortex pair in excellent qualitative agreement with the theoretical predictions. Quantitative analysis of the ejection velocity demonstrated the expected linear scaling with the square of the applied electric field. In addition both low and high frequency dispersions, associated with the relaxation times of the induced EDLs at the active electrodes and polarized corner, respectively, were observed. Numerical calculations yielded the same high frequency dispersion trend but with an order of magnitude larger RC time. Curiously, we obtained an opposite trend from what is commonly reported for ICEO over conducting (ideally -1.0E-052.0E-04 100 1000 10000 100000 E j e c t i on v e l o c i t y [ m / s ] f [Hz] Exp: 2uS projectedSimulation: 2uS projected ETExp: 50uS projectedSimulation: 50uS projected
0 polarizable) medium. Specifically, the numerical solution of the induced velocities underestimated the experimental data by an order of magnitude in the case of a dielectric wall, whereas in the case of conducting surfaces, the numerics tend to overestimate the experiments by at least an order of magnitude. Non-linear physical phenomenon corrections, such as surface conduction, non-linear capacitance and steric effects which tend to decrease the predicted velocity cannot explain this trend and thus the resolution of this discrepancy is left for future study. The strong divergence of the hydrodynamic flow pattern (from the classical ICEO ejection) with increasing electrolyte conductivity is proven to be a result of electro-thermal effects due to Joule heating by several means. Experimentally, two-dimensional mapping of the temperature distribution, using Rhodamine B, is used to illustrate an overall increase in solution temperature with increasing conductivity. Additionally, it is demonstrated that since the electric field within the narrow channel is larger than that in the wide channel, the former exhibits higher temperatures. These results were in excellent agreement with the numerical solution of the 2D heat equation we have formulated (Eq.(10)) which accounts for heat losses through all electrolyte-solid wall interfaces. Furthermore, accounting for the electro-thermal body force in the numerical simulations resulted in quantitative agreement with the experimental data in terms of both the flow patterns and the velocity dependence on the voltage including the shift from a quadratic to biquadratic scaling with increasing dominance of the electrothermal effects. Since the electro-thermal effects are arise from temperature gradients, which in turn stem from the electric field gradients, suppression of the latter (e.g. uniform channels and/or rounding of the corner) can minimize these electro-thermal effects. Thus, the current study besides providing, an in-depth description of the competition between ICEO and electrothermal effects, can also be used to guide the design of planar microfluidic systems.
ACKNOWLEDGMENTS
This work was supported by ISF Grant No. 1078/10. We thank the Technion RBNI (Russell Berrie Nanotechnology Institute) for their financial support. BIBLIOGRAPHY [1] F. Nadal, F. Argoul, P. Kestener, B. Pouligny, C. Ybert, and A. Ajdari, “Electrically induced flows in the vicinity of a dielectric stripe on a conducting plane,”
Eur. Phys. J. E , vol. 9, no. 4, pp. 387–399, Nov. 2002. [2] S. K. Thamida and H.-C. Chang, “Nonlinear electrokinetic ejection and entrainment due to polarization at nearly insulated wedges,”
Phys. Fluids , vol. 14, no. 12, pp. 4315–4328, Nov. 2002. [3] P. Takhistov, K. Duginova, and H.-C. Chang, “Electrokinetic mixing vortices due to electrolyte depletion at microchannel junctions,”
J. Colloid Interface Sci. , vol. 263, no. 1, pp. 133–143, Jul. 2003. [4] F. Zhang and D. Li, “Induced-charge electroosmotic flow around dielectric particles in uniform electric field,”
J. Colloid Interface Sci. , vol. 410, pp. 102–110, Nov. 2013. [5] G. Yossifon, I. Frankel, and T. Miloh, “On electro-osmotic flows through microchannel junctions,”
Phys. Fluids , vol. 18, no. 11, pp. 117108–117108–9, Nov. 2006. [6] Y. Eckstein, G. Yossifon, A. Seifert, and T. Miloh, “Nonlinear electrokinetic phenomena around nearly insulated sharp tips in microflows,”
J. Colloid Interface Sci. , vol. 338, no. 1, pp. 243–249, Oct. 2009. [7] M. Zehavi and G. Yossifon, “Particle dynamics and rapid trapping in electro-osmotic flow around a sharp microchannel corner,”
Phys. Fluids 1994-Present , vol. 26, no. 8, p. 082002, Aug. 2014. [8] T. M. Squires and M. Z. Bazant, “Induced-charge electro-osmosis,”
J. Fluid Mech. , vol. 509, pp. 217–252, Jun. 2004. [9] G. Yossifon, I. Frankel, and T. Miloh, “Macro-scale description of transient electro-kinetic phenomena over polarizable dielectric solids,”
J. Fluid Mech. , vol. 620, pp. 241–262, 2009. [10] C. Zhao and C. Yang, “ac electrokinetic phenomena over semiconductive surfaces: Effective electric boundary conditions and their applications,”
Phys. Rev. E , vol. 83, no. 6, p. 066304, Jun. 2011. [11] A. M. Boymelgreen and T. Miloh, “Induced-charge electrophoresis of uncharged dielectric spherical Janus particles,”
ELECTROPHORESIS , vol. 33, no. 5, pp. 870–879, 2012. [12] A. González, A. Ramos, N. Green, A. Castellanos, and H. Morgan, “Fluid flow induced by nonuniform ac electric fields in electrolytes on microelectrodes. II. A linear double-layer analysis,”
Phys. Rev. E , vol. 61, no. 4, pp. 4019–4028, Apr. 2000. [13] N. G. Green, A. Ramos, A. González, A. Castellanos, and H. Morgan, “Electrothermally induced fluid flow on microelectrodes,”
J. Electrost. , vol. 53, no. 2, pp. 71–87, Aug. 2001. [14] J. R. Anderson, D. T. Chiu, R. J. Jackman, O. Cherniavskaya, J. C. McDonald, H. Wu, S. H. Whitesides, and G. M. Whitesides, “Fabrication of Topologically Complex Three-Dimensional Microfluidic Systems in PDMS by Rapid Prototyping,”
Anal. Chem. , vol. 72, no. 14, pp. 3158–3164, Jul. 2000. [15] K. Haubert, T. Drier, and D. Beebe, “PDMS bonding by means of a portable, low-cost corona system,”
Lab. Chip , vol. 6, no. 12, pp. 1548–1549, Nov. 2006. 2 [16] S. T. Wereley and C. D. Meinhart, “Recent Advances in Micro-Particle Image Velocimetry,”
Annu. Rev. Fluid Mech. , vol. 42, no. 1, pp. 557–576, 2010. [17] D. Ross, M. Gaitan, and L. E. Locascio, “Temperature Measurement in Microfluidic Systems Using a Temperature-Dependent Fluorescent Dye,”
Anal. Chem. , vol. 73, no. 17, pp. 4117–4123, Sep. 2001. [18] J. E. Mark,
Polymer Data Handbook , 2 edition. Oxford ; New York: Oxford University Press, 2009. [19] D. Erickson, D. Sinton, and D. Li, “Joule heating and heat transfer in poly(dimethylsiloxane) microfluidic systems,”
Lab. Chip , vol. 3, no. 3, pp. 141–149, Aug. 2003. [20] A. Ramos, H. Morgan, N. G. Green, and A. Castellanos, “Ac electrokinetics: a review of forces in microelectrode structures,”
J. Phys. Appl. Phys. , vol. 31, no. 18, p. 2338, Sep. 1998. [21] F. P. Incropera, D. P. DeWitt, T. L. Bergman, and A. S. Lavine,
Introduction to Heat Transfer , 5 edition. Hobokenm NJ: Wiley, 2006. [22] J. Lyklema,
Fundamentals of Interface and Colloid Science: Solid-Liquid Interfaces , 1st ed. Academic Press, 1995. [23] Supporting Materials [24] G. Y. Tang, D. G. Yan, C. Yang, H. Q. Gong, C. J. Chai, and Y. C. Lam, “Joule heating and its effects on electroosmotic flow in microfluidic channels,”
J. Phys. Conf. Ser. , vol. 34, no. 1, p. 925, Apr. 2006. [25] M. Z. Bazant, M. S. Kilic, B. D. Storey, and A. Ajdari, “Towards an understanding of induced-charge electrokinetics at large applied voltages in concentrated solutions,”
Adv. Colloid Interface Sci. , vol. 152, no. 1–2, pp. 48–88, Nov. 2009. [26] L. Graftieaux, M. Michard, and N. Grosjean, “Combining PIV, POD and vortex identification algorithms for the study of unsteady turbulent swirling flows,”
Meas. Sci. Technol. , vol. 12, no. 9, p. 1422, Sep. 2001. [27] G. Soni, T. M. Squires, and C. D. Meinhart, “Nonlinear Phenomena in Induced Charge Electroosmosis,” pp. 761–771, Jan. 2007. [28] C. K. Harnett, J. Templeton, K. A. Dunphy-Guzman, Y. M. Senousy, and M. P. Kanouff, “Model based design of a microfluidic mixer driven by induced charge electroosmosis,”
Lab. Chip , vol. 8, no. 4, pp. 565–572, Mar. 2008. [29] M. Z. Bazant, M. S. Kilic, B. D. Storey, and A. Ajdari, “Nonlinear electrokinetics at large voltages,”
New J. Phys. , vol. 11, no. 7, p. 075016, Jul. 2009. [30] L. H. Olesen, H. Bruus, and A. Ajdari, “ac electrokinetic micropumps: The effect of geometrical confinement, Faradaic current injection, and nonlinear surface capacitance,”
Phys. Rev. E , vol. 73, no. 5, p. 056313, May 2006. [31] O. Schnitzer and E. Yariv, “Strong electro-osmotic flows about dielectric surfaces of zero surface charge,”
Phys. Rev. E , vol. 89, no. 4, p. 043005, Apr. 2014. [32] A. Castellanos, A. Ramos, A. González, N. G. Green, and H. Morgan, “Electrohydrodynamics and dielectrophoresis in microsystems: scaling laws,”