Analytical Solution of Time-Periodic Electroosmotic Flow through Cylindrical Microchannel with Non-Uniform Surface Potential

Abstract

Time-periodic electroosmotic flow (EOF) with heterogeneous surface charges on channel walls can potentially be used to mix species or reagent molecules in microfluidic devices. Although significant research efforts have been placed to understand different aspects of EOF, its role in the mixing process is still poorly understood, especially for non-homogeneous surface charge cases. In this work, dynamic aspects of EOF in a cylindrical capillary are analyzed for heterogeneous surface charges. Closed form analytical solutions for time-periodic EOF are obtained by solving the Navier–Stokes equation. An analytical expression of induced pressure is also obtained from the velocity field solution. The results show that several vortices can be formed inside the microchannel with sinusoidal surface charge distribution. These vortices change their pattern and direction as the electric field change its strength and direction with time. In addition, the structure and strength of the vorticity depend on the frequency of the external electric field and the size of the channel. As the electric field frequency or channel diameter increases, vortices are shifted towards the channel surface and the perturbed flow region becomes smaller, which is not desired for effective mixing. Moreover, the number of vorticities depends on the periodicity of the surface charge.

1. Introduction

Lab-on-a-chip microfluidic devices are utilized for numerous applications such as DNA sequencing, synthesis, crystallization, polymerization, and drug discovery. These applications require a number of unit operations such as pumping [1], mixing [2], etc. Among various functionalities, mixing in microfluidic devices is very challenging due to (ultra-low Reynolds number) creeping flow [3]. The quick mixing of species is key to reducing analysis time in a microfluidic device for its widespread application in the biomedical field [4]. To date, several mixing strategies have been proposed for microfluidic systems with different degrees of successes (see reviews [5,6,7]). Based on the methods used to achieve mixing, micromixers are generally classified as being passive or active. In passive micromixers, the mixing process depends entirely on diffusion; meanwhile, in active micromixers, the mixing process is accelerated by an external disturbance [6,8]. Although the passive mixing offers some benefits such as no external power requirements, it is not very effective in microfluidic devices because of the diffusion-based very slow mixing process. Active mixing has the potential to eliminate the inherent drawbacks of slow mixing using an acoustic, magnetic, or electrokinetic forcing [9]. Among various active micromixers (see reviews [2,5]), electrokinetic-based micromixers are preferred owing to their numerous advantages including ease of fabrication, no moving parts, ease of control, high reliability and repeatability, and quiet operation.

Electroosmotic mixing can be enhanced with an alternating electric field. Alternating electric field-based time-periodic electroosmotic flow (EOF) has been extensively studied with uniform surface charges for various types of channels such as parallel plate [10], circular [11], rectangular [12], and annulus [13] channels, etc. Oddy et al. [14] found that the rapid stretching and folding of fluid interfaces, induced by an alternating electric field, are able to stir fluid streams at a very low Reynolds number (Re < 1). Glasgow et al. [15] found that better mixing can be achieved in a T-junction channel with two out-of-phase alternating electric fields. The out-of-phase alternating electric field induces the oscillation of the fluid interface at the junction of the two inlet channels, which increases the contact surface area between the two fluid streams for significantly better mixing efficiency [16]. However, the alternating electric field-based EOF fails to take advantage of the vortex to enhance the mixing due to the unidirectional flow at any particular time.

The use of a nonhomogeneous channel surface is another way to improve electroosmotic mixing performance. The effects of a heterogeneous surface charge have also been extensively studied in direct current (DC) EOF. For example, Ajdari [17] theoretically showed that a surface with a sinusoidal charge can create vorticities within the bulk flow. Vortices formation within the bulk flow was also shown by Horiuchi et al. [18] for a step change in the surface charge in a DC EOF. Through a numerical investigation, Erickson and Li [19] showed that surface heterogeneity increases the mixing efficiency and reduces the required mixing length. Qian and Bau [20] presented a theoretical study of electroosmotic flows driven by a uniform electric field in a two-dimensional conduit with non-uniform surface potential distribution.

