1.Introduction

Electroosmotic flows have found wide applications in the development of a great variety of microfluidic systems consisting of valves, pumps, and mixers to be utilized as an efficient method for transporting micro volumes of fluids. Examples of such applications are drug delivery, DNA analysis, and biological/chemical agent detection sensors on microchips ^1-3. In this context, electroosmosis enables fluid pumping and flow control using electric fields, eliminating the need for mechanical pumps or valves with moving components. However, because of these applications, some times there is a need of modifying the cross-section area of the microchannels, for speeding up samples transported by electrokinetic effects ⁴, for transporting particle in micro/nano nozzles and diffusers ⁵, for determining the translocation speed of DNA in nanopores ⁶.

One of the first researchers that analyzed the EOF in microchannels with the simultaneous presence of charge and shape modulations on the surface in the direction of the applied electric field was Adajari and coworkers ^{7 8}. Recently, Alexander et al.⁹ studied the electroosmotic flow in wavy channels by expanding the solution into a double series in terms of the dimensionless amplitudes and of the dimensionless zeta potential for a binary dilute electrolyte. They demonstrated the importance of the varicose or sinuous character of the channels as well as the role of high frequency roughness. Meanwhile, Xia et al. ¹⁰ analyzed the electroosmotic flow in a channel bounded by a plane wall and a sinusoidal wall. They obtained an exact solution by using a complex function formulation together with the boundary integral method. The effects of the channel width and wave amplitude on the electric field, streamline pattern, and flow field were presented. More recently, Martínez et al. ¹¹ analyzed the electroosmotic flow of a viscoelastic fluid in a wavy wall microchannel. They obtained an asymptotic solution based on the domain perturbation method, which can be applied only for very small amplitude waviness.

On the other hand, some numerical works where EOF in wavy microchannels have been analyzed, are the following: ¹² investigated the mixing characteristics of electrokinetically-driven flow in microchannels with different wavy surface configurations. They conducted numerical simulations to analyze the influence of the wave amplitude and the length of the wavy section on the mixing efficiency within the microchannel. ¹³ carried out a numerical investigation to study the flow characteristics of non-Newtonian fluids (described by a power-law fluid) in rough microchannels defined by a complex-wavy surface, composed by the superimposition of two sinusoidal functions. The effects of the wave amplitude, geometry of the wave, Debye-Hückel parameter, and flow behavior index of the non-Newtonian fluids on the local velocity profiles, volumetric flow rate, and electric field distribution were examined. ¹⁴ analyzed the electroosmotic flow in a parallel-plate microchannel with sinusoidal surface roughness, and the solution of the governing equations was obtained by using the finite element method. Their simulation results indicate that the bulk flow velocity and the volumetric flow rate decrease slowly with the roughness height when the relative roughness is very small or very large, while decrease quickly when the relative roughness is moderate.

The most of theoretical studies concerning EOF consider that the walls of microchannels are perfectly smooth; however, a careful analysis reveals in reality that the walls are rugous, and therefore it geometrical factor plays a very important role in microscales which must be taken into account. Evaluating the perturbations on the flow due to such cause is a fundamental fluid mechanics problem of considerable interest in this kind of EOF, as demonstrated by ^{15 16}. Many of the mentioned works use the Helmholtz-Smoluchowski slip approximation under the assumption of infinitely thin electric double layers (EDLs) to simplify the determination of the flow field. Because of the assumption of thin EDLs used in the aforementioned works, the obtained velocity profiles are uniform in the transverse direction. Conversely to that, and given different applications that can be found in EOF, in the present work, we consider EDLs with finite thickness, which yields an EOF with a nonuniform velocity profile. We also determine the induced pressure along the channel, which cannot be obtained by an analysis of infinitely thin EDLs.

The mentioned works enable the research community to obtain significant physical insights on EOF in microfluidics channels; however, most of them have been performed by using numerical methods or complex mathematical formulations. Therefore, this work aims to present a simple analytical formulation, based on the lubrication theory that permits to determine the flow field of an electroosmotic flow in a microchannel with wavy walls in an easier way. The solution was asymptotically obtained in two limits: in the first case, we solve the mathematical problem by considering thin EDLs with values of the waviness amplitude of the same order as the microchannel height; while that in the second case, we conduct the analysis in the limit of very small waviness amplitude but the Debye layer being of order as the microchannel height. Besides a numerical solution was conducted to verify the asymptotic solution.

