Joule Heating Induced Transient Temperature Field and Its Effects on

7598

Langmuir 2005, 21, 7598-7607

Joule Heating Induced Transient Temperature Field and Its Effects on Electroosmosis in a Microcapillary Packed with Microspheres Y. Kang, C. Yang,* and X. Huang School of Mechanical and Aerospace Engineering, Nanyang Technological University, 50 Nanyang Avenue, Republic of Singapore 639798 Received January 8, 2005. In Final Form: June 2, 2005 The Joule heating induced transient temperature field and its effect on the electroosmotic flow in a capillary packed with microspheres is analyzed numerically using the control-volume-based finite difference method. The model incorporates the coupled momentum equation for the electroosmotic velocity, the energy equations for the Joule heating induced temperature distributions in both the packed column and the capillary wall, and the mass and electric current continuity equations. The temperature-dependent physical properties of the electrolyte solution are taken into consideration. The characteristics of the Joule heating induced transient development of temperature and electroosmotic flow fields are studied. Specifically, the simulation shows that the presence of Joule heating causes a noticeable axial temperature gradient in the thermal entrance region and elevates a significant temperature increment inside the microcapillary. The temperature changes in turn greatly affect the electroosmotic velocity by means of the temperaturedependent fluid viscosity, dielectric constant, and local electric field strength. Furthermore, the model predicts an induced pressure gradient to counterbalance the axial variation of the electroosmotic velocity so as to maintain the fluid mass continuity. In addition, under specific conditions, the present model is validated by comparing with the existing analytical model and experimental data from the literature.

1. Introduction Electroosmosis is originated from the electric charge separation that occurs at the fluid-solid interface because of complex electrochemical reactions. When in contact with a charged solid surface, the free ions in the polar medium will be preferentially distributed within an extremely thin region close to the interface. A compact layer of immobile counterions are attracted and adhere to the charged solid surface. Next to the compact layer is a diffuse layer where the counterions dominate over co-ions to neutralize the surface charge. Therefore, this thin region is called the electric double layer (EDL), and apparently, the local net charge density in the EDL is not zero.1,2 Upon application of a tangential electric field, the Coulombic body force causes the net migration of the mobile ions in the EDL. The momentum is transported to the adjacent and bulk liquid by viscosity, resulting in an electroosmotic flow (EOF). Because of its pressure-building capacity, electroosmosis has been extensively applied in variety of microfluidic devices for biochemical analyses such as liquid pumping,3,4 sample injection5-7 and focusing,8,9 and micromixing.10-12 Compared to the traditional mechanical * Author to whom correspondence should be addressed. Email: [email protected]. Telephone: (+65) 6790-4883. Fax: (+65) 67911859. (1) Probstein, R. F. Physicochemical Hydrodynamics: An Introduction; John Wiley and Sons: New York, 1994. (2) Masliyah, J. H. Electrokinetic Transport Phenomena; AOSTRA Technical Publication Series No. 12; AOSTRA: Edmonton, 1994. (3) Harrison, D. J.; Fluri, K.; Seiler, K.; Fan, Z.; Effenhauser, C. S.; Manz, A. Micromachining a miniaturized capillary electrophoresis-based chemical analysis system on a chip. Science 1993, 261, 895-897. (4) Conlisk, A. T.; McFerran, J.; Zheng, Z.; Hansford, D. Mass transfer and flow in electrically charged micro- and nanochannels. Anal. Chem. 2002, 74, 2139-2150. (5) Ermakov, S. V.; Jacobson, S. C.; Ramsey, J. M. Computer simulations of electrokinetic injection techniques in microfluidic devices. Anal. Chem. 2000, 72, 3512-3517. (6) Sinton, D.; Ren, L.; Xuan, X.; Li, D. Effect of liquid conductivity difference on multicomponent sample injection, pumping, and stacking in microfluidic chips. Lab Chip 2003, 3, 173-179.

pumping mechanisms, the electroosmotic pumping system has no mechanical friction, which greatly increases its durability and reduces the fabrication difficulty. A precise amount of liquid can be maneuvered by easily controlling the external electric field. Moreover, the electroosmotic flow velocity has a uniform pluglike profile, which can greatly reduce the solution dispersion and facilitate the sample detection.13 However, because the electroosmosis is electric fielddriven flow, the Joule (resistive) heat generation is inevitable because of the electric current incurred. Previous studies have shown that there are some major negative effects associated with the Joule heating phenomenon. The Joule heating induced temperature changes will in turn affect the EOF and thus the electrokinetic mass transfer by influencing the temperature-dependent electrical/transport properties such as viscosity, dielectric constant, and electric conductivity of the working fluid.14-18 In the application of electrokinetic separation, the tem(7) Zholkovskij, E. K.; Masliyah, J. H.; Czarnecki, J. Electroosmotic dispersion in microchannels with a thin double layer. Anal. Chem. 2003, 75, 901-909. (8) Schrum, D. P.; Culbertson, C. T.; Jacobson, S. C.; Ramsey, J. M. Microchip flow cytometry using electrokinetic focusing. Anal. Chem. 1999, 71, 4173-4179. (9) Fu, L.-M.; Yang, R. J.; Lee, G.-B. Electrokinetic focusing injection methods on microfluidic devices. Anal. Chem. 2003, 75, 1905-1910. (10) Oddy, M. H.; Santiago, J. G.; Mikkelsen, J. C. Electrokinetic instability micromixing. Anal. Chem. 2001, 73, 5822-5832. (11) Stroock, A. D.; Dertinger, S. K.; Ajdari, A.; Mezic, I.; Stone, H. A.; Whiteside, G. M. Chaotic mixer for microchannels. Science 2002, 295, 647-651. (12) Erickson, D.; Li, D. Influence of surface heterogeneity on electrokinetically driven microfluidic mixing. Langmuir 2002, 18, 18831892. (13) Bousse, L.; Cohen, C.; Nikiforov, T.; Chow, A.; Kopf-Sill, A. R.; Dubrow, R.; Parce, J. W. Eelectrokinetically controlled microfluidic analysis systems. Annu. Rev. Biophys. Biomol. Struct. 2000, 29, 155181. (14) Tang, G. Y.; Yang, C.; Chai, C. J.; Gong, H. Q. Modeling of electroosmotic flow and capillary electrophoresis with the Joule heating effect: The Nernst-Planck equation versus the Boltzmann distribution. Langmuir 2003, 19, 10975-10984.

