Analytical solution for transient electroosmotic ﬂ ow in a rotating microchannel

An analytical solution is developed for the unsteady ﬂ ow of ﬂ uid through a parallel rotating plate microchannel, under the in ﬂ uence of electrokinetic force using the Debye – H¨uckel (DH) approximation. Transient Navier – Stokes equations are solved exactly in terms of the cosine Fourier series using the separation of variables method. The e ﬀ ects of frame rotation frequency and electroosmotic force on the ﬂ uid velocity and the ﬂ ow rate distributions are investigated. The rotating system is found to have a damped oscillatory behavior. It is found that the period and the decay rate of the oscillations are independent of the DH parameter ( k ). A time dependent structure of the boundary layer is observed at higher rotational frequencies. Furthermore, the rotation is shown to generate a secondary ﬂ ow and a parameter is de ﬁ ned ( b ( t )) to examine the ratio of the ﬂ ow in the y and x directions. It showed that both the angular velocity and the Debye – H¨uckel parameters are in ﬂ uential on the induced transient secondary ﬂ ow in the y direction. At high values of the Debye – H¨uckel parameter and the rotation parameter the ﬂ ow rates in the x and y directions are found to be identical. The analytical solution results are found to be in good agreement with the numerical method results and previously published work in this ﬁ eld.