2.Formulation

For mathematical modeling, we consider the EOF in a wavy-wall microchannel as shown in Fig. 1. The length of the microchannel is L, and it is assumed that the walls are located at

Figure 1 Sketch of the EOF in a wavy-wall microchannel. 

Here, H represents the mean semi-channel height, and ε is a dimensionless geometrical parameter, which characterizes the waviness of the microchannel walls, and it can take values of 0≤ε≤0.8; we anticipate that this last value is the upper limit for which the solution obtained in this work exists. In this case, the cross-sectional area of the microchannel is varying periodically in the flow direction. Although the method used for solving this kind of EOF can be applied to other geometrical designs of the walls, in this work is assumed that the crest of the upper wall corresponds to the trough of the other.

We consider that the flow is only driven by the effect of electro-osmotic forces caused by imposing an external electric field of strength E₀ in the axial direction. The microchannel is supported at the two ends by two liquid reservoirs, which are found at a pressure P₀. According to typical applications of EOF, L≫H. Owing to the geometry, a two-dimensional rectangular coordinate system x, y is adopted with the origin at the microchannel inlet, x denoting the flow direction, and 𝑦 represents the transverse direction in the microchannel. Here, the flow will be 2D; therefore, the velocity component perpendicular to the plane (x,y) is zero.

On the other hand, and considering wavy-walls of the microchannel, the assumption of unidirectional flow is no longer valid, except in the case of ε = 0. Because of the fluid is incompressible, the velocity component in the x- direction should be increased for decreasing values of the dimension H in order to satisfy the mass conservation; a consequence of the above, is the appearance of an induced pressure the gradient along the microchannel. Besides other assumptions on which the mathematical model is obtained are the following: (i) constant physical properties and steady-state conditions; (ii) the electrolytic solution is assumed symmetric (z:z); (iii) the wall potentials ζ are axially invariant and low enough (≤25 mV), such that the Debye-Hückel linearization approximation is valid; (iv) the electric double layers (EDLs) on the inner surface of the microchannel do not overlap; (v) it is assumed that the electric field and current vectors make only a small angle to the x-axis, which is valid in the limit |h'(x)|≪1 along the x, where the prime denotes derivative , and therefore the lubrication theory can be used for solving this EOF ¹⁵.

3.Governing equations

Under the assumptions presented previously, the equations which govern the steady two-dimensional EOF in a wavy-wall microchannel are the continuity, momentum, and Poisson-Boltzmann equations. We start with the dimensionless governing equations, based on the lubrication theory, which are given by :

and

Here, u¯ and v¯ are the dimensionless velocities in the dimensionless x- longitudinal and Y- transverse coordinates, respectively; P¯ represents the dimensionless pressure. E¯x and ψ¯ denote the electric field in the x-direction and the electric potential in the EDL in dimensionless form, respectively. α=ζ/LE0 and κ¯=Hκ, with κ denoting the Debye-Hückel parameter, which is related to the Debye length λ_D, defined as κ=λD-1=2n∞z2e2/ϵkBT1/2 ¹⁸, where n∞, z, e, kB, and T are the bulk concentration of ions, the valence, the magnitude of the fundamental (elementary) charge on an electron, the Boltzmann constant, and the absolute temperature, respectively.

The dimensionless variables in Eqs. (2)-(5) are defined as follows:

where u and v represent the velocity components in the x- and y- directions, respectively; U_HS is the Helmholtz-Smoluchowski velocity, defined as UHS=-ϵζE0/μ, being μ the viscosity and ϵ the permittivity of the fluid. E_x represents the electrical field strength in the x-direction, and ψ is the electric potential in the EDL. Also, p̃=P-ϵκ2ψ2/2 (7)(8)(17), with P representing the hydrodynamic pressure. From the current continuity through any cross section orthogonal to x-axis, E¯χ can be written as E¯χ=1/h¯(χ) , where h¯(χ)=1+εsin(2πχ).

