Electrochemical MHD for Microfluidic Applications

48th AIAA Aerospace Sciences Meeting Including the New Horizons Forum and Aerospace Exposition 4 - 7 January 2010, Orlando, Florida

AIAA 2010-908

Electrochemical MHD for Microfluidics Applications K. M. Isaac* and Debamoy Sen† and Nicholas Leventis‡ Missouri University of Science & Technology, Rolla, MO, 65409, USA To investigate the feasibility of using magnetic fields for flow control in microfluidic devices, we employed potential sweep and potential step voltammetry in CFD simulations of the electrochemical magnetohydrodynamics (ECMHD) of a redox system. For all cases, the faradaic current at the working electrode was computed. Using time-varying boundary conditions based on the Butler-Volmer electrode kinetics model, we studied the interplay of Lorentz force, convection and redox species concentration distribution. The concentration contours obtained show strong effect of magnetoconvection due to Lorentz force. The evolution of the flow field in a two-dimensional electrochemical cell shows gradual development of a vortex structure with time. The potential step simulations have been validated by comparing with the approximate one-dimensional analytical solution for the diffusion only case. The cyclic voltammograms obtained from the potential sweep simulation have been compared with published experimental data. The results agree well. Simulations of a channelless microfluidic cell show that ECMHD is suitable for pumping, mixing and flow control.

I.

Introduction

E

lectrochemical MHD (ECMHD) based on redox electrode reactions have advantages for introducing and controlling fluid flow in microfluidics and lab-on-a-chip (LOAC) applications.1-18 Pumping, mixing, chemical reactions, sensing, and flow control are among the well known applications. ECMHD pumping is one of the few methods suitable for weakly conductive liquids such as buffer solutions. Microsyringe pumps have moving mechanical parts, and are prone to more frequent malfunctioning and failure. They also don’t fit the microfluidics paradigm because of their large size. Motion in a closed loop in situations in which fresh fluid need not be introduced, is another advantage of ECMHD pumping. MHD pumps with AC and DC electric fields have been proposed. AC causes inductive heating. DC-based MHD pumping has the disadvantage of bubble formation due to electrolysis if the applied voltage is high. Deterioration of electrode surface is another drawback of DC based MHD. ECMHD-based experiments and CFD simulations using redox reactions have been reported recently to avoid bubble formation and electrode deterioration. In redox-based systems, the applied voltage is low, typically in the mV-~1V range. Bubble formation is thus avoided. Electrode deterioration is not a serious issue because the electrodes simply facilitate electron transfer to and from the reduced species (R) and oxidized species (O) at the electrode surface, in contrast to methods in which the electrode material is consumed during the electrochemical reaction (as in electroplating). CFD simulation of ECMHD has been done by Qian and Bau10 and Kabbani et al.17 A rectangular channel with an inlet and outlet flow and equal size electrodes mounted on opposite walls were simulated. The magnetic field was applied parallel to the electrode surfaces and perpendicular to the incoming flow direction so that Lorentz force was in the direction of the main flow resulting in pumping action. Qian and Bau10 considered multiple electrode pairs and presented their results using variables in non-dimensional form. Results are presented for steady state conditions in the above studies. Henely et al.18 performed CFD simulations of a micro electrochemical reactor cell. They used the steady state form of the species conservation equations, and the time-dependent forms for the momentum equations. Information on how the species boundary conditions were implemented at the electrode surface is not given.

*

Professor, Mechanical & Aerospace Engineering, 400 W. 13th St, 271 Toomey Hall, AIAA Associate Fellow, [email protected]. † Graduate research assistant, Mechanical & Aerospace Engineering. ‡ Curators’ Professor, Chemistry Department. 1 American Institute of Aeronautics and Astronautics Copyright © 2010 by K. M. Isaac. Published by the American Institute of Aeronautics and Astronautics, Inc., with permission.

II.

Theory

The model to simulate the MHD phenomena comprises the oxidation and reduction of an electroactive species at the electrodes and the convection induced by Lorentz force through interaction of faradaic current and external magnetic field. The redox couple can be represented by the electron-transfer reaction O + ne ↔ R (1) where O and R represent the oxidized and reduced species, respectively, e is the electronic charge and n is the number of electrons exchanged in the reaction. For such systems, the faradaic current density is represented by the equation (2) j = nFDO (∇CO ) electrode