In microfluidic devices, surface heterogeneity can be introduced intentionally through micromanufacturing technology, such as microcontact printing or rapid prototyping. For instance, Biddiss et al. [4] manufactured several micromixers with heterogenous surface charge configurations by rapid prototyping for the quick mixing of species. Their analysis showed that heterogeneously charged micromixers significantly improve the mixing efficiency as compared to homogeneously charged micromixers. Strook et al. [21] also used microcontact printing to vary the surface charge in directions parallel and perpendicular to the applied electric field. They found that surface charge variation in the parallel direction generates recirculating flow, while surface charge variation in the perpendicular direction creates multidirectional flow. Electrokinetic characteristics of a microchannel surface can also be modulated through the dynamic or static coating of proteins, DNA, colloids, and/or nanosized particles on the surface. For example, Norde and Rouwendal [22] studied the change of surface charge by the protein adsorption on a glass surface. Wei et al. [23] generated non-uniform zeta potential distribution by periodically attaching DNA molecules in the microfluidic channel surface for increased DNA–DNA hybridization. The aforementioned studies showed that DC EOF with a heterogeneous surface charge can increase the mixing efficiency while reducing the mixing length by creating vorticity. However, a DC EOF fails to take advantages offered by an alternating electric field.

Thus, the primary focus of this work is to further increase the efficiency of electroosmotic mixing by introducing an alternating electric field in a heterogeneously charged channel. Through a numerical investigation, Lee et al. [24] demonstrated various kinds of flow circulation in a circular slit with time-periodic EOF for non-uniform zeta potential distribution along the microchannel walls. Luo [25] numerically investigated the transient electroosmotic flow induced by an alternating electric field with patch-wise surface heterogeneities in a planner microchannel. Tang et al. [26] studied the pulsating electroosmotic flow in a planar microchannel with non-uniform surface charges by the lattice Boltzmann method. These simulation results showed that time-periodic EOF with a heterogeneous surface can combine the advantages of the alternating electric field and the heterogeneous surface charge. However, no generalized analytical model exists for time-periodic EOF with non-uniform zeta potential in cylindrical systems. An analytical solution could help to find the optimum micromixer design. Motivated by the growing interest in time-periodic EOF as a reliable non-mechanical strategy to mix species in microfluidic devices, in this paper, a theoretical approach is presented. Specifically, we studied a time-periodic EOF field in a cylindrical microchannel with heterogeneous surface charge distribution. A closed-form analytical solution was obtained for the time-periodic electroosmotic flow velocity by solving the Navier–Stokes equations. In addition, an expression for pressure distribution was also obtained for time-periodic electroosmotic flow.

2. Mathematical Model

2.1. Governing Equations for Time-Periodic EOF in a Cylindrical Microchannel

Figure 1 shows a schematic for time-periodic EOF in a cylindrical microchannel, where an alternating electric field is applied along the channel length. For generality, the model was developed such that the surface zeta potential can take any periodic distribution (e.g., sinusoidal, square) along the length of the channel, but there is no variation in zeta potential along the circumferential direction.

The application of a time-periodic external electric field results in a net body force on the free ions within the electric double layer (EDL), inducing a time-periodic bulk fluid motion. The motion of the fluid inside the channel can be governed by the incompressible Navier–Stokes equation [10]:

$$\rho \left[\frac{\partial \overrightarrow{V}}{\partial t}+\left(\overrightarrow{V}\cdot \nabla \right)\overrightarrow{V}\right]=-\nabla P+\mu {\nabla }^2\overrightarrow{V}+{\rho }_e\overrightarrow{E}$$

(1)

where $\rho$ and $\mu$ are the density and viscosity of the fluid, respectively; $P$ is the pressure; ${\rho }_e$ is the charge density; and $\overrightarrow{E}$ is the applied electric field. The applied electric field can be expressed as $\overrightarrow{E}=E_z\widehat{z}=E_0\sin \left(\omega t\right)\widehat{z}$, where $\widehat{z}$ is the unit vector in the axial (z) direction and $E_0$ is the reference electric field. In addition, the mass conservation equation can be expressed as:

$$\frac{\partial \rho }{\partial t}+\nabla \cdot \left(\rho \overrightarrow{V}\right)=0.$$

(2)

In general, for an incompressible fluid, density is assumed to be a constant. Thus, for the axisymmetric flow scenario ($v_{\theta }=0$) presented here, the mass conservation equation can be reduced to:

$$\frac{1}{r}\frac{\partial }{\partial r}\left(r{v}_r\right)+\frac{\partial v_z}{\partial z}=0.$$

(3)

where $v_r$ and $v_z$ are the velocity components in the radial and axial direction, respectively. For the axisymmetric case, the radial and axial velocity can be related to the stream function, $\psi$ as follows:

