Investigation of the effects of time periodic pressure and potential gradients on viscoelastic fluid flow in circular narrow confinements

In this paper we present an in-depth analysis and analytical solution for time periodic hydrodynamic flow (driven by a time-dependent pressure gradient and electric field) of viscoelastic fluid through cylindrical micro- and nanochannels. Particularly, we solve the linearized Poisson–Boltzmann equation, together with the incompressible Cauchy momentum equation under no-slip boundary conditions for viscoelastic fluid in the case of a combination of time periodic pressure-driven and electro-osmotic flow. The resulting solutions allow us to predict the electrical current and solution flow rate. As expected from the assumption of linear viscoelasticity, the results satisfy the Onsager reciprocal relation, which is important since it enables an analogy between fluidic networks in this flow configuration and electric circuits. The results especially are of interest for micro- and nanofluidic energy conversion applications. We also found that time periodic electro-osmotic flow in many cases is much stronger enhanced than time periodic pressure-driven flow when comparing the flow profiles of oscillating PDF and EOF in micro- and nanochannels. The findings advance our understanding of time periodic electrokinetic phenomena of viscoelastic fluids and provide insight into flow characteristic as well as assist the design of devices for lab-on-chip applications.


Introduction
Micro-and nanofluidic applications (e.g., on-chip bioanalysis, on-chip diagnostic devices, DNA molecules separation, energy harvesting, and so on) require the transportation of fluids to be driven by an external driving force, which can be either a pressure gradient [pressuredriven flow (PDF)] or an external electric field [electroosmotic flow (EOF)] or the combination of these two driving forces. Force application results in the coupled flow of matter and ionic current, so-called electrokinetic flow. Based on the physical problem of interest, these driving forces can be steady or time-dependent. The application of steady driving forces for Newtonian fluids, like aqueous electrolyte solutions, whose viscosity is constant, was extensively investigated in the past (Masliyah and Bhattacharjee 2006;Bruus 2008). Recently, the necessity of manipulation of biofluids (for example blood, DNA solutions) and polymeric liquids in small confinements has triggered a renewed interest in the dynamics of non-Newtonian fluid. Berli theoretically studied the utilization of steady PDF (Berli 2010a), steady EOF (Olivares et al. 2009;Berli 2010b), and steady combined PDF-EOF (Berli and Olivares 2008) for inelastic non-Newtonian fluids using a power law constitutive equation in both rectangular and cylindrical microchannels. Experiments carried out for steady PDF non-Newtonian flow in a rectangular microchannel inspired by Berli's theory were also reported Abstract In this paper we present an in-depth analysis and analytical solution for time periodic hydrodynamic flow (driven by a time-dependent pressure gradient and electric field) of viscoelastic fluid through cylindrical micro-and nanochannels. Particularly, we solve the linearized Poisson-Boltzmann equation, together with the incompressible Cauchy momentum equation under no-slip boundary conditions for viscoelastic fluid in the case of a combination of time periodic pressure-driven and electro-osmotic flow. The resulting solutions allow us to predict the electrical current and solution flow rate. As expected from the assumption of linear viscoelasticity, the results satisfy the Onsager reciprocal relation, which is important since it enables an analogy between fluidic networks in this flow configuration and electric circuits. The results especially are of interest for micro-and nanofluidic energy conversion applications. We also found that time periodic electro-osmotic flow in many cases is much stronger enhanced than time periodic pressure-driven flow when comparing the flow profiles of oscillating PDF and EOF in micro-and nanochannels. The findings advance our understanding of time periodic electrokinetic phenomena of viscoelastic fluids and provide (Nguyen et al. 2013). Chakraborty and colleagues have theoretically studied transport of non-Newtonian fluid (inelastic power law fluids and recently viscoelastic constitutions) using separately steady PDF (Bandopadhyay and Chakraborty 2011), steady EOF (Chakraborty 2007;Ghosh and Chakraborty 2015), time periodic PDF (Bandopadhyay and Chakraborty 2012a, b; Bandopadhyay et al. 2014) and time periodic EOF (Bandopadhyay et al. 2013) in rectangular narrow confinements. Afonso et al. studied the combined steady PDF and EOF using two different viscoelastic fluid models, namely the Phan-Thien-Tanner (PTT) model and the finitely extensible nonlinear elastic with a Peterlin approximation (FENE-P) model (Afonso et al. 2009). Dhinakaran et al. (2010) studied the steady EOF for viscoelastic fluids using the PTT model and nonlinearity of the Poisson-Boltzmann equation. Liu et al. studied time periodic EOF of viscoelastic fluid in rectangular (Liu et al. 2011a), cylindrical (Liu et al. 2012) and semicircular microchannels (Bao et al. 2013). However, so far no author discussed on the time-dependent combined PDF-EOF of viscoelastic flow in a narrow confinement (micro-and nanochannels). In this context, our work aims to fill this gap by attempting to investigate the theoretical relations between fluxes and forces for time periodic electrokinetic (mixed PDF-EOF) flow of viscoelastic fluid in narrow confinement. It is important to note that knowing the relationships between driving forces and conjugate fluxes in electrokinetics [which for simple Newtonian fluid and steady mixed PDF-EOF can be described by transport equations and the Onsager relations of non-equilibrium thermodynamics (Masliyah and Bhattacharjee 2006)] is a crucial aspect for miniaturization and integration. It is thus relevant for the design and operation of micro-and nanochannels in fluidic networks (lab-on-chip platforms) as well as for understanding the underlying fundamental physics of fluids. The results are also of interest for energy conversion in micro-and nanofluidic systems.

