Unsteady electroosmotic flow of power-law fluid in a rectangular microchannel

Mechanics Research Communications 39 (2012) 9–14

Contents lists available at SciVerse ScienceDirect

Mechanics Research Communications journal homepage: www.elsevier.com/locate/mechrescom

Unsteady electroosmotic flow of power-law fluid in a rectangular microchannel S.Y. Deng a , Y.J. Jian a,∗ , Y.H. Bi b , L. Chang b , H.J. Wang c , Q.S. Liu a a

School of Mathematical Science, Inner Mongolia University, Hohhot 010021, China School of Statistics and Mathematics, Inner Mongolia Finance and Economics College, Huhhot 010051, China c Mathematical Science College, Inner Mongolia Normal University, Hohhot 010022, China b

a r t i c l e

i n f o

Article history: Received 26 August 2011 Received in revised form 19 September 2011 Available online 29 September 2011 Keywords: Electroosmotic flow Power-law fluid Flow behavior index Electric double layer thickness Finite difference method

a b s t r a c t This paper studies the unsteady electroosmotic flow of the power-law fluid in a rectangular microchannal. The electric potential distribution is given by linearized Poisson–Boltzmann equation. By using finite difference method, we solved Cauchy momentum equation and obtained the velocity profile of the power-law fluid. The influences of the flow behavior index and the ratio of the width of the rectangular microchannel to the double layer thickness on the velocity profiles are investigated. In addition, the evolution of the velocity profiles with the time is analyzed. © 2011 Elsevier Ltd. All rights reserved.

1. Introduction Microfluidic devices have been significantly applied in microelectric mechanic systems and bio-sensor areas, such as lab-ona-chip (Stone et al., 2004; Bayraktar and Pidugu, 2006). Since the electrolyte and microchannel wall contact with each other, it leads to the exchange of charges in the channel wall and electrolyte. It depends on the chemical composition of microchannel and the surface chemical reaction process of the electrolyte. Because surface charge of the fluid affects the distribution of anions and cations, which results in the formation of electric double layer (EDL) (Hunter, 1981). Positively charged ions quickly gather near the wall, forming the stern layer whose thickness is only one diameter of ions. Close to the stern layer, there is the diffuse layer which contains cations and anions, and the distribution of the ion density obeys the Boltzmann distribution (Karniadakis et al., 2005). For the applied electric field along the charged surface tangent direction, the free ions of the EDL migrate due to the effect of Coulomb force. Because of the viscosity of fluid, the free charged ions will lead to the whole fluid movement, resulting in an electroosmotic flow. Now the EOF is widely used in biological, chemical and pharmaceutical fields. At the macro scale, it is the mechanical pump often used to generate the pressure and then using the pressure to drive fluid in the length scale of the channel flow. However, in microchannels, in particular, when the length scale of channel reduced to nanometer

∗ Corresponding author. Tel.: +86 471 4992946x8313; fax: +86 471 4991650. E-mail address: [email protected] (Y.J. Jian). 0093-6413/$ – see front matter © 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.mechrescom.2011.09.003

scale, the model of the pressure-driven flow will become very difficult due to the reduction of the diameter of channel. In a variety of geometric shapes of microchannel, such as parallel plate channels (Burgreen and Nakache, 1964), cylindrical-shaped channel (Levine et al., 1975), circular channel (Tsao, 2000; Kang et al., 2002), ellipse-shaped channel (Hsu et al., 2002), rectangular channel (Yang et al., 1998; Arulanandam and Li, 2000; Marcos et al., 2004), T-shaped channel (Bianchi et al., 2000) and semi-circular channel (Wang et al., 2008), fully developed Newtonian fluid in the theoretical and experimental aspects of EOF has been greatly studied. All of the above mentioned studies are about Newtonian fluids. However, microfluidic devices are often used to analyze biological fluids, and these fluids are often long-chain molecules in solution. These long-chain molecules made the fluids showing some of the fluid nature of non-Newtonian fluid, for example, changes in viscosity, memory effect, normal stress effects, yield stress and hysteresis fluid properties. In order to describe non-Newtonian fluid flow characteristics inside the microchannel, the constitutive relation should be given to represent the connection between kinetic viscosity and the shear stress. The study of the EOF of non-Newtonian fluids is few, mainly because the complexity of the constitutive relation of fluid itself brings difficulties. Therefore, research in this area mainly limited to simple non-Newtonian fluid model. Das and Chakraborty (2006) and Chakraborty (2007) first studied the steady EOF of the inelastic power-law fluid. The same fluid model is used to the study of the steady EOF in slit microchannel by Zhao et al. (2008). The analytical expression of the velocity profile depending on the flow behavior index n has been given. Vasu