10.1021/la050061g CCC: $30.25 © 2005 American Chemical Society Published on Web 06/29/2005

Electroosmosis in Microcapillary with Microspheres

perature elevation due to Joule heating gives rise to band spreading or peak broadening in capillary electrophoresis and therefore reduces the separation efficiency.19,20 The significant temperature rise may also lead to the decomposition of the thermal labile samples21 and even cause vapor bubbles. In addition, thermal convection characteristics in combined EOF (in the presence of Joule heating) and pressure-driven flow were analytically and numerically examined using classic idealized thermal boundary conditions.22-25 The aforementioned Joule heating studies, however, exclusively focus on the electroosmotic flow in unpacked capillaries. Because the electroosmosis is an interfacial phenomenon, use of porous structures can greatly enlarge the interfacial area and thus significantly enhance the pressure-building capacity of electroosmotic flows. Taking advantages of such high pressure-building capacity of the electroosmotic flow in packed capillaries, two technologies have been developed, including the capillary electrochromatography for chemical analyses and molecular detection26-28 and the eletrokinetic micropumps.29-33 The first effort made to design an electrokinetic micropump may be attributed to Paul et al.,29 who reported the ability to generate high pressures using an electrokinetic pumping of liquid through porous media. Pressures in excess of (15) Tang, G. Y.; Yang, C.; Chai, C. J.; Gong, H. Q. Joule heating effect on electroosmotic flow and mass species transport in a microcapillary. Int. J. Heat Mass Transfer 2004, 47, 215-227. (16) Tang, G. Y.; Yang, C.; Chai, C. J.; Gong, H. Q. Numerical analysis of the thermal effect on electroosmotic flow and electrokinetic mass transport in microchannels. Anal. Chim. Acta 2004, 507, 27-37. (17) Xuan, X.; Sinton, D.; Li, D. Thermal end effects on electroosmotic flow in a capillary. Int. J. Heat Mass Transfer 2004, 47, 3145-3157. (18) Xuan, X.; Xu, B.; Sinton, D.; Li, D. Electroosmotic flow with Joule heating effects. Lab Chip 2004, 4, 230-236. (19) Knox, J. H. Thermal effects and band spreading in capillary electro-separation. Chromatographia 1988, 26, 329-337. (20) Xuan, X.; Li, D. Joule heating effects on peak broadening in capillary zone electrophoresis. J. Micromech. Microeng. 2004, 14, 11711180. (21) Rush, R. S.; Cohen, A. S.; Karger, B. L. Influence of column temperature on the Electrophoretic behavior of myoglobin and R-lactalbumin in high-performance capillary electrophoresis. Anal. Chem. 1991, 63, 1346-1350. (22) Zhao, T. S.; Liao, Q. Thermal effects on electroosmotic pumping of liquids in microchannels. J. Micromech. Microeng. 2002, 12, 962970. (23) Maynes, D.; Webb, B. W. Fully developed electroosmotic heat transfer in microchannels. Int. J. Heat Mass Transfer 2003, 46, 13591369. (24) Maynes, D.; Webb, B. W. Fully developed thermal transport in combined pressure and electroosmotically driven flow in microchannels. ASME J. Heat Transfer 2003, 125, 889-895. (25) Horiuchi, K.; Dutta, P. Joule heating effects in electroosmotically driven microchannel flows. Int. J. Heat Mass Transfer 2004, 47, 30853095. (26) Ivory, C.; Gobie, W. Continuous counteracting chromatographic electrophoresis. Biotechnol. Prog. 1990, 6, 21-32. (27) Rudeg, S. R.; Basak, S. K.; Ladisch, M. R. Solute retention in electrochromatograph by electrically induced sorption. AIChE J. 1993, 39, 797-808. (28) Cole, K. Preparative concentration and size fractionation of DNA by porous media using combination of flow and low electric field strength. Biotechnol. Prog. 1997, 13, 289-295. (29) Paul, P. H.; Arnold, D. W.; Rakestraw, D. J. Electrokinetic generation of high pressures using porous microstructures. In Micro Total Anal. Syst. 1998, Proc. µTAS 1998 Symp. 1998. (30) Zeng, S.; Chen, C. H.; Mikkelsen, J. C., Jr.; Santiago, J. G. Fabrication and characterization of electroosmotic micropumps. Sens. Actuators, B 2001, 79, 107-114. (31) Zeng, S.; Chen, C. H.; Santiago, J. G.; Chen, J. R.; Zare, R. N.; Tripp, J. A.; Svec, F.; Fre´chet, J. M. J. Electroosmotic flow pump with polymer frits. Sens. Actuators, B 2002, 82, 209-212. (32) Chen, L.; Ma, J.; Tan, F.; Guan, Y. Generating high-pressure submicroliter flow rate in packed microchannel by electroosmotic force: potential application in microfluidic systems. Sens. Actuators B 2003, 88, 260-265. (33) Yao, S.; Hertzog, D. E.; Zeng, S.; Mikkelsen, J. C., Jr.; Santiago, J. G. Porous glass electroosmotic pumps: design and experiments. J. Colloid Interface Sci. 2003, 268, 143-153.

Langmuir, Vol. 21, No. 16, 2005 7599