Theoretical model
We consider the flow of a linearized Maxwell fluid in an infinity long circular micro-or nanochannel (with channel radius R) under application of both an oscillating pressure and electric field using a cylindrical coordinate system (Fig. 1).

Potential distribution
When the charged channel surface is in contact with the fluid with dissolved ions, electrical double layers (EDL) are formed at the channel walls. The electrical potential (ψ) in the EDL is a function of r in cylindrical coordinate system and has the non-dimensional form as: in which the non-dimensional quantities are as follows: r = r R , ψ = ψ ζ , R = R or when converted back to dimensional quantities, in which κ = 1 and = ǫk B T 2n 0 z 2 e 2 · ζ is the zeta potential. n 0 is the bulk ionic density, k B is the Boltzmann constant, T is the operational temperature, e is the elementary charge, ϵ is the permittivity of the fluid, and z is the valency of the positively and negatively charged species (for a symmetric electrolyte, z + = −z − = z).
This model is classical for electrical double layers when we do not consider finite ionic size effects. A detailed model description on the effects of finite ionic size and solvent polarization for electrical double layers is beyond the scope of this work but can be found in . It is noticed that the electrical potential causes by EDL is normal to the wall and the convection is parallel to the wall, so there is no disturbance of the EDL potential.

Fluid velocity
The flow is governed by the incompressible Cauchy's momentum equation. Considering the flow in z direction (unidimensional flow), the scalar momentum equation can be expressed as: with ρ, the fluid density; u(r, t), the fluid velocity; − ∂ ∂z p(z, t), the applied pressure gradient; τ(r, t), the stress tensor; and E(z, t) the externally applied electric field.
(1) It is important to note that E(z, t) in Eq. (5) includes two components: (1) the induced electric field by the applied pressure gradient E S e −iωt (the streaming potential field) and (2) the applied electric field E A e −i(ωt+ϕ) . Here, ϕ is the phase difference between the applied pressure gradient and the applied electric field. We now define E 0 as: Viscoelastic behavior is presented using the linear Maxwell model.
where t n is the liquid relaxation time and η is the liquid viscosity.