10

S.Y. Deng et al. / Mechanics Research Communications 39 (2012) 9–14

is the Boltzmann constant, and T is the absolute temperature. The local net charge density can be expressed as e = e





ni∞ zi exp

−



zi e kb T



= −2en∞ sinh

e kb T



,

(3)

where the electrolyte is assumed to be z+ :z− = 1:1, and n+∞ = n−∞ = n∞ , the ionic number concentration in the bulk solution where (x, y) = 0. Substituting Eq. (3) into Eq. (1), we can obtain two-dimensional nonlinear Poisson–Boltzmann equation (Zhang et al., 2007) ∂2 ∂2 2en∞ + = sinh 2 εε0 ∂x ∂y2

Fig. 1. Schematic of the rectangular microchannel and the coordinate system used for modeling.

and De (2010) studied the EOF of power-law fluid with high zeta potential in rectangular microchannel. Zhao and Yang (2010) carried out Smoluchowski velocity of the EOF of power-law fluid in the microchannel with high zeta potential. Tang et al. (2009) using lattice Boltzmann method simulated the EOF of the power-law fluid in a microchannel. However, no one seems to have discussed, to the authors’ knowledge, the unsteady EOF of power-law fluid in two-dimensional rectangular microchannel. A mathematical model is developed to describe the EOF in a rectangular microchannel. The model includes solving the Poisson–Boltzmann equation for EDL potential distribution and the modified Cauchy momentum equation for EOF field. Based on the finite difference method, we obtained the velocity profile for various flow behavior index n and the value of the ratio of channel width to Debye length. The theoretical results of Newtonian fluid are compared with those of numerical computation by the present finite difference method, and good agreement is obtained. Furthermore, we compute the velocity profiles of powerlaw fluid with different flow behavior index n.

(4)

(5) 1/2

where k = (2e2 n∞ /εε0 kb T ) , 1/k is the EDL thickness. Right next to the stern layer, there is the diffuse layer containing cations and anions. Stern and diffuse layer form the EDL. Due to the symmetry of the potential and velocity fields, the solution domain can be reduced to a quarter cross section of the channel. The boundary conditions are x=0:

∂ = 0, ∂x

∂ = 0, ∂y

y=0:

x=H:

= ,

= .

(x, y) = 4

∞  (−1)n+1 cosh

(6)

·cos

2n − 1 x 2H

 ,

(2)

where ni∞ is the ionic number concentration at the neutral state where (x, y) = 0 and ni is the ionic number concentration of the ith ionic species at the state where the electric potential is (x, y), zi is the valence of ith ionic species, e is the elementary charge, kb

∞  (−1)m+1 cosh m=1

·cos



(((2n − 1)2 2 + 4k2 H 2 )/4H 2 )y



(2n − 1) cosh

n=1

(1)

where ε is the relative permittivity, ε0 is the permittivity in vacuum, and e is the local net charge density. It is assumed that the ionic concentration per volume in electrolyte solution obeys the Boltzmann distribution ni = ni∞ exp

.

∂2 ∂2 + = k2 , ∂x2 ∂y2

+ 4

∂2 ∂2 e + =− , εε0 ∂x2 ∂y2

ze − i kb T



Using the separation of variables method to Eq. (5) and the boundary conditions (6), the solution of (5) is given by (Yang et al., 1998)

The unsteady EOF of power-law fluid through a twodimensional rectangular microchannel with the z-axis being in the axial velocity direction is sketched in Fig. 1. The Cartesian axes are placed at the middle of the channel. According to the electrostatic theory, the connection between the electric potential (x, y) of the EDL and the local net charge density e can be expressed as Poisson equation