8000 psi have been achieved using silica capillaries packed with micron-size silica beads. Since then, a variety of electroosmotic micropumps have been designed and fabricated.30-33 In accordance with the experimental development, the theoretical studies on the underlying mechanism of electroosmotic flow in porous media have been widely reported.34-41 However, the investigation of the Joule heating effects on the electroosmotic flow in porous media is limited in the literature known to the authors. Knox19 pointed out that Joule heating imposes a main limitation of the performance for capillary electrochromatography. He derived a simple analytical expression for the cross-stream temperature difference within the core region and the capillary wall under thermally developed, steady-state conditions, but he did not consider the temperature dependence of the liquid electrical/transport properties. Recently, Keim and Ladisch42 developed a two-dimensional transient temperature model for capillary electrochromatography. However, their model does not address the full coupling among the temperature field, flow field, and electrical potential field. Optimal design and precisely controllable operation of electrokinetic micropumps and electrochromatography systems demand that the temporary development of the temperature and velocity fields can be modeled and predicted. This paper, therefore, presents a rigorous model describing the Joule heating induced temporary and spatial distributions of the temperature and its effect on the electroosmotic flow in a microcapillary packed with microspheres. The model includes the energy equations for the solid and fluid phases, the continuity and momentum equations for the EOF in porous media, and the equation for conservation of the electric current density. The coupled governing equations are solved numerically using the finite difference method. 2. Problem Formulation Consider a cylindrical microcapillary densely packed with electrically charged nonconducting microspheres of diameter, dp. The ζ potential of microparticles is ζp. The liquid flowing through the porous structure is assumed to be an incompressible, Newtonian, monovalence (z+ ) z- ) 1) electrolyte of density, Ff, and viscosity, µf. In normal applications, the diameter ratio of the capillary to the packing particles is larger than 50, i.e., dc/dp g 50.43 Thus, (34) Rathore, A. S.; Wen, E.; Horva´th Cs. Electroosmotic mobility and conductivity in columns for capillary electrochromatography. Anal. Chem. 1999, 71, 2633-2641. (35) Liapis, A. I.; Grimes, B. A. Modeling the velocity field of the electroosmotic flow in charged capillaries and in capillary columns packed with charged particles: interstitial and intraparticle velocities in capillary electrochromatography systems. J. Chromatogr., A 2000, 877, 181-215. (36) Tallarek, U.; Scheenen, T. W. J.; Van, As H. Macroscopic heterogeneities in electroosmotic and pressure-driven flow through fixed beds at low column-to-particle diameter ratio J. Phys. Chem. B 2001, 105, 8591-8599. (37) Coelho, D.; Shapiro, M.; Thovert, J. F.; Adler, P. M. Electroosmotic phenomena in porous media. J. Colloid Interface Sci. 1996, 181, 169190. (38) Marino, S.; Coelho, D.; Be´kri, S.; Adler, P. M. Electroosmotic phenomena in fractures. J. Colloid Interface Sci. 2000, 223, 292-304. (39) Yao, S.; Santiago, J. G. Porous glass electroosmotic pumps: theory. J. Colloid Interface Sci. 2003, 268, 133-142. (40) Kang, Y.; Yang, C.; Huang, X. Analysis of the electroosmotic flow in a microchannel packed with homogeneous microspheres under electrokinetic wall effect. Int. J. Eng. Sci. 2004, 42, 2011-2027. (41) Kang, Y.; Yang, C.; Huang, X. AC electroosmosis in microchannels packed with a porous medium J. Micromech. Microeng. 2004, 14, 12491257. (42) Keim, C.; Ladisch, M. Model for temperature profiles in large diameter electrochromatography columns AIChE J. 2003, 49, 402410.

7600

Langmuir, Vol. 21, No. 16, 2005

Kang et al.

(Fcp)s

∂Ts ) ∇‚(ks∇Ts) + q˘ s ∂t

(1)

and, for the fluid phase,

(

(Fcp)f

Figure 1. Schematic illustration of the two modeling subsystems of the packed microcapillary.

the entire porous medium can be assumed to be homogeneous, and the capillary wall effects can be safely neglected.40 The present analysis is based on the volume averaging method,40,41 in which macroscopic variables, such as temperature or Darcy velocity, are defined as an appropriate mean over a sufficiently large representative elementary volume (REV). It is assumed that the resultant macroscopic variables are independent of the size of the REV. The length scale of the REV is much larger than that of the pore scale, but considerably smaller than that of the length scale of the macroscopic flow domain. 2.1. Equations of Energy. In the following, we will derive the governing equation that expresses the first law of thermodynamics for an electroosmotic flow in a microcapillary packed with microspheres. For modeling purposes, we divide the packed capillary into two distinct subsystems (as shown in Figure 1). The fluid and stationary phases (packing particles) form the first subsystem, with axial boundaries consisting of the capillary inlet and outlet. The outer radial boundary of the first subsystem is the capillary inner wall. The second subsystem contains the capillary wall itself and is subject to the free convection at the outer surface and heat conduction through the solid wall. Subsystem I: Homogeneous Porous Media. We start with a simple situation in which the medium is assumed homogeneous. We also assume that there is local thermal equilibrium44,45 so that Ts ) Tf ) T, where Ts and Tf are the temperatures of the solid and fluid phases, respectively. When the system dimension L is much longer than the particle dp (L/dp . 1), and when the variation of temperature across dp is negligible compared to that across L for both the solid and fluid phases, we then can assume that, within a distance dp, both phases are in the local thermal equilibrium (LTE).44 This suggests that heat transfer in the solid and fluid phases takes place in parallel so that there is no net heat exchange from one phase to the other. In a representative elementary volume (REV) of the medium, we have, for the solid phase, (43) Kang, Y. Electroosmotic flow in a microchannel packed with microspheres. Ph.D. Thesis, Nanyang Technological University, 2004. (44) Kaviany, M. Principles of Heat Transfer in Porous Media; Springer: New York, 1995. (45) Hsu, C. T.; Cheng, P. Thermal dispersion in a porous medium. Int. J. Heat Mass Transfer 1990, 33, 1587-1597.

)

∂Tf +U B p‚∇Tf ) ∇‚(kf∇Tf) + µf Ω + q˘ f ∂t

(2)

where the subscripts s and f refer to the solid and fluid phases, respectively, (Fcp)s and (Fcp)f are the heat capacities of the solid and fluid phases, respectively, k is the thermal conductivity, q˘ is the heat production per unit volume, µf is the fluid viscosity, and U B p is the local pore velocity vector. The term Ω refers to the viscous dissipation of the fluid, which is proportional to the square of the fluid velocity gradients. Because of the low Reynolds number nature of electroosmosis in porous media, the viscous dissipation term is neglected in the present study. Taking averages over the REV of the porous medium, we have, for the solid phase,

(1 - φ)(Fcp)s

∂Ts ) (1 - φ)∇‚(ks∇Ts) + (1 - φ)q˘ s ∂t

(3)

and, for the fluid phase,

φ(Fcp)f