$$v_r=\frac{1}{r}\frac{\partial \psi }{\partial z},v_z=-\frac{1}{r}\frac{\partial \psi }{\partial r}.$$

(4)

2.2. Analysis of Time-Periodic EOF in a Cylindrical Microchannel

Our next objective is to simplify the Navier–Stokes equations. In most microfluidic applications, the buffer solutions have a concentration of the order of $mM$, which results in a very thin EDL. The bulk fluid flow outside of the EDL region can be modeled by dropping the electroosmotic body force, ${\rho }_e\overrightarrow{E}$, and by introducing a Helmholtz–Smoluchowski slip boundary condition at the channel wall [24]:

$${\overrightarrow{V}}_{r=r_0}=-\frac{\epsilon \zeta \left(z\right)}{\mu }\overrightarrow{E}$$

(5)

where $\epsilon$ is the permittivity and $\zeta \left(z\right)$ is the surface potential distribution along the length of the conduit. This slip boundary condition was originally developed for steady-state EOF, but it can still be used for time-periodic EOF because the applied electric field frequency (${10}^2\sim {10}^5\mathrm{Hz}$) is less than the charge relaxation frequency (${10}^6\sim {10}^{8}Hz$) [27]. Moreover, the Helmholtz–Smoluchowski formulation-based slip boundary condition is widely used for nonhomogeneous surface charges as well as time-periodic EOF analysis [3,24,28,29,30]. Thus by dropping the body force, the component form of equations of motion are:

$$\rho \left(\frac{\partial v_r}{\partial t}+v_r\frac{\partial v_r}{\partial r}+v_z\frac{\partial v_r}{\partial z}\right)=-\frac{\partial P}{\partial r}+\mu \left[\frac{\partial }{\partial r}\left\{\frac{1}{r}\frac{\partial }{\partial r}\left(r{v}_r\right)\right\}+\frac{{\partial }^2{v}_r}{\partial z^2}\right],$$

(6a)

$$\rho \left(\frac{\partial v_z}{\partial t}+v_r\frac{\partial v_z}{\partial r}+v_z\frac{\partial v_z}{\partial z}\right)=-\frac{\partial P}{\partial z}+\mu \left[\frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial v_z}{\partial r}\right)+\frac{{\partial }^2{v}_z}{\partial z^2}\right].$$

(6b)

Here, the $\theta$-component of the momentum equation is discarded due to axisymmetric flow, as there is no driving force in the $\theta$-direction. In general, the advection effects are negligible in EOF due to the low Reynolds number [31,32]. For example, if water ($\rho =1000$ $\mathrm{kg}·{\mathrm{m}}^{-3}$ and $\mu =0.001$ $\mathrm{kg}·{\mathrm{m}}^{-1}·{\mathrm{s}}^{-1}$) flows through a $100$ $\mathsf{\mu }\mathrm{m}$ diameter cylindrical microchannel with an axial velocity of $1$ $\mathrm{mm}·{\mathrm{s}}^{-1},$ the resulting Reynolds number would be $Re=0.1$. Thus, by dropping the advective term, eliminating $P$ between Equation (6), and introducing the stream function $\psi$, we get:

$$\left(D-\frac{1}{\vartheta }\frac{\partial }{\partial t}\right)D\psi =0$$

(7)

where $\vartheta$ is the kinetic viscosity and $D$ denotes the following operator:

$$D\equiv \frac{{\partial }^2}{\partial r^2}-\frac{1}{r}\frac{\partial }{\partial r}+\frac{{\partial }^2}{\partial z^2}.$$

(8)

The commutative properties of operators $D$ and $D-\frac{1}{\vartheta }\frac{\partial }{\partial t}$ lead us to separate the stream function, $\psi$, into two parts [32,33] as $\psi ={\psi }_1+{\psi }_2$, where ${\psi }_1$ satisfies the equation

$$D{\psi }_1=0$$

(9)

and ${\psi }_2$ satisfies

$$\left(D-\frac{1}{\vartheta }\frac{\partial }{\partial t}\right){\psi }_2=0.$$

(10)

Considering periodic boundary conditions in the upstream ($z=0$) and downstream ($z=L$) regions, we assume wave-like solutions of elliptic Equations (9) and (10), which can be given as:

$${\psi }_1={\displaystyle \sum _{n=-\infty }^{\infty }{\phi }_{1,n}\left(r\right)e^{i{k}_{n}z}{e}^{i\omega t}},$$

(11a)