e kb T

The inner boundary of the EDL with one hydrated ion radius is denoted as the Stern layer. The ions attached to the surface are located within the Stern layer and they are considered to be immobile. Ions whose centers are located beyond the Stern layer form the diffuse mobile part of the EDL, which is located between one to two radii away from the surface. This boundary is referred to as the shear plane. It is on this plane where the no-slip fluid flow boundary condition is assumed to apply. The potential at the shear plane is commonly known as the zeta potential. Assuming a small zeta potential, which means physically that the electrical potential is small compared with the thermal energy of the charged species, the electric potential due to charged wall can be described by linearizing hyperbolic sine function sinh(x) with x when x is small

y=H: 2. Model description



((2n − 1)2 2 + 4k2 H 2 )/4



(((2m − 1)2 2 + 4k2 H 2 )/4H 2 )x



(2m − 1) cosh

((2m − 1)2 2 + 4k2 H 2 )/4

2m − 1 y. 2H

(7)

The vector formation EOF of the Cauchy momentum equation can be expressed as



m

∂V + V · ∇ V ∂t



 e, = −∇ p + ∇ ·  + E

(8)

where m is the mass density, V is EOF velocity, t is the time, p is the pressure, E is the strength of electric field,  = 2 ˙ n−1 e is the stress tensor defined by with ˙ = (2ekl ekl )1/2 being the second invariant of the strain rate tensor eij = (1/2)((∂ui /∂xj ) + (∂uj /∂xi )),  is the dynamic viscosity coefficient of dimension [Nm−2 sn ]. For a single direction EOF in a rectangular microchannel, if the axial velocity

S.Y. Deng et al. / Mechanics Research Communications 39 (2012) 9–14

11

is w = w(x, y, t), then the constitutive equation of power-law fluid in axial direction is given as follow:



∂w = − ∂x

n−1



∂w ∂w + − ∂x ∂y

n−1

∂w , ∂y

(9)

where  is the shear stress and n is the flow behavior index. The liquid in the microchannel is assumed to be an incompressible fluid, and the flow is only along the z direction. Thus there is no pressure gradient along the microchannel. Substituting Eq. (9) into Eq. (8), it is reduced only to axial direction momentum equation



∂ ∂w = (−1)n−1

∂t ∂x



∂w ∂x

n

∂ + ∂y



∂w ∂y

n 

+

1 Ee , m

(10)

n = 1 represents a Newtonian fluid with constant kinematic coefficient of viscosity = /m . While n < 1 and n > 1 correspond to the case of pseudoplastic (shear-thinning) and dilatant (shearthickening) fluids, respectively. 3. Numerical algorithm Firstly, from Eq. (10), when the flow behavior index n = 1, the fluid is Newtonian fluid, so an analytical solution will be obtained, and Eq. (10) becomes ∂w =

∂t



∂2 w ∂2 w + 2 ∂x ∂y2



+

1 Ee , m

Fig. 2. The comparison of analytical and numerical solution for n = 1 (kH = 5 and y/H = 0.75).

We replace the differential Eq. (14) by a discretized version of it, using the forward difference approximation in terms of time and the central difference approximation in the space direction both first-order derivative and second-order derivative i+1 i wj,k − wj,k ∂w = , t ∂t

(11)

i i wj,k+1 − wj,k−1 ∂w = , 2 H ∂y

and the initial condition and symmetry, no-slip boundary conditions are ∂w = 0, ∂x

x=0:

y=H:

y=0:

w = 0,

∂w = 0, ∂y

t=0:

x=H:

w(x, y, t) =

n=1 m=1

Knm

where Cnm

E = −8k εε0 Hm

(1 − exp(−Knm t)) · cos ˛n x cos ˇm y,



+



i+1 wj,k



(2m − 1)˛n (1 + (2m /˛2n ))

,

˛n =

m =

 n =

×

∂2 w ∂x2



∂w ∂x

n−1

∂2 w + ∂y2



∂w ∂y

n−1

+

