1 Introduction

Microfluidic devices, e.g., lab-on-a-chip (LOC), intravenous drug delivery systems, and biochemical reactive platforms, are typically polymer-based devices possessing micro and nanoscale geometric features to hold and manipulate small volumes of biofluids^[1]. Such devices serve as tools for accurately controlling small volumes of liquids for a wide variety of chemical, medical, environmental, and biological applications^[2]. Microfluidic devices require the ability to pump, control, and manipulate samples. For these types of tasks, electroosmosis has been widely used to manipulate the fluid flows in the channels with the lengths of microns, and electroosmosis has been achieved through the electrostatic interaction between an applied external electric field and an electrical double layer (EDL)^[3]. An EDL is created when an electrolyte comes into contact with a dielectric material. This interaction generates an electric force near the wall, thereby driving the fluid motion subsequently transmitted to the bulk fluid through viscous forces^[4]. For transporting biofluids in the mentioned devices, there is need to mathematically characterize the transport mechanism for the efficient design of analytic systems. However, in several circumstances, the mathematical modeling in micro-conducts of LOC is conducted by assuming Newtonian fluids, which does not describe correctly the transport phenomena existing in such applications. This may be attributed to the fact that, the biofluids in their most general mathematical form are non-linear and strain-rate dependent, as governed by the material characteristics needed to analyze the response of a particular fluid against any disturbance being imposed^[5].

A considerable amount of research has been conducted on modeling the electroosmotic flows (EOFs) for Newtonian fluids under transient regimes. Peralta et al.^[6] studied the start-up of an oscillatory electroosmotic flow in a parallel-plate micro-channel, and showed that, with proper adjustment of the zeta potentials at the channel walls together with a dimensionless frequency, the velocity profile could be tuned to induce flow recirculation. This result should be useful in the design of microfluidic mixers. However, microfluidic devices are always used to analyze biofluids, whose rheological behaviors are very different from those of Newtonian fluids.

In the specialized literature, the transient EOFs of non-Newtonian fluids^[7-11] have received less attention compared with Newtonian fluids. Li et al.^[7] studied the transient EOFs of Maxwell fluids in a micro-parallel channel and in a micro-tube, and showed that, when the normalized relaxation time increased, it took longer for the flow to reach the steady state. Wang et al.^[8] investigated the transient EOFs of the generalized Maxwell fluids with a fractional derivative in a straight pipe with a circular cross section, and showed the effects of the relaxation time, fractional derivative parameter, and the Debye-Hückel parameter on the flow. Jian et al.^[9] and Liu et al.^[10] obtained the results regarding the transient EOFs driven by alternating current (AC) electric fields, and achieved an analytical solution of the time-periodic EOF for the generalized Maxwell fluids through a rectangular microchannel under the Debye-Hückel approximation. Bandopadhyay et al.^[11] characterized the electroosmotically driven flows for a linearized Maxwell fluid in presence of modulated surface charge, and studied the non-intuitive interactions between the patterned interfacial electrokinetics and the flow rheology.

The use of pulsatile flows has also been widely studied^[12-13]. However, this concept is scarce in the specialized literature on the EOFs. Chakraborty et al.^[14] and Chakraborty and Ray^[15] primarily conducted the related studies. They analyzed the EOFs driven by pulsating electric fields in microchannels to characterize and control the periodic mass flow rate. Rojas et al.^[16] recently studied a pulsatile EOF (PEOF) in a circular microchannel with high and low zeta potentials and slippage at the inner surface of the microchannel, where the fluid motion was caused by a pulsatile electric field. Some applications of PEOFs correspond to the active micromixers based on the disturbance induced by pulsatile external fields^[17]. In most of the cases, pressure-driven disturbance can be created by an external actuator. However, instabilities can also appear due to the variations in the electric fields^[18-19]. The instabilities in EOFs can also be generated due to the variations in the zeta potentials at the microchannel walls, which can be controlled by the shielding electrodes yielding the field effect^[20], by fabricating microchannels whose wall materials are different^[21], or by the variations in the pH of the solution along the microchannel^[22]. Therefore, the techniques used in the aforementioned papers to induce the changes at the walls, either in the surface or the zeta potentials of the walls, can modify the electrical potential distribution and thus the electrical body force in the flow field, allowing significant effects on the characteristics of the flow^[23-25].