where j is the faradaic current density, F is the Faraday constant, A is the area of the electrode, DO is the diffusion coefficient of species O. The current i can be calculated by integrating the current density over the electrode area. If we assume that current density is a constant over the electrode area, the current is given by i = nFADO (∇CO ) electrode (3) The term in parenthesis in Eqs. 2 and 3 is the concentration gradient of species O at the electrode surface. Thus, we find that the faradaic current is related to the concentration gradient of the redox species. In order to obtain the concentration gradient, we need to solve the species equation ∂Ci (4) + ∇ ⋅ (CiV ) = ∇ ⋅ ( Di ∇Ci ) ∂t where subscript i represents species O and R in the medium and Ci is the concentration of species i. Appropriate initial and boundary conditions must be implemented for the solution of the system of equations. The evolution of the concentration profiles of the electroactive species is governed by Eq. 4. For the redox system of Eq.1, Eq. 4 represents two species conservation equations for O and R. If initially both oxidized and reduced species are present, then the initial conditions become (5) CO (t = 0) = CO∗ (6) C R (t = 0) = C R∗ where the superscript * denotes the initial concentrations. We assume that the concentrations in the bulk remain at their initial values, not affected by the heterogeneous electrode reactions. To apply the appropriate electrode surface boundary conditions, we employ heterogeneous surface reactions under the Butler-Volmer (B–V) model which describes the kinetics of the electrode reactions.19,20 Using the B–V formulation we calculate flux from the surface reaction in Arrhenius form, given by DR

∂CR ∂ξ

= − DO electrode

∂CO ∂ξ

= CO electrode

electrode

k 0 exp(−α

nF η ) − CR RT

electrode

k 0 exp((1 − α )

nF η) RT

(7)

Where α is the charge transfer coefficient (ranging from 0.0 to 1.0), R is the universal gas constant, T is the temperature, η is the overpotential and k0 is standard reaction rate constant. We consider isothermal (no joule heating) conditions of operation, and hence do not include the temperature dependence in the Arrhenius form of the equation. Thus T is a constant in the above equation. The first and second terms on the right hand side of Eq. 7 represent the forward and backward reaction rates, respectively. Thus, the reaction rates depend on the concentrations of the electroactive species at the electrode surface and on the overpotential η. For the redox reaction represented by Eq. 1, the following condition also applies ∂C ∂C DR R = − DO O (8) ∂ξ electrode ∂ξ electrode where ξ is the direction normal to the electrode surface. The solution of the species conservation equations allows calculating the faradaic current and from it the Lorentz force, FL , which is obtained by taking the r cross product of ion flux vector, j , and magnetic flux vector, B . ( FL = j × B ) . The direction in which this Lorentz force acts is given by the right hand rule. This force acting on a volume element, ΔV, of the 2 American Institute of Aeronautics and Astronautics

electrolytic conductor, accelerates it and displaces in the direction of the force. A void is created, which gets replenished by fresh electrolyte rushing in from the bulk, setting up a convective motion. The Lorentz force appears as a body force term in the Navier–Stokes (N–S) momentum equations. Solution of the N–S equations gives the flow field and thereby the velocity in the medium, which in turn influences how the concentration profiles of the species develop. This in turn affects the faradaic current generated in the system and ultimately affects the Lorentz force. The flow and species conservation equations need to be solved in a coupled manner. Figure 1 illustrates the interaction between the electrochemistry and hydrodynamics. The conservation of mass and the momentum equations are as follows. ∇ ⋅V = 0 (9)

ρ

DV = −∇p + μ∇ 2 V + j × B + ρ g Dt

(10)

where, V = u eˆ x + v eˆ y + weˆ z is the fluid velocity, u, v and w are the velocity components in the x, y and z directions, respectively, ρ is the density, p is the pressure, μ is the dynamic viscosity, and g is the gravity. For a two–dimensional problem, Eq. 10 has two components in the x and y directions.

Figure 1. Simulation flow chart showing the interaction between the electrochemistry and hydrodynamics