∂Tf + (Fcp)f U B D‚∇Tf ) φ∇‚(kf∇Tf) + φq˘ f (4) ∂t

Here, we use the space-averaged velocity, i.e., Darcy velocity vector which is defined as

U BD )

1 Vp

∫V ∫ UB pdVp p

(5)

where Vp is the void space (which is also called pore where the fluid phase can pass through) of the elementary volume, and Ve is the total space of the elementary volume. The porosity is defined as φ ) Vp/Ve ) Vf/V. By setting Ts ) Tf ) T and adding eqs 3 and 4, we have

∂T h D‚∇T ) ∇‚(km∇T) + q˘ m (Fcp)m + (Fcp)fU ∂t

(6)

where we define the thermal capacity, thermal conductivity, and heat generation of the porous medium, respectively, as

(Fcp)m ) (1 - φ)(Fcp)s + φ(Fcp)f

(7a)

km ) (1 - φ)ks + φkf

(7b)

q˘ m ) (1 - φ)q˘ s + φq˘ f

(7c)

with further arrangement and using q˘ s ) 0 (as solid particles are nonconducting), and assuming constant heat capacities, we have

φq˘ f ∂T 1 + U B ‚∇T ) ∇‚(Rm∇T) + ∂t σ D (Fcp)m

(8)

where we introduce the capacity ratio

σ)

(Fcp)m (Fcp)f

and the thermal diffusivity of the porous medium

(9)


Rm )

km

Langmuir, Vol. 21, No. 16, 2005 7601

(Fcp)m

The Joule heat generation in the fluid phase can be obtained using Ohm’s law

I2 q˘ f ) ) λfE2p λf

(11)

where λf is the electric conductivity of the electrolyte which is given by

λf(T) ) λ+ (T)C+ + λ- (T)C-

(12)

Here, λ+ (T) ) λ+0 + 0.025λ+0(T - 298.13) and λ- (T) ) λ-0 + 0.025λ-0(T - 298.13) are the electrical conductivity of the cation and anion of the electrolyte, respectively; C+ and C- respectively denote the mole concentration of cations and anions of the electrolyte. For NaCl solution, λ+0 ) 50.08 × 10-4 m2 S/mol and λ+0 ) 76.31 × 10-4 m2 S/mol.14-16 Considering the tortuosity factor of the porous media, the local electric field strength is given by Ep ) Ez/xτ.40,41 Therefore, the Joule heat generation is expressed as

q˘ f )

E2z λf λf ∂Φ ) τ τ ∂z

( )

2

∂T uz ∂T ur ∂T + + ) ∂t σ ∂z σ ∂r ∂2T 1 ∂ ∂T φ λ ∂Φ Rm 2 + r + r ∂r ∂r τ (Fcp)m ∂z ∂z

[

( )

2

(14)

Subsystem II: Capillary Wall. Generally, the temperature boundary condition at the inner capillary wall is unknown. Because the heat generated by Joule heating in the electrolyte solution is primarily dissipated through the capillary wall to the surrounding environment, a conjugate heat transfer problem has to be solved to simultaneously account for heat transfer in both the solution and the capillary wall. The governing equation for heat conduction in the capillary wall is expressed as

[

( )]

∂2T 1 ∂ ∂T ∂T ) Rw 2 + r ∂t r ∂r ∂r ∂z

(15)

with the thermal diffusivity of the capillary wall defined as

Rw )

kw (Fcp)w

(17a)

|

∂T )0 ∂r r)0

∂T -kw ) h(T - T∞) ∂r r)Rw+d

(17b)

|

∂T )0 ∂z z)L

where h is the convection heat transfer coefficient at the outer capillary wall, Rw is the inner radius of the capillary, and d is the thickness of the capillary wall. 2.2. Continuity and Momentum Equations. Here, we apply the so-called capillary model, also known as the Carman-Kozeny theory,44 in which the porous medium is assumed to be equivalent to a series of parallel tortuous tubules with an effective pore size dpore ) 2Rpore and with the same charge of the particle ζ potential, ζp. According to the definition of the REV, the scale of pore size is smaller than that of REV. The Stokes equations and the Poisson-Boltzmann equation governing the EOF and interstitial EDL potential field inside each imaginary tubule should be solved before one can obtain the macroscopic Darcy velocity by taking an average over the REV, accounting for the porosity φ and tortuosity factor τ of the porous medium. The detailed derivation can be found in the authors’ recently published works.40,41 The developed analytical results will directly be applied here. The equation of continuity of the fluid flow in porous media is44

(13)

Ez ) -∂Φ/∂z is the axial electric field strength. Here Φ is the applied electric field potential. Further if the thermal conductivity km is assumed as constant, eq 8 can be expanded as

( )]

|

T(r,z)|z)0 ) T0

(10)

(16)

where kw is the thermal conductivity of the capillary wall, and (Fcp)w is the heat capacity of the capillary wall. The governing eqs 14-15 are subject to following boundary conditions

∂F + ∇‚(FU h D) ) 0 ∂t

(18)

For an incompressible electrolyte, eq 18 is reduced to

1 ∂(rur) ∂uz + )0 r ∂r ∂z

(19)

In the presence of Joule heating, the electroosmotic flow is strongly coupled the with temperature field through the temperature-dependent viscosity, dielectric constant, and electrical conductivity (local electric field). As a result, the velocity field deviates from its situation of no Joule heating effects because a pressure field will be induced to satisfy the mass conservation, i.e., constant flow rate. Therefore, a model for mixed electroosmotically and pressure-driven flow through porous media is introduced. According to the well-known Darcy’s law that describes the pressure-driven flow in porous media44 and the derivation of the electroosmotic flow inside porous media by Kang et al.,40-41 the Darcy velocity due to electroosmosis and induced pressure can be determined by

U BD )

φ r0 K ζ (1 - G)E B- P B τ µf p µf

(20)

where the electric field strength Ε B ) -∇Φ and the pressure gradient P B ) ∇p. K is the permeability of the porous medium defined as K ) φ3d2p/[a(1 - φ)2] with dp and a ) 150 being the packing particle diameter and the Ergun constant, respectively.45,46 Neglecting the radial gradient of the electric potential field, ∂Φ/∂r, eq 20 can be rewritten as (46) Patankar, S. V. Numerical Heat Transfer and Fluid Flow; McGraw-Hill: New York, 1980.