$${\psi }_2={\displaystyle \sum _{n=-\infty }^{\infty }{\phi }_{2,n}\left(r\right)e^{i{k}_{n}z}{e}^{i\omega t}}.$$

(11b)

Here, ${\varphi }_{1,n}$ and ${\varphi }_{2,n}$ are eigenfunctions, $\omega$ is the temporal angular frequency, and $k_n$ is the eigenvalue for the spatially periodic process. Thus, Equation (11) provides a general solution for axially periodic boundary conditions [32,33]. The values of $k_n$ can be given as $2\pi n/L$, where $n$ is an integer and $L$ is the length of the microfluidic conduit. The eigenfunctions ${\varphi }_{1,n}$ and ${\varphi }_{2,n}$ are solutions of the following differential equations:

$$\frac{d^2{\phi }_{1,n}}{d{r}^2}-\frac{1}{r}\frac{d{\phi }_{1,n}}{dr}-k_n^2{\phi }_{1,n}=0,$$

(12a)

$$\frac{d^2{\varphi }_{2,n}}{d{r}^2}-\frac{1}{r}\frac{d{\varphi }_{2,n}}{dr}-l_n^2{\varphi }_{2,n}=0,$$

(12b)

satisfying the boundary conditions at the tube wall. In Equation (12b), $l_n^2=k_n^2+i\omega /\vartheta$. The general solutions of Equation (12a,b) can be obtained as

$${\varphi }_{1,n}=r\left[A_{1,n}{I}_1\left(k_{n}r\right)+B_{1,n}{K}_1\left(k_{n}r\right)\right],$$

(13a)

$${\varphi }_{2,n}=r\left[A_{2,n}{I}_1\left(l_{n}r\right)+B_{2,n}{K}_1\left(l_{n}r\right)\right],$$

(13b)

respectively, where $I_1\left(r\right)$ and $K_1\left(r\right)$ are modified Bessel functions of the first and second kind of order $1$; $A_{1,n}$, $A_{2,n}$, $B_{1,n}$, $B_{2,n}$ are coefficients. Thus, the stream function, $\psi$, can be given as follows:

$$\psi ={\displaystyle \sum _{n=-\infty }^{\infty }r\left[A_{1,n}{I}_1\left(k_{n}r\right)+B_{1,n}{K}_1\left(k_{n}r\right)+A_{2,n}{I}_1\left(l_{n}r\right)+B_{2,n}{K}_1\left(l_{n}r\right)\right]e^{i{k}_{n}z}{e}^{i\omega t}}.$$

(14)

From the stream function solution (Equation (14)), the radial and axial velocity components can be obtained using Equation (4):

$$v_r={\displaystyle \sum _{n=-\infty }^{\infty }\left[A_{1,n}{I}_1\left(k_{n}r\right)+B_{1,n}{K}_1\left(k_{n}r\right)+A_{2,n}{I}_1\left(l_{n}r\right)+B_{2,n}{K}_1\left(l_{n}r\right)\right]}i{k}_n{e}^{i{k}_{n}z}{e}^{i\omega t},$$

(15a)

$$v_z={\displaystyle \sum _{n=-\infty }^{\infty }\left[-k_n{A}_{1,n}{I}_0\left(k_{n}r\right)+k_n{B}_{1,n}{K}_0\left(k_{n}r\right)-l_n{A}_{2,n}{I}_0\left(l_{n}r\right)+l_n{B}_{2,n}{K}_0\left(l_{n}r\right)\right]}{e}^{i{k}_{n}z}{e}^{i\omega t}.$$

(15b)

Our next objective is to find coefficients for velocity distribution. The velocity components $v_r$ and $v_z$ must be finite everywhere. This condition requires that the velocity components do not contain the function $K_1\left(k_{n}r\right)$; hence, $B_{1,n}=B_{2,n}=0$. At the channel wall ($r=r_0=d/2$), the boundary conditions for radial and axial velocity become no penetration to the wall ($v_r=0$) and Helmholtz–Smoluchowski slip velocity ($v_z=-\epsilon \zeta \left(z\right)E_0{e}^{i\omega t}/\mu$), respectively. The first boundary condition requires that $A_{2,n}=-A_{1,n}{I}_1\left(k_n{r}_0\right)/I_1\left(l_n{r}_0\right)$. Using the second Helmholtz–Smoluchowski boundary condition and applying complex Fourier analysis, one can find the remaining coefficient using the equation $A_{1,n}=-\frac{\epsilon E_0{I}_1\left(l_n{r}_0\right)}{\mu \left[-k_n{I}_0\left(k_n{r}_0\right)I_1\left(l_n{r}_0\right)+l_n{I}_1\left(k_n{r}_0\right)I_0\left(l_n{r}_0\right)\right]}\frac{1}{L}{\displaystyle \underset{0}{\overset{L}{\int }}\zeta \left(z\right)e^{-i{k}_{n}z}dz}$. Thus, the final forms of radial and axial velocity can be given as