Numerical Modeling and Simulation We use the commercially available software package Fluent22 for numerical solution of the governing equations. Fluent employs a finite volume method in which the nodes are at the center of the finite volumes, and for these nodes the conservation equations are written and discretized in their integral form. An unstructured grid technology is used, by which the shape of the elements can be quadrilaterals and triangles for the two–dimensional simulations, and hexahedrals and tetrahedral for three-dimensional simulations. Using unstructured mesh allows accurate simulation of many problems of practical interest. For our study, we chose the Fluent preprocessor Gambit as the meshing tool. We used a rectangular electrochemical cell model for the present simulations. The cell extends to +/-

∞ in the z direction (Figure 2a), thus allowing a two–dimensional formulation. The solution is identical in planes normal to the z axis, i.e. the current is only in the xy plane and there is no flow in the zdirection. The working electrode is a band electrode, placed at the center of the bottom surface as shown in Figure 2a. The entire top surface is the counter electrode. The meshed two–dimensional computational domain is shown in Figure 2b. The mesh has 20,000 quadrilateral cells, distributed non-uniformly with higher grid density near the working electrode. This is because the ionic flux drops off quickly away from the working electrode, and hence Lorentz force is significant only near the working electrode. To capture the fine features of the processes taking place near the working electrode, we used a densely clustered grid near the working electrode. 3 American Institute of Aeronautics and Astronautics

A time-accurate unsteady formulation is used to simulate the time–varying phenomena. A flow diagram of the major steps in the simulation is shown in Figure 3. The coupled mass, momentum and species conservation equations are solved at each time level, till the user–specified convergence criteria are met. The solution then advances to the next time level. The conservation equations differ from the basic fluid flow equations because of the electrochemistry and the presence of the magnetic field. The Lorentz force terms present as body forces in the momentum equations and the coupling of the species boundary condition with the applied transient voltage/current is handled by code written by the user, known in Fluent as User– Defined Functions (UDFs). UDFs allow great flexibility in solving a wide range of phenomena. UDFs are basically C/C++ functions written by the user and linked dynamically with the Fluent flow solver to enhance standard features of the code. These UDFs are compiled and linked during problem set up. Then they are invoked by the solver during the solution process. UDFs allow accessing solver data and perform tasks such as using user-specified source terms and special boundary conditions. In our case, using UDFs, we modify the governing equations to include source terms such as Lorentz force, and define dynamic species boundary conditions in terms of the time-varying applied potential/current. UDFs also allow the user’s model to include solution–dependent properties such as viscosity of non–Newtonian fluids. Auxiliary equations such as Laplace’s and Poisson’s, that govern electric and magnetic fields, can also be solved using this framework in a coupled manner, and made part of the overall solution. Using this approach we have successfully solved a wide range of problems including adsorption,23 fluid-structure interaction,24 multiphase flow,25,26 and thermophoresis.27

(a) (b) Figure 2. Schematic diagram of model two-dimensional electrolytic cell

As shown in the flow diagram in Figure 3, the major steps in the solution procedure involve generating the mesh file and reading it into the flow solver. This is followed by the case setup for the particular problem. This step involves appropriate selection of models and specification of the boundary conditions. Next, we initialize the problem by specifying the initial molar concentrations of the species. III. Results and Discussion The geometry of the computational domain used for the present simulations is shown in Figure 2. The domain size is 2x1 cm2. We used a 2 mm wide band working electrode, which is centrally placed at the bottom surface. The entire top surface acts as counter electrode. For the redox electrolytic solution, we used the redox species –tetramethyl–phenylenediamine (TMPD). The supporting electrolyte is tetrabutylammoniumperchlorate i.e., TBAP in acetonitrile. When setting up the problem, we define three species, the oxidized (O) and reduced (R) forms of TMPD and a bulk species for the supporting electrolyte. Molecular weight of TMPD is 164.25 kg / kmol and that of TBAP is 341.92 kg / kmol. The dynamic

4 American Institute of Aeronautics and Astronautics