7602


Kang et al.

K ∂p ur ) µf ∂r

(21a)

φ r0 ∂Φ K ∂p uz ) ζ (1 - G) τ µf p ∂z µf ∂z

(21b)

where r is the dielectric constant of the electrolyte and is considered as a function of temperature, expressed by2 r(T) ) 305.7 exp(-(T/219)). 0 is the permittivity of vacuum. µf is the viscosity of the electrolyte solution, and its variation with the temperature is given by19 µf(T) ) 2.761 × 10-6 exp(1713/T). ζp is the ζ potential at the surface of the packing particles. G is defined as

G)

2 ζpR2pore

∫0R

pore

rψi(r)dr

[

]

[ ]

∇‚(λf∇Φ) ) 0

∂ ∂Φ λ )0 ∂z f ∂z

( )

|

r)0

)0

{

(24)

where e0 is the elementary charge, n0 is the ionic concentration of the electrolyte in the bulk phase, and kB is the Boltzmann constant. The pore size is defined as40,41

Rpore ) dpore/2 )

φ d 3(1 - φ) p

(25)

It should be noted that the solutions of the Darcy velocity distribution in eq 20 are based on the volume-averaging method over a representative elementary volume (REV). As the solutions were derived without considering the temperature change due to the Joule heating effect, discussion of its validity will be provided as following. The length scale of the REV is considerably smaller than the system length scale. According to the definition, the REV is the basic element, accounting for the macroscopic spatial variations such as temperature and Darcy velocity.40 Under the volume-averaging method, the temperature and Darcy velocity within the REV are treated as constant, and so are the temperature-dependent physical properties, such as viscosity. Likewise, the P-B eq 23 is based on the interstitial EDL field inside each imaginary tubule (i.e., the pore size scale), which is even smaller than the length scale of REV. For the same reason, the interstitial temperature variation within each imaginary

Φ|z)0 ) Φ0

Φ|z)L ) 0

(28)

2.4. Normalized Governing Equations. The above derived governing equations can be nondimensionlized as

(0 e r e Rpore) (23)

ψi|r)Rpore ) ζp

(27)

where the electric potential field Φ is subject to the boundary conditions

subject to the following boundary conditions

dψi dr

(26)

Neglecting the radial gradient of the electric potential field, ∂Φ/∂r, eq 26 is reduced to

(22)

From its definition in above equation, G is the normalized average interstitial EDL potential inside each imaginary tubule. Physically, 1-G denotes the EDL-related correction factor to the classical Smoluchowski equation on the EOF velocity. Normally, the EDL thickness ranges from tens to several hundred nanometers. When the EDL is relatively thin compared with the capillary size, G is very small and negligible. However, when the EDL thickness is comparable to the capillary size, such as in the present condition where the pore size is only several microns, the EDL-related correction should be taken into account. According to the capillary model, ψi(r) is the interstitial EDL field inside the imaginary tubules, and it is governed by the Poisson-Boltzmann equation1,2

1 d dψi(r) r ) r dr dr e0ψi(r) 2n0e0 sinh r0 kBT

tubule is neglected. This means the temperature at the cross-section in each tubule is also treated as constant, and so are the physical properties. Therefore, although eq 23 is derived for constant dielectric constant r, the resulting correction factor G in eq 22 and the Darcy velocity in eq 20 are still valid in the present model. 2.3. Continuity of Electric Current. Because the capillary wall is nonconducting, the conservation of the electric current density gives

1 ∂T h ∂T h ∂T h + Peu j z + Peu jr ) Fom ∂th ∂zj ∂rj h ∂2T h 1 ∂ ∂T + rj + Bi Joq h˘ jh 2 j r ∂r j ∂r j ∂zj h 1 ∂T h ∂2T h 1 ∂ ∂T ) + rj Fow ∂th rj ∂rj ∂rj ∂zj2 jz j r) ∂u 1 ∂(rju + )0 rj ∂rj ∂zj ∂p j u j r ) -DaRe ∂rj φ ∂p j h z - DaRe u j z ) (1 - G)E τ ∂zj ∂ (λh E h ))0 ∂zj f z

[

[

( )] ( )]

(29)

The dimensionless parameters are defined as

ht )

t t0

rj )

r Rw

zj )

z Rw

T h)

T - T0 Tb - T0

p j)

p FfU2s

u Us Ez E hz ) (30) E0 u j)

where Rw is the radius of the packed capillary, t0 ) R2w/RmBi is the reference time. T0 is the reference temperature that can be the initial temperature; Tb ) 100 °C is the reference boiling temperature. The reference velocity is the Smoluchowski velocity defined as Us ) r0ζpE0/µf. The dimensionless Joule heat generated in a porous h˘ jh ) φλhfE medium is given by q h 2z /τ. The dimensionless electrical conductivity is defined as λhf ) λf(T)/λf0(T0). The Biot number is defined as Bi ) hRw/km. The Fourier numbers for the packed capillary and the capillary wall are defined as Fom ) Rmt0/R2w and Fow ) Rwt0/R2w, respectively. Thermal Peclet number is defined as Pe ) (FCp)fUsRw/km. The Joule number is defined as Jo ) λf0E20Rw/h(Tb - T0), which characterizes the ratio of the Joule heat to the convection-diffusion heat. It is noted that Jo is proportional to λf0, E20, and Rw, and inversely


Langmuir, Vol. 21, No. 16, 2005 7603

Table 1. Material Properties Used in Numerical Simulation6-10 electrolyte solution (NaCl)

properties density, F [kg m-3] heat capacity, Cp [J kg-1 K-1] thermal conductivity, k [W m-1 K-1] molar conductivity, λ0 [m2 S mol-1] dynamic viscosity, µ [kg m-1 s-1] dielectric constant, r

packing particles (ODS)

capillary wall (fused silica)

coating

1000

2600

2202

1420

4180

750

743

1100

0.61

1.2

1.38

0.15

126.4[1 + 0.025(T - 298)] × 10-4 2.761 exp(1713/T) × 10-6 305.7 exp(-T/219)