Flow rate
The flow rate q = ℜ(Qe (−iωt) ) in which the flow rate amplitude has the form: By integrating and taking − d dz P = P L and E 0 (z) = d dz Φ = �Φ L , the complex flow rate Q amplitude has the form: The flow rate amplitude as shown in Eq. (11) is composed of two parts. The first part is driven by the applied oscillating pressure, and the second part is driven by the applied oscillating electrical field.

Ionic current
The ionic current i cur = ℜ(Ie (−iωt) ) , in which the current amplitude has the form: Here, f is the Stokes-Einstein friction factor, f = k B T D and D is the diffusion coefficient, σ s is the conductivity of the Stern layer (Masliyah and Bhattacharjee 2006;Lee et al. 2012;Davidson and Xuan 2008). It is important to note that since we use a linear viscoelastic model, the f factor presented here is not dependent on the power law exponent [denoted as β in ] which is solely used for a power law (inelastic) fluid. For more discussion on the f factor in case of using an inelastic fluid, please refer to .
Changing to the non-dimensional variable r, we obtain: By substituting the velocity given by Eq. (9) into Eq. (12) and integrating, we obtain the complex amplitude current: At this point, the velocity profiles expressed in Eq. (6) for both oscillating pressure-driven and electro-osmotic flows are fully determined. Equation (16) is used for plotting velocity amplitudes as shown in the following sections.

Onsager's reciprocal relations
The Maxwell model for viscoelastic fluid is restricted to small deformations so that the fluid responds linearly. This phenomenon is known as linear viscoelasticity.
Because of this linear relation, the Onsager relations are expected to be obeyed (Onsager 1931a, b;Lebon et al. 2008;Rajagopal 2008) and indeed, we find that the complex flow rate amplitude and complex ionic current amplitude in Eqs. (11) and (13) can be re-written as follows: The transport Eq. (17) shows that flow rate amplitude Q and ionic current amplitude I are linear with applied pressure and electric potential amplitudes. L ij in Eq. (17) are phenomenological coefficients. In particularly, L 11 characterizes the hydraulic conductance and L 22 characterizes the electric conductance. L 12 characterizes the electro-osmosis and L 21 characterizes the streaming potential effect. Onsager's reciprocal relation is complied with if L 12 = L 21 . We see that this relation is indeed Here, Du is the Dukhin number and Du = σ s Rσ b , σ b is the conductivity of the bulk solution.
As with the flow rate amplitude, the current response of the system is caused by the oscillating pressure (the first term) and the oscillating electrical field (the second term).

Consideration of streaming potential and applied electric field
By substituting Eq. (7) into Eq. (9), the complex velocity amplitude can be written as: The velocity field therefore can be viewed as the superposition of the velocity fields caused by (1) the pressure gradient coupling with its streaming potential field [the first and the second terms on the right-hand side of Eq. (14)] and (2) the applied electric field [the third term on the right-hand side of Eq. (14)]. In this context, if one considers solely pressure-driven system, where no electric field is applied E A e −iϕ = 0, the streaming potential E S = E 0 . Since the total ionic current at maximal streaming potential is zero, this gives us the opportunity to extract the relation between U refP and U refE S from Eq. (12) as following (by taking I = 0): The velocity amplitude U(r) can therefore be expressed solely as a function of U refP and U refE A as: satisfied because from Eqs. (11), (13) and (17) it is obvious that: Equation (17) can be used to construct an analogy between micro-and nanofluidic channel networks and electric circuits because it describes the electrokinetic phenomenon as a generalization of Ohm's law where linear relations between currents (of mass or charges) and applied gradients (voltage or pressure) occur (Ajdari 2004;Campisi et al. 2006). In this context, it is interesting to apply our calculation results to examine the energy conversion efficiency of the streaming potential energy harvesting system in a manner comparable to the work of Bandopadhyay and Chakraborty (2012a).