and initial and boundary conditions are ∂w = 0, ∂x

y=H:

w = 0,

y=0: t=0:

∂w = 0, ∂y

x = H,

w(x, y) = g(x, y).

t

n

i i wj,k+1 − wj,k−1

n−1

2 H

i i + wi wj,k+1 − 2wj,k j,k−1



H 2

i i = w2,k , w1,k

k2 εε0 E t m

j,k .

(18)

i i wj,1 = wj,2 ,

i wN,k = 0,

i wj,N = 0,

(19)

where gj,k is initial condition, i = 1, 2, . . ., Nt , j, k = 1, 2, . . ., N. Consider the Eq. (18), when the flow behavior index n < 1, the denominator of the first-order derivative in the space direction is equal to zero, therefore, it is assumed that the initial condition is g(x, y) =

(15)

t +

The corresponding boundary conditions and the initial condition are

1 wj,k = gj,k ,

E e , m (14)

x=0:

n−1

n−1

2 H

H 2



(2n − 1)2 2 + 4k2 H 2 . 4H 2

n

i i wj+1,k − wj−1,k

i i + wi wj+1,k − 2wj,k j−1,k

+ (−1)

(2m − 1)2 2 + 4k2 H 2 , 4H 2



n−1



When the flow behavior index n = / 1, the fluid is non-Newtonian. The equation is ∂w = (−1)n−1

∂t

+ (−1)

×