$$v_r={\displaystyle \sum _{n=-\infty }^{\infty }{A}_{1,n}\left[I_1\left(k_{n}r\right)-\frac{I_1\left(k_n{r}_0\right)}{I_1\left(l_n{r}_0\right)}{I}_1\left(l_{n}r\right)\right]}i{k}_n{e}^{i{k}_{n}z}{e}^{i\omega t},$$

(16a)

$$v_z={\displaystyle \sum _{n=-\infty }^{\infty }{A}_{1,n}\left[-k_n{I}_0\left(k_{n}r\right)+l_n\frac{I_1\left(k_n{r}_0\right)}{I_1\left(l_n{r}_0\right)}{I}_0\left(l_{n}r\right)\right]}{e}^{i{k}_{n}z}{e}^{i\omega t},$$

(16b)

respectively. The imaginary part of these equations provides flow velocity corresponds to the external electric field, $\overrightarrow{E}=E_{0}sin\left(\omega t\right)=Im\left(E_0{e}^{i\omega t}\right)$, which we considered throughout the paper. Even though Equation (16) provides the analytical solution for time-periodic electroosmotic flow in a heterogeneously charged circular shaped channel, it cannot be used for time-periodic electroosmotic flow in a homogeneously charged channel since we started our analysis for a general solution (Equation (11)) with periodicity in the axial (z) direction. Thus, to obtain time-periodic electroosmotic velocity in a homogenously charged microchannel, one has to drop the z-dependency in the general solution and governing equations (Equation (6)) and follow the method described in this work.

3. Results and Discussion

In this section, results are provided for a sinusoidal surface charge distribution, although the general solution presented in the aforementioned section can be applied for any periodic surface charge distribution along the tube. The sinusoidal surface charge distribution can be given as:

$$\zeta \left(z\right)={\zeta }_0\sin \left(qz\right)$$

(17)

where ${\zeta }_0$ is the amplitude (reference) of the surface charge and $q$ is the angular frequency of the surface charge distribution. With this charge distribution, the coefficient for Equation (16) can be found as follows:

$$A_{1,n}=-\frac{\epsilon E_0{I}_1\left(l_n{r}_0\right)}{\mu \left[-k_n{I}_0\left(k_n{r}_0\right)I_1\left(l_n{r}_0\right)+l_n{I}_1\left(k_n{r}_0\right)I_0\left(l_n{r}_0\right)\right]}\frac{{\zeta }_0}{L}\left[\frac{\left\{i{k}_n\sin \left(qL\right)+q\cos \left(qL\right)\right\}{e}^{-i{k}_{n}L}-q}{k_n^2-q^2}\right]$$

(18)

Throughout this section, unless otherwise stated, the value of $q$ is considered as $2\pi /L$, which corresponds to a single period in the surface charge distribution as shown in Figure 2.

3.1. Velocity Profiles at Various Nondimensional Times