The dimensionless boundary conditions are:

and

From Eq. (4), P¯ is exclusively a function of the coordinate X; thus, the hydrodynamic problem consists of determining the solution of Eqs. (2) and (3). As can be appreciated, these equations are similar to that of unidirectional flow. However, in this case, due to variations of the cross-section of the microchannel, u¯ will depend on the dimensional coordinate X and v¯ will not be zero in general. Besides, because of u¯ is a function of the dimensionless coordinate X, and so is P¯.

In Eq. (3),

which satisfies Eq. (5), together with the third and second boundary conditions of Eqs. (7) and (8), respectively.

In Eq. (3), the last term on the right hand side can be neglected because of some values of the parameters involved in the definition of a can assume values such as ζ≤0.25mV, E0=1-100 V/mm, and L∼O(10-2) m. Therefore, α∼10-4. Thus, the solution of Eq. (3), subject to the no slip boundary condition u¯=0 at Y=h¯(χ); and the symmetry boundary condition ∂u¯/∂Y=0 at Y = 0 given, respectively, in (7) and (8), is given by

where the pressure gradient, dP¯/dχ, is unknown and can be obtained with the aid of the continuity Eq. (2). Substitution of Eq. (12) into Eq. (2), integrating with respect to Y, and applying the impermeability condition v¯=0 at Y=-h¯χ given in Eq. (9), we obtain

At Y=h¯(χ), v¯=0, thus, from Eq. (13) we obtain the differential equation that allows determining the pressure field and is given by:

Equation (14) must be solved subject to boundary conditions in Eq. (10), where P¯=0 at X = 0,1. In the following subsection, we obtain the solution for thin EDLs ( κ¯≫1).

4.Results and discussion

We have complemented the asymptotic analysis carried out in this work by obtaining the numerical solution of Eq. (14), and therefore the flow field. For solving Eq. (14), the well-known Shooting method was used ²¹. For the numerical integration, Δχ=ΔY steps of 1 X10^-2 have been used in all numerical runs. In the following figures, symbols represent the numerical solution, while lines correspond to the asymptotic solutions, showing an excellent agreement between them.

Figures 2 and 3 show the dimensionless velocity distribution as function of the dimensionless transversal coordinate Y for κ¯=100, different values of ε(=0.1,0.3,0.5,0.7,0.8), evaluated at χ=0.2 and 0.8, respectively. By comparing both figures, it can be appreciated that the velocity increases and decreases when the microchannel cross-section decreases and increases, respectively, to guarantee the mass conservation. Besides the velocity profiles, evaluated at X = 0.2, are weakly concave and are weakly convex when are evaluated at X = 0.8. These behaviors are present always that the cross-section increases and decreases, respectively. Also, from these figures, the flow is strongly accelerated (in regions where the cross-section of the channel diminishes) at X = 0.8, in comparison when is evaluated at X = 0.2.

Figure 2 The dimensionless velocity in the χ direction as function of the dimensionless transverse coordinate Y , evaluated at χ = 0:2, for different values of ε (= 0:1; 0:3; 0:5; 0:7; 0:8), with κ- = 100. 

Figure 3 The dimensionless velocity in the χ direction as function of the dimensionless transverse coordinate Y , evaluated at χ = 0:8, for various values of ε(= 0:1; 0:3; 0:5; 0:7; 0:8), with κ- · = 100. 

In Fig. 4 we plot the dimensionless velocity profile v¯, as a function of the dimensionless longitudinal coordinate X, evaluated at various values of the dimensionless transverse coordinate. In this figure, it is shown that the movement of the fluid in the direction Y is alternating from positive to negative values as the fluid flows in the axial direction X, according to the geometry of the walls. When the velocity v¯ is evaluated through the thin EDL solution at Y=h¯χ, and compared against the numerical solution, there is a small discrepancy, which would violate the impermeability condition; however, the above occurs because we have neglected the terms containing exp-2κ¯h¯ in Eq. (17) for obtaining Eq. (18).

Figure 4 The Y -component velocity, Eq. (13), plotted against χ for various values of Y . 