Streaming potential energy harvesting
The electrokinetic energy conversion efficiency (Eff) in a microchannel for a Newtonian fluid under steady pressuredriven flow was theoretically predicted to be less than 1% (Morrison and Osterle 1965), while for an inelastic polymer it was predicted to be about 1% (Berli 2010a). In a nanochannel, for a Newtonian fluid under no-slip boundary conditions and based on a Poisson-Boltzmann charge distribution, the theoretical prediction of energy conversion efficiency is up to 12% (van der Heyden et al. 2006). Recently, Bandopadhyay and Chakraborty (2012a) gave a valuable contribution to the theory of electrokinetic energy conversion by taking into account the utilization of Maxwell viscoelastic fluid and oscillating pressure-driven flow in slit micro-and nanochannels. Bandopadhyay et al. showed that for a slit-type microchannel ( H = 500, with H the half channel height and λ the Debye length), the conversion efficiency can be in the order of 10%, and that for a nanochannel ( H = 10) without taking into account surface conductance, the conversion efficiency can be even larger than 95% [see Fig. 1 and S3 in ref. Bandopadhyay and Chakraborty (2012a)]. Our calculation results for a cylindrical geometry show that an efficiency can be obtained in the same order for the case of a microchannel and that the maximum efficiency can be larger than 95% for a nanochannel (Fig. 2). For the purpose of comparison, plots are constructed using the same input data as provided by the work of Bandopadhyay and Chakraborty (2012a) (i.e., ϑ = 10 −4 , ζ = −1, Ω = −10, Du = 0).
It must be remarked that the maximal efficiencies shown in Fig. 2 and those predicted by Bandopadhyay et al. are thermodynamic efficiencies [Eff = I S �φ Q(�p) ], i.e., in the case no power is delivered by the system. For practical purposes, the maximal conversion efficiency under the condition of maximal output power at a load resistor is more relevant (Olthuis et al. 2005), Eff max = 1 4 I S �φ Q(�p) . Figure 3 shows that the maximum efficiencies at maximal output power are 24.3 and 7.7% for a cylindrical nanochannel (R = 15) and microchannel (R = 500), respectively. These values though much lower than the thermodynamic efficiencies are still much higher than the predictions for conventional systems using DC actuations and Newtonian fluids cited above, especially for microchannels.

Understanding the mechanism
In the work of Bandopadhyay and Chakraborty (2012a), the mechanism behind the massive enhancement of the energy conversion efficiency using viscoelastic fluid was not in detail described. Herewith, we will provide a description of the mechanism that enhances the efficiency. Figure 4 shows the maximal thermodynamic energy conversion efficiency following ω * and the inverse Deborah number ϑ[here, ϑ = ρR 2 ηt n , Bandopadhyay and Chakraborty (2012a)] for a nanochannel at R = 5 [in this context, for the comparison with the work of Bandopadhyay et al., the Deborah number is defined as De = ηt n ρR 2 . It is noticed that some other authors have also defined De = ωt n (Bao et al. 2013)]. It is obvious from Fig. 4 that in the limit ϑ → 0(high relaxation time, elastic dominant zone), the efficiencies are high, while at high ϑ, (low relaxation time, viscous dominant zone), no efficiency peaks appear. This behavior can be explained from the linear Maxwell viscoelastic model that presents fluid as a serial connection between a spring (elastic behavior) and a dashpot (viscous behavior). The closer to the dominantly elastic zone (lower ϑ), the more the fluid behaves as a Hookean solid in responding (large relaxation time), resulting in a shift of the resonant peak toward the higher ω * values. When ϑ → 0, at resonant frequencies, the fluid inside the channel exhibits an entirely elastic character and hence moves frictionlessly, as a result providing high conversion efficiencies. The peak locations at which maximal efficiencies are observed depend on the oscillation frequencies that are also determined by the channel dimension. This can be seen when ϑ is constant (10 −4 ), the maximal efficiency peaks shift to smaller frequencies at an increase in channel dimensionless radius, R (shown in Fig. 5). This frequency shift was also observed in the work of Bandopadhyay and Chakraborty (2012a, b). Furthermore, the peaks also split into two separate peaks so that they can be shifted to smaller frequencies when increasing the channel radii (for example the peak at ω * approximate 550 and R = 100 in Fig. 5).