proportional to h. The Darcy number is defined as Da ) K/R2w. The Reynolds number is defined as Re ) FfUsRw/µf. 3. Numerical Method Because the physical properties such as electric conductivity λf, dielectric constant r, and viscosity µf are temperature-dependent, the formulated governing equations: the energy eqs 14-15, the mass continuity eq 19, the momentum eq 21, the P-B eq 23, and the electric current continuity eq 27, are strongly coupled, and hence, they must be solved iteratively using a numerical method. An in-house code is developed to solve the proposed mathematical models using a finite-volume-based numerical method.46 In brief, the governing equations are discretized using the control volume integration. The nonlinear source terms are linearized using the Taylor series expansion, and the resulting discrete algebraic equations are solved iteratively using the line-by-line tridiagonal-matrix algorithm (TDMA). For this transient problem, we use iterations to update the equation coefficients and source terms so that the nonlinearity and coupling of the equations are taken into account. The physical properties are evaluated by using the temperature field from the previous iteration (initially at room temperature) to determine the EDL potential field from the P-B equation eq 23, the EDL-related correction factor G from eq 22, the electric field strength from eq 27, and then the velocity and pressure gradient from eqs 19-21. These results are then substituted into the energy equation eqs 14-15 to compute the temperature field for the present iteration. Within each iteration, we solve the discrete algebraic equation system by an iterative lineby-line method, which is a convenient combination of TDMA and the Gauss-Seidel iteration. The computations within every time step are iterated until the temperature change between two consecutive iterations is less than a given tolerance; as such, the computation for the temperature at the current time step is convergent, and can move on to the next time step. Finally, the overall computation stops when the temperature change between two consecutive time steps is less than the given tolerance. In calculation, a grid-independent study was carried out, and the results show that a system of 100 × 100 control volumes in the radial and axial directions produces grid-independent solutions. A relative tolerance of 10-4 for the convergence of iteration within every time step and between time steps was found to be able to greatly reduce the time cost without losing accuracy. The length of time step used in the computation is based on the reference time scale t0 for a specific case study. Because of the low Reynolds number, the diffusion (conduction) effect predominates over the convection effect in the present heat transfer process. Thus, the time scale for the temperature to reach its steady state can be roughly

estimated by letting the product of Biot number and Fourier number of the system equal to 1, i.e., Bi‚Fo ) Bi(Rmt0/R2w) ) 1. Hence, we have t0 ) R2w/(RmBi). A timestep independent study indicates that a time-step length ∆t ≈ 0.007 t0 is appropriate to give time-step independent results satisfactorily. To verify the code, we compared results under the present model with two published works on the Joule heating. One is the modified Knox’s formulas19 on the thermally developed, steady-state temperature excess across both the porous and the capillary wall subsystems. The other is the numerical and experimental study by Keim and Ladisch.42 It is found that the numerical predictions using the present model are in reasonable agreement with the analytical solution by Knox19 and the experimental data reported by Keim and Ladisch.42 The details will be shown in the following section. 4. Results and Discussion In this section, the transient joule heating effects on the electroosmotic flow in a microcapillary packed with microspheres are presented. In simulation, NaCl solution of concentration 10-4-10-3 M is chosen as the fluid phase; the stationary phase is octyldecyl silica (ODS) particles. The packing column is a fused silica capillary with polyimide coating. According to the product specifications (Polymicro Technologies, USA), the ratio of the capillary radius to the wall thickness is chosen as 3:1, and the ratio of the wall thickness to the coating layer thickness is chosen as 3:1. The physical properties of the capillary and particles are listed in Table 1. Because the electrokinetic boundary wall effect is neglected because of a large columnparticle ratio, the porosity and tortuosity factors of the packed capillary are fixed as 0.4 and 1.5, respectively. 4.1. Transient Developing Temperature Field and its Effect on EOF. The temporal development of the temperature field induced by Joule heating is shown in Figure 2. It is demonstrated that there exists a thermal entrance region where the temperature distribution varies not only along the capillary radius but also along the capillary axis. Figure 2a shows the axial temperature distributions at the capillary centerline. The temperature of the fluid starts to increase from the room temperature (25 °C), then gradually increases, and finally reaches a steady temperature at a downstream distance. The length of such spatially developing distance is traditionally called the thermal entrance length.47 The axial temperature gradient is highest at the inlet and decreases to zero at the end of the thermal entrance region. It is noted that the thermal entrance length for the electroosmotic flow in a packed microcapillary extends about 250 times the channel diameter, which is exceptionally longer than the (47) Bejan, A. Convection Heat Transfer; Wiley: New York, 1995.

7604


Kang et al.

Figure 2. Transient development of the temperature field: (a) along the axis, (b) radial distribution at downstream. Working parameters: capillary inner diameter dw ) 530 µm, capillary length L ) 15 cm, packing particle size dp ) 6 µm, porosity φ ) 0.4, tortuosity τ ) 1.5, electrolyte concentration C ) 10-3 M, ζ potential at particle surface ζp ) 50 mV, applied electric field Φ ) 3000 V, and convection heat transfer coefficient at the outer capillary surface h ) 5 W/m2 K. Reference time t0 ) 150 s, Tb ) 100 °C, T0 ) 25 °C.

conventional thermal entrance of a pressure-driven flow in a nonpacked channel (5-10 Dh).47,48 According to the references,47,48 the thermal entrance length is proportional to the thermal Peclet number, i.e., xfd,t/D ≈ 0.05 Pe. However, there is still no literature known to the authors that has studied the thermal entrance length for the EOF in a packed channel. Obviously, the thermal entrance length in this case does not follow this simple relationship. Because the thermal entrance length is significant, it merits a further study to include the thermal entrance effect when designing an electrokinetic chromatography and micropump system. Figure 2b shows the temporal development of radial temperature distributions at a downstream distance longer than the thermal entrance length. It is shown that the radial temperature difference inside the porous packing and the channel wall is insignificant compared with the transient temperature rise. This phenomenon is due to the low Biot number of the whole system (as a rough estimation, Bi ) O(10-3)). Because the Reynolds number of the macroscopic electroosmotic flow is very small (Re ) O(10-1)), the heat generated in the packed region is dissipated mainly by radial conduction through the channel wall and then by convection with the surrounding air rather than by axial advection of the fluid flow. If Bi , 1, the resistance to heat conduction within the capillary system is much less than the resistance to thermal convection across the channel wall-surrounding air boundary layer.48 Hence, virtually, the temperature variation is only the difference between the outer wall and the surrounding air, whereas the temperature within the whole capillary system remains nearly uniform. However, with a close examination of the steady-state radial temperature distribution (th > 4) as in the inset of Figure 2b, the temperature distribution in the porous packing shows a parabolic profile. The highest temperature occurs at the channel centerline. This confirms the findings by Knox,19 who studied the thermal effects on a capillary electrochromatography. A further quantitative comparison between the current results and the Knox’s formulas will be presented in the following section. One of the effects of the temperature rise is that it causes variation of the fluid electric conductivity. The local electric conductivity increases with temperature elevation. Because the total electric current along the capillary is