viscosity of the redox species is 0.0015 kg/(ms) and that of the supporting electrolyte is 0.0005179 kg/(ms). The densities of the redox species and supporting electrolyte are 992 and 826.1 kg/m3, respectively. The diffusivities of species O and R are 2.4x10 -9 m2/s and 2.3x10 -9 m2/s, respectively. The initial concentration of reduced species was 10.89 mM (milli-moles/liter) and that for supporting electrolyte was 0.5 M (moles/liter). The initial concentration of oxidized species was set to 0 mM, i.e., initially we have only reduced species.

Figure 3. Flow diagram showing main solution steps We apply the magnetic field in the positive z-direction, with a magnitude of 1.75 Tesla. Since we are employing the Butler–Volmer model of electrode kinetics, in which we assume it to have an Arrhenius form, we take the standard reaction rate constant as 0.55 cm/s.21 For the potential step mode, we apply an overpotential (η) = 0.7V, and for the potential sweep mode, the overpotential is swept in the -0.25V to 0.4V range. Referring to Eq. 7, we see that η and k0 are the two controlling parameters that govern the working electrode surface flux boundary conditions. Since we neglect joule heating, the effect of temperature is not considered in this study; we ran the simulations at the standard temperature, T = 25oC. At high negative η, the forward reaction dominates and as we move towards a more positive η the backward reaction dominates. Note for the potential step case, in the beginning, the backward reaction dominates as the stepped overpotential is positive. However, for the potential sweep case, though the overpotential is swept from a negative to a positive value, the forward reaction cannot dominate because the initial concentration of oxidized species is zero. While the swept overpotential becomes more and more positive, the backward reaction becomes more and more significant. Hence, we can conclude that as time progresses, based on the applied overpotential and the initial concentrations of the electroactive species, the contributions of the forward and the backward reactions will vary. We have used zero flux species boundary condition at the counter electrode, because of the larger area of the counter electrode compared to the working electrode (ACE >> AWE), and the current density at the counter electrode is negligible. A flux boundary condition similar to the one used for the working electrode can be used when the working electrode and the counter electrode are similar in size. We have used an integration time step Δt = 0.01 s for our simulations, which is well below the time constant of the process. Other values of Δt were studied to establish Δt-independence of the solution. The results obtained from the simulations were post-processed to get velocity vectors, velocity magnitude contours, and concentrations of redox species as they evolved with time. We also obtained profile plots of the species concentration at different times along a line perpendicular to the working electrode surface for comparison to the one-dimensional analytical solution in order to establish the accuracy of the numerical solution. First, we solve the potential step mode with and without the magnetic field. In the absence of magnetic field, i.e., for the diffusion–only problem, we allow the solution to proceed for a short time, 6 s, and 5 American Institute of Aeronautics and Astronautics

compare the concentration profiles with the analytical solution of the one-dimensional diffusion problem, given in terms of the error function. Finally, we sweep the potential, both in the absence and presence of the magnetic field. We then compare the corresponding cyclic voltammograms (i-E curves) to the available experimental results from literature. We also note the effects of sweep rates and standard reaction rate constant k0 on the nature of the i-E curve.

(a)

(b)

Figure 4. Profile plot of species R normal to working electrode at (a) t = 2 s, and (b) 6 s. B = 0 a. Potential Step Voltammetry Analytical solution of the one-dimensional (1D) diffusion problem is obtained as follows. The 1D diffusion equation is a second order PDE

∂CR ( x, t ) ∂ 2CR ( x, t ) = DR (11) ∂t ∂x 2 which requires one initial and two boundary conditions to be solved. In the present case, initially only * species R is present with concentration CR . The first boundary condition is that the bulk electrolyte is unperturbed by the redox reaction at the working electrode (far field condition). For the second boundary condition, we consider the overpotential to be stepped to the mass transfer–controlled regime, i. e., the current is diffusion limited, and concentration of species becomes effectively zero at the working electrode. In such cases, using Laplace transform, the concentration profile of species can be expressed in terms of the following error function ⎡ ⎤ x (12) CR ( x, t ) = CR* erf ⎢ 1/ 2 ⎥ 2( ) D t ⎣ ⎦ R We stepped the overpotential to a large value of 0.7 V, so that the concentration at the working electrode approaches zero and we approach mass–transfer controlled regime. We obtained concentration profiles at different times from 0 to 6s, and compared them to the above equation. Figure 4a shows the profile plot comparisons that we have obtained at time t = 2 s. Here the numerical and the 1D analytical solutions agree well. The differences between the two close to the electrode can be attributed to the approximate boundary condition used for the analytical solution. In the simulations, we don’t assume the concentration to go to zero at the electrode. The agreement becomes better as time proceeds, as shown in Figure 4b, at 6 s.