In the present study, we will analyze the PEOF of a biofluid, e.g., blood, under certain rheological conditions, e.g., circulating in arteries, which has a behavior similar to that of a non-Newtonian fluid. The motion of the biofluid is driven by an external pulsatile electric field, and asymmetric zeta potentials are considered at the microchannel walls. To understand the fundamental physical aspects of this phenomenon, which can be used for enhancing the dispersion and the mass transport when the concentration gradients of miscible mass species are introduced into the flow field^[26-28], the basic mechanism is investigated.

2 Theoretical modeling 2.1 Physical model

Figure 1 shows the physical model analyzed in this work. We consider the PEOF of a symmetric (z:z) electrolyte solution, where z is the valence of the electrolyte in a parallel-plate microchannel with the height H and the length L (L» H). The flow is driven by a pulsatile electroosmotic force induced by the simultaneous effect of the EDL formed at the interface between the liquid and the microchannel surface and an external time-dependent electric field given by

where t, ω, and ε are the time, the angular frequency, and a dimensionless constant determining the amplitude of the electric field fluctuations, respectively. E₀ and εE₀ sin(ωt) represent the steady and oscillating components of the external electric field, respectively. The origin of the coordinate system is located at the lower surface of the microchannel. We consider that the walls of the microchannel have different (asymmetric) zeta potentials, i.e., ζ₁≠ ζ₂, and the high concentration of the electric charges is localized near the channel walls within the EDLs, whose screening lengths are represented by κ^-1, where κ ≡ (2 e²z²n_∞/(ϵk_BT))^1/2[29]. The net charge density in the EDLs follows the well-known Boltzmann distribution, which remains valid if the frequency of the external electric field is not very high (e.g., less than 1 MHz)^[30]. Moreover, it is assumed that the EDLs on the inner surfaces of the microchannel do not overlap, i.e., H » κ^-1, and the ends of the microchannel have the same pressure P₀.

2.2 Governing equations 2.2.1 Electric potential

According to Hsu et al.^[31], the characteristic time scale of the electro-migration in the EDL is on the order ranging from 10^-8 s to 10^-7 s, whereas the characteristic time scale associated with the evolution of the EOF is on the order ranging from 10^-5 s to 10^-3 s^[32]. Therefore, the latter time scale is at least two orders of the magnitude of the first time scale. Thus, the temporal evolution of the hydrodynamic field is considerably slower than that corresponding to the electrical phenomenon. Under this quasi-steady-state assumption, the distribution of the electric potential near the inner surface of the microchannel is governed by Poisson's equation^[4] as follows:

(1)

where ϵ is the dielectric permittivity of the liquid. Φ(x, y, t)=φ(x, t)+ψ(y) is the total electric potential in the microchannel, where φ(x, t) is the local electric potential, which depends parametrically on the time t, due to the electric field E_x. ψ(y) denotes the electric potential due to the EDLs. The charge density follows the following Boltzmann distribution:

(2)

where e, k_B, n_∞, T, and ψ are the electron charge, the Boltzmann constant, the ionic number concentration, the absolute temperature, and the electric potential due to the EDL in the equilibrium state, respectively. Because we are considering a very long microchannel, i.e., L » H, the term in Eq.(1) may be neglected^[4]. Then, we have

(3)

To solve Eq.(3), the following boundary conditions are required:

(4)

Due to the asymmetry of the wall zeta potentials, the plug-like velocity profile of pure EOFs will be no longer present, which will lead to a velocity gradient across the transversal section of the microchannel.

2.2.2 Flow field

To determine the dynamics of the PEOF, we have assumed that the microchannel is very long, and the analysis focuses on the central region, which is far away from the entry and the exit of the channel, such that the flow can be assumed to be unidirectional^[33]. In the absence of a pressure gradient, we use the modified Cauchy equation and the Maxwell constitutive equations^[34] given by

(5)

(6)