Introduction
2][3][4][5][6] The thickness of the EDL is comparable with the size of the microchannel so the length scales of the electrostatic interactions enable the development of many highimpact technologies in the area of electrokinetic transport in microchannels.In such technologies, electroosmotic ow (EOF) is the prevalent phenomenon for the transport of small volumes of aqueous solutions when an electric eld is applied across a microchannel with charged walls. 7The electrical forces on the EDL near the interface of a solid surface and an electrolyte solution generate a uid ow. 8,9he characteristics of EOF in microchannels are typically examined by solving continuity and momentum equations in the presence of an electrical body force given by the Poisson-Boltzmann (PB) equation.1][12][13][14][15][16][17][18] This approximation can be used for small magnitudes of wall zeta-potential (<0.025V) which makes the PB equation linear, thus leading to an analytical and explicit expression of the potential prole.
Levine et al. 11 analytically solved the steady state problem of electrokinetic ow in a ne cylindrical capillary considering the DH approximation.This study provided a more precise calculation of the power requirements for pumping electrolytes through ne capillaries.Burgreen and Nakache 10 analytically studied electrokinetic laminar steady ow in a very thin rectangular channel without the DH approximation.Their analysis provided a solution for cases of high ionic energy and a small electrokinetic radius.Hsu et al. 12 theoretically studied the ow of an electrolyte solution through an elliptical microchannel to simulate the ow of a uid in a vein.They considered a time independent spatial variation for the liquid velocity with different types of boundary conditions on the channel wall.In their analysis, they didn't consider any assumption regarding the thickness of the double layer or the level of the electrical potential in order to present a more realistic description for biological systems.They also reported that the effects of boundary conditions are more apparent when the size of a microchannel decreases.EOF in a human meridian was studied by Sheu et al. 13 through a numerical model considering a steady state creeping ow for tissue uid motion.The tissue uids contain ions and nutrients and they interact with the blood in the capillary vessel through a complex nonlinear electro-osmosis transport process.EOF in microchannels in the presence of thermal effects has also attracted much attention.Chen 15 presented analytical expressions for the velocity and the temperature distributions of a fully-developed non-Newtonian uid ow in a microchannel subjected to heat ux at the boundary conditions.Das and Chakraborty 16 obtained closed form expressions for the velocity and temperature distributions of steady and fully developed ow of a power-law uid.They examined the effect of non-Newtonian parameters on ow behaviour.In all the above mentioned studies, electrokinetic ow (more specically EOF) was applied as the only means for controlling the ow of liquids in microuidic channels. 19A plug-like velocity prole along with easy controllability are the main characteristics of EOF in microchannels.However, electrokinetic pumping suffers from drawbacks so researchers have been directed toward techniques for generating a secondary ow.First, physicochemical properties, such as ionic strength and pH, may affect the pumping performance.For example, excessive Joule heating (the process by which the passage of an electric current through a conductor releases heat) prohibits pumping of liquids with high ionic strengths. 20Therefore, it may be difficult or impossible to use EOF to pump biological uids, such as blood and urine. 21Second, electrokinetic pumping requires high-voltage supplies which may have adverse safety effects.Third, electrokinetic pumping only works for uids without any trapped bubbles in the microchannels.Finally, due to Joule heating, EOF cannot be used for high ow rates (>1 mL s À1 ) in wide channels (>1 mm). 21A simple and safe alternative for controlling the ow of liquids in microuidic channels is based on centrifugal force which is not sensitive to the chemical properties of the liquid.This system can also be applied for liquids with trapped bubbles and a wide range of channel sizes.In a centrifugal system, the uid spatial and temporal controls have to be incorporated over the position of the uids in the micro-channels.
Duffy et al. 21conducted an experimental investigation on the effect of centrifugal force on the micro-ow transport through microscopic channels fabricated in a plastic disk.The uid was pumped outward using centrifugal force due to the rotation of the disk at 60-3000 rpm.They applied this method to pump biological uids, such as blood and urine, since the driving force for the uid ow was no longer sensitive to the physicochemical properties of the uids such as pH and ionic strength.Chang and Wang 22 investigated steady state rotating EOF over an innite plate and inside a rectangular channel using the DH approximation for the charge distributions.They have shown that, the speed of rotation and the electrokinetic width are the most important parameters affecting the steady ow in a microchannel.In addition, they have demonstrated that when the speed of rotation increases, axial and transverse ows vanish and the entire system forms a rigid body rotation.Xie and Jian 23 analysed the rotating EOF of power-law uids at high zeta potentials in a microchannel.They applied the nonlinear PB equation for the EDL potential distribution.Using the nite difference method, they have numerically computed the rotating EOF velocity proles for non-Newtonian uids.However, the velocity proles in their study are only computed for the steady-state conditions.Li et al. 24 continued this study for the third grade uids.They used the perturbation method to approximately obtain the solutions by assuming a slight non-Newtonian behaviour.Recently, Ng and Qi 25 presented a more detailed study on electro-osmotic ow in rotating rectangular microchannels.They considered the same steady state problem as in ref. 22 to represent the effect of the channel width and different combination of wall charges on the Ekman-EDL, ow structure and ow rate.
In this study an analytical solution is developed for transient electroosmotic ow in a rotating microchannel.To the best of the authors' knowledge, this time dependent behaviour has not been reported in an exact closed-form series solution.First, a general solution of the problem is obtained via the superposition and separation of variables techniques.By using a cosine Fourier series expansion, the governing complex partial differential equation (PDE) was solved to yield a closed-form solution.Next, the analytical solutions for the velocity and ow rate were applied to some examples.Structures of the boundary layer at higher rotational frequencies were studied and the rotational induced secondary ows are discussed.Moreover, transient secondary ow in the y direction is shown to be dependent on angular velocity and the Debye-Hückel parameters.Finally, the analytical solution results were compared with numerical solutions from Mathematica and previously published results.