For the numerical solution of the diffusion–only problem, we also calculated the current density at the working electrode. In the current density vs. time curve in response to the potential step, we see that the current density drops with time as expected. To verify the current density obtained from numerical simulation, we compared it with the diffusion–limited current given by the Cottrell equation19 (13) j (t ) = jd (t ) = nFDR1/ 2 C R* π 1/ 2 t 1/ 2 We find that the numerically obtained current density curve follows the same trend as in Eq. 13. However, the current density from numerical solution is slightly higher. This can be attributed to the fact that, for the numerical simulation we have used the two–dimensional model. Hence, it is expected that the species concentration would fall off faster at the working electrode because of the additional degree of freedom in space. With the higher concentration gradient, the current density obtained would also be higher. These results help validate the simulations, and the solution can now proceed for cases with the magnetic field applied. b. Cyclic Voltammetry

Figure 5. Triangular potential sweep waveform for cyclic voltammetry Cyclic voltammetry (CV) is the most widely used technique to characterize redox systems, providing qualitative and quantitative information about electrochemical reactions. Cyclic voltammetry (CV) provides considerable information on the thermodynamics of redox processes and the kinetics of heterogeneous electron–transfer reactions. CV consists of linearly scanning the overpotential of a working electrode using a triangular waveform, as shown in Figure 5. Present simulations are for one CV cycle. Figure 6 shows the cyclic voltammogram obtained by simulation, at a sweep rate of 20 mV/s, in the absence and presence of a magnetic field of intensity 1.75 Tesla (T), and a standard reaction rate constant k0 = 0.55 cm/s. In the absence of magnetic field, as the overpotential sweeps from a negative to a positive value, the current density rises, reaches its peak and then begins to drop. Under magnetic field, the current density behavior is that of the steady state mass transport–limited case. This cyclic voltammogram generated numerically is compared with that obtained experimentally by Leventis et al.3,4 shown in Figure 6b. Considering the numerical and experimental voltammograms in the absence of magnetic field, we find that they have comparable peak current densities. The experimental anodic and cathodic peak currents are 38x10-6 A and -11.44 23 x10-6 A, respectively. The above current values correspond to current densities of 18.9 A/m2 and A/m2, respectively, for an electrode area of 0.0201 cm2. The numerical peak current densities obtained from our simulations are 16.25 and -11.58 A/m2. The two compare well. The differences in the current density values can be attributed to the fact that the experimental data correspond to a three–dimensional domain with a circular working electrode, whereas the simulations are for a two–dimensional band electrode geometry.


(a) (b) Figure 6. Cyclic voltammograms obtained from simulations (a) and experiments (b). Sweep rate = 20 mV/s.

(a) (b) Figure 7. Cyclic voltammograms obtained for different sweep rates (a) and after scaling using square root of sweep rate.

In addition to the baseline simulations at 20 mV/s sweep rate and k0 = 0.55 cm/s, we have also varied both the sweep rate and rate constant values, in the absence of magnetic field, and observed expected current response in the simulations. When different sweep rates are employed, we see that each curve has the same shape, but the current density increases with increasing sweep rate (Figure 7a). This behavior can be explained by the size of the diffusion layer and the time taken to complete the sweep. As current density is proportional to the flux at the electrode surface, we conclude that the flux at the electrode surface is considerably smaller at slow sweep rates than it is at faster sweep rates. We also scaled the CV 8 American Institute of Aeronautics and Astronautics

by dividing the current density by the square root of sweep rate. We see that the scaled CVs collapse nearly on top of each other (Figure 7b). This suggests that the scaled voltammograms are nearly sweep rate independent. By varying the standard reaction rate constant in the simulations, we have been able to observe its effect on the voltammogram. The maximum current density locations are different for different values of k0. At smaller values of k0, peak–to–peak separation is larger. Also, decreasing k0 lowers the magnitude of the maximum current density (Figure 8). Keeping other parameters unchanged, we have run other simulation cases with different initial concentrations of the reduced species, and observed that the maximum current density can be increased with increase in initial species R concentration.