Oscillating pressure-driven flow profile
For the sake of generality, all the plots are presented using the non-dimensional quantity: 4η , see Eq. (16) for U(r). Figure 6 shows the oscillating pressure-driven flow profile of viscoelastic fluid following ω * and channel radius r at R = 20, ϑ = 10 −4 , ζ = −1, Ω = −10, Du = 0. In order to compare with the case of oscillating electro-osmotic flow, the velocity amplitude is also plotted and shown in Fig. 7.
It is important to stress that while the pressure gradient − ∂p ∂z and the velocity u(r, t) appear to have the same oscillatory form in the time variable t [see Eqs. (4) and (6)], this does not mean that they actually are in phase. The reason for this is that the other part of the velocity, namely, the U(r) or U(r) is a complex quantity. The product of this complex quantity with e −iωt as shown in Eq. (6) causes changes in the phases of the real and imaginary parts of the U(r) or U(r) and hence of the velocity u(r, t) so that a phase shift will occur with respect to the pressure gradient − ∂p ∂z .

Complex and real velocity amplitude
The velocity u has the form Since U is a complex number, we can express it as: Substituting Eq. (19) into Eq. (20) and isolating the Real part, we have: is the phase shift, and hence U c = |U| is the (real) velocity amplitude (Moyers-Gonzalez et al. 2009), see Fig. 7. Figure 8 shows the phase shift of the velocity following the dimensionless pressure frequency (ω * ) with two different values of ϑ. It can be seen that depending on the values of ϑ, the phases pass from negative (viscous zone) to positive (elastic zone) (Moyers-Gonzalez et al. 2009). The green line represents the phase for Newtonian dominant fluid (ϑ = 10 10 ) and stays in viscous zone (negative). As for ϑ = 10 −4 (the blue curve), the phase is in the elastic zone (positive) at low frequency. As the frequency increases, the fluid responds viscously indicating by the changing of the blue curve from positive to negative zone. When the frequency further increases and reaches the resonant frequency, the phase shifts back to the elastic zone (positive). At resonant frequency, the fluid behaves elastically and hence moves frictionlessly, as a result providing high energy conversion efficiencies as mentioned in previous section.