A. Electro-osmotic rotating microchannel ow equations and boundary conditions
The geometry for the present analysis is shown in Fig. 1.It consists of two parallel plates separated by a distance of 2h.The bottom and top plate are located at z ¼ Àh and z ¼ h, respectively.
The basic eld equations governing the ow of an incompressible uid between two parallel walls of a microchannel are the continuity equation and the modied Cauchy equation given by where u* ¼ (u*, v*, 0) is the velocity vector, r is the uid density, t* is the time, U ¼ (0, 0, U) is the angular velocity vector, p is the modied pressure including centrifugal force (p ¼ P À r/2|U Â r| 2 ), s is the extra stress tensor, and b is the body force per unit volume.The pressure gradient along the microchannel is assumed to be zero except for the centrifugal force that is included in the modied pressure.In EOF (except for in electrophoresis applications with high electrical elds 26 ) the Reynolds number is less than unity which results in a negligible convective term in the NS equation (u*$Vu* ¼ 0). 27A list of the constants and parameters used here is presented in Table 1.
To nd the electric potential distribution within the EDL the Poisson-Nernst-Planck (PNP) model is considered.The species conservation equation can be written as: where n * þ and n * À are the concentration numbers of the free charges in the uid, e is the elementary charge, m AE is the mobility, j* is the electric potential in the uid, k B is the Boltzmann constant, and T is the absolute temperature.By subtracting the positive and negative number concentrations the following equation is obtained: The electrical conductivity of the electrolyte (s) is found based on the molar conductivity of the solution, 8 s ¼ 1000Ms M .Substituting eqn (3) into eqn (2) yields: 8,28-30 where 3 is the electric permittivity and s c ¼ 3/s and s p ¼ mh 2 /3j* 2 are the charge relaxation time scale and the process time scale, respectively.Here, based on the parameters presented in Table 1, the process time is higher than the charge relaxation time (s p [ s c ), the dynamics of the charges become negligible compared to the uid motion and the PNP equation simplies to the Poisson-Boltzmann equation.The net electric charge density in the uid and the electric potential eld j*(z) within the EDL are related as follows: The term, r e is the Boltzmann distribution of the net electric charge density near to the charged walls described as 31 where n 0 is the bulk concentration of the electrolytes in the liquid (n 0 ¼ 1000N A M), e is the electron charge, and x (¼1) is the valence.The body force per unit volume, dened in eqn (1.b) can be written as b ¼ (r e E, 0, 0).
In the present study, it is assumed that the wall electrical potential is small (j ¼ 10 # 25 mV) compared to the thermal potential of the charged species (|exj| # k B T). Hence, the DH linearization principle (sinh(A) z A) 8,32 can be applied to obtain the following simpler linear relation: where k* ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2n 0 e 2 x 2 =3k B T p is the DH parameter and it is the inverse of the Debye length (l D ).It should be noted that, the DH approximation can be used when the Debye thickness is small but nite, i.e., for h/l D # 10 3 . 32caling factors are dened to simplify the analysis procedure as follows: Applying these scaling factors to eqn (1.b) and ( 7) along with considering the no-slip (zero velocities) boundary conditions on the channel walls, the following dimensionless governing equations are obtained in the x and y directions: where u ¼ rU/mk 2 is the dimensionless angular velocity.In the electrokinetic part, one can solve the dimensionless form of eqn (7) subjected to the boundary conditions of j(AE1) ¼ 1 as follows: To facilitate solutions of eqn (9.a) and (9.b), a complex function of c(z, t) ¼ u(z, t) + iv(z, t) with i ¼ À1 was used.Hence, the hydrodynamic equations (eqn (9.a) and (9.b) can be written as a single complex partial differential equation as follows:

B. Solution procedure
To nd the 'analytical' solution of the complex PDE described in eqn (11), the superposition property is used: where, c 1 (z) is the steady state component of the velocity prole and c 2 (z, t) is the transient component.Introducing eqn (12)  into eqn (11) while imposing the boundary and initial conditions, the following set of differential equations are obtained for c 1 (z): and for c 2 (z, t) Applying the corresponding boundary conditions to the differential equation eqn (13), yields the steady state component of the velocity prole as 22 follows: Eqn ( 14) is solved using the separation of variables method by assuming c 2 (z, t) ¼ Z(z)T(t), with Z(z) and T(t) being the spatial and temporal functions, respectively.
The new form for c 2 (z, t) is substituted into eqn ( 14) to obtain the following separated ordinary differential equations (ODEs) for z and t: where, l 2 is the separation constant.
A solution of eqn ( 16) for Z(z) satisfying its boundary conditions Z(1) ¼ Z(À1) ¼ 0 can be assumed to take the form Z(z) ¼ cos(l n z) with l n ¼ (2n + 1)p/2.Also, the exact solution for T(t) in eqn ( 16) is as follows: where C n is a constant.Hence the expression for c 2 (z, t) becomes: Next, the initial condition requires c 2 (z, 0) ¼ Àc 1 (z) and the coefficient C n is determined as Finally, the solution for c(z, t) is: Integrating eqn (19) with respect to z leads to the following expression for the normalized complex volume ow rate per width: Terms, Q x (t) and Q y (t) are the transient ow rates in the x and y direction, respectively.Finally, b(t) ¼ tan À1 [Q y (t)/Q x (t)] is dened as a parameter representing the rotationally induced transient secondary ow in the y direction.
The solution represented in eqn ( 19) goes through three stages of development.Firstly, when ut ( 1, the Coriolis force in the momentum eqn ( 9) is small compared to the stress gradient.Using the relationship e iut ¼ cos(ut) + i sin(ut) at this stage, eqn ( 19) is re-written as follows: At this stage, the solution is related to a non-rotating transport accelerated by a constant stress.Next, at longer times, the Coriolis force deects the ow towards steady-state conditions.The solution undergoes damped inertial oscillations near the steady-state conditions.Finally, aer a long time, for ut [ 1 momentum completely diffuses to the microchannel ow.At this stage, the velocity eld in the microchannel is at the nal phase of inertial oscillations, which are damped by the viscous stress gradient.Therefore, the inertial oscillations are dissipated and one may use the rst unsteady term of solution in eqn (19) for a long time solution, as represented in the literature. 22

Results and discussion
To illustrate the nature and general behaviour of the analytical solution, some numerical examples are considered in this section.The variation of the dimensionless velocity eld (u(z, t)) as a function of time for different values of the DH parameter (k) and dimensionless angular velocity (u) are discussed using contour plots.Then, the effects of k and u on the evolution of rotationally induced secondary ow in the y direction are investigated by examining b(t).The uid experiences a force imparted due to EOF and the frame rotation inuences the ow velocity eld with centrifugal force.In order to nd the overall transient signature of the ow both of these two factors and their relative strengths must be included.However, we note that the oscillation characteristic of the ow, is appropriately described using the rotation parameter (u), apart from the DH parameter (k).In all the numerical results, a truncation constant of n ¼ 5 was found to yield satisfactory results.
2-D dimensionless contour plots of velocity evolution (xcomponent, u(z, t)) and the centre-line velocity prole versus time are shown in Fig. 2. To investigate the inuence of rotation, different values for the dimensionless rotation parameter u ¼ 1.0, 2.5, 5.0, 10 and 20 are used.The transient centre-line velocity is also plotted to have a better comparison of the velocity magnitude and behaviour at the centre of the channel at various u values.In all these gures, the horizontal axis is the dimensionless time parameter (0 < t < 3).
Results show that, as the frequency of rotation increases, the oscillatory behaviour of the ow eld increases too.Generally, as is observed in Fig. 2a-e, the values and ranges of the velocity in the x-direction are reduced as the angular velocity increases.This is attributed to the effects of u on the initiation of the secondary ow or the y-component of the ow. 22It should be noted that, the velocity is positive in the whole cross-section of the channel when u ¼ 1, whereas when increasing the u, negative values are also observed in the velocity contour plots.This could be related to the fact that higher rotational centrifugal force pushes the ow back in the x-direction.The reduction of the centre-line velocity, as shown in Fig. 2f, conrms these observations.As is seen, the velocity at the centre is reduced when the frequency of rotation increases and becomes almost zero at u ¼ 20.Furthermore, the position of the maximum velocity moves from the centre of the microchannel toward the walls when u increases.The same trend was also reported by Chang and Wang. 22n Fig. 3, the transient 2-D dimensionless velocity u(z, t) and the centre-line velocity prole are plotted for k ¼ 5 for the same range of angular velocities as in Fig. 2. As is observed, the oscillatory behaviour of the uid ow didn't change when increasing the DH parameter.According to eqn (17), the period of the velocity oscillations is equal to p/u, regardless of the DH parameter.At higher rotational frequencies, the structure of the boundary layer is observed.The boundary layer structure is found to be time dependent in the present work which was not noticed in earlier studies. 22,25owever, it is clear that the velocity magnitude increased when increasing the DH parameter.At higher k values the effect of the induced electroosmotic forces (i.e., body force per unit volume on the uid) becomes more inuential, thus increasing the uid velocity in the x-direction.Also, the effect of u on u(z, t) is dampened at higher k values.It becomes more evident when we compare the ratios of the maximum velocity represented on the range bars at different angular velocities.For example, u max      The time interval is 0 < t < 3. The velocity in the y-direction is also observed to follow an oscillatory behaviour.However, the most important distinction between u and v is that for all the positive values of u and k, the sign of the velocity in the ydirection is always negative.This means that the direction of the induced secondary ow is always Ày.
The variations of the y-component velocity and centre line velocities for k ¼ 5 and k ¼ 1000 were also investigated and the results are presented in Fig. 6 and 7, respectively.The most important distinctions are as follow.
Higher rotation frequency leads to more oscillatory behaviour of the ow velocity.Also, velocity gets some positive values when the DH parameter increases.This is directly linked to the effect of the driving force of the electroosmotic forces at high k values on the y-component of the velocity eld.Furthermore, as the rotation frequency increases the probability of boundary layer initiation increases too which is in full agreement with previous reports. 22,25ig. 8 shows the evolution of b over time (0 < t < 3) for k ¼ 1000 with different dimensionless angular velocities, u ¼ (a) 0.1, (b) 1.0, (c) p, and (d) 5.0.It is indicated that u has a noticeable effect on the oscillatory behaviour of the ow.Smaller values of u led to an almost non-oscillatory behaviour for b(t) which was seen before for the velocity behaviour.As shown, the value of b(t) increases signicantly with the increase of u.This means that, physically, the secondary induced ow is enhanced as the frequency of rotation increases.Also, the secondary induced ow parameter varies most rapidly with smaller values of u which is in full agreement with the ndings of Chang and Wang. 22ig. 9 shows the evolution of b over time (0 < t < 3) for u ¼ 5 with different DH parameters, k ¼ (a) 1.0, (b) 5.0, (c) 50, and (d) 1000.It is clear that increasing the DH parameter (k) leads to a decrease in the value of b(t).For the large values of k, b(t) gradually approaches 45.This means that for large k, equal ows in both the x and y directions are obtained which matches well with the Chang and Wang calculations to generate Fig. 8 in their study. 22The governing equations are also solved numerically using Mathematica.As can be seen in Fig. 8 and 9, the exact analytical solution results are in very good agreement with those obtained through numerical computation.Moreover, as indicated in Fig. 8 and 9, it should be pointed out that the effect of the angular velocity gradually fades as the momentum diffusion increases which is also observed in Fig. 1.As a result, b(t) approaches its steady state value, as reported by Chang and Wang, 22 when the time is long enough.

Conclusions
In the present paper, the dynamics of rotationally induced secondary ow in the EOF between two parallel plates have been analytically studied.A general series solution expression for the non-dimensional velocity eld and ow rate has been derived through the separation of variables technique.The effect of the Debye-Hückel parameter (k) and the dimensionless angular velocity (u) on the velocity contour plots and parameter b(t) has been examined.The main contribution of this study is providing a closed-form solution for the transient EOF in a rotating microchannel.It is found that the oscillatory behaviour of the uid ow is only related to u while the DH parameter only changes the amplitude of the velocity eld.
A time dependent structure of the boundary layer was observed at higher rotational frequencies.The effects of u on the initiation of the secondary ow were examined.Furthermore, the b(t) parameter was investigated to show the effect of the angular velocity and the Debye-Hückel parameter on the induced transient secondary ow in the y direction.It is found that when the Debye-Hückel parameter and the rotation parameter are high, b(t) oscillates near 45 degrees implying identical ow rates in the x and y directions.Moreover, excellent agreement was found between the analytical results and the numerical results obtained using Mathematica.

Fig. 1
Fig.13D view of the rotating microchannel with EOF.
mentioned in the problem formulation section.The ratio of u max at u ¼ 1.0 over u max at u ¼ 2.5 in this case is about 1.0.From Fig.2-4 and eqn (17) it can be concluded that the period of the velocity oscillations is constant (p/u) and the rate of decay of oscillations is l n 2 which is independent of the DH parameter.The inertial oscillations are dissipated (i.e., third stage of development) through momentum diffusion.The results related to the transient electroosmotic ow in a nonrotating microchannel are plotted for three DH parameters (k ¼ 1, 5 and 1000) in Fig. 2f and 4f.As is observed, there is an

Table 1
Constants and parameters