Figure 8. Effect of standard state reaction rate constant on voltammogram. We have also run various potential sweep cases in the presence of magnetic field. The cyclic voltammograms obtained under magnetic field are found to be sweep rate independent. The effect of concentration on the limiting current density from our simulations are in agreement with experimental results–increasing with increase in concentration of the electroactive species. Also, we have seen that the limiting current density depends strongly on the magnetic field intensity–increasing the intensity increases the value of the limiting current density. Preliminary results have been obtained for a channelless microfluidic cell. The results (not presented here) show that redox-based ECMHD is viable for pumping, mixing and flow control in microfluidic devices.

IV.

Conclusions and Future Work

We have successfully performed numerical simulation of the transient electrochemical magnetohydrodynamics (ECMHD) in a redox electrochemical cell. The model used to solve the transient ECMHD phenomena has been presented. The present work not only aids our understanding of the ECMHD phenomena, but also gives us a picture of the flow field and concentration field as it develops with time. Lorentz force is found to exert great influence on the process. The good agreement between the numerical solution and the approximate analytical solution, and experimental results from literature validates the present approach. Using non–rectangular mesh with clustering near the working electrode gives highly accurate numerical solutions. The small size of the working electrode compared to the counter electrode is more similar to experimental geometries. The present simulations have no externally-imposed flow, which helps better observation of the ECMHD phenomena. Three–dimensional simulations for other geometries such as circular electrodes and band electrodes are being planned in order to make quantitative comparisons with the magnetic field coupled electrochemistry experiments. We also plan to extend this work on ECMHD to labon-a-chip using micro–electrode arrays in microfluidic systems. The work will also be extended to flow control in microfluidic devices. The numerical model developed here provides the framework for the above studies.


Acknowledgments We acknowledge support from the National Science Foundation Grant, CHE-0719097. We also acknowledge fruitful discussions with Ingrid Fritsch, Melissa Weston and Matthew Gerner.

References 1.

Sen, D.; Isaac, K. M.; Leventis, N., "Modeling and Numerical Simulation of Magnetic Field Coupled Electrochemical Processes". ECS Transactions 2008, 13, (16), 33-43.

2.

Sen, D.; Isaac, K. M.; Leventis, N., "Investigation of Electrochemical Magnetohydrodynamic Phenomena using Numerical Simulations". J. Electrochem. Soc. 2009 (in preparation).

3.

Leventis, N.; Chen, M.; Gao, X.; Canalas, M.; Zhang, P., "Electrochemistry with Stationary Disk and Ring-Disk Millielectrodes in Magnetic Fields". J. Phys. Chem. B 1998, 102, 3512-3522.

4.

Leventis, N.; Gao, X., "Steady-State Voltammetry with Stationary Disk Millielectrodes in Magnetic Fields: Nonlinear Dependence of the Mass-Transfer Limited Current on the Electron Balance of the Faradaic Process". J. Phys. Chem. B 1999, 103, 5832-5840.

5.

Fritsch, I.; Carter, C. S.; Aguilar, Z. P. Microfluidics and Small Volume Mixing Based on Redox Magnetohydrodynamics Methods. 10/317,777 (Patent App. No.); 7,147,441.

6.

Fritsch, I.; Carter, C. S.; Aguilar, Z. P. Microfluidics and Small Volume Mixing based on Redox Magnetohydrodynamics Methods. 10/026,748 (Patent App. No.); 6,733,244.

7.

Henry, C. S.; Vandaveer IV, W. R.; Mubarak, I.; Gray, S. R.; Fritsch, I. In Self-Contained Microelectrochemical Detectors for Analysis in Small Volumes of Static and Flowing Fluids, Micro Total Analysis Systems 2001, Monterey, CA, 2001; Ramsey, J. M.; Berg, A. v. d., Eds. Kluwer Academic Publishers: Monterey, CA, 2001; pp 321-322.

8.

Aguilar, Z. P.; Vandaveer, W. R.; Fritsch, I., "Self-Contained Microelectrochemical Immunoassay for Small Volumes using Mouse IgG as a Model System". Anal. Chem. 2002, 74, 3321-3329.