where u(y, t) represents the velocity in the x-direction, and ρ, µ, λ₁, and τ_xy are the mass density, the viscosity of the fluid, the relaxation time, and the shear stress, respectively. In the present study, the only component of the Maxwell model required to solve the problem is given by Eq.(5)^[35].

Equation (5) is subject to the following initial and boundary conditions:

(7)

To obtain the momentum equation exclusively in terms of the velocity u(y, t), Eq.(6) is derived once with respect to the variable y, and the result is substituted in Eq.(5). Then, we have

(8)

It is a hyperbolic partial differential equation, where the second derivative of u(y, t) with respect to t requires that Eq.(8) necessarily has a wave-like solution.

2.2.3 Dimensionless governing equations

We rescale the governing equations by introducing the dimensionless variables as follows:

where ψ_c=k_BT/(ze), u_HS=-ϵζ₁E₀/η₀ is the Helmholtz-Smoluchowski equation^[29], =τ_xyH/(u_HSη₀), and t_c=1 /ω. Because the present analysis is performed for low zeta potentials, i.e., zeψ/k_BT ≤ 25~mV, under the Debye-Hückel linearization, the dimensionless form of the Poisson-Boltzmann equation (3) becomes

(9)

where ≡ κH. The dimensionless forms of the boundary conditions associated with Eq.(9) are given by

(10)

where

Then, the dimensionless version of the momentum equation (8) can be rewritten as follows:

(11)

where R_ω=ωρH² /η₀ is the angular Reynolds number^[36]. It represents the ratio of the characteristic diffusive time to the characteristic time associated with the oscillatory electric field, and determines the importance of the acceleration effects in the fluid relative to the viscous effects (momentum diffusion). =λ₁η₀/(ρH²) is the dimensionless relaxation time representing the competition between the elastic and viscous effects of the fluid, and it is also called the elasticity number^[37].

The dimensionless initial and boundary conditions to solve Eq.(11) are as follows:

(12)

The equations given by Eqs.(9)-(12) describe the transient PEOF for low zeta potentials. However, in the present work, we focus only on the determination of the time-periodic PEOF, which means that the estimated solution corresponds to the stage after the transient component has died out.

3 Solution methodology

The solution for the electric field distribution of Eq.(9) subject to the boundary conditions defined in Eq.(12) is given by

(13)

where

(14)

Because Eq.(11) and the associated initial and boundary conditions are linear, we can write the corresponding solution as the sum of two terms, i.e., the response to the steady component of the electric field and the response to the time-dependent component of the oscillatory electric field, as follows:

(15)

where and represent the steady and oscillating solutions of the velocity field, respectively. Substituting Eq.(15) into Eq.(11), we have

(16)

with

(17)

and

(18)

with

(19)

3.1 Steady-state solution for

The solution of Eq.(16), together with the boundary conditions given by Eq.(17), is given by

(20)

where

3.2 Periodic solution for

The periodic solution for the flow field will be a repetitive oscillation, which can be found from the set of Eqs.(18) and (19) and ignoring the initial condition. To solve this problem, we define the complementary complex velocity in a manner such that

where Im denotes the imaginary part of the complementary complex velocity. Therefore, in terms of the complex velocity, Eq.(18) satisfies the following problem:

(21)

subject to the following boundary conditions:

(22)

To solve Eq.(21), a solution is assumed for as follows:

(23)

where denotes the imaginary number. Substituting Eq.(23) into Eq.(21) yields

(24)

subject to

(25)

Therefore, the solution , which is obtained from Eq.(24), is given by

(26)

where

By substituting into Eq. (23), we have

(27)

Therefore, from Eqs.(27), (20), and (15), the periodic solution for the dimensionless velocity profile of the PEOF with asymmetric zeta potentials at the wall can be given by

(28)

The procedure to obtain Eq.(28) is straightforward. However, to evaluate this equation for any value of R_ω, MATHEMATICA software^[38] is used. To obtain more physical insights about the nature of the exact solution, in the following paragraphs, we obtain an asymptotic solution of Eq.(24) for R_ω« 1 (low-frequency limit).

3.2.1 Asymptotic solution in the limit R_ω« 1