2m − 1 2 = ), , Knm = (˛2n + ˇm 2H



=

i wj,k



2n − 1 , 2H

(17)

where t is the step for t direction and H is the step for x direction and y direction. Finally the corresponding differential equation is

2 )) (2n − 1)ˇm (1 + (2n /ˇm

(−1)n+m

ˇm

(13)

(−1)n+m

2

i i + wi wj,k+1 − 2wj,k ∂2 w j,k−1 = , 2 ∂y H 2

(12)

Using the separation of variable method (Marcos et al., 2004), the solution of Eq. (11) can be written as ∞ ∞   Cnm

(16)

i + wi i wj+1,k − 2wj,k ∂2 w j−1,k = , ∂x2 H 2

w = 0,

w(x, y) = 0.

i i − wj−1,k wj+1,k ∂w = , 2 H ∂x

∞ ∞   Cnm n=1 m=1

Knm

cos ˛n x cos ˇm y.

(20)

12

S.Y. Deng et al. / Mechanics Research Communications 39 (2012) 9–14

Fig. 3. The changes in velocity profile for the different time, when the value of flow behavior index n = 0.8 (kH = 30 and 2H = 200 ␮m): (a) t = 1 × 10−4 s; (b) t = 0.001 s; (c) t = 0.003 s; (d) t = 0.005 s; (e) t = 0.008 s.

When n > 1, there is no above situation. The initial condition is g(x, y) = 0.

(21)

4. Results and discussion The parametric values are taken as following. The strength of electric field is E = 10,000 V/m, and the reference velocity is chosen as W = 0.001 m/s, and the concentration of electrolyte solution n∞ = NA c, c is the molar concentration and c = 10−5 mol/l, NA = 6.02 × 1023 /mol, e = 1.6 × 10−19 C, kb = 1.38 × 10−23 J/K, m = 1000 kg/m3 , = 9 × 10−7 m2 s(n−2) , ε = 80, ε0 = 8.85 × 10−12 C2 /J m, T = 293 K,  = −0.025 V. When the power-law index n = 1, kH = 5, y/H = 0.75, we compared the numerical solution of the velocity profile given by the differential equation (18) with the analytical solution (13) for long time evolution in Fig. 2. It can be seen from Fig. 2 that the numerical result is agreed well with the analytical result. Therefore, the above numerical method is feasible. Using the same differential method, we computed the velocity profile of power-law fluid when n = / 1. Figs. 3 and 4 show the changes of velocity profile at different time when flow behavior index n = 0.8 and 1.2, respectively. The

initial values of Figs. 3 and 4 are given by Eqs. (20) and (21). From Figs. 3 and 4, when the time is very small, the velocity in the middle region of microchannel is at rest, and the velocity in the region close to the channel wall is not equal to zero. With the evolution of the time, the velocity far away from EDL increases gradually due to the application of the electric field. This scenario reveals the unique feature of the EOF. As time elapse, the liquid inside the EDL region exist hydrodynamic shear stress on its adjacent liquid, and the liquid outside the EDL region is set in motion layer by layer, finally extending to the entire channel and the velocity reaches its maximum value and tends to steady state. The shapes of velocity profiles in Figs. 3 and 4 are plug-like (n < 1) and bowl-like (n > 1) respectively. Plug-like velocity profiles refer to the shape of the velocity section likes a plug, while bowl-like velocity profiles means the velocity profiles like a bowl. It still can be seen that the time reaching steady status is very short. Fig. 5 described the velocity profile at different time in the section of y/H = 0.65. Fig. 6 show the comparison of the velocity profile for various flow behavior index n when t = 0.006 s. At this time and after, the velocity profile reached the steady status. From Fig. 6, the velocity profiles depend greatly on the fluid behavior index. For pseudoelastic fluids, namely n < 1, the dimensionless velocity increases as the

S.Y. Deng et al. / Mechanics Research Communications 39 (2012) 9–14

13

Fig. 4. When the flow behavior n = 1.2, the changes velocity profile for different time (kH = 3 and 2H = 200 ␮m): (a) t = 1 × 10−4 s; (b) t = 0.001 s; (c) t = 0.002 s; (d) t = 0.003 s.

Fig. 5. The comparison of the velocity profile for different time (y/H = 0.65): (a) n = 0.8 and (b) n = 1.2.

Fig. 6. The comparison of the velocity profile for various value of flow behavior index n (t = 0.006 s, 2H = 200 ␮m and kH = 10).

value of n decreases and the velocity profile becomes more pluglike. It is explained that the shear stress fast increases as n decreases. For the dilatant fluid, namely n > 1, the velocity profile changes from parabolic to plug-like shapes as n approaches to 1. This is because the viscosity increases gradually as the value of n increases. Fig. 7 described the variation of the velocity profile with kH for different flow behavior index n [(a) n = 0.8; (b) n = 1.1; (c) n = 1]. It can be seen from Fig. 7 that regardless of the power-law index n, the velocity profiles become more plug-like as the value of kH transits from 10 to 50 when keeping other parameters unchanged. Increasing the value of kH indicates that the ratio of the length scale H to the double layer thickness increases. Namely, the density of free charged ions decrease. The free ions only exist in the region close to the channel wall, so that the shear stress changes slowly almost in the entire region of the channel except for the channel wall where the shear stress rapidly jumps to a certain magnitude and the changes in the magnitude of dimensionless velocity is rapidly.

14

S.Y. Deng et al. / Mechanics Research Communications 39 (2012) 9–14

Fig. 7. The comparison of the velocity profile for various width ratio K = kH under the situation of various flow behavior index: (a) n = 0.8; (b) n = 1.1; (c) n = 1.

5. Summary and conclusions This paper numerically studied the EOF of power-law fluid in a rectangular microchannel. The electric potential distribution in the electrolyte and the flow field is given by linearized Poisson–Boltzmann equation. Using the forward difference approximation in terms of time and the central difference approximation in the space direction both first-order derivative and second-order derivative, we solved the Cauchy momentum equation by the finite difference method. The numerical velocity for power-law fluid at different times is obtained. It is found that the flow behavior index n plays an important role and the value of kH is capable of changing the EOF pattern significantly. The shapes of velocity profiles are plug-like (n < 1) and bowel-like (n > 1) respectively. Regardless of the power-law index n, the velocity profiles become more plug-like shapes as the value of kH transits from 10 to 50 when keeping other parameters unchanged. In addition, we found the time reaching steady status is very short. Acknowledgements The work was supported by the National Natural Science Foundation of China (No. 11062005), the support of Knowledge Innovation Programs of the Chinese Academy of Sciences under Contract Nos. KZCX2-YW-201 and KZCX1-YW-12, the Natural Science Foundation of Inner Mongolia (Grant No. 2010BS0107), the research start up fund for excellent talents at Inner Mongolia University (Grant No. Z20080211) and the support of Natural Science Key Fund of Inner Mongolia (Grant No. 2009ZD01). References Arulanandam, S., Li, D., 2000. Liquid transport in rectangular microchannels by electroosmotic pumping. Colloids Surf. A 161, 89–102. Bayraktar, T., Pidugu, S.B., 2006. Characterization of liquid flows in microfluidic systems. Int. J. Heat Mass Transfer 49, 815–824.

Bianchi, F., Ferrigno, R., Girault, H.H., 2000. Finite element simulation of an electroosmotic driven flow division at a t-junction of microscale dimensions. Anal. Chem. 72, 1987–1993. Burgreen, D., Nakache, F.R., 1964. Electrokinetic flow in ultrafine capillary slits. J. Phys. Chem. 68, 1084–1091. Chakraborty, S., 2007. Electroosmotically driven capillary transport of typical non-Newtonian biofluids in rectangular microchannels. Anal. Chim. Acta 605, 175–184. Das, S., Chakraborty, S., 2006. Analytical solutions for velocity, temperature and concentration distribution in electroosmotic microchannel flows in a nonNewtonian bio-fluid. Anal. Chim. Acta 559, 15–24. Hsu, J.P., Kao, C.Y., Tseng, S.J., Chen, C.J., 2002. Electrokinetic flow through an elliptical microchannel: effects of aspect ratio and electrical boundary conditions. J. Colloid Interface Sci. 248, 176–184. Hunter, R.J., 1981. Zeta Potential in Colloid Science. Academic Press, New York. Kang, Y.J., Yang, C., Huang, X.Y., 2002. Electroosmotic flow in a capillary annulus with high zeta potentials. J. Colloid Interface Sci. 253, 285–294. Karniadakis, G., Beskok, A., Aluru, N., 2005. Microflows and Nanoflows: Fundamentals and Simulation. Springer Press, New York. Levine, S., Marriott, J.R., Neale, G., Epstein, N., 1975. Theory of electrokinetic flow in fine cylindrical capillaries at high zeta potentials. J. Colloid Interface Sci. 52, 136–149. Marcos, Yang, C., Teck, N.W., Kim, T.O., 2004. Dynamic aspects of electroosmotic flow in rectangular microchannels. Int. J. Eng. Sci. 42, 1459–1481. Stone, H.A., Stroock, A.D., Ajdari, A., 2004. Engineering flows in small devices: microfluidics toward a lab-on-a-chip. Annu. Rev. Fluid Mech. 36, 381–411. Tang, G.H., Li, X.F., He, Y.L., Tao, W.Q., 2009. Electroosmotic flow of non-Newtonian fluid in microchannels. J. Non-Newtonian Fluid Mech. 157, 133–137. Tsao, H.K., 2000. Electroosmotic flow through an annulus. J. Colloid Interface Sci. 225, 247–250. Vasu, N., De, S., 2010. Electroosmotic flow of power-law fluids at high zeta potentials. Colloids Surf. A: Physicochem. Eng. Aspects 368, 44–52. Wang, C.Y., Liu, Y.H., Chang, C.C., 2008. Analytical solution of electro-osmotic flow in a semicircular microchannel. Phys. Fluids 20, 063105. Yang, C., Li, D., Masliyah, J.H., 1998. Modeling forced liquid convection in rectangular microchannels with electrokinetic effects. Int. J. Heat Mass Transfer 41, 4229–4249. Zhang, P., Zuo, C.C., Zhou, D.Y., 2007. Numerical simulation on influence factors of electroosmotic flow in rectangular microchannel. J. Chin. J. Mech. Eng. 43 (3), 38–40. Zhao, C., Yang, C., 2010. Nonlinear Smoluchowski velocity for electroosmosis of power-law fluids over a surface with arbitrary zeta potentials. Electrophoresis 31, 973–979. Zhao, C., Zholkovskij, E., Masliyah, J.H., Yang, C., 2008. Analysis of electroosmotic flow of power-law fluids in a slit microchannel. J. Colloid Interface Sci. 326, 503–510.