Figure 3. Transient development of the electric field strength. Working parameters are identical to those in Figure 2.

constant, as determined by the current continuity in eq 27, a higher local electric conductivity leads to a lower local electric field strength. Therefore the electric field strength becomes much higher at the thermal entrance region where the temperature is under development, and it becomes lower and finally becomes constant in the spatially developed region, as shown in Figure 3. According to eq 21, the EOF velocity field is coupled with the temperature field through the temperaturedependent dielectric constant, r, viscosity, µf, and electric field strength, ∂Φ/∂z. The transient development of the axial and radial velocity distributions is shown in Figure 4. Although the viscosity and dielectric constant decrease with increasing temperature (see Table 1), the effect on viscosity is much more significant, thus resulting in an increase in EOF velocity. Further, as shown in Figure 4a, the high local electric field strength greatly increases the local velocity at the thermal entrance region. In Figure 4b, it is noted that the EOF velocity begins to increase gradually over time until it reaches a steady state. The time for the velocity development is the same as that for the temperature field development. Actually, (48) Incropera, F. P.; DeWitt, D. P. Fundamentals of Heat and Mass Transfer; Wiley: New York, 2002.


Langmuir, Vol. 21, No. 16, 2005 7605

Figure 4. Transient development of the electroosmotic velocity field: (a) axial distribution, (b) radial distribution at downstream. Working parameters are identical to those in Figure 2. Reference velocity Us ) 2.03 × 10-4 m/s.

if there is no Joule heating effect, the velocity field will remain uniform in the porous domain, corresponding to the curve of t/t0 ) 0. According to the reference,43 without the presence of the Joule heating effect, it can be inferred that the time scale for the EOF velocity in porous media to reach its steady state after the application of an electric field is t* ) FR2pore/µ ) O(10-6 s). This time scale is very small because of the small interparticulate pore size, which is about 1 µm in this case. However, when the Joule heating effect is present, the time for the velocity to be fully developed becomes almost eight orders higher (O(102s)), depending on the size of the capillary. Also the steadystate velocity is increased to about 60% higher than its initial magnitude (which is also the EOF velocity without Joule heating effects), suggesting significant Joule heating effects. Furthermore, the steady-state velocity distribution in Figure 4b also shows a parabolic profile, with the highest velocity in the channel center. This is because the highest temperature rise occurs at the center, but the radial temperature difference is so small that the radial velocity is virtually uniform in its transient development. In the literature, the distorted velocity profiles in an unpacked capillary under the Joule heating effect were experimentally verified by Xuan et al.18 using a caged-fluorescent dye-based microfluidic visualization technique. It is interesting to note that the present model predicts an induced pressure gradient along the flow direction, as shown in Figure 5. This phenomenon is not due to the effect of packing microspheres in the packed capillary, but due to the Joule heating-induced temperature gradients and their effect on electrical/transport properties. The presence of the induced pressure gradient inside the capillary is to maintain the fluid mass conservation, i.e., constant flow rate. This is due to the fact that, for scenarios where the electrical and transport properties vary with temperature and where Joule heating effects are simulated, there is local imbalance between the electroosmotic body force and the shear force, which results in the induced pressure gradient. Similar to the findings reported by Xuan et al.17 on the Joule heating in nonpacked capillaries, the induced pressure gradient is positive at the entrance region, whereas it becomes negative at the spatially developed region. However, the velocity profile at the entrance in this study is in the convex shape (lower close to the wall and higher near the center), which is contrary to the concave (higher close to the wall and lower near the center) profile in a nonpacked capillary reported by Xuan

Figure 5. Transient development of the induced pressure distribution. Working parameters are identical to those in Figure 2.

et al.17 This discrepancy is due to the different flow characteristics in packed and nonpacked capillaries. In nonpacked capillaries, the sole electroosmotic flow profile for thin EDL thickness is pluglike without consideration of Joule heating effects. In the presence of the Joule heating effect, the combination of electroosmosis and the induced positive-pressure gradient (back pressure) at the entrance can generate a concave velocity profile in case of a nonpacked capillary. However, in packed capillaries, according to Darcy’s law, the pressure-induced (macroscopic) flow profile is also pluglike. Thus the induced back pressure only reduces the magnitude of the Darcy velocity uniformly in a transverse direction. However, the net velocity profile is slightly concave because of the presence of a Joule heating-induced parabolic temperature profile as shown in Figure 2b. Nonetheless, the prediction of the presence of an induced pressure gradient due to the Joule heating effect suggests that the velocity profile apparently exhibits a mixed characteristic with a combination of electroosmotic flow and pressure-driven flow in a packed microcapillary. 4.2. Comparison with Published Works. To verify the numerical model presented in this work, case studies compared with Knox’s formulas19 and the numerical and experimental work by Keim and Ladisch42 will be shown in the following.

7606


Kang et al.

Figure 6. Comparison of the present model predictions with the published work by Keim and Ladisch:42 (a) experiment data of transient temperature rise at the center of the column outlet, (b) modeling of the fully developed radial temperature profile at column outlet. Working parameters: capillary inner diameter dw ) 3.81 cm, capillary length L ) 38.1 cm, packing particle size dp ) 5 µm, porosity φ ) 0.36, tortuosity τ ) 1.5, electrolyte concentration C ) 3.9 × 10-3 M, ζ potential at particle surface ζp ) 60 mV, applied electric voltage Φ ) 419 V, and convection heat transfer coefficient at the outer capillary surface h ) 5 W/m2 K. Reference time t0 ) 2.7 h, Tb ) 100 °C, T0 ) 25 °C.

