- Email: [email protected]
Startup electroosmotic flow of a viscoelastic fluid characterized by Oldroyd-B model in a rectangular microchannel with symmetric and asymmetric wall zeta potentials
Accepted Manuscript
Startup electroosmotic flow of a viscoelastic fluid characterized by Oldroyd-B model in a rectangular microchannel with symmetric and asymmetric wall zeta potentials P Kaushik , Suman Chakraborty PII: DOI: Reference:
S0377-0257(16)30113-6 10.1016/j.jnnfm.2017.06.003 JNNFM 3902
To appear in:
Journal of Non-Newtonian Fluid Mechanics
Received date: Revised date: Accepted date:
28 July 2016 9 June 2017 10 June 2017
Please cite this article as: P Kaushik , Suman Chakraborty , Startup electroosmotic flow of a viscoelastic fluid characterized by Oldroyd-B model in a rectangular microchannel with symmetric and asymmetric wall zeta potentials, Journal of Non-Newtonian Fluid Mechanics (2017), doi: 10.1016/j.jnnfm.2017.06.003
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT
AC
CE
PT
ED
M
AN US
CR IP T
Highlights Transient electroosmotic flow of a viscoelastic fluid in a rectangular microchannel Rheology of viscoelastic fluid is described by Oldroyd-B fluid model Complex nature of the fluid leads to enhancement of flowrates for symmetric zeta potentials Complex velocity profiles are observed for asymmetric zeta potentials
1
ACCEPTED MANUSCRIPT
Startup electroosmotic flow of a viscoelastic fluid characterized by Oldroyd-B model in a rectangular microchannel with symmetric and asymmetric wall zeta potentials P Kaushik and Suman Chakraborty* Department of Mechanical Engineering,
CR IP T
Indian Institute of Technology Kharagpur, Kharagpur, 721302, India
AN US
___________________________________ *Corresponding Author Email address: [email protected]
ABSTRACT
M
We study the startup startup flow of a viscoelastic fluid in a microfluidic channel, by considering the Oldroyd-B constitutive model. In particular, we aim to capture the confluence
ED
of the interfacial charge and the combined viscous and elastic properties of the fluid towards dictating the resulting startup flow characteristics. We further highlight the implications in the
PT
asymmetric nature of the wall zeta potentials on the flow behavior. We observe an increase in flow rate with increase in elastic nature of the fluid for studies with symmetric wall zeta
CE
potentials. For cases with asymmetric zeta potentials, we observe large movement of the velocity peaks, creating crests and troughs. Results from the present study may turn out to be useful in designing lab on a chip based microfluidic devices transporting bio-fluids such as
AC
blood, saliva and mucus by means of electroosmotic mechanisms.
Keywords: Viscoelastic Fluid; Rectangular microchannel; Startup flow; Transient analysis; Oldroyd-B model I. Introduction Recent advances in microfluidic applications, especially in Lab-on-chip devices, have triggered interests in exploring the electroosmotic method of flow activation. Lab-on-chip devices can be used for the analysis of chemical and biological samples which find many 2
ACCEPTED MANUSCRIPT applications, especially in medicine [1,2]. Many of the bio-fluids usually involved in such applications, such as blood, saliva, mucus, etc., are non-Newtonian in nature [3–5]. Implications of studying electrosmotically actuated flows of viscoelastic fluids may therefore be far-reaching, ranging from the development of medical devices with artificial flow actuation to portable diagnostic kits. The basic idea in electroosmotic flow is that the existence of a small potential at the wall called zeta potential causes the ions of the liquid to move. The ions charged oppositely
CR IP T
to the surface charge tend to move close to the wall surface, creating a thin layer (in which ionic species of charge of sign opposite to that of the wall are more abundant than ionic species of charge of same sign as that of the wall), so as neutralize the surface charge. The charged surface, together with the charged interfacial layer, is called the electrical double layer (EDL). The EDL formation creates an imbalance of charge in the fluid, which can drive
AN US
a fluid flow by exploiting the viscous drag between the ions and the fluid, on applying an external electrical potential gradient [6]. It is well known that the lack of moving parts makes it easier for fabrication and pumping of a fluid in microchannels by this method. Further control over fluid motion in electroosmosis can also be well accomplished by methods such
M
as surface patterning [7]. Velocity profiles in a typical electroosmotic flow with uniform wall zeta potential are plug-like, unlike pressure driven flows where the velocity profiles are
ED
usually parabolic. The plug-like velocity profiles can be used for separation of charged species such as DNA [8]. This tailored zeta potential at the walls makes controllability of the flow much better within the microchannels. And thus makes the study of electroosmotic flow
PT
with tailored within microchannels very important. Therefore, we will further look at the important developments in the study of electroosmotic flows existing in literature.
CE
Electroosmotic flow has been extensively studied for Newtonian fluids via theory [9– 11], experiments [12,13] and numerical methods [14–16]. However, the startup
AC
characteristics of electroosmotic flows of complex fluids have not been well explored. Some recent works on startup flows of viscoelastic fuids include studies in idealized geometries such as parallel plate channel [17], microtubes [18,19] and rotating channels [20,21], where Maxwell and Oldroyd-B models have been considered. Zhao and Yang [22] could find the exact solution to some special cases of generalized Oldroyd-B fluid which includes results for a Maxwell fluid and a generalized second grade fluid for small symmetric wall zeta potentials. The authors in their study also mentioned that full analytical solution may not be possible for Oldroyd-B fluid. Additionally the above studies have considered only symmetric zeta potentials at the walls. It is important to note that surface modifications can be used to 3
ACCEPTED MANUSCRIPT vary the wall zeta potentials and can be used for better control of the fluid [23,24]. Taking the possibility of asymmetry in zeta potential into consideration, some studies have been reported for Newtonian fluids [25], Phan-Thien-Tanner(PTT) fluids [26] and Maxwell fluids [27] between two parallel plates. Experimental study with Newtonian fluids in a channel with asymmetric zeta potentials was conducted and reported by Miller et al. [25]. In a more recent study, Jiménez et al. [28,29] studied the startup electroosmotic flow of Maxwell fluid in a rectangular microchannel with high assymetric zeta potentials by a semi-analytical method.
CR IP T
In this present study, we attempt to numerically analyze the electroosmotically driven transport of Oldroyd-B fluids in a rectangular microfluidic channel. We also compare and validate the results from our numerical method with Jiménez et al. [28,29] for a special case in which the retardation time is set to zero, for which the fluid behaves asymptotically like an idealized Maxwell fluid. First we focus on the cases where there is symmetry in the wall zeta
AN US
potentials; here we also compare how the complex nature of the fluid affects the flow rates. Then we focus on the hydrodynamics involved when there is asymmetry in the wall zeta potentials. The results are then concluded.
M
II. Problem description and mathematical formulation.
ED
A. Physical model
Here, we consider the transient startup flow of viscoelastic fluid characterized by Oldroyd-B
PT
constitutive relation in a rectangular microfluidic channel as shown in Fig. 1. The coordinate axes are chosen to be x, y and z , which are along the length, width and height of the channel
CE
respectively. The following assumptions have been considered for the present study: (i) The length ( L ) of the channel is much larger than its width and height. (ii) The fluid is assumed to be initially at rest and unstressed.
AC
(iii) The channel cross section is taken to be of height 2H and width 2W and the origin is placed at the channel centerline, as shown in Fig 1.
(iv) The buffer solution is a symmetric
z : z
electrolyte with constant physical
properties. (v) The wall zeta potentials are taken as shown in Fig. 1.
u H (vi) The flow Reynolds number Re c
is assumed to be very small for
4
ACCEPTED MANUSCRIPT electroosmotic flows Re
1 [28].
(vii)
The flow is considered to be unidirectional.
(viii)
No external pressure gradient is applied.
k BTref Ex Here, uc is the characteristic velocity scale, is the permittivity of z e the fluid, k B is the Boltzmann constant, Ex is the externally applied electric field parallel
CR IP T
to the flow direction. Tref is the reference temperature of the fluid (temperature at which the fluid is flowing), z is the valence of the electrolyte, e is the fundamental electron
AN US
charge and is the sum of base fluid and polymer viscosities.
M
Figure 1. (Color online) Schematic depicting the problem considered in the present study.
ED
B. Governing equations of the flow field
The governing equations of the flow dynamics in the channel considered are given by
CE
PT
[22,30,31]:
u 0
u u u p T e E t
(1) (2)
AC
In Eq. (2), is the fluid density, T is the stress tensor, e is the ion density, E is the total electric field and u is the velocity vector. Following the Oldroyd-B constitutive model for the stress tensor equation, we get, [32]: T 1 T 2 D 2 D
(3)
where, 1 and 2 are the fluid relaxation and retardation times respectively and the upper convected derivative of T is defined as [33]:
T
T T u T u T T u t
(4)
5
ACCEPTED MANUSCRIPT Here, D is defined as
D
1 u uT 2
(5)
In order to simplify the above equation, we split T as T T1 T2 , where
T2 2
2 D 1
(6)
Therefore, the momentum equation further reduces to:
where,
AN US
T1 1 T1 2 1 2 D 1
CR IP T
u u u p T1 2 2 D e E t 1
(7)
(8)
The components of T1 matrix appearing in Eq. (7) are as follows: Txx T1 Tyx Tzx
Txy Tyy Tzy
Txz Tyz Tzz
(9)
M
As a consequence of the approximately unidirectional flow assumption, we try to find
ED
the flow solution for u t , y, z at t 0 assuming v w 0 . In the subsequent section, we show the method of resolving the electroosmotic body force term in Eq. (7).
PT
C. Resolution of electroosmotic body force term It is important to mention beforehand that the time scale related to the ionic transport effects
CE
within the EDL is very short and of the order
107 108 s , which is at least 2-3 orders
smaller than characteristic time associated with the electroosmotic flow [34]. Hence the
AC
transient term in the equation that governs the electrical body force term e E in Eq. (7) can be neglected and we can consider the Poisson-Boltzmann equation as the equation which governs the electrical potential distribution within the EDL, as [20,35,36]: e z 2 2 2 zen Sinh y 2 z 2 k BTref
(10)
Here, e is the total ionic charge density, which for a z : z symmetric electrolyte
6
ACCEPTED MANUSCRIPT e z using Boltzmann distribution is given by e 2 zen Sinh k T B ref
and is the inverse of the
EDL thickness. Here, z is the valence of the ions, n is the concentration of the bulk and
Tref is the reference temperature (ambient temperature of the fluid).
D. Simplification and Non-dimensionalization of the governing equations
y*
CR IP T
The dimensionless parameters used are given as follows [28,36]:
e z Wi y z u TH t , z* , * , u* , t * , T* , El 1 2 , 2, 2 k BTref H H uc uc 1 H Re H
* kH , b* W / H
(11)
Where, El is the elasticity number which considers the effect of fluid elasticity and inertia. It is
AN US
defined as the ratio of the relaxation time scale of the fluid to the diffusion time scale [37]. The ratio of the fluid retardation to relaxation time is represented by . We would like to mention here that since the velocities in the y and z directions are much smaller, we neglect them and consider our flow to be approximately unidirectional. This then reduces the equations of stress in the Oldroyd-B model [See Eq. (8)] to their simple form given below in Eqs. (14) and (15). The simplified governing
M
equations are given as:
ED
2 * 2 * *2 *2 Sinh * *2 y z
AC
CE
PT
* 2u* 2u* *2 u* Txy Txz* * *2 *2 Sinh t * y* z * y z
Txy*
(13)
u* y*
(14)
Txz* u* Txz El * 1 * t z
(15)
* * u* Txx* * u 2 Wi T T xy xz t * y* z *
(16)
Txy* El
t *
1
*
Txx* El
(12)
In order to simplify the equations further we combine equations (13)-(15) and we get,
El
2u* u* 2u* 2u* *2 1 El *2 Sinh * *2 * * *2 t t t y z
(17)
In the present study, equations (12) and (17) have been solved numerically along with
7
ACCEPTED MANUSCRIPT no-slip condition at the four walls, to obtain the solution till steady state is reached. As assumed before, the fluid is considered to be free of any stresses at the start of the flow with
u* zero velocity. Mathematically the initial conditions used are u 0, y , z * 0, y* , z* 0 t *
*
*
and boundary conditions used are u* t* , y* , 1 u* t* , b* , z* 0 . It is worth noting that for a special case of 0 , the fluid model is reduced to the Maxwell model and the equation
El
CR IP T
(17) reduces to
2u* u* 2u* 2u* * *2 *2 *2 Sinh * *2 t t y z
(18)
From here on, we will drop the superscript „*‟ in the equations and non-dimensional
AN US
terms during further discussion of the problem for the convenience of representation. III. Numerical Method and Model Benchmarking
A finite difference code has been developed for solving the transport equations along with their respective boundary conditions as given in the preceding section. The Poisson-
M
Boltzmann equation (12) is solved by a simple second order in space discretization using SOR method. The equation (17) is hyperbolic and the code has been developed based on
ED
numerical method suggested by Dehghan and Shokri [38] for a hyperbolic telegraph equation. The method uses a second order scheme for both the temporal discretization as well as the discretization of the spatial terms. The discretized equation is as shown below:
u t t , y, z 2u t , y , z u t t , y , z u t t , y , z u t t , y , z 2 2 t t
PT
El
(19)
AC
CE
El 2 El 2 2 = u t t , y, z u t t , y, z 1 u t , y, z 2 t 2 t 2 Sinh y, z
The value of chosen is 0.5 as suggested by Dehghan and Shokri [38], which is the
well known Crank-Nicholson scheme [39] to track the steady shear diffusion term. For the solution, we use the well-known Gauss-Seidel method and the code is run till the steady state is reached. Validation study – In order to ensure that the results obtained from present numerical method
8
ACCEPTED MANUSCRIPT are correct, we try to validate our method with the results available in the literature. Semianalytical solution for the electroosmotic flow of a Maxwell fluid in a rectangular channel is available in literature [28]; using this we validate our model for special case of 0 as shown in equation (18). Fig 2(a) shows depicts the benchmarking of our numerical code. The circular markers are used to represent the results reported by Jiménez et al. [28], while solid line represents the results from the present method. There is a good match obtained between the present study
CR IP T
and reported results. We have also benchmarked our code with semi-analytical solution of Zhao and Yang [22] for a generalized Oldroyd-B fluid in electroosmotic flow with symmetric zeta potential case considering Boltzmann distribution at the center of the channel (0,0) as depicted in Fig. 2(b). We observe very close match with the results reported in the literature.
PT
ED
M
AN US
We have mentioned other parameters used in plotting the figure in the caption.
(a)
(b)
CE
Fig 2 (color online) Model Validation: (a) Plot showing the temporal variation of dimensionless velocity u t, y, z at two different locations which were studied by Jiménez et al. [28]. Location A represents
AC
y, z 3.12, 0.91
y, z 2.75, 0.81
and location B represents
for an aspect ratio of W : H 2 :1 i.e. b 2 . The other parameters
considered here are: 10 , 1 0.2 , 2 0.2 , 3 0.2 , 4 0.2 and El 0.5 . (b) Plot showing the temporal variation of dimensionless velocity u t, 0, 0 at channel center which was studied by Zhao and Yang [22] considering Boltzmann distribution. The other parameters considered here are: 10 , 1 2 3 4 1 , 0.1 and El 2 .
9
CR IP T
ACCEPTED MANUSCRIPT
(b)
AN US
(a)
M
(c) Fig 3 (color online) Plot showing (a) the variation of dimensionless velocity u t, y, z vs z
ED
at y 0 with various grid size and (b) the variation of dimensionless velocity u t, y, z vs z at y 0 with various time step size (c) the variation of normal stress Txx y, z vs z at y 0
CE
PT
with various grid size. Here we choose the following parameters: 10 , b 1 , 0.1 and El 2.0 , considering asymmetric zeta potentials at walls 1 1 , 2 1 , 3 1 , 4 1 (symmetric zeta potentials at walls in the inset 1 2 3 4 1 ). Grid independence and convergence study – To make sure that the results obtained in the study are independent of the grid size and time step size, we present the grid and time step
AC
independence studies as well. We depict, in Figs. 3(a)-(b), the variation of axial velocity
u t, y, z as a function of z for a case with asymmetric wall zeta potentials (cases with
symmetric wall zeta potentials are shown in the in inset) for different grid resolutions
y z and time step t variations. It is worth noting from the figures that the variation of velocity becomes insignificant as changes from 1/96 to 1/128 and while a change in t from 102 to 103 does not show any appreciable change in velocity as well. This is true for both the symmetric and asymmetric wall zeta potential cases. We also show
10
ACCEPTED MANUSCRIPT
the variation of normal stress Txx for different grid resolutions (cases with symmetric wall zeta potentials are shown in the in inset). Similar to the variation of axial velocity the variation of Txx becomes insignificant as changes from 1/96 to 1/128. Accordingly, we have considered t 103 and 1 96 throughout this study. In order to understand the important findings of the present study, we would like to show the typical values of different model parameters chosen for the numerical computation. 106 M
(the bulk concentration
CR IP T
A typical concentration of the ionic solution
n 6.023 1020 m3 ) [40,41] at temperature Tref 298K with Boltzmann constant
kB 1.38 1023 and the permittivity of the solvent 80 8.85 1012 we find that the EDL thickness 1 is ~300 nm . Here we assume the characteristic channel dimension as
AN US
~ 10 m . Using these values, we find the dimensionless EDL thickness as: * H 30 . The typical values of Elasticity number used in the literature are varies between 0 and 1 [20].
M
IV. Results and Discussions
1 [37] and usually
The results and discussions will be divided into two parts. In the first part we present the
ED
results from the studies of symmetric zeta potentials at the walls. In the second part we attempt to highlight the results for asymmetric zeta potentials at the walls.
PT
A. Studies with symmetric wall zeta potentials
CE
We will begin our discussion with the studies on the temporal variations on the velocity profiles when two fluid parameters namely El and are varied. This is depicted in figs. 4(a)-
AC
(d), where the other parameters chosen are mentioned in the figure caption. It can be seen from Figs. 4(a)-(d), that the velocities tend to increase in magnitude as the time progresses. For smaller values of and larger values of El, the elastic nature of the fluid becomes more pronounced. This can be attributed to the fact that fluid relaxation time increases with increase in El as well as decrease in . Larger relaxation time means that the fluid takes more time to settle down to a steady state value. It is evident from fig. 4(c) that for a case with larger value of El and smaller value of , the velocity field exhibits a spring-like behavior, with the fluid centerline velocity going to a very large value at an intermediate time
11
ACCEPTED MANUSCRIPT of t 1 and then going to a smaller value at t 2 and eventually settling down slowly into and intermediate value which is the steady state profile after a certain number of oscillations. Considering the opposite case as depicted in fig. 4(b), where the value of El is small and the value of is large, the viscous effect tends to dominate over the elastic effects and we do not see much spring-like behavior in the velocity profile. Elastic nature of viscoelastic fluids may be attributed to the stretching of the polymer chains within the solution. For the case of
CR IP T
1 , the fluid will tend to behave like a Newtonian fluid with no elastic behavior. Inherently from equation (17) it is evident that the steady state solution of the equation is same as the solution for a Newtonian fluid irrespective of the value of El and . This fact is clearly evident in the solution at steady state as well as depicted by t 1 in figs. 4(a) and (b) and by t 10 in figs. 4(c) and (d). It can also be observed from Figs. 4(a)-(d) that the time
ED
M
AN US
taken to reach steady state is strongly affected by the variation of El and .
(b)
AC
CE
PT
(a)
(c) (d) Fig 4 (color online) Plot showing the temporal variation of dimensionless velocity u t, y, z vs z at y 0 for (a) El 0.1 and 0.1 , (b) El 0.1 and 0.8 , (c) El 1 and 0.1 , (d) El 1 and 0.8 . Other parameters chosen are as follows: 10 , b 1 and
12
ACCEPTED MANUSCRIPT 1 2 3 4 1 . Before discussing the results of figs. 5(a)-(d), it is important to define the flow rate of the fluid in the cross-section. In study of flows in microchannels, flow rates tend to play an important role as they give the overall picture of how effective our forcing/pumping mechanism is. Augmentation of flow rates in microchannels is sought by various means. In our present study, we define the dimensionless flow rate as:
u t , y, z dzdy
W 1
CR IP T
Q t
W 1
(19)
The effect of various fluid, flow and geometric parameters on the flow rate has been considered in figs. 5(a)-(d). The other parameters which have been used in the plots have been described in the figure caption. The symmetry of the wall-zeta potentials has an aiding
AN US
effect on the flow. In fig 5(a), the effect of has been considered on the flow rate. It is clear from the figure that decreasing increases the flow rate initially, before reaching steady state eventually. The time taken by the fluid to reach steady state also drastically increases with decrease in . This may be attributed to the fact that the fluid retardation time is very
M
short compared to the relaxation time for smaller . This means that the fluid tends to relax slowly to any external forcing or changes in the nearby fluid. Also we see that for smaller
ED
value of the oscillatory behavior in the fluid becomes more prominent, with larger number oscillations in the flow rate magnitude, due to the dominance of elastic effects as compared to
PT
the case with larger where viscous effects dominate. Fig 5(b) shows the effect of El on the flow rate. Increase in El essentially signifies the increase in relaxation time of the fluid.
CE
Therefore, as evident from fig. 5(b), the increase in El effectively increases the relaxation time and therefore increases the elastic behavior of the fluid. It can also be noted that with
AC
increase in El, the time taken by the fluid to reach the steady state also drastically increases along with increase in the magnitude of oscillations. However, the number of oscillations tends to remain more or less same. In order to have large flow rates it may be desirable to use fluids with larger values of El and smaller values of . Fig. 5(c) depicts the effect of or the inverse of EDL thickness on the flow rate. It may be observed clearly that the effect of
on the flow rate is quite less and we may consider the flow rate as a weak function of . Clearly, the flow rate increases in magnitude with increase in , both during transience and at steady state. Increase in represents decrease in the EDL thickness, which leads to more
13
ACCEPTED MANUSCRIPT ion mobility in the bulk flow. This could be attributed as the reason for increase in flow rate with . Also, we observe for larger values of such as 30 and 50 , the flow rate becomes almost independent of . This goes to show that beyond a certain reduction in the EDL thickness, the mobility of ions cannot increase drastically. Fig. 5(d) describes the effect of aspect ratio of the channel b on the flow rate. It is quite evident that for larger aspect ratios, the wall effects become less prominent, however the total area available for flow increases and therefore the flow rate increases. However, it is important to note that when the
CR IP T
aspect ratio is doubled, the flow rate does not double. This may be attributed to the fact that in electroosmotic flows, the fluid forcing essentially occurs close to the walls, which makes the aspect ratio study very important. Having small aspect-ratios may enhance the forcing of the fluid close to the wall. This causes the flow rate not to have proportional increase with
PT
ED
M
AN US
increase in aspect rate as seen from fig. 5(d).
(b)
AC
CE
(a)
(c)
(d)
Fig 5 (color online) Plot showing the temporal variation of dimensionless flow rate Q t for (a) various taking El 2.0 , b 1 and 10 , (b) various El taking 0.1 , b 1 and 14
ACCEPTED MANUSCRIPT 10 , (c) various taking El 1.0 , b 1 and 0.1 , (d) ) various b taking El 1.0 , 10 and 0.1 . Other parameters chosen are: and 1 2 3 4 1 . In flows of viscoelastic fluids, it becomes important to consider how the flow velocity varies with location. This is addressed in fig 6(a)-(c), where the temporal variations of the velocities are shown at certain locations within the cross-section in case with symmetric wall zeta potentials. The other parameters used for creating the figure have been given in the
CR IP T
caption. From fig. 6(a), we observe that the magnitude of the velocity is the highest at (0,0) and lower closer to the walls (0.9,-0.9). We also see that the increase in velocity magnitude takes place closer to the walls (0.9,-0.9) earlier compared to the case with center of the channel (0,0). This may be attributed to the fact that in electroosmotic flow, the forcing essentially happens closer to the walls in the area of the EDL. We also observe that the elastic
AN US
behavior of the fluid is more prominent away from the walls (0.9,-0.9) and closer to the center (0,0). This may be attributed to the fact that the wall tends to aid the viscous effects within the fluid close to it. We also see from fig. 6(a) that the transient behavior of the fluid velocity at (0,0) and (0.5,-0.5) differs, however, the steady state value is almost the same indicating a plug-like flow at steady state with values same as that of a Newtonian fluid. In
M
fig. 6(b), we see the temporal variation of velocity with different at two different locations. It is observed that the oscillations during transience tends to die out faster in case larger
ED
values and reaches a common value at each location at steady state. It can also be seen that the value of the velocity is higher with more number of oscillations for smaller ,
PT
irrespective of the location. We also see that the time at which the peaks occur at the different locations does not vary much with change in . This may be attributed to the fact that for
CE
smaller , the retardation time is considered to be smaller, which causes the viscous effects to be weaker. Hence, we see a small delay in occurrence of the peak at smaller values of .
AC
This effect, however, is seen to be quite weak. Fig. 6(c) depicts the temporal variation of the velocity with El for two different locations. We again see that the larger value of velocity occurs as El increases and also we see the enhancement in the elastic behavior of the fluid. This is because the fluid relaxation time increases and it takes more time to adjust to any change in the nearby zones. This also causes the velocity peaks to strongly shift further downstream in time with larger El.
15
CR IP T
ACCEPTED MANUSCRIPT
(b)
M
AN US
(a)
(c) Fig 6 (color online) Plot showing the temporal variation of dimensionless velocity u t, y, z
PT
ED
for (a) various locations within the cross-section taking El 2.0 , 0.1 , b 1 and 10 , (b) various locations within the cross-section and varying while taking El 2.0 , b 1 and 10 , (c) various locations within the cross-section and varying El while taking 0.1 , b 1 and 10 . Other parameters chosen are: and 1 2 3 4 1 .
CE
The above analysis concludes the cases studied for symmetric wall zeta potentials. We now move on to cases where each of the walls has different zeta potentials, which may be
AC
either positive or negative. B. Studies with asymmetric wall zeta potentials For studies considering asymmetric wall zeta potentials, we start the discussion with temporal evolution of velocity profiles for a case where all zeta potentials are positive while being different, and another case with zeta potentials at certain walls are negative. Fig. 7(a) and (b) show the temporal evolution of the velocity profile at the middle of the y direction, along zdirection. The other parameters which have been chosen are mentioned in the figure caption.
16
ACCEPTED MANUSCRIPT The commonality in both fig. 7(a) and 7(b) is that the at early time the flow essentially is restricted close to the EDL zone and as time progresses this flow velocity gets diffused into the entire domain. The important observation in this case is that unlike Maxwell fluid, where the diffusive terms only start dominating at larger times, for the Oldroyd-B fluid model the diffusive terms play an important role throughout the flow time. This essentially means that the velocity profiles are quite essentially smooth for an Oldroyd-B fluid unlike a Maxwell fluid. This phenomenon is clearly observed in Fig. 4(a)-(d) as well as Fig 7(a)-(b). In fig.
CR IP T
7(a), we observe that despite the skewed nature of the wall zeta potentials, at time t 1.5 , the maximum value of u t , y, z occurs almost at the centerline i.e. at z 0 . On either side of t 1.5 , i.e. t 1 and t 2 , we see the maximum of u t , y, z occurring on either side of z 0 . As the fluid is pulled towards z 0 side at t 1 , the fluid relaxes by the maximum
AN US
moving across to z 0 side at t 2 . The fluid then eventually attains a steady state after some oscillation. In fig. 7(b), the zeta potentials across z direction are equal in magnitude and opposite in direction. Here we observe that the spring like nature causes the fluid to have
u t , y, z 0 at some locations in z 0 and u t , y, z 0 at some locations in z 0 at t 2 .
M
This is converse to what happens at steady state, where, u y, z 0 for z 0 and u y, z 0
AC
CE
PT
ED
for z 0 .
(a) (b) Fig 7 (color online) Plot showing the temporal variation of dimensionless velocity u t, y, z vs z at y 0 for (a) Asymmetric non-negative zeta potentials: 1 1.0 , 2 2.0 , 3 2.0 and 4 1.0 , (b) Asymmetric with negative zeta potentials: 1 2.0 , 2 2.0 , 3 2.0 and 4 2.0 . Other parameters chosen are as follows: El 2 , 0.1 , 10 and b 1 . Fig 8. shows the temporal variation of the dimensionless velocity distribution
17
ACCEPTED MANUSCRIPT u t, y, z for asymmetric wall zeta potentials. The parameters used for generating the figure are mentioned in the caption. In the figure, dark red colour indicates large positive value of
u t, y, z and dark blue colour represents large negative value of u t, y, z , while the green colour represents values close to zero. We observe from figs. 8(a) and (b) that for very early times the flow of fluid is essentially restricted close to walls near the EDL zone. As time progresses, the flow spreads into the complete domain by diffusion as seen in figs. 8(c)-(f). It
CR IP T
is also important to note that the velocity profile u t, y, z is essentially smooth unlike the case of Maxwell fluid, which has been studied in detail by Jiménez et al. [28]. This is due to the fact that the viscous effects in the flow are prominent even at small times in an Oldroyd-B fluid unlike Maxwell fluid. We see that during the intermediate part of the transience, the combination of elastic and viscous effects cause a smooth maximum which moves away from
AN US
the walls and closer to the center of the domain considered. At t 0.5 in fig. 8(c), where the maximum is near the bottom edge of z 1 , the shift in the maximum takes place towards the center at t 1.5 , after which the maximum again shifts back towards zone near the
z 1 wall, where, 4 4.0 as the flow reaches steady state. At steady state, as depicted in
M
figs. 8(e)-(h), the u t, y, z profile becomes independent of El and , which also is evident from figs. 7(a)-(b). Elastic behavior dominates viscous behavior in the case described by fig.
ED
8 as 0.1 is small and El chosen is large. The stretching and relaxing nature of the fluid is typical of dilute polymer solutions, where the polymer chains act similar to springs. The
PT
stretching causes the potential energy to get stored within the polymer chains and it is released when the fluid starts relaxing. Despite this, the viscous nature is quite evident
AC
CE
because of strong diffusion within the fluid leading to smooth velocity profiles.
18
ACCEPTED MANUSCRIPT
(b)
M ED
(g)
(e)
AN US
(d)
(c)
CR IP T
(a)
(h)
(f)
(i)
Fig 8 (color online) Plot showing the temporal variation of dimensionless velocity u t, y, z vs for asymmetric zeta potentials: 1 1.0 , 2 2.0 , 3 3.0 and 4 4.0 . Other parameters
PT
y, z
AC
CE
chosen are as follows: El 2 , 0.1 , 10 and b 1 . The times displayed are (a) t 0.1 , (b) t 0.2 , (c) t 0.5 , (d) t 1 , (e) t 1.5 , (f) t 2 , (g) t 4 , (h) t 10 and (i) t 20 .
19
ACCEPTED MANUSCRIPT
(b)
(c)
CR IP T
(a)
y, z
AN US
(d) (e) (f) Fig 9 (color online) Plot showing the temporal variation of dimensionless velocity u t, y, z vs for asymmetric zeta potentials: 1 1.0 , 2 2.0 , 3 2.0 and 4 1.0 . Other
M
parameters chosen are as follows: El 2 , b 1 , 10 . The times and chosen are (a) 0.1 and t 0.5 , (b) 0.1 and t 1 , (c) 0.1 and t 1.5 , (d) 0.8 and t 0.5 , (e) 0.8 and t 1 , (f) 0.8 and t 1.5 .
ED
In order to highlight the effect of on asymmetric wall zeta potentials, we consider
0.1 and 0.8 which are used to generate Fig 9, while keeping all other parameters
PT
same. The other parameters chosen are mentioned clearly in the figure caption. The top line of fig. 9, i.e. figs. 9(a)-(c) show the temporal variation of the dimensionless velocity
CE
distribution u t, y, z for 0.1 with times t 0.5, 1, 1.5 and the bottom line i.e. figs. 9(d)(f) show the temporal variation for 0.8 . In the figure, dark red colour indicates large
AC
positive value of u t, y, z and dark blue colour represents large negative value of u t, y, z , while green being an intermediate colour representing values close to zero. We observe from figs. 9(a)-(c) that the elastic nature of the fluid dominates over the viscous nature for 0.1 and vice versa for 0.8 , as depicted by 9(d)-(f). We also observe that the flow reaches steady state quickly for the case with 0.8 . This may be observed by seeing figs. 9(e) and (f), which are quite similar; this matches closely with the steady state case shown in fig. 8(i). For 0.8 case, the relaxation time is smaller and therefore the fluid relaxes quickly to any disturbance that may occur in the fluid and hence steady state is achieved more quickly than
20
ACCEPTED MANUSCRIPT in the case of 0.1 . Effect of El on the asymmetric wall zeta potentials is highlighted in fig. 10, where a small El value of El 0.1 and a large value of El 2 have been considered for study. Large value of El implies that the relaxation time of the fluid is large and it takes time to adjust to any changes in flow. The other parameters used in generating fig. 10 are displayed in the figure caption. Magnitudes of velocity represented by the various colours are as followed in
CR IP T
the previous figures. It is clearly evident from figs. 10(a)-(f) that for larger value of El in case of figs. 10(d)-(f), the elastic nature of the fluid is more prominent, as the fluid takes longer time to adjust to changes in the flow field. For very small value of El 0.1 , which is chosen in 10(a)-(c), the fluid behavior is close to that of a Newtonian fluid because of small relaxation times. The flow reaches steady state quickly for El 0.1 case as evident from
(b)
(c)
AC
CE
PT
(a)
ED
M
AN US
similarity in flow profiles shown in figs. 10(b) and (c).
(d)
(e)
(f)
Fig 10 (color online) Plot showing the temporal variation of dimensionless velocity u t, y, z vs
y, z
for asymmetric zeta potentials: 1 1.0 , 2 2.0 , 3 2.0 and 4 1.0 . Other
parameters chosen are as follows: 0.1 , b 1 , 10 . The times and El chosen are (a) El 0.1 and t 0.5 , (b) El 0.1and t 1 , (c) El 0.1 and t 1.5 , (d) El 2 and t 0.5 , (e) El 2 and t 1 , (f) El 2 and t 1.5 . The effect of aspect ratio has been shown in fig. 11. The top line i.e. figs. 11(a)-(c)
represents the case where aspect ratio defined by b is taken as b 1 and the bottom line i.e.
21
ACCEPTED MANUSCRIPT figs 11(d)-(f) represents the case with b 2 . The effect of aspect ratio becomes quite important for asymmetric wall zeta potentials because of inherent dimensional asymmetry as well. The other parameters used are displayed in the figure caption. We see from the two figures that the corner effect is more prominent for the b 2 case at early times, as displayed in figs. 11(a) and 11(d) for t 0.5 . We observe from fig. 11(d) that there is a maximum at the corner near (-2,1), whereas, such a maximum is not there in the case of fig 11(a). We also observe that in case of larger aspect ratio for b 2 , the amount of time required for the
CR IP T
momentum to diffuse into the fluid is more because of larger length scale in the y-direction. This causes a delay in the flow velocity profiles as compared to the case with smaller aspect ratio of b 1 . This fact is evident from fig. 11(f), where we see that the flow is still in the main transience stage, where as in case of fig. 11(c), the flow has started to approach steady state. It can be concluded that larger aspect ratio case takes longer to reach steady state as
(b)
(c)
CE
PT
(a)
ED
M
AN US
compared to the case with smaller aspect ratio.
(e)
(f)
AC
(d)
Fig 11 (color online) Plot showing the temporal variation of dimensionless velocity u t, y, z vs
y, z
for asymmetric zeta potentials: 1 1.0 , 2 2.0 , 3 2.0 and 4 1.0 . Other
parameters chosen are as follows: El 2 , 0.1 , 10 . The times and aspect ratio chosen are (a) b 1 and t 0.5 , (b) b 1 and t 1 , (c) b 1 and t 1.5 , (d) b 2 and t 0.5 , (e) b 2 and t 1 , (f) b 2 and t 1.5 . In electroosmotic flows of viscoelastic fluids, it becomes important to consider how
the flow velocity varies with location, especially in the case with asymmetric wall zeta
22
ACCEPTED MANUSCRIPT potentials. Hence, we plot fig 12(a)-(c) with the intention of showing how the flow velocity
u t, y, z varies with at certain locations within the domain. All other parameters used are given clearly in the figure caption. We essentially observe from fig. 12 that the velocity variation with time is smooth for all locations. This is quite unlike fig. 2, where Jiménez et al. [28], studied a similar case for a Maxwell fluid. In case of Maxwell fluid, due to weak diffusive effects at small times, the velocity profile tends to be coarser with sharp oscillation.
CR IP T
Here, for an Oldroyd-B fluid, we have strong diffusive nature right from small times which makes the oscillations smoother. Also, the strong diffusive effects tend to damp and kill the oscillation faster. Fig 12(a), shows how the flow velocity varies for a case where equal and opposite zeta potentials are taken across from each other. We observe that at the center of the channel, there is almost no velocity developed through all times. We also observe that close
AN US
to the walls i.e. (0.9,-0.9), the flow velocity magnitude is lower, with almost no oscillatory effects. The flow velocity magnitude also virtually reaches the steady state value at a quicker time. This may be attributed to the viscous effects of the nearby wall. The flow velocity magnitude midway between the walls and the center i.e. at (0.5,-0.5) is large with more oscillatory nature and takes larger time to reach steady state. In fig 12(b), we see the temporal
M
variation at two different locations with varying . We observe from the figure that the magnitude of the velocity increases with decrease in , irrespective of the location.
ED
However, the damping of the oscillatory nature is more pronounced at (0.9,-0.9) as compared to at (0.5,-0.5). This may be due to the fact that (0.9,-0.9) is closer to the wall, where
PT
diffusive effects are stronger. Also, for larger values of , oscillations are damped faster because of strong diffusive nature of the fluid as compared to the elastic nature. The effect of
CE
El on the velocity magnitudes at two different locations has been considered in fig. 12(c). Here we observe that the velocity magnitude increases at all locations with increase in El, and
AC
also that the oscillatory nature becomes more damped closer to the wall at (0.5,-0.5) as compared to away from the wall at (0.9,-0.9). We also observe that the time taken for the first peak to appear and time taken to reach steady state also increases with increase in El, irrespective of the location.
23
CR IP T
ACCEPTED MANUSCRIPT
(b)
M
AN US
(a)
ED
(c) Fig 12 (color online) Plot showing the temporal variation of dimensionless velocity u t, y, z
CE
3 4 2.0 .
PT
for (a) various locations within the cross-section taking El 2.0 , 0.1 , b 1 and 10 , (b) various locations within the cross-section and varying while taking El 2.0 , b 1 and 10 , (c) various locations within the cross-section and varying El while taking 0.1 , b 1 and 10 . Other parameters chosen are: and 1 2 2.0 and
The number of parametric variations, considering the immense number of
AC
combinations of wall zeta potentials that is possible is very large. Hence, we would like to restrict our studies to the above considered cases, as we feel they give a broad picture of how an Oldroyd-B fluid would behave in an electroosmotic flow through a rectangular channel.
V. Conclusions The effect of combining electroosmotic flow with rheological behavior of an Oldryod-B fluid in a rectangular channel was studied thoroughly, with focus being on symmetric and asymmetric zeta potentials. The following are some of the major conclusions from this study:
24
ACCEPTED MANUSCRIPT
The mathematical method and the numerical model can be used to study cases with both symmetric and asymmetric wall zeta potentials.
The flow startup is strongly dependent on El and , where large El and small tends to make the fluid more elastic and less viscous in nature.
We observe that the aspect ratio of the channel plays a strong role in determining the flow transience behavior. The flow velocity profiles of an Oldroyd-B fluid is essentially smooth as compared to sharp and coarse profiles for a Maxwell fluid.
CR IP T
We strongly believe that the results from this study will enhance the understanding of the behavior of a viscoelastic fluid within microchannels actuated by eletroosmosis, with farreaching implications in various applications including microfluidics based medical
AN US
diagnostics.
References
N.T. Nguyen, S.A.M. Shaegh, N. Kashaninejad, D.T. Phan, Design, fabrication and
M
[1]
characterization of drug delivery systems based on lab-on-a-chip technology, Adv.
[2]
ED
Drug Deliv. Rev. 65 (2013) 1403–1419.
Á. Ríos, M. Zougagh, M. Avila, Miniaturization through lab-on-a-chip: Utopia or reality for routine laboratories? A review, Anal. Chim. Acta. 740 (2012) 1–11. P.C. Braga, M. Moretti, A. Piacenza, C.C. Montoli, E.E. Guffanti, Effects of seaprose
PT
[3]
on the rheology of bronchial mucus in patients with chronic bronchitis. A double-blind
[4]
CE
study vs placebo., Int. J. Clin. Pharmacol. Res. 13 (1993) 179–85. J.R. Stokes, G.A. Davies, Viscoelasticity of human whole saliva collected after acid
AC
and mechanical stimulation, Biorheology. 44 (2007) 141–160.
[5]
Y.-J. Kang, S.-J. Lee, Blood viscoelasticity measurement using steady and transient flow controls of blood in a microfluidic analogue of Wheastone-bridge channel, Biomicrofluidics. 7 (2013) 54122.
[6]
S. Chakraborty, Microfluidics and Microscale Transport Processes, CRC Press, 2012.
[7]
A.D. Stroock, M. Weck, D.T. Chiu, W.T.S. Huck, P.J.A. Kenis, R.F. Ismagilov, G.M. Whitesides, Patterning Electro-osmotic Flow with Patterned Surface Charge, Phys. Rev. Lett. 84 (2000) 3314–3317.
25
ACCEPTED MANUSCRIPT [8]
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.
[9]
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(5) (2000) 1053–1057.
Phys. Chem. 69(11) (2002) 4017–4024. [11]
CR IP T
[10] C.L. Rice, R. Whitehead, Electrokinetic Flow in a Narrow Cylindrical Capillary, J.
P. Dutta, A. Beskok, Analytical Solution of Time Periodic Electroosmotic Flows: Analogies to Stokes‟ Second Problem, Anal. Chem. 73(21) (2001) 5097–5102.
[12] R. Sadr, M. Yoda, Z. Zheng, A.T. Conlisk, An experimental study of electro-osmotic flow in rectangular microchannels, J. Fluid Mech. 506 (2004) 357–367.
AN US
[13] S. Devasenathipathy, J.G. Santiago, K. Takehara, Particle Tracking Techniques for Electrokinetic Microchannel Flows, Anal. Chem. 74(15) (2002) 3704–3713. [14] N.A. Patankar, H.H. Hu, Numerical Simulation of Electroosmotic Flow, Anal. Chem. 70(9) (1998) 1870–1881.
M
[15] L. Ren, D. Li, Electroosmotic Flow in Heterogeneous Microchannels, J. Colloid Interface Sci. 243 (2001) 255–261. Dutta,
A.
Beskok,
T.C.
Warburton,
Numerical
simulation
of
mixed
ED
[16] P.
electroosmotic/pressure driven microflows, Numer. Heat Transf. Part A Appl. 41(2) (2002) 131–148.
PT
[17] X.X. Li, Z. Yin, Y.J. Jian, L. Chang, J. Su, Q.S. Liu, Transient electro-osmotic flow of generalized Maxwell fluids through a microchannel, J. Nonnewton. Fluid Mech. 187
CE
(2012) 43–47.
[18] M. Zhao, S. Wang, S. Wei, Transient electro-osmotic flow of Oldroyd-B fluids in a
AC
straight pipe of circular cross section, J. Nonnewton. Fluid Mech. 201 (2013) 135–139. [19] S. Wang, M. Zhao, X. Li, Transient electro-osmotic flow of generalized Maxwell fluids in a straight pipe of circular cross section, Open Phys. 12 (2014) 445–451.
[20] P. Abhimanyu, P. Kaushik, P.K. Mondal, S. Chakraborty, Transiences in rotational electro-hydrodynamics microflows of a viscoelastic fluid under electrical double layer phenomena, J. Nonnewton. Fluid Mech. 231 (2016) 56–67. [21] P. Kaushik, P. Abhimanyu, P.K. Mondal, S. Chakraborty, Confinement effects on the rotational microflows of a viscoelastic fluid under Electrical double layer phenomenon, J. Nonnewton. Fluid Mech. 244 (2017) 123–137. 26
ACCEPTED MANUSCRIPT [22] C. Zhao, C. Yang, Exact solutions for electro-osmotic flow of viscoelastic fluids in rectangular micro-channels, Appl. Math. Comput. 211 (2009) 502–509. [23] G. Bin Lee, L.M. Fu, C.H. Lin, C.Y. Lee, R.J. Yang, Dispersion control in microfluidic chips by localized zeta potential variation using the field effect, Electrophoresis. 25 (2004) 1879–1887. [24] C.Y. Lee, G. Bin Lee, L.M. Fu, K.H. Lee, R.J. Yang, Electrokinetically driven active micro-mixers utilizing zeta potential variation induced by field effect, J.
CR IP T
Micromechanics Microengineering. 14 (2004) 1390–1398.
[25] A. Miller, A. Villegas, F. Javier Diez, Characterization of the startup transient electrokinetic flow in rectangular channels of arbitrary dimensions, zeta potential distribution, and time-varying pressure gradient, Electrophoresis. 36 (2015) 692–702. [26] A.M. Afonso, M.A. Alves, F.T. Pinho, Electro-osmotic flow of viscoelastic fluids in
AN US
microchannels under asymmetric zeta potentials, J. Eng. Math. 71 (2011) 15–30.
[27] J. Escandón, E. Jiménez, C. Hernández, O. Bautista, F. Méndez, Transient electroosmotic flow of Maxwell fluids in a slit microchannel with asymmetric zeta potentials, Eur. J. Mech. - B/Fluids. 53 (2015) 180–189.
M
[28] E. Jiménez, J. Escandón, O. Bautista, F. Méndez, Start-up electroosmotic flow of Maxwell fluids in a rectangular microchannel with high zeta potentials, J. Nonnewton. Fluid Mech. 227 (2016) 17–29.
ED
[29] E. Jiménez, J. Escandón, O. Bautista, F. Méndez, A note on “Start-up electroosmotic flow of Maxwell fluids in a rectangular microchannel with high zeta potentials” [J.
PT
Non-Newton Fluid Mech. 227 (2016) 17–29], 2016. [30] P. Kaushik, P.K. Mondal, S. Chakraborty, Flow dynamics of a viscoelastic fluid
CE
squeezed and extruded between two parallel plates, J. Nonnewton. Fluid Mech. 227 (2016) 56–64.
AC
[31] J.G. Oldroyd, On the Formulation of Rheological Equations of State, Proc. R. Soc. London A Math. Phys. Eng. Sci. 200 (1950).
[32] C. Fetecau, C. Fetecau, M. Kamran, D. Vieru, Exact solutions for the flow of a generalized Oldroyd-B fluid induced by a constantly accelerating plate between two side walls perpendicular to the plate, J. Nonnewton. Fluid Mech. 156 (2009) 189–201. [33] H. Qi, M. Xu, Stokes‟ first problem for a viscoelastic fluid with the generalized Oldroyd-B model, Acta Mech. Sin. 23 (2007) 463–469. [34] Y. Kang, C. Yang, X. Huang, Dynamic aspects of electroosmotic flow in a cylindrical microcapillary, Int. J. Eng. Sci. 40 (2002) 2203–2221. 27
ACCEPTED MANUSCRIPT [35] R.J. Huter, Zeta Potential in Colloid Science, Academic Press, London, 1981. [36] J.H. Masliyah, S. Bhattacharjee, Electrokinetic and colloid transport phenomena, John Wiley & Sons, 2006. [37] L.E. Rodd, J.J. Cooper-White, D.V. Boger, G.H. McKinley, Role of the elasticity number in the entry flow of dilute polymer solutions in micro-fabricated contraction geometries, J. Nonnewton. Fluid Mech. 143 (2007) 170–191. [38] M. Dehghan, A. Shokri, A numerical method for solving the hyperbolic telegraph
CR IP T
equation, Numer. Methods Partial Differ. Equ. 24 (2008) 1080–1093.
[39] S.A. Teukolsky, Stability of the iterated Crank-Nicholson method in numerical relativity, Phys. Rev. D. 61 (2000) 87501.
[40] K. Macounova, C.R. Cabrera, M.R. Holl, P. Yager, Generation of natural pH gradients in microfluidic channels for use in isoelectric focusing, Anal. Chem. 72 (2000) 3745–
AN US
3751.
[41] Y. Gao, T.N. Wong, C. Yang, K.T. Ooi, Two-fluid electroosmotic flow in
AC
CE
PT
ED
M
microchannels, J. Colloid Interface Sci. 284 (2005) 306–314.
28
Startup electroosmotic flow of a viscoelastic fluid characterized by Oldroyd-B model in a rectangular microchannel with symmetric and asymmetric wall zeta potentials P Kaushik , Suman Chakraborty PII: DOI: Reference:
S0377-0257(16)30113-6 10.1016/j.jnnfm.2017.06.003 JNNFM 3902
To appear in:
Journal of Non-Newtonian Fluid Mechanics
Received date: Revised date: Accepted date:
28 July 2016 9 June 2017 10 June 2017
Please cite this article as: P Kaushik , Suman Chakraborty , Startup electroosmotic flow of a viscoelastic fluid characterized by Oldroyd-B model in a rectangular microchannel with symmetric and asymmetric wall zeta potentials, Journal of Non-Newtonian Fluid Mechanics (2017), doi: 10.1016/j.jnnfm.2017.06.003
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT
AC
CE
PT
ED
M
AN US
CR IP T
Highlights Transient electroosmotic flow of a viscoelastic fluid in a rectangular microchannel Rheology of viscoelastic fluid is described by Oldroyd-B fluid model Complex nature of the fluid leads to enhancement of flowrates for symmetric zeta potentials Complex velocity profiles are observed for asymmetric zeta potentials
1
ACCEPTED MANUSCRIPT
Startup electroosmotic flow of a viscoelastic fluid characterized by Oldroyd-B model in a rectangular microchannel with symmetric and asymmetric wall zeta potentials P Kaushik and Suman Chakraborty* Department of Mechanical Engineering,
CR IP T
Indian Institute of Technology Kharagpur, Kharagpur, 721302, India
AN US
___________________________________ *Corresponding Author Email address: [email protected]
ABSTRACT
M
We study the startup startup flow of a viscoelastic fluid in a microfluidic channel, by considering the Oldroyd-B constitutive model. In particular, we aim to capture the confluence
ED
of the interfacial charge and the combined viscous and elastic properties of the fluid towards dictating the resulting startup flow characteristics. We further highlight the implications in the
PT
asymmetric nature of the wall zeta potentials on the flow behavior. We observe an increase in flow rate with increase in elastic nature of the fluid for studies with symmetric wall zeta
CE
potentials. For cases with asymmetric zeta potentials, we observe large movement of the velocity peaks, creating crests and troughs. Results from the present study may turn out to be useful in designing lab on a chip based microfluidic devices transporting bio-fluids such as
AC
blood, saliva and mucus by means of electroosmotic mechanisms.
Keywords: Viscoelastic Fluid; Rectangular microchannel; Startup flow; Transient analysis; Oldroyd-B model I. Introduction Recent advances in microfluidic applications, especially in Lab-on-chip devices, have triggered interests in exploring the electroosmotic method of flow activation. Lab-on-chip devices can be used for the analysis of chemical and biological samples which find many 2
ACCEPTED MANUSCRIPT applications, especially in medicine [1,2]. Many of the bio-fluids usually involved in such applications, such as blood, saliva, mucus, etc., are non-Newtonian in nature [3–5]. Implications of studying electrosmotically actuated flows of viscoelastic fluids may therefore be far-reaching, ranging from the development of medical devices with artificial flow actuation to portable diagnostic kits. The basic idea in electroosmotic flow is that the existence of a small potential at the wall called zeta potential causes the ions of the liquid to move. The ions charged oppositely
CR IP T
to the surface charge tend to move close to the wall surface, creating a thin layer (in which ionic species of charge of sign opposite to that of the wall are more abundant than ionic species of charge of same sign as that of the wall), so as neutralize the surface charge. The charged surface, together with the charged interfacial layer, is called the electrical double layer (EDL). The EDL formation creates an imbalance of charge in the fluid, which can drive
AN US
a fluid flow by exploiting the viscous drag between the ions and the fluid, on applying an external electrical potential gradient [6]. It is well known that the lack of moving parts makes it easier for fabrication and pumping of a fluid in microchannels by this method. Further control over fluid motion in electroosmosis can also be well accomplished by methods such
M
as surface patterning [7]. Velocity profiles in a typical electroosmotic flow with uniform wall zeta potential are plug-like, unlike pressure driven flows where the velocity profiles are
ED
usually parabolic. The plug-like velocity profiles can be used for separation of charged species such as DNA [8]. This tailored zeta potential at the walls makes controllability of the flow much better within the microchannels. And thus makes the study of electroosmotic flow
PT
with tailored within microchannels very important. Therefore, we will further look at the important developments in the study of electroosmotic flows existing in literature.
CE
Electroosmotic flow has been extensively studied for Newtonian fluids via theory [9– 11], experiments [12,13] and numerical methods [14–16]. However, the startup
AC
characteristics of electroosmotic flows of complex fluids have not been well explored. Some recent works on startup flows of viscoelastic fuids include studies in idealized geometries such as parallel plate channel [17], microtubes [18,19] and rotating channels [20,21], where Maxwell and Oldroyd-B models have been considered. Zhao and Yang [22] could find the exact solution to some special cases of generalized Oldroyd-B fluid which includes results for a Maxwell fluid and a generalized second grade fluid for small symmetric wall zeta potentials. The authors in their study also mentioned that full analytical solution may not be possible for Oldroyd-B fluid. Additionally the above studies have considered only symmetric zeta potentials at the walls. It is important to note that surface modifications can be used to 3
ACCEPTED MANUSCRIPT vary the wall zeta potentials and can be used for better control of the fluid [23,24]. Taking the possibility of asymmetry in zeta potential into consideration, some studies have been reported for Newtonian fluids [25], Phan-Thien-Tanner(PTT) fluids [26] and Maxwell fluids [27] between two parallel plates. Experimental study with Newtonian fluids in a channel with asymmetric zeta potentials was conducted and reported by Miller et al. [25]. In a more recent study, Jiménez et al. [28,29] studied the startup electroosmotic flow of Maxwell fluid in a rectangular microchannel with high assymetric zeta potentials by a semi-analytical method.
CR IP T
In this present study, we attempt to numerically analyze the electroosmotically driven transport of Oldroyd-B fluids in a rectangular microfluidic channel. We also compare and validate the results from our numerical method with Jiménez et al. [28,29] for a special case in which the retardation time is set to zero, for which the fluid behaves asymptotically like an idealized Maxwell fluid. First we focus on the cases where there is symmetry in the wall zeta
AN US
potentials; here we also compare how the complex nature of the fluid affects the flow rates. Then we focus on the hydrodynamics involved when there is asymmetry in the wall zeta potentials. The results are then concluded.
M
II. Problem description and mathematical formulation.
ED
A. Physical model
Here, we consider the transient startup flow of viscoelastic fluid characterized by Oldroyd-B
PT
constitutive relation in a rectangular microfluidic channel as shown in Fig. 1. The coordinate axes are chosen to be x, y and z , which are along the length, width and height of the channel
CE
respectively. The following assumptions have been considered for the present study: (i) The length ( L ) of the channel is much larger than its width and height. (ii) The fluid is assumed to be initially at rest and unstressed.
AC
(iii) The channel cross section is taken to be of height 2H and width 2W and the origin is placed at the channel centerline, as shown in Fig 1.
(iv) The buffer solution is a symmetric
z : z
electrolyte with constant physical
properties. (v) The wall zeta potentials are taken as shown in Fig. 1.
u H (vi) The flow Reynolds number Re c
is assumed to be very small for
4
ACCEPTED MANUSCRIPT electroosmotic flows Re
1 [28].
(vii)
The flow is considered to be unidirectional.
(viii)
No external pressure gradient is applied.
k BTref Ex Here, uc is the characteristic velocity scale, is the permittivity of z e the fluid, k B is the Boltzmann constant, Ex is the externally applied electric field parallel
CR IP T
to the flow direction. Tref is the reference temperature of the fluid (temperature at which the fluid is flowing), z is the valence of the electrolyte, e is the fundamental electron
AN US
charge and is the sum of base fluid and polymer viscosities.
M
Figure 1. (Color online) Schematic depicting the problem considered in the present study.
ED
B. Governing equations of the flow field
The governing equations of the flow dynamics in the channel considered are given by
CE
PT
[22,30,31]:
u 0
u u u p T e E t
(1) (2)
AC
In Eq. (2), is the fluid density, T is the stress tensor, e is the ion density, E is the total electric field and u is the velocity vector. Following the Oldroyd-B constitutive model for the stress tensor equation, we get, [32]: T 1 T 2 D 2 D
(3)
where, 1 and 2 are the fluid relaxation and retardation times respectively and the upper convected derivative of T is defined as [33]:
T
T T u T u T T u t
(4)
5
ACCEPTED MANUSCRIPT Here, D is defined as
D
1 u uT 2
(5)
In order to simplify the above equation, we split T as T T1 T2 , where
T2 2
2 D 1
(6)
Therefore, the momentum equation further reduces to:
where,
AN US
T1 1 T1 2 1 2 D 1
CR IP T
u u u p T1 2 2 D e E t 1
(7)
(8)
The components of T1 matrix appearing in Eq. (7) are as follows: Txx T1 Tyx Tzx
Txy Tyy Tzy
Txz Tyz Tzz
(9)
M
As a consequence of the approximately unidirectional flow assumption, we try to find
ED
the flow solution for u t , y, z at t 0 assuming v w 0 . In the subsequent section, we show the method of resolving the electroosmotic body force term in Eq. (7).
PT
C. Resolution of electroosmotic body force term It is important to mention beforehand that the time scale related to the ionic transport effects
CE
within the EDL is very short and of the order
107 108 s , which is at least 2-3 orders
smaller than characteristic time associated with the electroosmotic flow [34]. Hence the
AC
transient term in the equation that governs the electrical body force term e E in Eq. (7) can be neglected and we can consider the Poisson-Boltzmann equation as the equation which governs the electrical potential distribution within the EDL, as [20,35,36]: e z 2 2 2 zen Sinh y 2 z 2 k BTref
(10)
Here, e is the total ionic charge density, which for a z : z symmetric electrolyte
6
ACCEPTED MANUSCRIPT e z using Boltzmann distribution is given by e 2 zen Sinh k T B ref
and is the inverse of the
EDL thickness. Here, z is the valence of the ions, n is the concentration of the bulk and
Tref is the reference temperature (ambient temperature of the fluid).
D. Simplification and Non-dimensionalization of the governing equations
y*
CR IP T
The dimensionless parameters used are given as follows [28,36]:
e z Wi y z u TH t , z* , * , u* , t * , T* , El 1 2 , 2, 2 k BTref H H uc uc 1 H Re H
* kH , b* W / H
(11)
Where, El is the elasticity number which considers the effect of fluid elasticity and inertia. It is
AN US
defined as the ratio of the relaxation time scale of the fluid to the diffusion time scale [37]. The ratio of the fluid retardation to relaxation time is represented by . We would like to mention here that since the velocities in the y and z directions are much smaller, we neglect them and consider our flow to be approximately unidirectional. This then reduces the equations of stress in the Oldroyd-B model [See Eq. (8)] to their simple form given below in Eqs. (14) and (15). The simplified governing
M
equations are given as:
ED
2 * 2 * *2 *2 Sinh * *2 y z
AC
CE
PT
* 2u* 2u* *2 u* Txy Txz* * *2 *2 Sinh t * y* z * y z
Txy*
(13)
u* y*
(14)
Txz* u* Txz El * 1 * t z
(15)
* * u* Txx* * u 2 Wi T T xy xz t * y* z *
(16)
Txy* El
t *
1
*
Txx* El
(12)
In order to simplify the equations further we combine equations (13)-(15) and we get,
El
2u* u* 2u* 2u* *2 1 El *2 Sinh * *2 * * *2 t t t y z
(17)
In the present study, equations (12) and (17) have been solved numerically along with
7
ACCEPTED MANUSCRIPT no-slip condition at the four walls, to obtain the solution till steady state is reached. As assumed before, the fluid is considered to be free of any stresses at the start of the flow with
u* zero velocity. Mathematically the initial conditions used are u 0, y , z * 0, y* , z* 0 t *
*
*
and boundary conditions used are u* t* , y* , 1 u* t* , b* , z* 0 . It is worth noting that for a special case of 0 , the fluid model is reduced to the Maxwell model and the equation
El
CR IP T
(17) reduces to
2u* u* 2u* 2u* * *2 *2 *2 Sinh * *2 t t y z
(18)
From here on, we will drop the superscript „*‟ in the equations and non-dimensional
AN US
terms during further discussion of the problem for the convenience of representation. III. Numerical Method and Model Benchmarking
A finite difference code has been developed for solving the transport equations along with their respective boundary conditions as given in the preceding section. The Poisson-
M
Boltzmann equation (12) is solved by a simple second order in space discretization using SOR method. The equation (17) is hyperbolic and the code has been developed based on
ED
numerical method suggested by Dehghan and Shokri [38] for a hyperbolic telegraph equation. The method uses a second order scheme for both the temporal discretization as well as the discretization of the spatial terms. The discretized equation is as shown below:
u t t , y, z 2u t , y , z u t t , y , z u t t , y , z u t t , y , z 2 2 t t
PT
El
(19)
AC
CE
El 2 El 2 2 = u t t , y, z u t t , y, z 1 u t , y, z 2 t 2 t 2 Sinh y, z
The value of chosen is 0.5 as suggested by Dehghan and Shokri [38], which is the
well known Crank-Nicholson scheme [39] to track the steady shear diffusion term. For the solution, we use the well-known Gauss-Seidel method and the code is run till the steady state is reached. Validation study – In order to ensure that the results obtained from present numerical method
8
ACCEPTED MANUSCRIPT are correct, we try to validate our method with the results available in the literature. Semianalytical solution for the electroosmotic flow of a Maxwell fluid in a rectangular channel is available in literature [28]; using this we validate our model for special case of 0 as shown in equation (18). Fig 2(a) shows depicts the benchmarking of our numerical code. The circular markers are used to represent the results reported by Jiménez et al. [28], while solid line represents the results from the present method. There is a good match obtained between the present study
CR IP T
and reported results. We have also benchmarked our code with semi-analytical solution of Zhao and Yang [22] for a generalized Oldroyd-B fluid in electroosmotic flow with symmetric zeta potential case considering Boltzmann distribution at the center of the channel (0,0) as depicted in Fig. 2(b). We observe very close match with the results reported in the literature.
PT
ED
M
AN US
We have mentioned other parameters used in plotting the figure in the caption.
(a)
(b)
CE
Fig 2 (color online) Model Validation: (a) Plot showing the temporal variation of dimensionless velocity u t, y, z at two different locations which were studied by Jiménez et al. [28]. Location A represents
AC
y, z 3.12, 0.91
y, z 2.75, 0.81
and location B represents
for an aspect ratio of W : H 2 :1 i.e. b 2 . The other parameters
considered here are: 10 , 1 0.2 , 2 0.2 , 3 0.2 , 4 0.2 and El 0.5 . (b) Plot showing the temporal variation of dimensionless velocity u t, 0, 0 at channel center which was studied by Zhao and Yang [22] considering Boltzmann distribution. The other parameters considered here are: 10 , 1 2 3 4 1 , 0.1 and El 2 .
9
CR IP T
ACCEPTED MANUSCRIPT
(b)
AN US
(a)
M
(c) Fig 3 (color online) Plot showing (a) the variation of dimensionless velocity u t, y, z vs z
ED
at y 0 with various grid size and (b) the variation of dimensionless velocity u t, y, z vs z at y 0 with various time step size (c) the variation of normal stress Txx y, z vs z at y 0
CE
PT
with various grid size. Here we choose the following parameters: 10 , b 1 , 0.1 and El 2.0 , considering asymmetric zeta potentials at walls 1 1 , 2 1 , 3 1 , 4 1 (symmetric zeta potentials at walls in the inset 1 2 3 4 1 ). Grid independence and convergence study – To make sure that the results obtained in the study are independent of the grid size and time step size, we present the grid and time step
AC
independence studies as well. We depict, in Figs. 3(a)-(b), the variation of axial velocity
u t, y, z as a function of z for a case with asymmetric wall zeta potentials (cases with
symmetric wall zeta potentials are shown in the in inset) for different grid resolutions
y z and time step t variations. It is worth noting from the figures that the variation of velocity becomes insignificant as changes from 1/96 to 1/128 and while a change in t from 102 to 103 does not show any appreciable change in velocity as well. This is true for both the symmetric and asymmetric wall zeta potential cases. We also show
10
ACCEPTED MANUSCRIPT
the variation of normal stress Txx for different grid resolutions (cases with symmetric wall zeta potentials are shown in the in inset). Similar to the variation of axial velocity the variation of Txx becomes insignificant as changes from 1/96 to 1/128. Accordingly, we have considered t 103 and 1 96 throughout this study. In order to understand the important findings of the present study, we would like to show the typical values of different model parameters chosen for the numerical computation. 106 M
(the bulk concentration
CR IP T
A typical concentration of the ionic solution
n 6.023 1020 m3 ) [40,41] at temperature Tref 298K with Boltzmann constant
kB 1.38 1023 and the permittivity of the solvent 80 8.85 1012 we find that the EDL thickness 1 is ~300 nm . Here we assume the characteristic channel dimension as
AN US
~ 10 m . Using these values, we find the dimensionless EDL thickness as: * H 30 . The typical values of Elasticity number used in the literature are varies between 0 and 1 [20].
M
IV. Results and Discussions
1 [37] and usually
The results and discussions will be divided into two parts. In the first part we present the
ED
results from the studies of symmetric zeta potentials at the walls. In the second part we attempt to highlight the results for asymmetric zeta potentials at the walls.
PT
A. Studies with symmetric wall zeta potentials
CE
We will begin our discussion with the studies on the temporal variations on the velocity profiles when two fluid parameters namely El and are varied. This is depicted in figs. 4(a)-
AC
(d), where the other parameters chosen are mentioned in the figure caption. It can be seen from Figs. 4(a)-(d), that the velocities tend to increase in magnitude as the time progresses. For smaller values of and larger values of El, the elastic nature of the fluid becomes more pronounced. This can be attributed to the fact that fluid relaxation time increases with increase in El as well as decrease in . Larger relaxation time means that the fluid takes more time to settle down to a steady state value. It is evident from fig. 4(c) that for a case with larger value of El and smaller value of , the velocity field exhibits a spring-like behavior, with the fluid centerline velocity going to a very large value at an intermediate time
11
ACCEPTED MANUSCRIPT of t 1 and then going to a smaller value at t 2 and eventually settling down slowly into and intermediate value which is the steady state profile after a certain number of oscillations. Considering the opposite case as depicted in fig. 4(b), where the value of El is small and the value of is large, the viscous effect tends to dominate over the elastic effects and we do not see much spring-like behavior in the velocity profile. Elastic nature of viscoelastic fluids may be attributed to the stretching of the polymer chains within the solution. For the case of
CR IP T
1 , the fluid will tend to behave like a Newtonian fluid with no elastic behavior. Inherently from equation (17) it is evident that the steady state solution of the equation is same as the solution for a Newtonian fluid irrespective of the value of El and . This fact is clearly evident in the solution at steady state as well as depicted by t 1 in figs. 4(a) and (b) and by t 10 in figs. 4(c) and (d). It can also be observed from Figs. 4(a)-(d) that the time
ED
M
AN US
taken to reach steady state is strongly affected by the variation of El and .
(b)
AC
CE
PT
(a)
(c) (d) Fig 4 (color online) Plot showing the temporal variation of dimensionless velocity u t, y, z vs z at y 0 for (a) El 0.1 and 0.1 , (b) El 0.1 and 0.8 , (c) El 1 and 0.1 , (d) El 1 and 0.8 . Other parameters chosen are as follows: 10 , b 1 and
12
ACCEPTED MANUSCRIPT 1 2 3 4 1 . Before discussing the results of figs. 5(a)-(d), it is important to define the flow rate of the fluid in the cross-section. In study of flows in microchannels, flow rates tend to play an important role as they give the overall picture of how effective our forcing/pumping mechanism is. Augmentation of flow rates in microchannels is sought by various means. In our present study, we define the dimensionless flow rate as:
u t , y, z dzdy
W 1
CR IP T
Q t
W 1
(19)
The effect of various fluid, flow and geometric parameters on the flow rate has been considered in figs. 5(a)-(d). The other parameters which have been used in the plots have been described in the figure caption. The symmetry of the wall-zeta potentials has an aiding
AN US
effect on the flow. In fig 5(a), the effect of has been considered on the flow rate. It is clear from the figure that decreasing increases the flow rate initially, before reaching steady state eventually. The time taken by the fluid to reach steady state also drastically increases with decrease in . This may be attributed to the fact that the fluid retardation time is very
M
short compared to the relaxation time for smaller . This means that the fluid tends to relax slowly to any external forcing or changes in the nearby fluid. Also we see that for smaller
ED
value of the oscillatory behavior in the fluid becomes more prominent, with larger number oscillations in the flow rate magnitude, due to the dominance of elastic effects as compared to
PT
the case with larger where viscous effects dominate. Fig 5(b) shows the effect of El on the flow rate. Increase in El essentially signifies the increase in relaxation time of the fluid.
CE
Therefore, as evident from fig. 5(b), the increase in El effectively increases the relaxation time and therefore increases the elastic behavior of the fluid. It can also be noted that with
AC
increase in El, the time taken by the fluid to reach the steady state also drastically increases along with increase in the magnitude of oscillations. However, the number of oscillations tends to remain more or less same. In order to have large flow rates it may be desirable to use fluids with larger values of El and smaller values of . Fig. 5(c) depicts the effect of or the inverse of EDL thickness on the flow rate. It may be observed clearly that the effect of
on the flow rate is quite less and we may consider the flow rate as a weak function of . Clearly, the flow rate increases in magnitude with increase in , both during transience and at steady state. Increase in represents decrease in the EDL thickness, which leads to more
13
ACCEPTED MANUSCRIPT ion mobility in the bulk flow. This could be attributed as the reason for increase in flow rate with . Also, we observe for larger values of such as 30 and 50 , the flow rate becomes almost independent of . This goes to show that beyond a certain reduction in the EDL thickness, the mobility of ions cannot increase drastically. Fig. 5(d) describes the effect of aspect ratio of the channel b on the flow rate. It is quite evident that for larger aspect ratios, the wall effects become less prominent, however the total area available for flow increases and therefore the flow rate increases. However, it is important to note that when the
CR IP T
aspect ratio is doubled, the flow rate does not double. This may be attributed to the fact that in electroosmotic flows, the fluid forcing essentially occurs close to the walls, which makes the aspect ratio study very important. Having small aspect-ratios may enhance the forcing of the fluid close to the wall. This causes the flow rate not to have proportional increase with
PT
ED
M
AN US
increase in aspect rate as seen from fig. 5(d).
(b)
AC
CE
(a)
(c)
(d)
Fig 5 (color online) Plot showing the temporal variation of dimensionless flow rate Q t for (a) various taking El 2.0 , b 1 and 10 , (b) various El taking 0.1 , b 1 and 14
ACCEPTED MANUSCRIPT 10 , (c) various taking El 1.0 , b 1 and 0.1 , (d) ) various b taking El 1.0 , 10 and 0.1 . Other parameters chosen are: and 1 2 3 4 1 . In flows of viscoelastic fluids, it becomes important to consider how the flow velocity varies with location. This is addressed in fig 6(a)-(c), where the temporal variations of the velocities are shown at certain locations within the cross-section in case with symmetric wall zeta potentials. The other parameters used for creating the figure have been given in the
CR IP T
caption. From fig. 6(a), we observe that the magnitude of the velocity is the highest at (0,0) and lower closer to the walls (0.9,-0.9). We also see that the increase in velocity magnitude takes place closer to the walls (0.9,-0.9) earlier compared to the case with center of the channel (0,0). This may be attributed to the fact that in electroosmotic flow, the forcing essentially happens closer to the walls in the area of the EDL. We also observe that the elastic
AN US
behavior of the fluid is more prominent away from the walls (0.9,-0.9) and closer to the center (0,0). This may be attributed to the fact that the wall tends to aid the viscous effects within the fluid close to it. We also see from fig. 6(a) that the transient behavior of the fluid velocity at (0,0) and (0.5,-0.5) differs, however, the steady state value is almost the same indicating a plug-like flow at steady state with values same as that of a Newtonian fluid. In
M
fig. 6(b), we see the temporal variation of velocity with different at two different locations. It is observed that the oscillations during transience tends to die out faster in case larger
ED
values and reaches a common value at each location at steady state. It can also be seen that the value of the velocity is higher with more number of oscillations for smaller ,
PT
irrespective of the location. We also see that the time at which the peaks occur at the different locations does not vary much with change in . This may be attributed to the fact that for
CE
smaller , the retardation time is considered to be smaller, which causes the viscous effects to be weaker. Hence, we see a small delay in occurrence of the peak at smaller values of .
AC
This effect, however, is seen to be quite weak. Fig. 6(c) depicts the temporal variation of the velocity with El for two different locations. We again see that the larger value of velocity occurs as El increases and also we see the enhancement in the elastic behavior of the fluid. This is because the fluid relaxation time increases and it takes more time to adjust to any change in the nearby zones. This also causes the velocity peaks to strongly shift further downstream in time with larger El.
15
CR IP T
ACCEPTED MANUSCRIPT
(b)
M
AN US
(a)
(c) Fig 6 (color online) Plot showing the temporal variation of dimensionless velocity u t, y, z
PT
ED
for (a) various locations within the cross-section taking El 2.0 , 0.1 , b 1 and 10 , (b) various locations within the cross-section and varying while taking El 2.0 , b 1 and 10 , (c) various locations within the cross-section and varying El while taking 0.1 , b 1 and 10 . Other parameters chosen are: and 1 2 3 4 1 .
CE
The above analysis concludes the cases studied for symmetric wall zeta potentials. We now move on to cases where each of the walls has different zeta potentials, which may be
AC
either positive or negative. B. Studies with asymmetric wall zeta potentials For studies considering asymmetric wall zeta potentials, we start the discussion with temporal evolution of velocity profiles for a case where all zeta potentials are positive while being different, and another case with zeta potentials at certain walls are negative. Fig. 7(a) and (b) show the temporal evolution of the velocity profile at the middle of the y direction, along zdirection. The other parameters which have been chosen are mentioned in the figure caption.
16
ACCEPTED MANUSCRIPT The commonality in both fig. 7(a) and 7(b) is that the at early time the flow essentially is restricted close to the EDL zone and as time progresses this flow velocity gets diffused into the entire domain. The important observation in this case is that unlike Maxwell fluid, where the diffusive terms only start dominating at larger times, for the Oldroyd-B fluid model the diffusive terms play an important role throughout the flow time. This essentially means that the velocity profiles are quite essentially smooth for an Oldroyd-B fluid unlike a Maxwell fluid. This phenomenon is clearly observed in Fig. 4(a)-(d) as well as Fig 7(a)-(b). In fig.
CR IP T
7(a), we observe that despite the skewed nature of the wall zeta potentials, at time t 1.5 , the maximum value of u t , y, z occurs almost at the centerline i.e. at z 0 . On either side of t 1.5 , i.e. t 1 and t 2 , we see the maximum of u t , y, z occurring on either side of z 0 . As the fluid is pulled towards z 0 side at t 1 , the fluid relaxes by the maximum
AN US
moving across to z 0 side at t 2 . The fluid then eventually attains a steady state after some oscillation. In fig. 7(b), the zeta potentials across z direction are equal in magnitude and opposite in direction. Here we observe that the spring like nature causes the fluid to have
u t , y, z 0 at some locations in z 0 and u t , y, z 0 at some locations in z 0 at t 2 .
M
This is converse to what happens at steady state, where, u y, z 0 for z 0 and u y, z 0
AC
CE
PT
ED
for z 0 .
(a) (b) Fig 7 (color online) Plot showing the temporal variation of dimensionless velocity u t, y, z vs z at y 0 for (a) Asymmetric non-negative zeta potentials: 1 1.0 , 2 2.0 , 3 2.0 and 4 1.0 , (b) Asymmetric with negative zeta potentials: 1 2.0 , 2 2.0 , 3 2.0 and 4 2.0 . Other parameters chosen are as follows: El 2 , 0.1 , 10 and b 1 . Fig 8. shows the temporal variation of the dimensionless velocity distribution
17
ACCEPTED MANUSCRIPT u t, y, z for asymmetric wall zeta potentials. The parameters used for generating the figure are mentioned in the caption. In the figure, dark red colour indicates large positive value of
u t, y, z and dark blue colour represents large negative value of u t, y, z , while the green colour represents values close to zero. We observe from figs. 8(a) and (b) that for very early times the flow of fluid is essentially restricted close to walls near the EDL zone. As time progresses, the flow spreads into the complete domain by diffusion as seen in figs. 8(c)-(f). It
CR IP T
is also important to note that the velocity profile u t, y, z is essentially smooth unlike the case of Maxwell fluid, which has been studied in detail by Jiménez et al. [28]. This is due to the fact that the viscous effects in the flow are prominent even at small times in an Oldroyd-B fluid unlike Maxwell fluid. We see that during the intermediate part of the transience, the combination of elastic and viscous effects cause a smooth maximum which moves away from
AN US
the walls and closer to the center of the domain considered. At t 0.5 in fig. 8(c), where the maximum is near the bottom edge of z 1 , the shift in the maximum takes place towards the center at t 1.5 , after which the maximum again shifts back towards zone near the
z 1 wall, where, 4 4.0 as the flow reaches steady state. At steady state, as depicted in
M
figs. 8(e)-(h), the u t, y, z profile becomes independent of El and , which also is evident from figs. 7(a)-(b). Elastic behavior dominates viscous behavior in the case described by fig.
ED
8 as 0.1 is small and El chosen is large. The stretching and relaxing nature of the fluid is typical of dilute polymer solutions, where the polymer chains act similar to springs. The
PT
stretching causes the potential energy to get stored within the polymer chains and it is released when the fluid starts relaxing. Despite this, the viscous nature is quite evident
AC
CE
because of strong diffusion within the fluid leading to smooth velocity profiles.
18
ACCEPTED MANUSCRIPT
(b)
M ED
(g)
(e)
AN US
(d)
(c)
CR IP T
(a)
(h)
(f)
(i)
Fig 8 (color online) Plot showing the temporal variation of dimensionless velocity u t, y, z vs for asymmetric zeta potentials: 1 1.0 , 2 2.0 , 3 3.0 and 4 4.0 . Other parameters
PT
y, z
AC
CE
chosen are as follows: El 2 , 0.1 , 10 and b 1 . The times displayed are (a) t 0.1 , (b) t 0.2 , (c) t 0.5 , (d) t 1 , (e) t 1.5 , (f) t 2 , (g) t 4 , (h) t 10 and (i) t 20 .
19
ACCEPTED MANUSCRIPT
(b)
(c)
CR IP T
(a)
y, z
AN US
(d) (e) (f) Fig 9 (color online) Plot showing the temporal variation of dimensionless velocity u t, y, z vs for asymmetric zeta potentials: 1 1.0 , 2 2.0 , 3 2.0 and 4 1.0 . Other
M
parameters chosen are as follows: El 2 , b 1 , 10 . The times and chosen are (a) 0.1 and t 0.5 , (b) 0.1 and t 1 , (c) 0.1 and t 1.5 , (d) 0.8 and t 0.5 , (e) 0.8 and t 1 , (f) 0.8 and t 1.5 .
ED
In order to highlight the effect of on asymmetric wall zeta potentials, we consider
0.1 and 0.8 which are used to generate Fig 9, while keeping all other parameters
PT
same. The other parameters chosen are mentioned clearly in the figure caption. The top line of fig. 9, i.e. figs. 9(a)-(c) show the temporal variation of the dimensionless velocity
CE
distribution u t, y, z for 0.1 with times t 0.5, 1, 1.5 and the bottom line i.e. figs. 9(d)(f) show the temporal variation for 0.8 . In the figure, dark red colour indicates large
AC
positive value of u t, y, z and dark blue colour represents large negative value of u t, y, z , while green being an intermediate colour representing values close to zero. We observe from figs. 9(a)-(c) that the elastic nature of the fluid dominates over the viscous nature for 0.1 and vice versa for 0.8 , as depicted by 9(d)-(f). We also observe that the flow reaches steady state quickly for the case with 0.8 . This may be observed by seeing figs. 9(e) and (f), which are quite similar; this matches closely with the steady state case shown in fig. 8(i). For 0.8 case, the relaxation time is smaller and therefore the fluid relaxes quickly to any disturbance that may occur in the fluid and hence steady state is achieved more quickly than
20
ACCEPTED MANUSCRIPT in the case of 0.1 . Effect of El on the asymmetric wall zeta potentials is highlighted in fig. 10, where a small El value of El 0.1 and a large value of El 2 have been considered for study. Large value of El implies that the relaxation time of the fluid is large and it takes time to adjust to any changes in flow. The other parameters used in generating fig. 10 are displayed in the figure caption. Magnitudes of velocity represented by the various colours are as followed in
CR IP T
the previous figures. It is clearly evident from figs. 10(a)-(f) that for larger value of El in case of figs. 10(d)-(f), the elastic nature of the fluid is more prominent, as the fluid takes longer time to adjust to changes in the flow field. For very small value of El 0.1 , which is chosen in 10(a)-(c), the fluid behavior is close to that of a Newtonian fluid because of small relaxation times. The flow reaches steady state quickly for El 0.1 case as evident from
(b)
(c)
AC
CE
PT
(a)
ED
M
AN US
similarity in flow profiles shown in figs. 10(b) and (c).
(d)
(e)
(f)
Fig 10 (color online) Plot showing the temporal variation of dimensionless velocity u t, y, z vs
y, z
for asymmetric zeta potentials: 1 1.0 , 2 2.0 , 3 2.0 and 4 1.0 . Other
parameters chosen are as follows: 0.1 , b 1 , 10 . The times and El chosen are (a) El 0.1 and t 0.5 , (b) El 0.1and t 1 , (c) El 0.1 and t 1.5 , (d) El 2 and t 0.5 , (e) El 2 and t 1 , (f) El 2 and t 1.5 . The effect of aspect ratio has been shown in fig. 11. The top line i.e. figs. 11(a)-(c)
represents the case where aspect ratio defined by b is taken as b 1 and the bottom line i.e.
21
ACCEPTED MANUSCRIPT figs 11(d)-(f) represents the case with b 2 . The effect of aspect ratio becomes quite important for asymmetric wall zeta potentials because of inherent dimensional asymmetry as well. The other parameters used are displayed in the figure caption. We see from the two figures that the corner effect is more prominent for the b 2 case at early times, as displayed in figs. 11(a) and 11(d) for t 0.5 . We observe from fig. 11(d) that there is a maximum at the corner near (-2,1), whereas, such a maximum is not there in the case of fig 11(a). We also observe that in case of larger aspect ratio for b 2 , the amount of time required for the
CR IP T
momentum to diffuse into the fluid is more because of larger length scale in the y-direction. This causes a delay in the flow velocity profiles as compared to the case with smaller aspect ratio of b 1 . This fact is evident from fig. 11(f), where we see that the flow is still in the main transience stage, where as in case of fig. 11(c), the flow has started to approach steady state. It can be concluded that larger aspect ratio case takes longer to reach steady state as
(b)
(c)
CE
PT
(a)
ED
M
AN US
compared to the case with smaller aspect ratio.
(e)
(f)
AC
(d)
Fig 11 (color online) Plot showing the temporal variation of dimensionless velocity u t, y, z vs
y, z
for asymmetric zeta potentials: 1 1.0 , 2 2.0 , 3 2.0 and 4 1.0 . Other
parameters chosen are as follows: El 2 , 0.1 , 10 . The times and aspect ratio chosen are (a) b 1 and t 0.5 , (b) b 1 and t 1 , (c) b 1 and t 1.5 , (d) b 2 and t 0.5 , (e) b 2 and t 1 , (f) b 2 and t 1.5 . In electroosmotic flows of viscoelastic fluids, it becomes important to consider how
the flow velocity varies with location, especially in the case with asymmetric wall zeta
22
ACCEPTED MANUSCRIPT potentials. Hence, we plot fig 12(a)-(c) with the intention of showing how the flow velocity
u t, y, z varies with at certain locations within the domain. All other parameters used are given clearly in the figure caption. We essentially observe from fig. 12 that the velocity variation with time is smooth for all locations. This is quite unlike fig. 2, where Jiménez et al. [28], studied a similar case for a Maxwell fluid. In case of Maxwell fluid, due to weak diffusive effects at small times, the velocity profile tends to be coarser with sharp oscillation.
CR IP T
Here, for an Oldroyd-B fluid, we have strong diffusive nature right from small times which makes the oscillations smoother. Also, the strong diffusive effects tend to damp and kill the oscillation faster. Fig 12(a), shows how the flow velocity varies for a case where equal and opposite zeta potentials are taken across from each other. We observe that at the center of the channel, there is almost no velocity developed through all times. We also observe that close
AN US
to the walls i.e. (0.9,-0.9), the flow velocity magnitude is lower, with almost no oscillatory effects. The flow velocity magnitude also virtually reaches the steady state value at a quicker time. This may be attributed to the viscous effects of the nearby wall. The flow velocity magnitude midway between the walls and the center i.e. at (0.5,-0.5) is large with more oscillatory nature and takes larger time to reach steady state. In fig 12(b), we see the temporal
M
variation at two different locations with varying . We observe from the figure that the magnitude of the velocity increases with decrease in , irrespective of the location.
ED
However, the damping of the oscillatory nature is more pronounced at (0.9,-0.9) as compared to at (0.5,-0.5). This may be due to the fact that (0.9,-0.9) is closer to the wall, where
PT
diffusive effects are stronger. Also, for larger values of , oscillations are damped faster because of strong diffusive nature of the fluid as compared to the elastic nature. The effect of
CE
El on the velocity magnitudes at two different locations has been considered in fig. 12(c). Here we observe that the velocity magnitude increases at all locations with increase in El, and
AC
also that the oscillatory nature becomes more damped closer to the wall at (0.5,-0.5) as compared to away from the wall at (0.9,-0.9). We also observe that the time taken for the first peak to appear and time taken to reach steady state also increases with increase in El, irrespective of the location.
23
CR IP T
ACCEPTED MANUSCRIPT
(b)
M
AN US
(a)
ED
(c) Fig 12 (color online) Plot showing the temporal variation of dimensionless velocity u t, y, z
CE
3 4 2.0 .
PT
for (a) various locations within the cross-section taking El 2.0 , 0.1 , b 1 and 10 , (b) various locations within the cross-section and varying while taking El 2.0 , b 1 and 10 , (c) various locations within the cross-section and varying El while taking 0.1 , b 1 and 10 . Other parameters chosen are: and 1 2 2.0 and
The number of parametric variations, considering the immense number of
AC
combinations of wall zeta potentials that is possible is very large. Hence, we would like to restrict our studies to the above considered cases, as we feel they give a broad picture of how an Oldroyd-B fluid would behave in an electroosmotic flow through a rectangular channel.
V. Conclusions The effect of combining electroosmotic flow with rheological behavior of an Oldryod-B fluid in a rectangular channel was studied thoroughly, with focus being on symmetric and asymmetric zeta potentials. The following are some of the major conclusions from this study:
24
ACCEPTED MANUSCRIPT
The mathematical method and the numerical model can be used to study cases with both symmetric and asymmetric wall zeta potentials.
The flow startup is strongly dependent on El and , where large El and small tends to make the fluid more elastic and less viscous in nature.
We observe that the aspect ratio of the channel plays a strong role in determining the flow transience behavior. The flow velocity profiles of an Oldroyd-B fluid is essentially smooth as compared to sharp and coarse profiles for a Maxwell fluid.
CR IP T
We strongly believe that the results from this study will enhance the understanding of the behavior of a viscoelastic fluid within microchannels actuated by eletroosmosis, with farreaching implications in various applications including microfluidics based medical
AN US
diagnostics.
References
N.T. Nguyen, S.A.M. Shaegh, N. Kashaninejad, D.T. Phan, Design, fabrication and
M
[1]
characterization of drug delivery systems based on lab-on-a-chip technology, Adv.
[2]
ED
Drug Deliv. Rev. 65 (2013) 1403–1419.
Á. Ríos, M. Zougagh, M. Avila, Miniaturization through lab-on-a-chip: Utopia or reality for routine laboratories? A review, Anal. Chim. Acta. 740 (2012) 1–11. P.C. Braga, M. Moretti, A. Piacenza, C.C. Montoli, E.E. Guffanti, Effects of seaprose
PT
[3]
on the rheology of bronchial mucus in patients with chronic bronchitis. A double-blind
[4]
CE
study vs placebo., Int. J. Clin. Pharmacol. Res. 13 (1993) 179–85. J.R. Stokes, G.A. Davies, Viscoelasticity of human whole saliva collected after acid
AC
and mechanical stimulation, Biorheology. 44 (2007) 141–160.
[5]
Y.-J. Kang, S.-J. Lee, Blood viscoelasticity measurement using steady and transient flow controls of blood in a microfluidic analogue of Wheastone-bridge channel, Biomicrofluidics. 7 (2013) 54122.
[6]
S. Chakraborty, Microfluidics and Microscale Transport Processes, CRC Press, 2012.
[7]
A.D. Stroock, M. Weck, D.T. Chiu, W.T.S. Huck, P.J.A. Kenis, R.F. Ismagilov, G.M. Whitesides, Patterning Electro-osmotic Flow with Patterned Surface Charge, Phys. Rev. Lett. 84 (2000) 3314–3317.
25
ACCEPTED MANUSCRIPT [8]
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.
[9]
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(5) (2000) 1053–1057.
Phys. Chem. 69(11) (2002) 4017–4024. [11]
CR IP T
[10] C.L. Rice, R. Whitehead, Electrokinetic Flow in a Narrow Cylindrical Capillary, J.
P. Dutta, A. Beskok, Analytical Solution of Time Periodic Electroosmotic Flows: Analogies to Stokes‟ Second Problem, Anal. Chem. 73(21) (2001) 5097–5102.
[12] R. Sadr, M. Yoda, Z. Zheng, A.T. Conlisk, An experimental study of electro-osmotic flow in rectangular microchannels, J. Fluid Mech. 506 (2004) 357–367.
AN US
[13] S. Devasenathipathy, J.G. Santiago, K. Takehara, Particle Tracking Techniques for Electrokinetic Microchannel Flows, Anal. Chem. 74(15) (2002) 3704–3713. [14] N.A. Patankar, H.H. Hu, Numerical Simulation of Electroosmotic Flow, Anal. Chem. 70(9) (1998) 1870–1881.
M
[15] L. Ren, D. Li, Electroosmotic Flow in Heterogeneous Microchannels, J. Colloid Interface Sci. 243 (2001) 255–261. Dutta,
A.
Beskok,
T.C.
Warburton,
Numerical
simulation
of
mixed
ED
[16] P.
electroosmotic/pressure driven microflows, Numer. Heat Transf. Part A Appl. 41(2) (2002) 131–148.
PT
[17] X.X. Li, Z. Yin, Y.J. Jian, L. Chang, J. Su, Q.S. Liu, Transient electro-osmotic flow of generalized Maxwell fluids through a microchannel, J. Nonnewton. Fluid Mech. 187
CE
(2012) 43–47.
[18] M. Zhao, S. Wang, S. Wei, Transient electro-osmotic flow of Oldroyd-B fluids in a
AC
straight pipe of circular cross section, J. Nonnewton. Fluid Mech. 201 (2013) 135–139. [19] S. Wang, M. Zhao, X. Li, Transient electro-osmotic flow of generalized Maxwell fluids in a straight pipe of circular cross section, Open Phys. 12 (2014) 445–451.
[20] P. Abhimanyu, P. Kaushik, P.K. Mondal, S. Chakraborty, Transiences in rotational electro-hydrodynamics microflows of a viscoelastic fluid under electrical double layer phenomena, J. Nonnewton. Fluid Mech. 231 (2016) 56–67. [21] P. Kaushik, P. Abhimanyu, P.K. Mondal, S. Chakraborty, Confinement effects on the rotational microflows of a viscoelastic fluid under Electrical double layer phenomenon, J. Nonnewton. Fluid Mech. 244 (2017) 123–137. 26
ACCEPTED MANUSCRIPT [22] C. Zhao, C. Yang, Exact solutions for electro-osmotic flow of viscoelastic fluids in rectangular micro-channels, Appl. Math. Comput. 211 (2009) 502–509. [23] G. Bin Lee, L.M. Fu, C.H. Lin, C.Y. Lee, R.J. Yang, Dispersion control in microfluidic chips by localized zeta potential variation using the field effect, Electrophoresis. 25 (2004) 1879–1887. [24] C.Y. Lee, G. Bin Lee, L.M. Fu, K.H. Lee, R.J. Yang, Electrokinetically driven active micro-mixers utilizing zeta potential variation induced by field effect, J.
CR IP T
Micromechanics Microengineering. 14 (2004) 1390–1398.
[25] A. Miller, A. Villegas, F. Javier Diez, Characterization of the startup transient electrokinetic flow in rectangular channels of arbitrary dimensions, zeta potential distribution, and time-varying pressure gradient, Electrophoresis. 36 (2015) 692–702. [26] A.M. Afonso, M.A. Alves, F.T. Pinho, Electro-osmotic flow of viscoelastic fluids in
AN US
microchannels under asymmetric zeta potentials, J. Eng. Math. 71 (2011) 15–30.
[27] J. Escandón, E. Jiménez, C. Hernández, O. Bautista, F. Méndez, Transient electroosmotic flow of Maxwell fluids in a slit microchannel with asymmetric zeta potentials, Eur. J. Mech. - B/Fluids. 53 (2015) 180–189.
M
[28] E. Jiménez, J. Escandón, O. Bautista, F. Méndez, Start-up electroosmotic flow of Maxwell fluids in a rectangular microchannel with high zeta potentials, J. Nonnewton. Fluid Mech. 227 (2016) 17–29.
ED
[29] E. Jiménez, J. Escandón, O. Bautista, F. Méndez, A note on “Start-up electroosmotic flow of Maxwell fluids in a rectangular microchannel with high zeta potentials” [J.
PT
Non-Newton Fluid Mech. 227 (2016) 17–29], 2016. [30] P. Kaushik, P.K. Mondal, S. Chakraborty, Flow dynamics of a viscoelastic fluid
CE
squeezed and extruded between two parallel plates, J. Nonnewton. Fluid Mech. 227 (2016) 56–64.
AC
[31] J.G. Oldroyd, On the Formulation of Rheological Equations of State, Proc. R. Soc. London A Math. Phys. Eng. Sci. 200 (1950).
[32] C. Fetecau, C. Fetecau, M. Kamran, D. Vieru, Exact solutions for the flow of a generalized Oldroyd-B fluid induced by a constantly accelerating plate between two side walls perpendicular to the plate, J. Nonnewton. Fluid Mech. 156 (2009) 189–201. [33] H. Qi, M. Xu, Stokes‟ first problem for a viscoelastic fluid with the generalized Oldroyd-B model, Acta Mech. Sin. 23 (2007) 463–469. [34] Y. Kang, C. Yang, X. Huang, Dynamic aspects of electroosmotic flow in a cylindrical microcapillary, Int. J. Eng. Sci. 40 (2002) 2203–2221. 27
ACCEPTED MANUSCRIPT [35] R.J. Huter, Zeta Potential in Colloid Science, Academic Press, London, 1981. [36] J.H. Masliyah, S. Bhattacharjee, Electrokinetic and colloid transport phenomena, John Wiley & Sons, 2006. [37] L.E. Rodd, J.J. Cooper-White, D.V. Boger, G.H. McKinley, Role of the elasticity number in the entry flow of dilute polymer solutions in micro-fabricated contraction geometries, J. Nonnewton. Fluid Mech. 143 (2007) 170–191. [38] M. Dehghan, A. Shokri, A numerical method for solving the hyperbolic telegraph
CR IP T
equation, Numer. Methods Partial Differ. Equ. 24 (2008) 1080–1093.
[39] S.A. Teukolsky, Stability of the iterated Crank-Nicholson method in numerical relativity, Phys. Rev. D. 61 (2000) 87501.
[40] K. Macounova, C.R. Cabrera, M.R. Holl, P. Yager, Generation of natural pH gradients in microfluidic channels for use in isoelectric focusing, Anal. Chem. 72 (2000) 3745–
AN US
3751.
[41] Y. Gao, T.N. Wong, C. Yang, K.T. Ooi, Two-fluid electroosmotic flow in
AC
CE
PT
ED
M
microchannels, J. Colloid Interface Sci. 284 (2005) 306–314.
28