At this point, we consider it important to briefly discuss the magnitudes of the physical parameters involved in this low-frequency limit. Some typical values of the parameters used in the EOFs are as follows: 10µm≤ H ≤ 100µm, the strengths of the electric fields are of a few kilovolts/centimeters, and 0 s^-1≤ ω≤ 10⁴s^-1 or higher^[39], depending on the configuration and of the physical phenomenon^[11]. According to Refs.[40] and [41], ρ ~ 10³ kg·m^-3, and η₀~ 10^-3 Pa·s. Therefore, with a suitable combination of the above physical parameters, the values of R_ω« 1 can be estimated.

Let us consider Eq.(24) in the limit R_ω« 1 together with the symmetric case of zeta potentials, i.e., . First, we propose an asymptotic expansion as follows:

(29)

Here, we have assumed up to the terms of O(R_ω²) in the expansion (29) because in the dimensionless momentum governing equation, Eq.(21), the asymptotic solution of this equation, will appear when the terms of O(R_ω²) are retained.

Substituting Eq.(29) into Eq.(24), we can obtain the following set of equations.

(ⅰ) The leading-order problem is given by

(30)

with the boundary conditions

(31)

The solution for this order is given by

(32)

(ⅱ) The O(R_ω¹) problem is

(33)

with the boundary conditions

(34)

The solution of F₁ is

(35)

where Π₁ and Π₂ are defined as follows:

(36)

(ⅲ) At the order O(R_ω²), the problem is given by

(37)

with the boundary conditions

(38)

The solution of order R_ω² is

(39)

where

(40)

4 Results and discussion

In this section, we present and discuss the results obtained for the PEOF of Maxwell fluids with asymmetric low zeta potentials at the microchannel walls. To describe the results, all variables and parameters are presented in dimensionless form. The hydrodynamic behavior of the PEOF is described in terms of the dimensionless parameters involved in the study. To estimate the values of the dimensionless parameters, we use the common values of the physical parameters reported in Ref.[42], some of which have already been presented in Subsection 3.2.1. For the numerical calculations, we select a suitable combination of values for the following physical parameters: 1 nm≤ κ^-1≤ 300 nm^[43], |ζ₁|≤25 mV, and |ζ₂|≤25 mV. The relaxation time can assume to be in the range of 10^-4 s≤ λ₁ ≤10³ s^[44]. However, according to Liu et al.^[10], in order to ensure the validity of the fundamental assumption of undisturbed EDL, the relaxation time λ₁ should be smaller than the oscillating period of the electric field 2π/ω. Thus, the product of the external electric field frequency and the relaxation time should be smaller than 2π, i.e., λ₁ω < 2π. In terms of the dimensionless parameters used in this analysis, the above restriction can be written as follows:

In addition, in the typical applications of AC-driven electroosmotic flows, the range of the angular frequency is very broad. For instance, in Ref.[45], 0≤ R_ω≤ 100, 0s^-1≤ ω≤ 10⁴ s^-1, ρ = 10³ kg·m^-3, the radius of the microchannel R = 100 µm, and η₀= 10^-3 Pa·s. In Ref.[42], 10³s^-1≤ ω≤ 10⁵ s^-1. In the following calculations, 0≤ R_ω≤ 7.

(ⅰ) Asymptotic solution

From examining Eq.(24), the qualitative nature of the relationship between and the external electric field E_x(t) can be observed by considering the limiting case of R_ω«1. It is evident that in such an asymptotic limit, the inertial and elastic terms in this equation can be neglected. Thus, the solution for will have the form This result is verified by the asymptotic solution conducted in Subsection 3.2.1. As can be appreciated, in Eq.(30), the velocity profiles are in-phase with the signals of the external electric field for all time, while the elastic effect is not appreciated in the leading and first-order solutions, i.e., O(1) and O(R_ω¹), respectively (see Eqs.(30) and (33)). We can observe that in this asymptotic limit, the elastic effects appear up to the terms of O(R_ω²), which are very small. In Fig. 2, we present a comparison between the asymptotic and the exact solutions for the dimensionless velocity profiles given by Eqs.(29) and (28), respectively, and find a very good agreement. We note that the asymptotic solution is valid up to values of R_ω=0.3 (see Fig. 2). For R_ω«1, a plug-like velocity Helmholtz-Smoluchowski EOF velocity profile is expected when symmetric zeta potentials are assumed. All the velocity variations are restricted to very narrow EDL regions close to the walls. The above can be explained from the definition of R_ω if we consider the values of the parameter of O(1), which means that the diffusion time ρH² /η₀ and the imposed time scale of the periodic electric field 1/ω are of the same order.