Knox19 presented simple formulas to estimate the spatially developed, steady-state temperature difference between the channel axis and its inner wall, θcore ) E2λCφD2c /16km, and temperature excess across the channel wall, θwall ) E2λCφD2c /8kw ln(Do/Dc). With simple modification, including the tortuosity factor of the porous packing, Knox’s formulas can give θcore ) (φ/τ)(E2λCD2c /16km) ) 0.046 °C and θwall ) (φ/τ)(E2λCD2c / 8kw) ln(Do/Dc) ) 0.057 °C. The numerical model in this study as shown in the inset of Figure 2b gives quantitatively similar predictions as those by Knox’s formulas. Keim and Ladisch42 developed a numerical model to investigate the temperature profiles in a large diameter (around 3 cm) electrochromatography. In addition, they measured the transient temperature elevations at the center of the column outlet. A case study was performed using the present numerical model to simulate the temperature changes under the same conditions as reported in Keim and Ladisch’s work. The working parameters are: capillary inner diameter dw ) 3.81 cm, capillary length L ) 38.1 cm, packing particle size dp ) 5 µm, porosity φ ) 0.36, tortuosity τ ) 1.5, electrolyte concentration C ) 3.9 × 10-3 M, ζ potential at particle surface ζp ) 60 mV, applied electric voltage Φ ) 419 V, and convection heat transfer coefficient at the outer capillary surface h ) 5 W/m2 K. In Figure 6a, the simulation results obtained from the present model are compared with the experimental data of Keim and Ladisch.42 When the electric field is switched on, the temperature rise in their experiment is slightly higher than that predicted by the present model. Whereas, during the cooling process, namely, when the electric field is turned off after 3 h and 35 min, the temperature drop in the experiment is slightly faster than that predicted by this model. This scenario is due to the special condition applied in their experiment, i.e., the flow rate was kept constant during the course of their experiment, although they did not mention how to realize it. However, in this study, the electroosmotic flow velocity is strongly coupled with the temperature field, i.e., the flow rate through the packed microcapillary may vary over time. As discussed in the previous section, during the heating session, the high axial temperature gradient at the entrance greatly

enhances the local electric field strength resulting in an increased flow field. Therefore, the increased axial flow convection has a better heat dissipation effect, which keeps the temperature rise at a relatively low level. During the cooling session, there is no driving force for the fluid flow because of the absence of electric field. Thus, the flow velocity will decrease until the fluid ends with a full stop. In this situation, the heat dissipation is mainly dependent on the heat conduction in radial direction. No wonder that the temperature drop is faster in their experiment in which the flow rate was not affected. Figure 6b shows a comparison of the radial temperature profile predicted by Keim and Ladisch’s numerical model and the model developed in this study. It can be seen that the two models give close results. The major discrepancy lies at the boundary between the porous packing (subsystem I) and the capillary wall (subsystem II), where their model presented a discontinuous temperature gradient. This is because, in their model, an empirical boundary condition between the two subsystems was imposed. This only enforces the continuity of the temperature itself, not the temperature gradient. In contrast, in the present study, the conjugated heat transfer equations within the two subsystems are solved simultaneously, and hence, the problem of the temperature gradient discontinuity is resolved. 4.3. Effect of the Joule Number. From the source term in the nondimensional energy equation in eq 29, it can be inferred that the Joule heating is dependent on the dimensionless Joule number, Jo. According to its definition, Jo ) λf0E20Rw/h(Tb - T0), Jo is proportional to λf0, E20, and Rw, and inversely proportional to h. A parametric study of the effect of Joule number on the steady-state axial temperature distribution is shown in Figure 7. It can be seen that the temperature elevation increases almost proportionally with the Joule number. For a nominal value of Jo ) 4.72, the solution temperature may increase to the boiling temperature. As aforementioned, such a high operation temperature has deleterious effects on the applications of the electroosmotic flow. Therefore, it is very important to control the Joule number to be as low as possible by controlling the working parameters,


Figure 7. Effect of the Joule number on the axial temperature distribution. Working parameters: capillary length L ) 15 cm, packing particle size dp ) 6 µm, porosity φ ) 0.4, tortuosity τ ) 1.5, Tb ) 100 °C, and T0 ) 25 °C.

such as the electric field strength, electrolyte concentration, channel size, and heat exchange at the outer wall. 5. Conclusions A numerical model for evaluating the transient Joule heating induced temperature field and its effect on the electroosmotic flow in a capillary packed with microspheres is developed. The model incorporates the coupled momentum equation, the energy equations, and the mass and electric current continuity equations. The coupled governing equations are solved numerically using the finite difference method. The following conclusions can be drawn from this study:

Langmuir, Vol. 21, No. 16, 2005 7607

(1) In the presence of Joule heating, there exists a significant axial temperature gradient in the thermal entrance region. This high temperature gradient strongly enhances the local electric field at the entrance, resulting in a nonuniform distributed electric field along the flow direction. (2) Compared to the conventional thermal entrance region in unpacked channels, the thermal entrance region in a packed microchannel is two orders longer. (3) Though the radial temperature profile virtually exhibits flat because of the small system Biot number, the presence of Joule heating significantly elevates the magnitude of temperature, which in turn greatly increases the EOF velocity by means of changing the local fluid viscosity, dielectric constant, and local electric field strength. (4) In the presence of the Joule heating effect, the electroosmotic velocity field is synchronized with the Joule heating induced temperature field. Thus, the duration of transient velocity development is several orders longer than that in the absence of Joule heating effects. (5) It is found that a pressure gradient is induced along the capillary. The velocity profile in fact exhibits a combined electroosmotically and pressure-driven flow in a packed microcapillary. (6) The numerical predictions using the present model are in reasonable agreement with the analytical solution by Knox19 and the experimental data reported by Keim and Ladisch.42 Acknowledgment. The financial support of the Academic Research Fund to C.Y. from Ministry of Education of Singapore (contract number RG07/02) and the Ph.D. scholarship to Y.J.K. from Nanyang Technological University are gratefully acknowledged. LA050061G

Joule Heating Induced Transient Temperature Field and Its Effects on

Recommend Documents