Figure 3 illustrates the velocity vector at different nondimensional times for a sinusoidal zeta potential distribution along the channel surface with a reference potential of −100 mV. When the zeta potential is changed along the channel wall, the driving force in electroosmotic flow is also changed. Consequently, a non-uniform flow field is developed along the channel (Figure 3). At $\omega t=\pi /2$, the electric field is positive. Thus, the fluids near the walls move in the positive x-direction in the left half of the channel, since the zeta potential is negative in the left half. In contrast, due to the positive surface charge in the right half of the channel, the fluids near the walls move in the negative x-direction. These two opposite directional flows create two counter-rotating vortices in the upper part ($0\le r/r_0\le 1$) of the channel (Figure 3a). In fact, two counter-rotating vortices are formed at any angle $\theta$ because of the axisymmetric conditions. In other words, for a single period of surface potential, two oppositely rotating vortices (each one extending one half of the total length) are formed within the tube. So, if the tube is filled with two fluids in such a way that each fluid fills the entire cross-section and one-fourth length of the channel in an alternative patch, this design will have strong potential to enhance the mixing efficiency. Figure 3b shows the vector plot of the flow field at a nondimensional time, $\omega t=3\pi /4$. This yields similar results to the previous case (Figure 3a) due to the similar conditions; however, the magnitude of the velocity components at the surface is reduced because of the reduced strength of the electric field. At $\omega t=\pi$, the slip velocity at the wall reaches zero due to the no applied electric field (Figure 3c). However, bulk fluid motion is still observed at the center part of the tube due to the phase lag between fluids in the electric double layer and bulk fluids. The phase lag occurs due to the finite time requirement for momentum diffusion from the surface to the center line of the tube. This out-of-phase behavior is clearly visible at time $\omega t=5\pi /4$ (Figure 3d), where velocity is negligible at the centerline of the channel but slip velocity at the channel surface is moderate. The centerline velocity approaches zero at $\omega t=1.26\pi \mathrm{and} 2.26\pi .$ Based on these results, one can calculate the time lag as ~0.26$\pi$ or 81.25 $\mathsf{\mu }\mathrm{s}$ for an applied electric field frequency of $1.6$ $\mathrm{kHz}.$

When the electric field changes its direction, the velocity vector also switches its course (Figure 3d,e). For example, the distribution of velocity vectors at nondimensional time $\omega t=3\pi /2$ appears to be opposite that at $\omega t=\pi /2$ (Figure 3e vs. Figure 3a). It has been found that in EOF, the net flow and its direction depend on the overall strength of the zeta potential:

$$\overline{\zeta }=\frac{{\displaystyle \sum _i{\zeta }_i\Delta L_i}}{{\displaystyle \sum _i\Delta L_i}}$$

(19)

where ${\zeta }_i$ and $\Delta L_i$ are the zeta potential and channel length at the $i$-th section, respectively [24]. Since the considered sinusoidal zeta potential has an overall strength of zero, the net flow in the current scenario is zero at any time.

3.2. Effect of Nondimensional Frequency

Previous works have shown that both electric field frequency and channel height have a similar effect on electroosmotic vortices [3]. Thus, it is wise to combine these two parameters to obtain a nondimensional frequency as $\Omega =\omega r_0^2/\vartheta$. This nondimensional frequency represents the ratio of the diffusion time scale, $t_{diff}=r_0^2/\vartheta$, to the period of the external electric field, $t_p=1/\omega$ [34]. The nondimensional frequency increases if either applied electric field frequency increases or the size of the tube increases.

The effect of nondimensional frequencies is shown in Figure 4 for the sinusoidal zeta potential distribution case presented in Figure 3. Velocity vectors are presented for four different nondimensional frequencies: (a) $\Omega =2.5,$ (b) $\Omega =25,$ (c) $\Omega =125$, and (d) $\Omega =250$ at $\omega t=\pi /2$. For a $100\mathsf{\mu }\mathrm{m}$ diameter tube, these normalized frequencies correspond to $f=0.16,$ $1.6,$ $8\mathrm{and}$ $16\mathrm{kHz}$. As the nondimensional frequency increases, the pattern of the flow field changes significantly (Figure 4). At a low value of $\Omega$ (e.g., 2.5), the diffusion time scale is on the same order of magnitude of the external electric field period. As a result, the flow has enough time to propagate from the surface to the center of the channel, resulting in large vortices which extend from the surface to the center region of the tube (Figure 4a). However, as $\Omega$ increases from 2.5 to 25, these vortices are slightly shifted toward the surface and the bulk fluid velocity diminishes slightly (Figure 4b). When $\Omega$ increases from 25 to 125, the vortices are mainly confined to near the surface and the bulk fluid motion reduces significantly, as shown in Figure 4c. At a very high value of $\Omega$ (250 or higher), the bulk fluids are virtually motionless, despite the very fast oscillating flow occurring near the channel surface (Figure 4d). Thus, it is obvious from Figure 4 that the perturbed flow regime becomes smaller at higher nondimensional frequencies. The effect of nondimensional frequency on the flow field can be represented by damped viscous waves traveling away from the wall.

3.3. Effect of Surface Potential