The mentioned concave and convex behaviors of the velocity profiles are due to the induced pressure gradient along the microchannel by the effects of varying the cross-section. In Fig. 5a), the dimensionless pressure is plotted, and in Fig. 5b), the corresponding pressure gradient is shown, for different values of the parameter ε, with κ¯=100. It is obvious that P¯→0 as ε→0, as expected, because this corresponds to the EOF in parallel flat plates microchannel case. Also, P¯∼Oκ¯-1, and therefore, the pressure in the microchannel can dissapear in the limit of κ¯→∞ for any value of 0≤ε<1. The above is evident from Eq. (20), where the pressure is inversely proportional to κ¯. On the other hand, when the amplitude of the wall waviness increases, i.e., ε→1, this condition yields negative pressures over extensive portions of the microchannel.

Figure 5 a) The dimensionless pressure distribution and b) the pressure gradient plotted as functions of the dimensionless coordinate χ, for values of ε(= 0:1; 0:2; 0:3; 0:4; 0:5), with k̄ · = 100. 

Referring to the geometry of the microchannel walls, and for the used values of the parameters shown in the figure, we see that there is a maximum positive of the dimensionless pressure P¯ at χ=0.64 and a corresponding minimum negative at χ≈0.8. On the other hand, when the pressure gradient is zero, the velocity is uniform (plug-like). Another important aspect is that the pressure has negative values in the region where the cross-section decreases in the direction of the flow of the microchannel and positive values where the cross-section is wider.

Although in principle the approximate solution for the flow field was obtained in the limit of k̄ ≫1, this remains valid even for relatively small values of κ¯, as shown in the following figures. In Fig. 6, the velocity profiles u¯ at X = 0.2 are plotted for various values of the parameter ε, and k̄ =5. It is clear that the flow behavior is similar to that presented in the previous figures. However, it is necessary to emphasize that decreasing the value of k̄ , the velocity profiles are not uniform across the transversal section of the microchannel, adopting a concave form. In this case, for ε = 0.1, the velocity profile has a paraboloid shape. For ε= 0.5, the velocity is more concave near to the center of the microchannel. This latter behavior is a consequence of the induced pressure is greater in comparison with the case of ε=0.1, as shown later in Fig. 8.

Figure 6 Dimensionless velocity profile u- as function of the dimensionless transverse coordinate Y , evaluated at χ = 0:2, for different values of ε (= 0:1; 0:2; 0:3; 0:4; 0:5) with κ- · = 5. 

Figure 7 Dimensionless velocity profile u- as function of the dimensionless transverse coordinate Y , evaluated at χ = 0:8, for different values of ε(= 0:1; 0:2; 0:3; 0:4; 0:5), and κ- · = 5. 

Figure 8 The dimensionless pressure distribution and the pressure gradient plotted as functions of the dimensionless coordinate χ, for values of ε(= 0:1; 0:2; 0:3; 0:4; 0:5), and κ- · = 5. 

Similar to the previous Fig. 6, in 7 we show the velocity profiles evaluated at χ≈0.8, where a concave behavior is appreciated. The coordinate χ≈0.8 corresponds to the region where the pressure gradient is negative (See Fig. 8), meaning that the induced pressure decreases in the direction of the flow.

Figure 9 shows the validity range of the thin EDL and numerical solutions for the flow rate as a function of κ¯ and ε. The percentage error in the flow rate is defined as

Figure 9 Contour plot of the relative percentage error Er(%) in the flow rate as a function of κ- · and ε. 

The numerical flow rate is Q¯num = ∫0h¯(χ)u¯dY, and the velocity field u¯ is given by Eq. ([e11]). The approximate solution Q¯edl is defined in Eq. ([thin7]). We can note from Fig. 9 that the blue region corresponds to Er (%) = 0 in the flow, which means a large range of the parameters κ¯ and ε, where Q¯edl is valid. The violet region corresponds to an error ≈0.1 and the pink region to 10%. The last region allows the determination of the flow rate for κ¯≈5 and ε≈0.7, with a relative error of 10%. Green and red regions have no physical sense for values of κ¯<4 when ε→0.8, due to negative values in the flow rate are obtained. In this context, good agreement in results is obtained with κ¯≥4 using 0<ε≤0.8.