(ⅱ) Periodic stage of the Newtonian fluid and asymmetric zeta potentials

For , the PEOF of a Newtonian fluid with asymmetric wall zeta potentials is recovered (see Fig. 3). Figure 3 shows the effects of the parameter R_ζ on the dimensionless velocity profiles, where the periodic solution of the dimensionless velocity profiles given by Eq.(28) is plotted as a function of the dimensionless transversal coordinate , ε=0.5, R_ω=5, and R_ζ=(1.0, 0.5, -1.0, -0.5). The results are obtained by use of the MATHEMATICA software^[38]. From the figure, we can see that, for R_ω=5, the velocity profiles are not uniform through the transversal coordinate due to the effect of the parameter R_ω (see Fig. 3(a)), the magnitude of the velocity near the upper wall is smaller than that at the lower wall due to the fact that the electroosmotic force is larger at (see Fig. 3(b)), which is caused by the corresponding zeta potential. The results agree well with those reported by Afonso et al.^[24]. Moreover, when the ratio of the zeta potentials R_ζ changes from symmetric (R_ζ=1.0) to anti-symmetric (R_ζ=-1.0), the corresponding dimensionless velocity profiles vary from fully symmetric to fully anti-symmetric. Figures 3(c) and 3(d) also show that any negative value of R_ζ will result in an inverse flow.

(ⅲ) Effects of the relaxation time on the PEOF with asymmetric zeta potentials

In Fig. 4, we plot the velocity profiles for the PEOF of a Maxwell fluid ( >0.00) as a function of the dimensionless transversal coordinate , where =50, ε=0.5, R_ω=5, R_ζ=(1.0, 0.5, -1.0, -0.5), and = 0.50. The fluid motion begins near the walls of the microchannel, where the electroosmotic force is acting, and propagates into the central region of the channel through momentum diffusion. This physical behaviors agree well with those reported by Escandón et al.^[35]. Besides, the behaviors of the asymmetric zeta potentials on the velocity profiles are similar to those presented in Fig. 3.