9.

Gerner, M. D. Factors affecting redox magnetohydrodynamics for flow in small volumes. M.S., University of Arkansas, Fayetteville, 2009.

10. Qian, S.; Bau, H. H., "Magneto-hydrodynamics based microfluidics". Mech. Res. Commun. 2009, 36, 10-21. 11. Jang, J.; Lee, S. S., "Theoretical and Experimental Study of MHD (Magnetohydrodynamic) Micropump". Sensors and Actuators A 2000, 80, 84-89. 12. Eijkel, J. C. T.; Dalton, C.; Hayden, C. J.; Burt, J. P. H.; Manz, A., "A circular ac magnetohydrodynamic micropump for chromatographic applications". Sens. Actuators B 2003, B92, 215-221. 13. Bau, H. H. In A case for mageto hydro dynamics (MHD), 2001 ASME International Mechanical Engineering Congress and Exposition, New York, NY, 2001; ASME: New York, NY, 2001; pp 1-8. 14. Zhong, J.; Yi, M.; Bau, H. H., "Mageto hydrodynamic (MHD) pump fabricated with ceramic tapes". Sens. Actuators A 2002, 96, 59-66. 15. Qian, S.; Zhu, J.; Bau, H. H., "A stirrer for magnetohydrodynamically controlled minute fluidic networks". Phys. Fluids 2002, 14, (10), 3584-3592. 16. Qian, S. and Bau, H. H., “Magnetohydrodynamic Flow of Redox Electrolyte,” Phys. Fluids, 17, 2005, 067105, 1-15. 17. Kabbani, H.; Wang, A.; Luo, X.; and Qian, A.; “Modeling Redox-based Magnetohydrodynamics in threedimensional microfluidic channels,” Phys. Fluids, 19, 083604, 1-12. 18. Henley, I. E.; Yunus, K.; and Fisher, A. C., “ Voltammetry under Microfluidic Control: Computer-Aided Design Development and Application of Novel Microelectrochemical Reactors,” J. Phys. Chem. B 2003, 107, 3878-3884. 19. A. J. Bard and L. R. Faulkner, Electrochemical Methods, p. 152, 2nd Edition, John Wiley and Sons, 2001.


20. J. Newman and K. E. Thomas – Alyea, Electrochemical Systems, 3rd Edition, John Wiley and Sons, 2004. 21. A. D. Clegg, N. V. Rees, O. V. Klymenko, B. A. Coles and R. G. Compton, “Experimental Validation of Marcus Theory for Outer-Sphere Electron-Transfer Reactions: The Oxidation of Substituted 1,4Phenylenediamines,” ChemPhysChem, 1234 – 1240, 5, 2004. 22. FLUENT v12, 2009, Ansys Inc. 23. Shivaram, P.; Bai, X.; Benne, K.; Isaac, K.M.; and Banerjee, R.,“Experimental and Numerical Investigation of n-butane Adsorption in a Carbon Canister,” Int. J. Heat Mass Transfer(under revision). 24. Isaac, K. M.; Rolwes, J.; and Colozza, A., “Aerodynamics of a Flapping and Pitching Wing using Simulations and Experiments,” AIAA Journal, Vol. 46, No. 6, Jun. 2008, 1505-1515. 25. Banerjee, R. and Isaac, K. M.,“A Study to Determine Vapor Generation from the Surface of Gasoline Flowing in an Inclined Channel Using a Continuous Thermodynamics Approach,” Numerical Heat Transfer: Part A: Applications, Vol. 50, Number 8, November 2006, pp. 705–729. 26. Banerjee, R. and Isaac, K. M., “An Algorithm to determine the Mass Transfer Rate from a Pure Liquid Surface using the Volume of Fluid Multiphase Model,” International Journal of Engine Research, Vol. 5, No. 1, 2004, pp. 23-37. 27. A. Janakiraman, K. M. Isaac, P. Whitefield and D. Hagen, “A Numerical Thermophoretic Model for Nanoparticle Deposition in an Aerosol Sampling Probe,” AIAA-2005-5351, 17th AIAA Computational Flow Dynamics Conference; Toronto, Canada, June 2005.