The complexity of obtaining the exact solution of the differential equation Eq. (14) for any value of κ¯ is difficult. However, in the limit of κ¯≫1, we have determined an approximate solution, for the velocity, pressure, and volumetric flow rate, which are given by Eqs. (15), (20) and (21), respectively. In Fig. 10, the approximate solution for the volumetric flow rate, in the limit of κ¯≫1, is shown and is compared against the numerical solution. As can be seen, the approximate solution provides an excellent agreement even for values of κ¯=5 when ε≤0.7 with an Er=0.1%, as shown in Fig. 9.

Figure 10 Comparison between the numerical and thin EDL solutions of the volumetric flow rate in the limit of κ-≫1. The thin EDL solution is given by Eq. (21). 

On the other hand, for values of ε≪1, we have applied the domain perturbation method (see Appendix A for details) for determining the flow field (u¯,v¯,P¯,Q¯), where the EDL solution of κ¯≫1 is no longer valid. In Fig. 11 a comparison between numerical and the perturbation domain method solutions is presented, and a good agreement between them is found in the range of 0<ε≤0.1 and 1≤κ¯≤100, with a maximum Er(%) = 5.

Figure 11 Contour plot of the relative percentage error Er(%) in the flow rate determined by the domain perturbation solution, as a function of κ- · and ε. 

For example, for the used value of ε = 0.03, the domain perturbation solution is indistinguishable from the numerical solution. It is worth noting that the obtained thin EDL solution in the limit of κ¯≫1 is still valid for values of κ¯>2.5. In Fig. 12 we show the thin EDL and the domain perturbation solutions for the dimensionless volumetric flow rate, given by Eqs. (21) and (A.17).

Figure 12 Numerical solution of the dimensionless volumetric flow rate as a function of the parameter κ- ·, for ε = 0:03. Also shown are the two asymptotes, (21) and (A:17). 

Finally, a plot of 1-Q¯κ¯ as a function of ε is given in Fig. 13. It is evident that Q¯ diminishes for increasing values of ε. Of course, in the limit of κ¯→∞, Q¯→1, as can be appreciated from Eq. (21). It should be mentioned that this figure is valid for values of κ¯≫2.5. In particular, for values of ε>0.8 (this value was found by the relative percentage error Er(%) in the flow rate, shown in Fig. 9), the thin EDL solution is no longer valid because of the limit of Eq. (21) as ε→1 is undetermined. A physical interpretation is that the wavy-microchannel is obstructed at χ=0.75, when ε→1. This condition has no physical sense, then ε must be less than unity. According to Fig. 9, the upper limit of ε in Fig. 13 can be established as 0.8 for values of 5≤κ¯≤100, and the value of κ¯(1-Q) is 4.12458 with ε = 0.8.

Figure 13 The dimensionless volumetric flow rate as a function of parameter ε, which characterizes the waviness amplitude. 

5.Conclusion

In this work, we have conducted a theoretical analysis of an EOF in a wavy wall microchannel of a Newtonian fluid, by using the lubrication approximation, under the Debye-Hückel approximation. Two approximate solutions were determined: in a first case we consider κ¯≫1 with ε∼O(1), and in the second case, ε≪1 with finite values of κ¯. These asymptotic solutions were compared against a numerical solution, and an excellent agreement was found. The effects of the waviness were investigated, and we have shown that the most important factors that affect the flow field are the parameter related to the amplitude of the microchannel waviness, ε, and the ratio of the mean semi-channel height to the Debye length, κ¯.

This theoretical analysis can predict general trends in the data and basic aspects of the observed flow field in EOF in wavy wall microchannels. Future work will involve the effect of nonuniform zeta potential. Although the present analysis considers that the geometry of the microchannel is such that the crest of a wall corresponds to the through of the other, this can be applied to analyze other configurations: (i) that the crest of a wall corresponds to the crest of the other wall of the microchannel; (ii) one of the walls has a phase-advance/lag respect to the other. Besides it is very important to complement this study with a thermal analysis, to understand the effect of the waviness on the heat transfer process, due to the inevitable Joule heating effect, which is present in EOF.