Next, we studied the effect of surface potential distribution in the flow profile. As shown in Figure 5, the number of vortices formed within the channel can be controlled by the periodicity of the surface zeta potential distribution. With a half period ($q=\pi /L$), a large vortex is formed within the tube (Figure 5a), while with two periods ($q=4\pi /L$), four smaller vortices (each one extending one-fourth of the total length) are formed within the length of the tube (Figure 5b). Thus, a reduction in the surface potential periodicity reduces the number of vortices but extends the size of the vortices. Even though both scenarios (Figure 5a,b) provide an opportunity to enhance mixing through vortex formation, the higher periodicity case (Figure 5b) may be better for effective mixing because more vortices may yield a shorter mixing length. However, one has to be mindful about the sample loading to realize any positive outcome. For instance, in the case of half periodicity (Figure 5a), each mixing constituent needs to be loaded in each half of the channel. For example, fluid 1 should be loaded at length $0\le z\le L/2$ and fluid 2 should be loaded at length $L/2\le z\le L$. In contrast, in the case of double periodicity (Figure 5b), each mixing constituent needs to be loaded at one-eighth of the channel length in an alternative patch. For example, fluid 1 should be loaded at $0\le z\le L/8$, $L/4\le z\le 3L/8$, $L/2\le z\le 5L/8$, and $3L/4\le z\le 7L/8$, whereas fluid 2 should be loaded at $L/8\le z\le L/4$, $3L/8\le z\le L/2$, $5L/8\le z\le 3L/4$, and $7L/8\le z\le L$ for effective mixing. Thus, our analytical results show that one can precisely control the flow field to achieve the desired level of mixing by modifying the zeta potential patterning on the tube surface.

3.4. Pressure Distribution

Although no external pressure was imposed for our proposed micromixer, pressure can be induced inside the channel to ensure the constant flow rate between different regions with different zeta potentials on the wall. This induced pressure can be quantified from the velocity profile. An analytical expression for pressure can be obtained from the modified Navier–Stokes equations. For instance, in the absence of an advection term, one can replace the axial velocity, i.e., the z-directional velocity component (Equation (16b)), in the z-component of the Navier–Stokes equation (Equation (6b)) to obtain the z-directional pressure gradient:

$$\frac{\partial P}{\partial z}={\displaystyle \sum _{n=-\infty }^{\infty }{A}_{1,n}\left[i\rho \omega k_n{I}_0\left(k_{n}r\right)+\left(\mu l_n^3-\mu l_n{k}_n^2-i\rho \omega l_n\right)\frac{I_1\left(k_n{r}_0\right)}{I_1\left(l_n{r}_0\right)}{I}_0\left(l_{n}r\right)\right]}{e}^{i{k}_{n}z}{e}^{i\omega t}.$$

(20)

The pressure distribution along the tube can be obtained by integrating with respect to $z$ as follows:

$$P={\displaystyle \sum _{n=-\infty }^{\infty }{A}_{1,n}\left[i\rho \omega k_n{I}_0\left(k_{n}r\right)+\left(\mu l_n^3-\mu l_n{k}_n^2-i\rho \omega l_n\right)\frac{I_1\left(k_n{r}_0\right)}{I_1\left(l_n{r}_0\right)}{I}_0\left(l_{n}r\right)\right]}\frac{1}{i{k}_n}{e}^{i{k}_{n}z}{e}^{i\omega t}+P_1(r,t)$$

(21)

where $P_1$ is an integration factor, which can be a fixed value or a function of $r$ and $t$. Now, by manipulating Equations (6a), (16a), and (21), we obtain:

$$\frac{\partial P_1}{\partial r}={\displaystyle \sum _{n=-\infty }^{\infty }{A}_{1,n}\left[\left(-\rho \omega k_n+i\mu \frac{l_n^4}{k_n}-2i\mu k_n{l}_n^2+\rho \omega \frac{l_n^2}{k_n}+i\mu k_n^3\right)\frac{I_1\left(k_n{r}_0\right)}{I_1\left(l_n{r}_0\right)}{I}_1\left(l_{n}r\right)\right]}{e}^{i{k}_{n}z}{e}^{i\omega t}.$$

(22)

The aforementioned equation indicates that $\partial P_1/\partial r$ might be a function of $z$, $r$, and $t$. However, following numerical evaluation with a sinusoidal surface charge distribution (Equation (17)) and applied electric field, $\overrightarrow{E}=E_{0}sin\left(\omega t\right)=Im\left(E_0{e}^{i\omega t}\right)$, it is found that $\partial P_1/\partial r$ is zero in space and time (data not shown). Thus, the integration of Equation (22) will yield $P_1$ equal to a constant. Since we are not seeking the absolute pressure, we can neglect that constant, and quantify the gauge pressure distribution from Equation (21):

