- Email: [email protected]
Exact solutions for electro-osmotic flow of viscoelastic fluids in rectangular micro-channels
Applied Mathematics and Computation 211 (2009) 502–509
Contents lists available at ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
Exact solutions for electro-osmotic flow of viscoelastic fluids in rectangular micro-channels Cunlu Zhao *, Chun Yang School of Mechanical and Aerospace Engineering, Nanyang Technological University, 50 Nanyang Avenue, 639798, Singapore
a r t i c l e
i n f o
Keywords: Electro-osmotic flow Lab-on-a-chip Viscoelastic fluids Generalized Oldroyd-B model
a b s t r a c t Transient electro-osmotic flow of viscoelastic fluids in rectangular micro-channels is investigated. The general twofold series solution for the velocity distribution of electro-osmotic flow of viscoelastic fluids with generalized fractional Oldroyd-B constitutive model is obtained by using finite Fourier and Laplace transforms. Under three limiting cases, the generalized Oldroyd-B model simplifies to Newtonian model, fractional Maxwell model and generalized second grade model, where all the explicit exact solutions for velocity distribution are found through the discrete Laplace transform of the sequential fractional derivatives. These exact solutions may be able to predict the flow behavior of viscoelastic biological fluids in BioMEMS and Lab-on-a-chip devices and thus could benefit the design of these devices. Ó 2009 Elsevier Inc. All rights reserved.
1. Introduction Accompanied by the recently fast and thriving development of BioMEMS and Lab-on-a-chip technologies, renewed interests are also being brought up in electro-osmosis due to the fact that electro-osmosis has been demonstrated as a powerful activation mechanism to manipulate liquid samples of nanovolumes in capillary network of these devices. In reviewing the general concept of electro-osmosis, liquids at the solid–liquid interface ionize because of the zeta potential at the wall. The electric double layer (EDL) is a very thin region close to the charged surface in which there is an excess of counterions over coions to neutralize the surface charge. Therefore, a net fluid flow in the desired direction is induced when imposing external electric fields [1]. Compared to conventional mechanical pumping modes, electro-osmosis-based pumping has numerous advantages, including ease of fabrication and control, no movable mechanical parts and high reliability, no noise, etc. More favorably, the fluid moves as a plug in the electro-osmosis, rather than with the parabolic-flow profile observed when pumping is accomplished by applying pressure gradient to the fluid. This feature helps to minimize the broadening of plugs of sample that occurs with many pressure-driven systems, and allows very high-resolution separations of ionic species. It is a key contributor to electrophoretic separations of DNA in micro-channels [2,3]. Extensive studies of electro-osmotic flow of Newtonian fluids in micro-channels have been reported in the literature theoretically [4–7], experimentally [8–11] and numerically [12–15]. However, forenamed studies are just concerned with Newtonian fluids. Actually, Lab-on-a-chip devices are frequently used to analyze biofluids which usually exhibit viscoelastic behavior. Thus, the more general Cauchy momentum equation, instead of the Navier–Stokes equation should be used to describe the flow characteristics of viscoelastic fluids provided that proper constitutive equations are available. Recently, fractional calculus has gained considerable interest due to its applications in different areas of physics and engineering, including complex fluid dynamics. The constitutive equations with fractional derivative have been proved to be valuable * Corresponding author. E-mail address: [email protected] (C. Zhao). 0096-3003/$ - see front matter Ó 2009 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2009.01.068
503
C. Zhao, C. Yang / Applied Mathematics and Computation 211 (2009) 502–509
tools to handle viscoelastic properties. In general, these equations are derived from known models via substituting time ordinary derivatives of stress and strain by derivatives of fractional order from which one can define the non-integer order integral and derivatives [16,17]. Up to date, various researchers [18–30] have introduced the fractional calculus approach into various rheological problems at macroscale and gained great success. To our best knowledge, no exaction solution for the electro-osmotic flow of viscoelastic fluids has been reported yet. In present study the generalized Oldroyd-B constitutive model [23–25,29] constructed from original Oldroyd-B model using the concept of fractional calculus is adopted to describe the elastic and memory effects exhibited by most polymer and biological liquids. By using finite Fourier transform method, one can show that the whole problem (two dimensional and time-dependent) is transformed to just determining the time-dependent component in a twofold infinite series. Subsequently, exact solutions of the velocity distribution for electro-osmotic flow of viscoelastic fluids in a rectangular micro-channel are obtained in terms of generalized Mittag–Leffler functions for some special cases by means of discrete inverse of Laplace transform of the sequential fractional derivatives.
2. Governing equations and generalized Oldroyd-B model The fundamental equations governing the unsteady motion of an incompressible viscoelastic fluid are continuity equation
rV ¼0
ð1Þ
and the general Cauchy momentum equation
q
oV þ ðV rÞV ¼ r r þ F; ot
ð2Þ
where V is the velocity vector, q is the fluid density, r is the Cauchy stress tensor, F is the external body force vector. The Cauchy stress tensor r for a viscoelastic fluid with the generalized Oldroyd-B model [23–25,29] is
r ¼ pI þ s;
! b Da b D A1 ; 1þk s¼l 1þh Dta Dt b a
ð3Þ
in which p is the pressure, I the identity tensor, s is the shear stress tensor, k and h is the relaxation time and the retardation time respectively, l the dynamic viscosity of the fluid, a and b are fractional calculus parameters such that 0 6 a 6 b 6 1. For a > b the relaxation fraction is increasing, which is generally not reasonable [17]. First Rivlin–Ericksen tensor A1 is given by
A1 ¼ rV þ ðrVÞT ;
ð4Þ a
a
where r is the gradient operator and the superscript T denotes the tensor transpose. D /Dt refers to the so-called upper convected time derivative of order a – a special time derivative devised for constitutive equations to meet the requirements of continuum mechanics [31]
Da s oa s ¼ þ ðV rÞs ðrVÞ s s ðrVÞT Dta ot a
ð5Þ
op/otp is the fractional differentiation operator of order p with respect to t, and may be defined as [16,32] p
op ½f ðtÞ d ½f ðtÞ 1 d ¼ ¼ p otp Cð1 pÞ dt dt
Z 0
t
f ðsÞ ds; ðt sÞp
0 6 p 6 1;
ð6Þ
where C() is the Gamma function and op/otp is equivalent to dp/dtp if the function f is just dependent on time t. Actually, generalized Oldroyd-B model encompasses a large class of fluids. When k = h = 0, it reduces to Newtonian fluid. When h = 0 and k – 0, it reduces to fractional Maxwell model and when k = 0 and h – 0, it reduces to generalized second grade fluid. 3. Exact solutions for electro-osmotic flow of generalized Oldroyd-B fluids Fig. 1 shows the dimensions of the micro-channel and the coordinate systems considered in this work. The channel is filled with a liquid solution of dielectric constant, e. It is assumed that the channel wall is uniformly charged with a zeta potential, ww, and the liquid solution is a typical viscoelastic fluid whose behavior can be described by the well-known generalized Oldroyd-B model. When an external electric field E0f(t) is imposed along the x-axis direction, the fluid in the microchannel sets in motion due to electro-osmosis. f(t) represents a time-dependent electric field and concrete forms will be given in Section 4. Because of symmetry, the analysis is restricted in the first quadrant of the y–z plane. For unidirectional flow, we consider the velocity of the form
V ¼ uðy; z; tÞi;
ð7Þ
504
C. Zhao, C. Yang / Applied Mathematics and Computation 211 (2009) 502–509
2H
y
x
o
z
2W Fig. 1. Configuration of the rectangular micro-channel.
where u is the x-component of velocity and i is the unit vector in the x-direction. Thus using Eq. (7), the continuity Eq. (1) is satisfied automatically. Furthermore, for electro-osmotic flow, the only driving force is due to the interaction of the applied electrical field E0f(t) and the net charge density in the EDL region near the channel wall. Such force acts only along x-direction, and is expressed by [1]
F x ¼ qe E0 f ðtÞ:
ð8Þ
For an open-end horizontally placed channel, no pressure gradient is induced and hence the pressure gradient term in the Cauchy momentum equation disappears. According to the theory of electrostatics, the net charge density qe in the diffuse layer of the wall EDL is given by the Poisson equation, which takes the form
o2 w o2 w q þ ¼ e: oy2 oz2 e
ð9Þ
¼ zv ew=ðkB TÞ and with the assumptions of the ¼ y=Dh ; z ¼ z=Dh ; K ¼ jDh , and w Introducing the dimensionless groups: y Boltzmann distribution and a small surface (zeta) potential, the electrical potential profile in the EDL is governed by the linearized Poisson–Boltzmann equation expressed by [1,33]
o2 w o2 w þ ¼ K 2 w; 2 oz2 oy
ð10Þ
which is subject to the following boundary conditions:
wj wj z¼W=Dh ¼ f; y¼H=Dh ¼ f; ow ow ¼ 0; ¼ 0: y¼0 oy oz z¼0
ð11aÞ ð11bÞ
Dh is the hydrodynamic diameter of the rectangular micro-channel and is defined as Dh = 4HW/(H + W), the dimensionless wall zeta potential is given by f ¼ zv eww =kB T, the Debye length j1 is defined as j1 ¼ ðekB T=2e2 z2v n1 Þ1=2 , where n1 and zv are the bulk number concentration and the valence of ions, respectively, e is the fundamental charge, kB is the Boltzmann constant, and T is the absolute temperature. The solution for the electrical potential distribution in the EDL region is found by using the separation of variable method and has the following form:
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi m 2 2 Þ ðK þ m z ð1Þ cosh ðK 2 þ c2m Þy 1 1 n X X Þ þ 2f ; zÞ ¼ 2f qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi cosðmn y qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi cosðcmzÞ; wðy n¼0 n þ 1 p cosh m¼0 m þ 1 p cosh ðK 2 þ m2n Þ DW ðK 2 þ c2m Þ DH 2 2 ð1Þn cosh
h
ð12Þ
h
where the two eigenvalues are given by
mn ¼
Dh 1 nþ p; 2 H
cm ¼
Dh 1 mþ p; 2 W
m; n ¼ 0; 1; 2; . . .
ð13Þ
For the sake of convenience, a geometrical function are defined as
Pðy; zÞ ¼ 2
1 X n¼0
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi m ðK 2 þ m2n Þz ð1Þ cosh ðK 2 þ c2m Þy 1 X Þ þ 2 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi cosðmn y qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi cosðcm zÞ: m¼0 m þ 1 p cosh ðK 2 þ m2n Þ DW ðK 2 þ c2m Þ DH n þ 12 p cosh 2 ð1Þn cosh
h
ð14Þ
h
And thus the electrical potential distribution in the EDL can be expressed concisely as
¼ fPðy ; zÞ: w
ð15Þ
C. Zhao, C. Yang / Applied Mathematics and Computation 211 (2009) 502–509
505
Meanwhile the net charge density qe can be expressed as a function of the EDL potential [1],
qe ðy; zÞ ¼ j2 ew:
ð16Þ
For this specific problem the general constitutive relationship Eq. (3) becomes
oa ob ou 1 þ ka a syx ¼ l 1 þ hb b ; ot ot oy oa ob ou 1 þ ka a szx ¼ l 1 þ hb b : ot ot oz
ð17aÞ ð17bÞ
The general Cauchy momentum equation given by Eq. (2) can be simplified to
q
ou osyx oszx ¼ þ j2 ewE0 f ðtÞ; ot oy oz
ð18Þ
which has the following initial and boundary conditions:
ujt¼0 ¼
ou ¼ 0; ot t¼0
ð19aÞ
ujy¼H ¼ ujz¼W ¼ 0; ou ou ¼ ¼ 0: oy y¼0 oz z¼0
ð19bÞ ð19cÞ
From Eqs. (17) and (18), we can get the governing equation of the velocity distribution
! 2 oa ou ob o u o2 u 1 þ ka a q þ j2 ewE0 f ðtÞ ¼ l 1 þ hb b þ ot oy2 oz2 ot ot
ð20Þ
Non-dimensionalizing Eq. (20) using the below listed non-dimensional parameters
¼ u
u ; us
t ¼
l l l ew E0 t; k ¼ k; h ¼ h; us ¼ w ; K ¼ jDh l qD2h qD2h qD2h
ð21Þ
we can show that the non-dimensional counterpart for Eq. (20) takes the form
! 2 o2 u oa ou ob o u 2 ; zÞf ðtÞ ¼ 1 þ hb : 1 þ ka a P ð y þ K 2 oz2 oy ot ot b ot
ð22Þ
Also, it is straightforward that the corresponding dimensionless initial and boundary conditions are
ou ujt¼0 ¼ ¼ 0; ot t¼0 jy¼H=D ¼ u jz¼W=D ¼ 0; u h h ou ou ¼ ¼ 0: y¼0 oz z¼0 oy
ð23aÞ ð23bÞ ð23cÞ
According to the fast finite Fourier transform (FFT) method [34], the exact solution of velocity determined by Eqs. (22) and (23) can be expressed in the form of an twofold infinite series
ðy ; z; tÞ ¼ u
1 X 1 X
ÞZ m ðzÞ; T nm ðtÞY n ðy
ð24Þ
n¼0 m¼0
where two orthogonal basis functions are, respectively,
rffiffiffiffiffiffiffiffiffi 2Dh ; cos mn y H rffiffiffiffiffiffiffiffiffi 2Dh Z m ðzÞ ¼ cos cm z: W Þ ¼ Y n ðy
ð25aÞ ð25bÞ
Here, the eigenvalues mn and cm are the same as given by Eq. (13). By using the basis functions, the original fractional partial differential equation (22) and its initial and boundary conditions Eq. (23) can be modified to a fractional ordinary equation with respect to Tnm(t) by using finite Fourier transform
! a a b dT nm 2 a d a d 2 2 b d K Anm 1 þ k a f ðtÞ ¼ ðmn þ cm Þ 1 þ h b T nm 1 þ k a dt dt dt dt
ð26Þ
506
C. Zhao, C. Yang / Applied Mathematics and Computation 211 (2009) 502–509
with a simple initial boundary condition
T nm jt¼0
dT nm ¼ ¼ 0; dt t¼0
ð27Þ
where
2Dh ð1Þnþm mn cm : Anm ¼ pffiffiffiffiffiffiffiffiffi 2 þ HW K þ m2n þ c2m cm mn
ð28Þ
Subsequently, the Laplace transform and inverse Laplace transforms
T snm ðsÞ ¼ L½T nm ðtÞ ¼
Z
þ1
T nm ðtÞest dt;
ð29aÞ
0
Z rþj1 1 T s ðsÞest ds T nm ðtÞ ¼ L1 ½T snm ðsÞ ¼ 2pj rj1 nm
ð29bÞ
are applied to solve Eq. (26), where r is a real number so that the contour path of integration is in the region of convergence of T snm ðsÞ normally requiring r > Re(sp) for every singularity sp of T snm ðsÞ and j2 = 1. Transforming both sides of Eq. (26) by taking into account the initial condition expressed by Eq. (27) and then we can get
ð1 þ ka sa ÞsT snm ðsÞ K 2 Anm ð1 þ ka sa ÞFðsÞ ¼ ðm2n þ c2m Þð1 þ hb sb ÞT snm ðsÞ:
ð30Þ
Hence, the Laplace transform of Tnm(t) can be shown to bear the form
T snm ðsÞ ¼
K 2 Anm ð1 þ ka sa ÞFðsÞ ; ð1 þ ka sa Þs þ ðm2n þ c2m Þð1 þ hb sb Þ
ð31Þ
where FðsÞ ¼ L½f ðtÞ. We can obtain the exact solution to the model if the inverse Laplace transform of T snm ðsÞ is known. At last, according to Eq. (24), the exact solution can be further shown as
ðy ; z; tÞ ¼ u
1 X 1 X
ÞZ m ðzÞ ¼ T nm ðtÞY n ðy
n¼0 m¼0
1 X 1 X
ÞZ m ðzÞ: L1 ½T snm ðsÞY n ðy
ð32Þ
n¼0 m¼0
4. Special cases Up to now, we only need to know Tnm(t) to construct exact solutions from Eq. (32). But for generalized Oldroyd-B fluids ( k – 0; h – 0Þ, it is found that Eq. (31) can not be inversely transformed into to the time domain analytically, one can resort to numerical methods, such as Durbin’s method [35], Stehfest’s method [36], etc. However, for some special types of generalized Oldroyd-B fluids, such as Newtonian fluids ð k¼ h ¼ 0Þ, Fractional Maxwell fluids ð k – 0; h ¼ 0Þ and generalized second grade fluids ð k ¼ 0; h – 0Þ, the analytical solution does exist and will be given below. In view of expressing formulae compactly, only results of Tnm(t) are provided. 4.1. Newtonian fluid model ð k¼ h ¼ 0Þ As mentioned in Section 2, the Newtonian fluid is just a special kind of generalized Oldroyd-B fluid with k ¼ 0 and h ¼ 0, from general formula Eq. (31), the Laplace transform of Tnm(t) is of the form
T snm ðsÞ ¼
K 2 Anm FðsÞ : s þ ðm2n þ c2m Þ
ð33Þ
Also, in this paper, two most widely used electric fields in practical applications are considered, which will be reflected in the different forms of f(t) and thus different F(s). Case I: the external electric field is applied and remains constant from the time t = 0 (i.e., the electric field follows stepchange):
f ðtÞ ¼ HðtÞ;
ð34Þ
where H(t) is the well-known Heaviside unit step function. Then, Eq. (33) becomes
T snm ðsÞ ¼
K 2 Anm s½s þ ðm2n þ c2m Þ
ð35Þ
507
C. Zhao, C. Yang / Applied Mathematics and Computation 211 (2009) 502–509
and its inverse Laplace transform reads
T nm ðtÞ ¼
K 2 Anm 2 2 ½1 eðmn þcm Þt : m2n þ c2m
ð36Þ
Case II: we have a unit impulsive electric field applied at t=0, which means that f(t) = d(t), where d(t) represents Dirac delta function. Upon knowing this, Eq. (33) reduces to
T snm ðsÞ ¼
K 2 Anm s þ ðm2n þ c2m Þ
ð37Þ
and the corresponding inverse Laplace transform is 2
2
T nm ðtÞ ¼ K 2 Anm eðmn þcm Þt :
ð38Þ
4.2. Fractional Maxwell fluid model ( h ¼ 0 and k – 0) With h ¼ 0 in Eq. (31), the Laplace transform of Tnm (t) for fractional Maxwell fluids takes the form
T snm ðsÞ ¼
K 2 Anm ð1 þ ka sa ÞFðsÞ : ð1 þ ka sa Þs þ ðm2n þ c2m Þ
ð39Þ
Just like in 4.1, we just consider two forms of electric field namely a step-change electric field and an impulsive electric field. In the rest of this work, without otherwise notifications, Case I stands for the case of step-change electric field and Case II represents the case of impulsive electric field. Case I: Eq. (39) is changed to
T snm ðsÞ ¼
K 2 Anm ð1 þ ka sa Þ Bnm 2 2 ka sa2 Bnm ; ¼ K A s þ nm 1 þ Bnm 1 þ Bnm s½ð1 þ ka sa Þs þ ðm2n þ c2m Þ
ð40Þ
where
Bnm ¼
ka s : s1þa þ ðm2n þ c2m Þka
ð41Þ
In order to avoid the burdensome calculation of residues and contour integrals in Eq. (29b), we will apply the discrete inverse Laplace transform to get Tnm(t). At first, one can find dk and make dk to satisfy when Re(s) > Re(dk), jBnmj < 1, then Eq. (40) can be expanded as
T snm ðsÞ ¼ K 2 Anm
1 X ð1Þk kaðkþ1Þ k¼0
sk1 2 Þ a kþ1 m k
½s1þa þ ðm2n þ c
þ K 2 Anm
1 X
ð1Þk kak
k¼0
sk1þa ½s1þa þ ðm2n þ c2m Þka kþ1
:
ð42Þ
From [16], we obtain an important Laplace transform of the Mittag–Leffler function
h i p L t pnþq1 EðnÞ p;q ðat Þ ¼
n!spq ðs aÞnþ1
;
ð43Þ
where Ep,q(z) is the generalized Mittag–Leffler function and its derivative of order n is expressed by n
EðnÞ p;q ðzÞ ¼
1 X d ðj þ nÞ!zj : n Ep;q ðzÞ ¼ j!Cðpj þ pn þ qÞ dz j¼0
ð44Þ
The term-by-term inversion of Eq. (42), based on the general expansion theorem for the Laplace transform, using Eq. (43), it can be shown that the solution in the time domain is
(
1 X
ð1Þk aðkþ1Þaðkþ1Þþ1 ðkÞ t E1þa;kþaþ2 ðm2n þ c2m Þkat1þa k k! k¼0 ) 1 X
ð1Þk akakþ1 ðkÞ þ E1þa;kþ2 ðm2n þ c2m Þkat1þa : k t k! k¼0
T nm ðtÞ ¼ K 2 Anm
ð45Þ
Case II: for impulsive electric field, following the similar procedure in Case I, the results are consecutively shown below:
T snm ðsÞ ¼
K 2 Anm ð1 þ ka sa Þ Bnm 2 1 ka sa1 Bnm ¼ K A s þ ; nm 1 þ Bnm 1 þ Bnm ð1 þ ka sa Þs þ ðm2n þ c2m Þ
ð46Þ
508
C. Zhao, C. Yang / Applied Mathematics and Computation 211 (2009) 502–509
where Bnm has been given by Eq. (41). The Taylor expansion produces 1 X ð1Þk kaðkþ1Þ
T snm ðsÞ ¼ K 2 Anm
k¼0
sk ½s1þa þ ðm2n þ c
2 Þ a kþ1 m k
þ K 2 Anm
1 X ð1Þk kak k¼0
skþa ½s1þa þ ðm2n þ c2m Þka kþ1
:
ð47Þ
Finally, the inverse Laplace transform of Eq. (47) result in 2
T nm ðtÞ ¼ K Anm
( 1 X ð1Þk k¼0
) 1
X
ð1Þk akak ðkÞ 2 2 a1þa 2 2 a1þa kaðkþ1Þt aðkþ1Þ EðkÞ : t t t E1þa;kþ1 ðmn þ cm Þk k þ 1þa;kþ1þa ðmn þ cm Þk k! k! k¼0 ð48Þ
4.3. Generalized second grade fluid model ( k ¼ 0 and h – 0) For the generalized second grade fluids, it is straightforward that Eq. (31) reduces to
T snm ðsÞ ¼
K 2 Anm FðsÞ : s þ ðm þ c2m Þð1 þ hb sb Þ
ð49Þ
2 n
For both Case I and Case II, implementing similar procedure in 4.2, we can show the results: Case I: the Laplace transform of the time-dependent part Tnm (t) in Eq. (24) is expressed as
T snm ðsÞ ¼
K 2 Anm K 2 Anm 1 C nm ¼ 2 s ; 2 2 b b m 1 þ C nm þ cm Þð1 þ h s Þ n þ cm
s½s þ ðm
2 n
ð50Þ
where
C nm ¼
ðm2n þ c2m Þsb : þ ðm2n þ c2m Þhb
ð51Þ
s1b
Expanding Eq. (51) in Taylor series form as 1 X ðm2n þ c2m Þk sbðkþ1Þ1 ð1Þk : ½s1b þ ðm2n þ c2m Þhb kþ1 k¼0
T snm ðsÞ ¼ K 2 Anm
ð52Þ
Then, using discrete Laplace transform of the sequential fractional derivatives, the solution in the time domain is shown as
T nm ðtÞ ¼ K 2 Anm
1 X ð1Þk 2 ðkÞ ðmn þ c2m Þktkþ1 E1b;2þkb ½ðm2n þ c2m Þhbt1b : k! k¼0
ð53Þ
Case II:
T snm ðsÞ ¼
K 2 Anm sþ
ðk2n
þ c2m Þð1 þ hb sb Þ
ð54Þ
:
Eq. (55) can be expanded as
T snm ðsÞ ¼ K 2 Anm
1 X ð1Þk k¼0
ðm2n þ c2m Þk sbðkþ1Þ : þ ðm2 þ c2 Þhb kþ1
½s1b
n
ð55Þ
m
Using discrete inverse Laplace transform, we can get Tnm(t) as
T nm ðtÞ ¼ K 2 Anm
1 X
ð1Þk 2 ðkÞ ðmn þ c2m Þkt k E1b;bkþ1 ðm2n þ c2m Þhbt 1b : k! k¼0
ð56Þ
5. Concluding remark The generalized Oldroyd-B constitutive model with fractional calculus is introduced to describe the viscoelastic behavior of biological fluids in BioMEMS and Lab-on-a-chip devices. It is noted that the generalized Oldroyd-B model has extensive adaptability as it is a generalization of Newtonian model, fractional Maxwell model and fractional second grade model. At first, we use finite Fourier transform method to construct a general twofold infinite series solution for the velocity distribution of electro-osmotic flow of viscoelastic fluids in a rectangular micro-channel and however the analytical solution is not completed known unless we can obtain the inverse Laplace transform of the time-dependent component T snm ðsÞ. Then, for three special generalized Oldroyd-B models (Newtonian model, fractional Maxwell model and generalized second grade
C. Zhao, C. Yang / Applied Mathematics and Computation 211 (2009) 502–509
509
model), explicit analytical solutions of velocity distribution are obtained by using the discrete Laplace transform of the sequential fractional derivatives and the generalized Mittag–Leffler function. These exaction solutions are of potential importance because it can describe the flow behavior of electro-osmotic flow of biological fluids in Lab-on-a-chip devices. Moreover, they may also serve as accuracy checks for experimental, numerical and asymptotic methods. Acknowledgements The Ph.D. scholarship awarded to Z.C.L. from Nanyang Technological University is sincerely appreciated. The authors also would like to thank the referees for their useful comments and suggestions regarding an earlier version of this paper. References [1] J.H. Masliyah, S. Bhattacharjee, Electrokinetic and Colloid Transport Phenomena, Wiley Interscience, Hoboken, NJ, 2006. [2] A. Wainright, U.T. Nguyen, T. Bjornson, T.D. Boone, Preconcentration and separation of double-stranded DNA fragments by electrophoresis in plastic microfluidic devices, Electrophoresis 24 (2003) 3784–3792. [3] G.M. Whitesides, The origins and the future of microfluidics, Nature 442 (2006) 368–373. [4] A.E. Herr, J.I. Molho, J.G. Santiago, M.G. Mungal, T.W. Kenny, M.G. Garguilo, Electroosmotic capillary flow with nonuniform zeta potential, Anal. Chem. 72 (2000) 1053–1057. [5] C.L. Rice, R. Whitehead, Electrokinetic flow in a narrow cylindrical capillary, J. Phys. Chem. 69 (1965) 4017–4024. [6] J.G. Santiago, Electroosmotic flows in microchannels with finite inertial and pressure forces, Anal. Chem. 73 (2001) 2353–2365. [7] R.J. Yang, L.M. Fu, Y.C. Lin, Electroosmotic flow in microchannels, J. Colloid Interface Sci. 239 (2001) 98–105. [8] S. Devasenathipathy, J.G. Santiago, K. Takehara, Particle tracking techniques for electrokinetic microchannel flows, Anal. Chem. 74 (2002) 3704–3713. [9] M.J. Pikal, S. Shah, Transport mechanisms in iontophoresis: III. An experimental study of the contributions of electroosmotic flow and permeability change in transport of low and high molecular weight solutes, Pharm. Res. 7 (1990) 222–229. [10] D. Ross, T.J. Johnson, L.E. Locascio, Imaging of electroosmotic flow in plastic microchannels, Anal. Chem. 73 (2001) 2509–2515. [11] D. Sinton, D. Li, Electroosmotic velocity profiles in microchannels, Colloids Surf., A 222 (2003) 273–283. [12] P. Dutta, A. Beskok, T.C. Warburton, Numerical simulation of mixed electroosmotic/pressure driven microflows, Numer. Heat Transfer, Part A 41 (2002) 131–148. [13] Y. Hu, C. Werner, D. Li, Electrokinetic transport through rough microchannels, Anal. Chem. 75 (2003) 5747–5758. [14] N.A. Patankar, H.H. Hu, Numerical simulation of electroosmotic flow, Anal. Chem. 70 (1998) 1870–1881. [15] L. Ren, D. Li, Electroosmotic flow in heterogeneous microchannels, J. Colloid Interface Sci. 243 (2001) 255–261. [16] I. Podlubny, Fractional Differential Equations, Academic Press, San Diego, 1999. [17] R. Hilfer, Applications of Fractional Calculus in Physics, World Scientific Press, Singapore, 2000. [18] W. Tan, F. Xian, L. Wei, An exact solution of unsteady Couette flow of generalized second grade fluid, Chin. Sci. Bull. 47 (2002) 1783–1785. [19] W. Tan, M. Xu, Plane surface suddenly set in motion in a viscoelastic fluid with fractional Maxwell model, Acta Mech. Sin./Lixue Xuebao 18 (2002) 342– 349. [20] W. Tan, W. Pan, M. Xu, A note on unsteady flows of a viscoelastic fluid with the fractional Maxwell model between two parallel plates, Int. J. Non Linear Mech. 38 (2003) 645–650. [21] T. Hayat, S. Nadeem, S. Asghar, Periodic unidirectional flows of a viscoelastic fluid with the fractional Maxwell model, Appl. Math. Comput. 151 (2004) 153–161. [22] R.L. Magin, Fractional calculus in bioengineering, part 3, Crit. Rev. Biomed. Eng. 32 (2004) 195–377. [23] S. Hyder Ali, M. Khan, H. Qi, Exact solutions for a viscoelastic fluid with the generalized Oldroyd-B model, Nonlinear Anal.: Real World Applications (2008), doi:10.1016/j.nonrwa.2008.03.012. [24] M. Khan, S. Hyder Ali, H. Qi, Some accelerated flows for a generalized Oldroyd-B fluid, Nonlinear Anal.: Real World Applications (2007), doi:10.1016/ j.nonrwa.2007.11.017. [25] M. Khan, T. Hayat, S. Asghar, Exact solution for MHD flow of a generalized Oldroyd-B fluid with modified Darcy’s law, Int. J. Eng. Sci. 44 (2006) 333–339. [26] D. Vieru, C. Fetecau, C. Fetecau, Flow of a viscoelastic fluid with the fractional Maxwell model between two side walls perpendicular to a plate, Appl. Math. Comput. 200 (2008) 459–464. [27] A. Hernandez-Jimenez, J. Hernandez-Santiago, A. Macias-Garca, J. Sanchez-Gonzalez, Relaxation modulus in PMMA and PTFE fitting by fractional Maxwell model, Polym. Test. 21 (2002) 325–331. [28] T. Hayat, M. Khan, S. Asghar, On the MHD flow of fractional generalized Burgers’ fluid with modified Darcy’s law, Acta Mech. Sin./Lixue Xuebao 23 (2007) 257–261. [29] M. Khan, K. Maqbool, T. Hayat, Influence of Hall current on the flows of a generalized Oldroyd-B fluid in a porous space, Acta Mech. 184 (2006) 1–13. [30] M. Khan, S. Nadeem, T. Hayat, A.M. Siddiqui, Unsteady motions of a generalized second-grade fluid, Math. Comput. Model. 41 (2005) 629–637. [31] D.D. Joseph, Fluid Dynamics of Viscoelastic Liquids, Springer-Verlag, New York, 1990. [32] L. Debnath, D. Bhatta, Integral Transforms and Their Applications, Taylor & Francis Group, New York, 2007. [33] D. Li, Electrokinetics in Microfluidics, Elsevier Academic Press, 2004. [34] W.M. Deen, Analysis of Transport Phenomena, Oxford University Press, New York, 1998. [35] F. Durbin, Numerical inversion of Laplace transforms: an efficient improvement to Dubner and Abate’s method, Comput. J. 17 (1974) 371–376. [36] H. Stehfest, Algorithm 368: numerical inversion of Laplace transforms [D5], Commun. ACM 13 (1970) 47–49.
Contents lists available at ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
Exact solutions for electro-osmotic flow of viscoelastic fluids in rectangular micro-channels Cunlu Zhao *, Chun Yang School of Mechanical and Aerospace Engineering, Nanyang Technological University, 50 Nanyang Avenue, 639798, Singapore
a r t i c l e
i n f o
Keywords: Electro-osmotic flow Lab-on-a-chip Viscoelastic fluids Generalized Oldroyd-B model
a b s t r a c t Transient electro-osmotic flow of viscoelastic fluids in rectangular micro-channels is investigated. The general twofold series solution for the velocity distribution of electro-osmotic flow of viscoelastic fluids with generalized fractional Oldroyd-B constitutive model is obtained by using finite Fourier and Laplace transforms. Under three limiting cases, the generalized Oldroyd-B model simplifies to Newtonian model, fractional Maxwell model and generalized second grade model, where all the explicit exact solutions for velocity distribution are found through the discrete Laplace transform of the sequential fractional derivatives. These exact solutions may be able to predict the flow behavior of viscoelastic biological fluids in BioMEMS and Lab-on-a-chip devices and thus could benefit the design of these devices. Ó 2009 Elsevier Inc. All rights reserved.
1. Introduction Accompanied by the recently fast and thriving development of BioMEMS and Lab-on-a-chip technologies, renewed interests are also being brought up in electro-osmosis due to the fact that electro-osmosis has been demonstrated as a powerful activation mechanism to manipulate liquid samples of nanovolumes in capillary network of these devices. In reviewing the general concept of electro-osmosis, liquids at the solid–liquid interface ionize because of the zeta potential at the wall. The electric double layer (EDL) is a very thin region close to the charged surface in which there is an excess of counterions over coions to neutralize the surface charge. Therefore, a net fluid flow in the desired direction is induced when imposing external electric fields [1]. Compared to conventional mechanical pumping modes, electro-osmosis-based pumping has numerous advantages, including ease of fabrication and control, no movable mechanical parts and high reliability, no noise, etc. More favorably, the fluid moves as a plug in the electro-osmosis, rather than with the parabolic-flow profile observed when pumping is accomplished by applying pressure gradient to the fluid. This feature helps to minimize the broadening of plugs of sample that occurs with many pressure-driven systems, and allows very high-resolution separations of ionic species. It is a key contributor to electrophoretic separations of DNA in micro-channels [2,3]. Extensive studies of electro-osmotic flow of Newtonian fluids in micro-channels have been reported in the literature theoretically [4–7], experimentally [8–11] and numerically [12–15]. However, forenamed studies are just concerned with Newtonian fluids. Actually, Lab-on-a-chip devices are frequently used to analyze biofluids which usually exhibit viscoelastic behavior. Thus, the more general Cauchy momentum equation, instead of the Navier–Stokes equation should be used to describe the flow characteristics of viscoelastic fluids provided that proper constitutive equations are available. Recently, fractional calculus has gained considerable interest due to its applications in different areas of physics and engineering, including complex fluid dynamics. The constitutive equations with fractional derivative have been proved to be valuable * Corresponding author. E-mail address: [email protected] (C. Zhao). 0096-3003/$ - see front matter Ó 2009 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2009.01.068
503
C. Zhao, C. Yang / Applied Mathematics and Computation 211 (2009) 502–509
tools to handle viscoelastic properties. In general, these equations are derived from known models via substituting time ordinary derivatives of stress and strain by derivatives of fractional order from which one can define the non-integer order integral and derivatives [16,17]. Up to date, various researchers [18–30] have introduced the fractional calculus approach into various rheological problems at macroscale and gained great success. To our best knowledge, no exaction solution for the electro-osmotic flow of viscoelastic fluids has been reported yet. In present study the generalized Oldroyd-B constitutive model [23–25,29] constructed from original Oldroyd-B model using the concept of fractional calculus is adopted to describe the elastic and memory effects exhibited by most polymer and biological liquids. By using finite Fourier transform method, one can show that the whole problem (two dimensional and time-dependent) is transformed to just determining the time-dependent component in a twofold infinite series. Subsequently, exact solutions of the velocity distribution for electro-osmotic flow of viscoelastic fluids in a rectangular micro-channel are obtained in terms of generalized Mittag–Leffler functions for some special cases by means of discrete inverse of Laplace transform of the sequential fractional derivatives.
2. Governing equations and generalized Oldroyd-B model The fundamental equations governing the unsteady motion of an incompressible viscoelastic fluid are continuity equation
rV ¼0
ð1Þ
and the general Cauchy momentum equation
q
oV þ ðV rÞV ¼ r r þ F; ot
ð2Þ
where V is the velocity vector, q is the fluid density, r is the Cauchy stress tensor, F is the external body force vector. The Cauchy stress tensor r for a viscoelastic fluid with the generalized Oldroyd-B model [23–25,29] is
r ¼ pI þ s;
! b Da b D A1 ; 1þk s¼l 1þh Dta Dt b a
ð3Þ
in which p is the pressure, I the identity tensor, s is the shear stress tensor, k and h is the relaxation time and the retardation time respectively, l the dynamic viscosity of the fluid, a and b are fractional calculus parameters such that 0 6 a 6 b 6 1. For a > b the relaxation fraction is increasing, which is generally not reasonable [17]. First Rivlin–Ericksen tensor A1 is given by
A1 ¼ rV þ ðrVÞT ;
ð4Þ a
a
where r is the gradient operator and the superscript T denotes the tensor transpose. D /Dt refers to the so-called upper convected time derivative of order a – a special time derivative devised for constitutive equations to meet the requirements of continuum mechanics [31]
Da s oa s ¼ þ ðV rÞs ðrVÞ s s ðrVÞT Dta ot a
ð5Þ
op/otp is the fractional differentiation operator of order p with respect to t, and may be defined as [16,32] p
op ½f ðtÞ d ½f ðtÞ 1 d ¼ ¼ p otp Cð1 pÞ dt dt
Z 0
t
f ðsÞ ds; ðt sÞp
0 6 p 6 1;
ð6Þ
where C() is the Gamma function and op/otp is equivalent to dp/dtp if the function f is just dependent on time t. Actually, generalized Oldroyd-B model encompasses a large class of fluids. When k = h = 0, it reduces to Newtonian fluid. When h = 0 and k – 0, it reduces to fractional Maxwell model and when k = 0 and h – 0, it reduces to generalized second grade fluid. 3. Exact solutions for electro-osmotic flow of generalized Oldroyd-B fluids Fig. 1 shows the dimensions of the micro-channel and the coordinate systems considered in this work. The channel is filled with a liquid solution of dielectric constant, e. It is assumed that the channel wall is uniformly charged with a zeta potential, ww, and the liquid solution is a typical viscoelastic fluid whose behavior can be described by the well-known generalized Oldroyd-B model. When an external electric field E0f(t) is imposed along the x-axis direction, the fluid in the microchannel sets in motion due to electro-osmosis. f(t) represents a time-dependent electric field and concrete forms will be given in Section 4. Because of symmetry, the analysis is restricted in the first quadrant of the y–z plane. For unidirectional flow, we consider the velocity of the form
V ¼ uðy; z; tÞi;
ð7Þ
504
C. Zhao, C. Yang / Applied Mathematics and Computation 211 (2009) 502–509
2H
y
x
o
z
2W Fig. 1. Configuration of the rectangular micro-channel.
where u is the x-component of velocity and i is the unit vector in the x-direction. Thus using Eq. (7), the continuity Eq. (1) is satisfied automatically. Furthermore, for electro-osmotic flow, the only driving force is due to the interaction of the applied electrical field E0f(t) and the net charge density in the EDL region near the channel wall. Such force acts only along x-direction, and is expressed by [1]
F x ¼ qe E0 f ðtÞ:
ð8Þ
For an open-end horizontally placed channel, no pressure gradient is induced and hence the pressure gradient term in the Cauchy momentum equation disappears. According to the theory of electrostatics, the net charge density qe in the diffuse layer of the wall EDL is given by the Poisson equation, which takes the form
o2 w o2 w q þ ¼ e: oy2 oz2 e
ð9Þ
¼ zv ew=ðkB TÞ and with the assumptions of the ¼ y=Dh ; z ¼ z=Dh ; K ¼ jDh , and w Introducing the dimensionless groups: y Boltzmann distribution and a small surface (zeta) potential, the electrical potential profile in the EDL is governed by the linearized Poisson–Boltzmann equation expressed by [1,33]
o2 w o2 w þ ¼ K 2 w; 2 oz2 oy
ð10Þ
which is subject to the following boundary conditions:
wj wj z¼W=Dh ¼ f; y¼H=Dh ¼ f; ow ow ¼ 0; ¼ 0: y¼0 oy oz z¼0
ð11aÞ ð11bÞ
Dh is the hydrodynamic diameter of the rectangular micro-channel and is defined as Dh = 4HW/(H + W), the dimensionless wall zeta potential is given by f ¼ zv eww =kB T, the Debye length j1 is defined as j1 ¼ ðekB T=2e2 z2v n1 Þ1=2 , where n1 and zv are the bulk number concentration and the valence of ions, respectively, e is the fundamental charge, kB is the Boltzmann constant, and T is the absolute temperature. The solution for the electrical potential distribution in the EDL region is found by using the separation of variable method and has the following form:
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi m 2 2 Þ ðK þ m z ð1Þ cosh ðK 2 þ c2m Þy 1 1 n X X Þ þ 2f ; zÞ ¼ 2f qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi cosðmn y qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi cosðcmzÞ; wðy n¼0 n þ 1 p cosh m¼0 m þ 1 p cosh ðK 2 þ m2n Þ DW ðK 2 þ c2m Þ DH 2 2 ð1Þn cosh
h
ð12Þ
h
where the two eigenvalues are given by
mn ¼
Dh 1 nþ p; 2 H
cm ¼
Dh 1 mþ p; 2 W
m; n ¼ 0; 1; 2; . . .
ð13Þ
For the sake of convenience, a geometrical function are defined as
Pðy; zÞ ¼ 2
1 X n¼0
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi m ðK 2 þ m2n Þz ð1Þ cosh ðK 2 þ c2m Þy 1 X Þ þ 2 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi cosðmn y qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi cosðcm zÞ: m¼0 m þ 1 p cosh ðK 2 þ m2n Þ DW ðK 2 þ c2m Þ DH n þ 12 p cosh 2 ð1Þn cosh
h
ð14Þ
h
And thus the electrical potential distribution in the EDL can be expressed concisely as
¼ fPðy ; zÞ: w
ð15Þ
C. Zhao, C. Yang / Applied Mathematics and Computation 211 (2009) 502–509
505
Meanwhile the net charge density qe can be expressed as a function of the EDL potential [1],
qe ðy; zÞ ¼ j2 ew:
ð16Þ
For this specific problem the general constitutive relationship Eq. (3) becomes
oa ob ou 1 þ ka a syx ¼ l 1 þ hb b ; ot ot oy oa ob ou 1 þ ka a szx ¼ l 1 þ hb b : ot ot oz
ð17aÞ ð17bÞ
The general Cauchy momentum equation given by Eq. (2) can be simplified to
q
ou osyx oszx ¼ þ j2 ewE0 f ðtÞ; ot oy oz
ð18Þ
which has the following initial and boundary conditions:
ujt¼0 ¼
ou ¼ 0; ot t¼0
ð19aÞ
ujy¼H ¼ ujz¼W ¼ 0; ou ou ¼ ¼ 0: oy y¼0 oz z¼0
ð19bÞ ð19cÞ
From Eqs. (17) and (18), we can get the governing equation of the velocity distribution
! 2 oa ou ob o u o2 u 1 þ ka a q þ j2 ewE0 f ðtÞ ¼ l 1 þ hb b þ ot oy2 oz2 ot ot
ð20Þ
Non-dimensionalizing Eq. (20) using the below listed non-dimensional parameters
¼ u
u ; us
t ¼
l l l ew E0 t; k ¼ k; h ¼ h; us ¼ w ; K ¼ jDh l qD2h qD2h qD2h
ð21Þ
we can show that the non-dimensional counterpart for Eq. (20) takes the form
! 2 o2 u oa ou ob o u 2 ; zÞf ðtÞ ¼ 1 þ hb : 1 þ ka a P ð y þ K 2 oz2 oy ot ot b ot
ð22Þ
Also, it is straightforward that the corresponding dimensionless initial and boundary conditions are
ou ujt¼0 ¼ ¼ 0; ot t¼0 jy¼H=D ¼ u jz¼W=D ¼ 0; u h h ou ou ¼ ¼ 0: y¼0 oz z¼0 oy
ð23aÞ ð23bÞ ð23cÞ
According to the fast finite Fourier transform (FFT) method [34], the exact solution of velocity determined by Eqs. (22) and (23) can be expressed in the form of an twofold infinite series
ðy ; z; tÞ ¼ u
1 X 1 X
ÞZ m ðzÞ; T nm ðtÞY n ðy
ð24Þ
n¼0 m¼0
where two orthogonal basis functions are, respectively,
rffiffiffiffiffiffiffiffiffi 2Dh ; cos mn y H rffiffiffiffiffiffiffiffiffi 2Dh Z m ðzÞ ¼ cos cm z: W Þ ¼ Y n ðy
ð25aÞ ð25bÞ
Here, the eigenvalues mn and cm are the same as given by Eq. (13). By using the basis functions, the original fractional partial differential equation (22) and its initial and boundary conditions Eq. (23) can be modified to a fractional ordinary equation with respect to Tnm(t) by using finite Fourier transform
! a a b dT nm 2 a d a d 2 2 b d K Anm 1 þ k a f ðtÞ ¼ ðmn þ cm Þ 1 þ h b T nm 1 þ k a dt dt dt dt
ð26Þ
506
C. Zhao, C. Yang / Applied Mathematics and Computation 211 (2009) 502–509
with a simple initial boundary condition
T nm jt¼0
dT nm ¼ ¼ 0; dt t¼0
ð27Þ
where
2Dh ð1Þnþm mn cm : Anm ¼ pffiffiffiffiffiffiffiffiffi 2 þ HW K þ m2n þ c2m cm mn
ð28Þ
Subsequently, the Laplace transform and inverse Laplace transforms
T snm ðsÞ ¼ L½T nm ðtÞ ¼
Z
þ1
T nm ðtÞest dt;
ð29aÞ
0
Z rþj1 1 T s ðsÞest ds T nm ðtÞ ¼ L1 ½T snm ðsÞ ¼ 2pj rj1 nm
ð29bÞ
are applied to solve Eq. (26), where r is a real number so that the contour path of integration is in the region of convergence of T snm ðsÞ normally requiring r > Re(sp) for every singularity sp of T snm ðsÞ and j2 = 1. Transforming both sides of Eq. (26) by taking into account the initial condition expressed by Eq. (27) and then we can get
ð1 þ ka sa ÞsT snm ðsÞ K 2 Anm ð1 þ ka sa ÞFðsÞ ¼ ðm2n þ c2m Þð1 þ hb sb ÞT snm ðsÞ:
ð30Þ
Hence, the Laplace transform of Tnm(t) can be shown to bear the form
T snm ðsÞ ¼
K 2 Anm ð1 þ ka sa ÞFðsÞ ; ð1 þ ka sa Þs þ ðm2n þ c2m Þð1 þ hb sb Þ
ð31Þ
where FðsÞ ¼ L½f ðtÞ. We can obtain the exact solution to the model if the inverse Laplace transform of T snm ðsÞ is known. At last, according to Eq. (24), the exact solution can be further shown as
ðy ; z; tÞ ¼ u
1 X 1 X
ÞZ m ðzÞ ¼ T nm ðtÞY n ðy
n¼0 m¼0
1 X 1 X
ÞZ m ðzÞ: L1 ½T snm ðsÞY n ðy
ð32Þ
n¼0 m¼0
4. Special cases Up to now, we only need to know Tnm(t) to construct exact solutions from Eq. (32). But for generalized Oldroyd-B fluids ( k – 0; h – 0Þ, it is found that Eq. (31) can not be inversely transformed into to the time domain analytically, one can resort to numerical methods, such as Durbin’s method [35], Stehfest’s method [36], etc. However, for some special types of generalized Oldroyd-B fluids, such as Newtonian fluids ð k¼ h ¼ 0Þ, Fractional Maxwell fluids ð k – 0; h ¼ 0Þ and generalized second grade fluids ð k ¼ 0; h – 0Þ, the analytical solution does exist and will be given below. In view of expressing formulae compactly, only results of Tnm(t) are provided. 4.1. Newtonian fluid model ð k¼ h ¼ 0Þ As mentioned in Section 2, the Newtonian fluid is just a special kind of generalized Oldroyd-B fluid with k ¼ 0 and h ¼ 0, from general formula Eq. (31), the Laplace transform of Tnm(t) is of the form
T snm ðsÞ ¼
K 2 Anm FðsÞ : s þ ðm2n þ c2m Þ
ð33Þ
Also, in this paper, two most widely used electric fields in practical applications are considered, which will be reflected in the different forms of f(t) and thus different F(s). Case I: the external electric field is applied and remains constant from the time t = 0 (i.e., the electric field follows stepchange):
f ðtÞ ¼ HðtÞ;
ð34Þ
where H(t) is the well-known Heaviside unit step function. Then, Eq. (33) becomes
T snm ðsÞ ¼
K 2 Anm s½s þ ðm2n þ c2m Þ
ð35Þ
507
C. Zhao, C. Yang / Applied Mathematics and Computation 211 (2009) 502–509
and its inverse Laplace transform reads
T nm ðtÞ ¼
K 2 Anm 2 2 ½1 eðmn þcm Þt : m2n þ c2m
ð36Þ
Case II: we have a unit impulsive electric field applied at t=0, which means that f(t) = d(t), where d(t) represents Dirac delta function. Upon knowing this, Eq. (33) reduces to
T snm ðsÞ ¼
K 2 Anm s þ ðm2n þ c2m Þ
ð37Þ
and the corresponding inverse Laplace transform is 2
2
T nm ðtÞ ¼ K 2 Anm eðmn þcm Þt :
ð38Þ
4.2. Fractional Maxwell fluid model ( h ¼ 0 and k – 0) With h ¼ 0 in Eq. (31), the Laplace transform of Tnm (t) for fractional Maxwell fluids takes the form
T snm ðsÞ ¼
K 2 Anm ð1 þ ka sa ÞFðsÞ : ð1 þ ka sa Þs þ ðm2n þ c2m Þ
ð39Þ
Just like in 4.1, we just consider two forms of electric field namely a step-change electric field and an impulsive electric field. In the rest of this work, without otherwise notifications, Case I stands for the case of step-change electric field and Case II represents the case of impulsive electric field. Case I: Eq. (39) is changed to
T snm ðsÞ ¼
K 2 Anm ð1 þ ka sa Þ Bnm 2 2 ka sa2 Bnm ; ¼ K A s þ nm 1 þ Bnm 1 þ Bnm s½ð1 þ ka sa Þs þ ðm2n þ c2m Þ
ð40Þ
where
Bnm ¼
ka s : s1þa þ ðm2n þ c2m Þka
ð41Þ
In order to avoid the burdensome calculation of residues and contour integrals in Eq. (29b), we will apply the discrete inverse Laplace transform to get Tnm(t). At first, one can find dk and make dk to satisfy when Re(s) > Re(dk), jBnmj < 1, then Eq. (40) can be expanded as
T snm ðsÞ ¼ K 2 Anm
1 X ð1Þk kaðkþ1Þ k¼0
sk1 2 Þ a kþ1 m k
½s1þa þ ðm2n þ c
þ K 2 Anm
1 X
ð1Þk kak
k¼0
sk1þa ½s1þa þ ðm2n þ c2m Þka kþ1
:
ð42Þ
From [16], we obtain an important Laplace transform of the Mittag–Leffler function
h i p L t pnþq1 EðnÞ p;q ðat Þ ¼
n!spq ðs aÞnþ1
;
ð43Þ
where Ep,q(z) is the generalized Mittag–Leffler function and its derivative of order n is expressed by n
EðnÞ p;q ðzÞ ¼
1 X d ðj þ nÞ!zj : n Ep;q ðzÞ ¼ j!Cðpj þ pn þ qÞ dz j¼0
ð44Þ
The term-by-term inversion of Eq. (42), based on the general expansion theorem for the Laplace transform, using Eq. (43), it can be shown that the solution in the time domain is
(
1 X
ð1Þk aðkþ1Þaðkþ1Þþ1 ðkÞ t E1þa;kþaþ2 ðm2n þ c2m Þkat1þa k k! k¼0 ) 1 X
ð1Þk akakþ1 ðkÞ þ E1þa;kþ2 ðm2n þ c2m Þkat1þa : k t k! k¼0
T nm ðtÞ ¼ K 2 Anm
ð45Þ
Case II: for impulsive electric field, following the similar procedure in Case I, the results are consecutively shown below:
T snm ðsÞ ¼
K 2 Anm ð1 þ ka sa Þ Bnm 2 1 ka sa1 Bnm ¼ K A s þ ; nm 1 þ Bnm 1 þ Bnm ð1 þ ka sa Þs þ ðm2n þ c2m Þ
ð46Þ
508
C. Zhao, C. Yang / Applied Mathematics and Computation 211 (2009) 502–509
where Bnm has been given by Eq. (41). The Taylor expansion produces 1 X ð1Þk kaðkþ1Þ
T snm ðsÞ ¼ K 2 Anm
k¼0
sk ½s1þa þ ðm2n þ c
2 Þ a kþ1 m k
þ K 2 Anm
1 X ð1Þk kak k¼0
skþa ½s1þa þ ðm2n þ c2m Þka kþ1
:
ð47Þ
Finally, the inverse Laplace transform of Eq. (47) result in 2
T nm ðtÞ ¼ K Anm
( 1 X ð1Þk k¼0
) 1
X
ð1Þk akak ðkÞ 2 2 a1þa 2 2 a1þa kaðkþ1Þt aðkþ1Þ EðkÞ : t t t E1þa;kþ1 ðmn þ cm Þk k þ 1þa;kþ1þa ðmn þ cm Þk k! k! k¼0 ð48Þ
4.3. Generalized second grade fluid model ( k ¼ 0 and h – 0) For the generalized second grade fluids, it is straightforward that Eq. (31) reduces to
T snm ðsÞ ¼
K 2 Anm FðsÞ : s þ ðm þ c2m Þð1 þ hb sb Þ
ð49Þ
2 n
For both Case I and Case II, implementing similar procedure in 4.2, we can show the results: Case I: the Laplace transform of the time-dependent part Tnm (t) in Eq. (24) is expressed as
T snm ðsÞ ¼
K 2 Anm K 2 Anm 1 C nm ¼ 2 s ; 2 2 b b m 1 þ C nm þ cm Þð1 þ h s Þ n þ cm
s½s þ ðm
2 n
ð50Þ
where
C nm ¼
ðm2n þ c2m Þsb : þ ðm2n þ c2m Þhb
ð51Þ
s1b
Expanding Eq. (51) in Taylor series form as 1 X ðm2n þ c2m Þk sbðkþ1Þ1 ð1Þk : ½s1b þ ðm2n þ c2m Þhb kþ1 k¼0
T snm ðsÞ ¼ K 2 Anm
ð52Þ
Then, using discrete Laplace transform of the sequential fractional derivatives, the solution in the time domain is shown as
T nm ðtÞ ¼ K 2 Anm
1 X ð1Þk 2 ðkÞ ðmn þ c2m Þktkþ1 E1b;2þkb ½ðm2n þ c2m Þhbt1b : k! k¼0
ð53Þ
Case II:
T snm ðsÞ ¼
K 2 Anm sþ
ðk2n
þ c2m Þð1 þ hb sb Þ
ð54Þ
:
Eq. (55) can be expanded as
T snm ðsÞ ¼ K 2 Anm
1 X ð1Þk k¼0
ðm2n þ c2m Þk sbðkþ1Þ : þ ðm2 þ c2 Þhb kþ1
½s1b
n
ð55Þ
m
Using discrete inverse Laplace transform, we can get Tnm(t) as
T nm ðtÞ ¼ K 2 Anm
1 X
ð1Þk 2 ðkÞ ðmn þ c2m Þkt k E1b;bkþ1 ðm2n þ c2m Þhbt 1b : k! k¼0
ð56Þ
5. Concluding remark The generalized Oldroyd-B constitutive model with fractional calculus is introduced to describe the viscoelastic behavior of biological fluids in BioMEMS and Lab-on-a-chip devices. It is noted that the generalized Oldroyd-B model has extensive adaptability as it is a generalization of Newtonian model, fractional Maxwell model and fractional second grade model. At first, we use finite Fourier transform method to construct a general twofold infinite series solution for the velocity distribution of electro-osmotic flow of viscoelastic fluids in a rectangular micro-channel and however the analytical solution is not completed known unless we can obtain the inverse Laplace transform of the time-dependent component T snm ðsÞ. Then, for three special generalized Oldroyd-B models (Newtonian model, fractional Maxwell model and generalized second grade
C. Zhao, C. Yang / Applied Mathematics and Computation 211 (2009) 502–509
509
model), explicit analytical solutions of velocity distribution are obtained by using the discrete Laplace transform of the sequential fractional derivatives and the generalized Mittag–Leffler function. These exaction solutions are of potential importance because it can describe the flow behavior of electro-osmotic flow of biological fluids in Lab-on-a-chip devices. Moreover, they may also serve as accuracy checks for experimental, numerical and asymptotic methods. Acknowledgements The Ph.D. scholarship awarded to Z.C.L. from Nanyang Technological University is sincerely appreciated. The authors also would like to thank the referees for their useful comments and suggestions regarding an earlier version of this paper. References [1] J.H. Masliyah, S. Bhattacharjee, Electrokinetic and Colloid Transport Phenomena, Wiley Interscience, Hoboken, NJ, 2006. [2] A. Wainright, U.T. Nguyen, T. Bjornson, T.D. Boone, Preconcentration and separation of double-stranded DNA fragments by electrophoresis in plastic microfluidic devices, Electrophoresis 24 (2003) 3784–3792. [3] G.M. Whitesides, The origins and the future of microfluidics, Nature 442 (2006) 368–373. [4] A.E. Herr, J.I. Molho, J.G. Santiago, M.G. Mungal, T.W. Kenny, M.G. Garguilo, Electroosmotic capillary flow with nonuniform zeta potential, Anal. Chem. 72 (2000) 1053–1057. [5] C.L. Rice, R. Whitehead, Electrokinetic flow in a narrow cylindrical capillary, J. Phys. Chem. 69 (1965) 4017–4024. [6] J.G. Santiago, Electroosmotic flows in microchannels with finite inertial and pressure forces, Anal. Chem. 73 (2001) 2353–2365. [7] R.J. Yang, L.M. Fu, Y.C. Lin, Electroosmotic flow in microchannels, J. Colloid Interface Sci. 239 (2001) 98–105. [8] S. Devasenathipathy, J.G. Santiago, K. Takehara, Particle tracking techniques for electrokinetic microchannel flows, Anal. Chem. 74 (2002) 3704–3713. [9] M.J. Pikal, S. Shah, Transport mechanisms in iontophoresis: III. An experimental study of the contributions of electroosmotic flow and permeability change in transport of low and high molecular weight solutes, Pharm. Res. 7 (1990) 222–229. [10] D. Ross, T.J. Johnson, L.E. Locascio, Imaging of electroosmotic flow in plastic microchannels, Anal. Chem. 73 (2001) 2509–2515. [11] D. Sinton, D. Li, Electroosmotic velocity profiles in microchannels, Colloids Surf., A 222 (2003) 273–283. [12] P. Dutta, A. Beskok, T.C. Warburton, Numerical simulation of mixed electroosmotic/pressure driven microflows, Numer. Heat Transfer, Part A 41 (2002) 131–148. [13] Y. Hu, C. Werner, D. Li, Electrokinetic transport through rough microchannels, Anal. Chem. 75 (2003) 5747–5758. [14] N.A. Patankar, H.H. Hu, Numerical simulation of electroosmotic flow, Anal. Chem. 70 (1998) 1870–1881. [15] L. Ren, D. Li, Electroosmotic flow in heterogeneous microchannels, J. Colloid Interface Sci. 243 (2001) 255–261. [16] I. Podlubny, Fractional Differential Equations, Academic Press, San Diego, 1999. [17] R. Hilfer, Applications of Fractional Calculus in Physics, World Scientific Press, Singapore, 2000. [18] W. Tan, F. Xian, L. Wei, An exact solution of unsteady Couette flow of generalized second grade fluid, Chin. Sci. Bull. 47 (2002) 1783–1785. [19] W. Tan, M. Xu, Plane surface suddenly set in motion in a viscoelastic fluid with fractional Maxwell model, Acta Mech. Sin./Lixue Xuebao 18 (2002) 342– 349. [20] W. Tan, W. Pan, M. Xu, A note on unsteady flows of a viscoelastic fluid with the fractional Maxwell model between two parallel plates, Int. J. Non Linear Mech. 38 (2003) 645–650. [21] T. Hayat, S. Nadeem, S. Asghar, Periodic unidirectional flows of a viscoelastic fluid with the fractional Maxwell model, Appl. Math. Comput. 151 (2004) 153–161. [22] R.L. Magin, Fractional calculus in bioengineering, part 3, Crit. Rev. Biomed. Eng. 32 (2004) 195–377. [23] S. Hyder Ali, M. Khan, H. Qi, Exact solutions for a viscoelastic fluid with the generalized Oldroyd-B model, Nonlinear Anal.: Real World Applications (2008), doi:10.1016/j.nonrwa.2008.03.012. [24] M. Khan, S. Hyder Ali, H. Qi, Some accelerated flows for a generalized Oldroyd-B fluid, Nonlinear Anal.: Real World Applications (2007), doi:10.1016/ j.nonrwa.2007.11.017. [25] M. Khan, T. Hayat, S. Asghar, Exact solution for MHD flow of a generalized Oldroyd-B fluid with modified Darcy’s law, Int. J. Eng. Sci. 44 (2006) 333–339. [26] D. Vieru, C. Fetecau, C. Fetecau, Flow of a viscoelastic fluid with the fractional Maxwell model between two side walls perpendicular to a plate, Appl. Math. Comput. 200 (2008) 459–464. [27] A. Hernandez-Jimenez, J. Hernandez-Santiago, A. Macias-Garca, J. Sanchez-Gonzalez, Relaxation modulus in PMMA and PTFE fitting by fractional Maxwell model, Polym. Test. 21 (2002) 325–331. [28] T. Hayat, M. Khan, S. Asghar, On the MHD flow of fractional generalized Burgers’ fluid with modified Darcy’s law, Acta Mech. Sin./Lixue Xuebao 23 (2007) 257–261. [29] M. Khan, K. Maqbool, T. Hayat, Influence of Hall current on the flows of a generalized Oldroyd-B fluid in a porous space, Acta Mech. 184 (2006) 1–13. [30] M. Khan, S. Nadeem, T. Hayat, A.M. Siddiqui, Unsteady motions of a generalized second-grade fluid, Math. Comput. Model. 41 (2005) 629–637. [31] D.D. Joseph, Fluid Dynamics of Viscoelastic Liquids, Springer-Verlag, New York, 1990. [32] L. Debnath, D. Bhatta, Integral Transforms and Their Applications, Taylor & Francis Group, New York, 2007. [33] D. Li, Electrokinetics in Microfluidics, Elsevier Academic Press, 2004. [34] W.M. Deen, Analysis of Transport Phenomena, Oxford University Press, New York, 1998. [35] F. Durbin, Numerical inversion of Laplace transforms: an efficient improvement to Dubner and Abate’s method, Comput. J. 17 (1974) 371–376. [36] H. Stehfest, Algorithm 368: numerical inversion of Laplace transforms [D5], Commun. ACM 13 (1970) 47–49.