In Fig. 5, the velocity field and the associated memory effects of the PEOF of a Maxwell fluid are plotted, where =50, R_ω=5, ε=0.5, and R_ζ=1.0 with low and symmetric zeta potentials are assumed. All the velocity profiles are evaluated at =π /2. For a given value of R_ω, when the dimensionless relaxation time increases, the velocity profiles will oscillate rapidly. At the same time, the amplitudes of the PEOF velocity decrease gradually. This is because that the diffusion time scale is much greater than the oscillation time period. Therefore, there is no sufficient time for the momentum to diffuse far into the central region of the microchannel, and the PEOF velocity variations are restricted only to a thin layer near the two microchannel walls. Note that for the Newtonian fluid case (=0.00), a uniform profile is present across the microchannel, whereas for ( > 0.00, the elastic effects will affect the velocity profiles, causing wave motions, which are due to the fading and elastic memory phenomena of the viscoelastic fluids^[10]. Conversely, increasing the values of the relaxation time more easily leads to the variation in the velocity profiles caused by the external pulsatile electric field because of the shear thinning effect of the Maxwell fluids. Shear thinning is an effect, in which the viscosity decreases with the increase in the rate of the shear stress. Thus, the amplitude of the waves in the velocity profiles becomes smaller with the gradual increase in the relaxation time. However, with the increase in the relaxation time, the elasticity of the fluid becomes notable, and the velocity variations can be expanded to the whole flow field since the elasticity is a physical property of the whole fluid, from which we can easily find that higher angular Reynolds number R_ω results in quicker oscillation of the velocity profiles.

Figure 6 shows the effects of the relative amplitude of the sinusoidal electric signal on the velocity profiles. Here, we have considered ε=(0.0, 0.2, 0.5, 1.0). All curves are plotted at the selected dimensionless time =π/2 for =50, R_ω=5, and R_ζ=1.0. The solid line corresponds to the case of a direct current (DC)-driven EOF of a Maxwell fluid, and shows a plug-like behavior. The above can be inferred from Eq.(28), where the contribution of the oscillating electric field is null by making ε=0.0. In such a case, in the periodic stage, the velocity profiles do not depend on the rheological characteristics of the Maxwell fluid. In contrast, when ε increases, the amplitude of the waves in the velocity profiles increases.

(ⅳ) Effects of R_ω and on the flow field

The effects of R_ω and on the velocity profiles of a Maxwell fluid are shown in Figs. 7(a), 7(b), 7(c), and 7(d). By comparing these figures against the case of a pure EOF (=0.00) and considering the values of both R_ω > 1 and » 1, we can see that the velocity profiles show a wave motion in contrast to a pure EOF, where plug-like velocity profiles are found. Evidently, this behavior could be used as a benchmark for designing microfluidic devices, where micromixing is needed, at higher frequencies, viscoelastic fluids are expected to give better mixing and dispersion as compared with a Newtonian fluid. Moreover, in this direction, the flow is in an oscillating stage, the fluid motion begins near the surface of the microchannel where the electroosmotic force is acting, and the effect of the electroosmotic force propagates into the central region of the microchannel due to momentum diffusion. However, elastic forces always exist, which will affect the aforementioned wave motion all time as long as R_ω>1.

Figure 8 compares the periodic evolution of the dimensionless velocity profiles as functions of the dimensionless coordinate for fixed values of R_ζ=-1.0, ε=0.5, and = 0.50 for two cases of R_ω and . In Fig. 8(a), R_ω=1.0 and =10 are used, and the second case in Fig. 8(d) uses R_ω=7 and =100. Here, the anti-symmetric behavior of the velocity can be observed. Similar to Fig. 7, oscillations in the velocity profiles exist in Figs. 7(c) and 7(d) for large values of .

(ⅴ) Volumetric flow rate

The instantaneous dimensionless volumetric flow rate Q and the time-averaged volumetric flow rate 〈Q 〉 can be evaluated as follows:

(41)

(42)

where

By substituting the velocity profile derived in the analysis and after integrating the above equations, Q is obtained and plotted in Fig. 9. The oscillatory behaviors of the volumetric flow rate as a function of the dimensionless time for low zeta potentials and for three different values of the dimensionless relaxation time =(0.00, 0.50, 1.25) are depicted, where =50, ε=0.5, and R_ζ=1.0, and R_ω(=5, 7). The cycles of Q depend on the frequency of the applied electric field, and as expected, when increases, there is a decrease in the oscillation between the frequency and the volumetric flow rate.

5 Conclusions

In this work, an analytical solution for the PEOF of Maxwell fluids in a parallel-plate microchannel under asymmetric zeta potentials is analyzed. Through the solution, the following conclusions can be drawn:

(ⅰ) In the asymptotic limit of R_ω« 1, the velocity profiles are in-phase with the signal of the external electric field for all time, and the elastic effect is not appreciable.

(ⅱ) The mathematical model allows the behavior of the PEOF with symmetric and asymmetric zeta potentials to be determined.

(ⅲ) The negative values of R_ζ produce an inverse flow, causing asymmetric velocity profiles.

(ⅳ) Compared with the PEOFs of Newtonian fluids, the wave motions in the velocity profiles of Maxwell fluids are observed, even if » 1.

The analytical model developed in this work is of importance because it can describe the flow behavior of the PEOF of biofluids and is useful as a benchmark for designing microfluidic devices. The solution derived in this work may be extended to analyzing the mass transport in EOFs^[40] and studing the dispersion of solutes in microchannels. In addition, the results could be used for studying the instability in channels where biofluids are being transported. Such PEOFs may also be useful for improving the performance of micromixing due to the transversal transport of solutes and Taylor dispersion.

Acknowledgement M. PERALTA acknowledges the support of the Consejo Nacional de Ciencia y Tecnología program for a postdoctoral fellowship at the Escuela Superior de Ingeniería Mecánica Azcpotzalco from the Instituto Politécnico Nacional of Mexico.