$$P(r,t)={\displaystyle \sum _{n=-\infty }^{\infty }{A}_{1,n}\left[i\rho \omega k_n{I}_0\left(k_{n}r\right)+\left(\mu l_n^3-\mu l_n{k}_n^2-i\rho \omega l_n\right)\frac{I_1\left(k_n{r}_0\right)}{I_1\left(l_n{r}_0\right)}{I}_0\left(l_{n}r\right)\right]}\frac{1}{i{k}_n}{e}^{i{k}_{n}z}{e}^{i\omega t}.$$

(23)

The distribution of pressure at various nondimensional times is shown in Figure 6. From this figure, it can be seen that the pressure varies mainly in the axial direction. The reason behind this is discussed shortly. Since the flow field varies with time, the induced pressure also varies accordingly. For any particular time, for instance, $\omega t=\pi /2$ (Figure 6a), an adverse (positive) pressure gradient is observed in the left half of the channel, while a favorable (negative) pressure gradient is formed in the right half of the channel. This induced positive pressure gradient tends to suppress electroosmotic flow, resulting in a mixed electroosmotic and pressure-driven flow, as shown in Figure 3a. This kind of velocity profile has also been obtained in an experimental study [35]. As the electric field direction switches, the flow field changes and, as a consequence, the pressure distribution is also switched, as shown in Figure 6b.

The axial variation of pressure is shown in Figure 7a for various nondimensional times, where pressure is obtained at the centerline of the tube. As seen from Figure 7a, for a sinusoidal variation of the zeta potential, the axial variation of pressure is also sinusoidal. However, there is a phase lag of 90° between the pressure distribution and zeta potential distribution. Thus, for a sine distribution of the zeta potential, the axial variation of induced pressure is cosine. At $\omega t=\pi /2$, an adverse pressure gradient is observed in the left half of the channel, whereas a favorable pressure gradient is observed in the right half of the channel. As the electric field changes, the pattern of pressure variation in the axial direction also changes. For example, in contrast to $\omega t=\pi /2$, at $\omega t=3\pi /2$ a favorable pressure gradient is observed in the left half of the channel while an adverse pressure gradient is observed in the right half of the channel. It should be noted that in this work, a pressure gradient was also observed in the radial direction, however, the variation was very small in comparison to the axial variation and, hence, it was not captured in the pressure contour plot (Figure 6). The radial directional variation of pressure is shown in Figure 7b by line plot for various nondimensional times, where pressure is obtained along the radial direction of the channel at the axial location $z/r_0=500$ ($z=L/2$). Since the radial variation of pressure is very small, instead of gauge pressure, its change from the center point ($r/r_0=0,z/r_0=500$) is plotted in Figure 7b. This figure shows that the radial variation of pressure is several orders magnitude lower than that of the axial variation due to the low radial velocity. As the electric field changes direction, the pattern of pressure variation in the radial direction is also changed. For example, at $\omega t=\pi /2$, the pressure variation in the radial direction follows a concave parabolic curve; whereas at $\omega t=3\pi /2$ the pressure variation in the radial direction becomes a convex parabolic curve (Figure 7b).

4. Conclusions

Motivated by the growing interest in electroosmosis as a reliable non-mechanical strategy to control fluid motion and mix species in microfluidic devices, we analytically studied the time-periodic electroosmotic flow in cylindrical microchannels with a heterogeneous surface charge distribution. The analytical solution was derived by solving a simplified Navier–Stokes equation with Helmholtz–Smoluchowski slip boundary conditions at the wall. Although the derived analytical solution is valid for any periodic surface charge distribution, only the sinusoidal surface charge case was analyzed in this paper. The induced pressure field was also obtained from the velocity profiles. The results show that several vortices are formed inside the microchannel with a sinusoidal surface charge. These vortices change their pattern and direction as the external electric field changes with time. In addition, the dominance of the vorticity depends on the frequency of the external electric field and the size of the channel. As the electric field frequency or channel diameter increases, vortices are shifted towards the channel surface and the perturbed flow regions become smaller. Unlike the homogenous surface charge case, velocity variation in the radial direction was also observed in this paper. The results also show that the number of vorticities can be changed easily by changing the periodicity of the surface charge. Thus, the flow field can be precisely controlled with a specific surface pattern, which will provide the desired mixing efficiency.