Effectiveness of electro-osmotic flow compared to pressure-driven flow
It can be seen from Figs. 6,7 and Figs. 9, 10 that at resonant frequencies, the maximal velocity in the case of oscillating EOF is much higher than in the case of oscillating PDF even though at low frequencies these flows have the same maximal velocities (see Fig. 11).
In the textbook, for DC electrokinetic flow, the concept of effectiveness (B) of electro-osmotic flow as compared to pressure-driven flow is given by the ratio of volume flow rate, see page 244 of ref. Masliyah and Bhattacharjee (2006).
In our case, for time periodic electrokinetics with a Maxwell fluid, the volume flow is expressed by Eq. (11). The effectiveness B therefore has the form: Figure 12 shows the frequency-dependent effectiveness of oscillating EOF over oscillating PDF. It is clear that at the resonant frequencies, the effectiveness of oscillating EOF is much higher than oscillating PDF, while at small frequencies, effectiveness is equal (as also evident from Fig. 11). Furthermore, in nanochannels, the effectiveness is much more strongly increased than in microchannels. This observation could be explained by noticing that we have the like-standing waves in the channel (see Figs. 6,9). For oscillating PDF, the applied pressure force is exerted over the entire cross section of the channel. This flow behavior allows all energy to be coupled into the actuation in one direction (for example first harmonic, the peak around ω * = 250, see Fig. 6). For the first harmonic of oscillating EOF (see Fig. 9), also all energy is coupled in one direction; hence, both have equal effectiveness at low ω * . However, with the third harmonic (the peak around ω * = 500), the situation is quite different. As with oscillating PDF, the pressure force in the center of the channel is directed against the direction of the movement; hence, the center velocity is  lower than in the first harmonic. With oscillating EOF, there is no force exerted in the center of the channel, but only in a thin layer at the wall. Hence the force exerted in the wide area close to the walls can be coupled to the much narrower area at the center. This concentration of energy in a small cross section (especially for nanochannel) causes strong increase in velocity in the center, hence much higher effectiveness than oscillating PDF. The question can be posed whether the high velocities generated will not disturb the electrical double layer composition. It is important to realize that our model concerns an infinitely long channel of constant fluid properties and homogeneous wall charge density. In this channel the potential and ionic composition in the electrical double layer only vary in the direction normal to the channel wall. Only when turbulence occurs, the double layer composition will hence be disturbed. The Reynolds number in our case is Re = ω * R2 2 ρ t n η (Jian et al. 2010;Liu et al. 2011b). For the optimal dimensionless parameter values as found in this work namely R = 10, ω * = 250 , and the practical values mentioned in the work of Bandopadhyay and Chakraborty (2012b), ρ = 10 3 kg/m 3 , t n = 10 −2 , η = 10 −3 Pa s, we find that Re = 2.5 × 10 11 2 . Since Debye length λ is always below 1 µm, turbulence is not expected. From practical point of view, in future experimental systems, the interfacing to an electrical system would need to be considered. This would involve electrode/solution interfaces with local storage and exchange of charge and possibly channel openings. At every interface where an inhomogeneity of flow or fixed charge concentration would occur, conservation of charge and matter would give rise to local gradients of electrical field, pressure and/or concentration. This would cause additional losses that would need to be considered in the design of such systems. One single aspect of the interfacing, namely the disturbances of the electrical double layer composition by advective fluxes can be estimated in isolation. By comparing the advective flux parallel to the wall, disturbing the electrical double layer composition, with the restoring diffusion flux normal to the wall, restoring equilibrium, we can estimate the severity of the disturbances in double layer composition. The ratio of the two fluxes provides a Péclet number, Pe = ω * R 2 Dt n . For R = 10, ω * = 250, D = 10 −9 m 2 /s and t n = 10 −2 s, we find Pe = 2.5 × 10 14 λ 2 . For λ < 60 nm, Pe < 1 and diffusional equilibration will be sufficiently rapid.

Conclusions
We report for the first time an analytical solution for time-dependent electrokinetic flow (mixed oscillating pressure gradient and electrical field) when using a linear Maxwell viscoelastic fluid in cylindrical microand nanochannels. The analytical solution is derived by solving the linearized Poisson-Boltzmann equation, together with the incompressible Cauchy's momentum equation in no-slip boundary conditions for the case of a combination of time periodic pressure-driven flow and electro-osmotic flow (PDF/EOF). The results show that the Onsager' reciprocal relations are complied with due to using the linear constitutive Maxwell fluid model. The validity of these Onsager's relations is important for practical implementation since it enables the analogy between fluidic networks in this flow configuration and electric circuits. We applied our calculation results for energy conversion systems in cylindrical micro-and nanochannels and compare the results with the work of Bandopadhyay and Chakraborty (2012a) which was performed in slit micro-nanochannels. It is shown that for both case the enhancement is in the same order. We furthermore provided a mechanism to understand the massive efficiency enhancement. We also found that time periodic electro-osmotic flow in many cases is much stronger enhanced than time periodic pressuredriven flow when comparing the flow profiles of oscillating PDF and EOF in micro-and nanochannels. The findings advance our understanding of time periodic electrokinetic phenomena of viscoelastic fluids and provide insight into flow characteristic as well as assist the design of devices for lab-on-chip applications.