1 Introduction

Fluid flow in microchannels has been recognised as a major role in the applications of microelectro-mechanical systems (MEMS) and microfluidics, for instance, microchannel heat sinks for cooling microchips and laser diode arrays^[1-2], lab-on-chip devices for chemical and biomedical analyses^[3], and microfluid pumps^[4]. Several theoretical and experimental studies have been performed towards understanding the fundamental transport mechanisms of microchannel flow. Among them, Eringen^[5] suggested a simple model in which he described the fluid flow in the microchannels by using the Navier-Stokes equations with jump conditions. Pfahler et al.^[6] observed that the predictions by the Navier-Stokes equations are in agreement with the experimental results to some extent for the flow in the large channels, but different from the experimental data of the flow within small ones. Mala et al.^[7-8] discovered that the influence of the electric double layer (EDL) at the solid-liquid interface on the flow in a microchannel has to be considered to make the theoretical results consistent with the experimental ones. Ren et al.^[9] found that it is necessary to address the additional flow resistance due to the electrokinetic effect in microchannels. Their theoretical analysis based on the electro-viscous flow models gave better agreement with the experimental results compared with those by the Navier-Stokes equations.

It is worth mentioning that the electro-viscous flow model used by Mala et al.^[7-8] employed an arbitrary boundary condition to force the electrostatic potential ϕ in the middle of the microchannel to be 0, which leads to the fact that both the electrostatic potential and the flow velocity profiles in the microchannel are not physical. Moreover, the mathematical model for flow field used by Mala et al.^[7-8] as well as that by Ren et al.^[9] contains a constant pressure parameter whose value was set to be known from the experimental data. Though it significantly simplified the procedure to obtain results with acceptable accuracies, it is far from being a correct way to develop mathematical models for the flow in microchannels. Theoretically, as has been discussed in detail in the literature^[10-19], this pressure parameter has to be taken as an unknown constant to be calculated based on the boundary conditions.

In this paper, we develop a theoretical model of the liquid flow in a microchannel driven by its upper moving wall together with the effects of the EDL. Compared with the models used by Mala et al.^[7-8] and Ren et al.^[9] for the flow in the microchannel due to the pressure difference, two drawbacks have been overcome in our model. Firstly, we remove the unphysical artificial zero electric potential condition in the middle of the microchannel. Secondly, we calculate the unknown pressure constant in the model by using the boundary conditions on the microchannel walls, which makes our model consistent with the models that are commonly employed in the field of fluid mechanics to describe the channel flow. After this short introduction section, we describe the flow model for the microchannel flow due to the moving wall together with the electric double layer effect in detail in Section 2 and Section 3. We then give the detailed procedure to find the analytical solution of the governing equation by using the homotopy analysis method (HAM)^[20] in Section 4 that is also validated with the numerical results obtained by the finite difference technology^[21]. The flow phenomena are studied under different flow conditions in Section 5, which is followed by a short conclusion in Section 6.

2 Analytical solution of electrostatic potential

Consider a fluid flow between two horizontally parallel flat plates separated by a small distance H, as shown in Fig. 1. In practice, the liquid usually is a dilute solution that contains positive and negative ions. Also, the walls of the microchannel are commonly slightly charged. When the flow channel is small, the electrostatic effects on the flow cannot be ignored. In our model, we set that the electric potentials on both surfaces are identical and reduce gradually into fluid. The concentration of the suspended ions is uniform in the channel. In the vicinity of the surfaces the electrostatic potential ϕ at any point is associated with the net number of electrical charges per unit volume ρ_e, which can be expressed, based on the theory of electrostatics, in the following form:

(1)

where ε and ε₀ are, respectively, the dielectric constant of the fluid and the permittivity of vacuum.

Invoking the equilibrium Boltzmann distribution, in which the uniform dielectric constant is applicable, the solution of the number of ion distribution in a symmetric electrolyte can be then written as

(2)

where n_0i is the bulk ionic concentration, is the valence of type i ions, e is the fundamental charge, k_b is the Boltzmann's constant, and is the absolute temperature. The net charge density in a unit volume of the fluid is defined by

(3)

Substituting Eq. (3) into the Poisson equation (1), we obtain

(4)

By using the dimensionless variables

(5)

Eq. (1) and Eq. (4) are simplified to

(6)

(7)

where κ=H k is a dimensionless constant with being the reciprocal of the Debye length. When the electrical potential is small compared with the thermal energy of the ions, , Eq. (6) can be linearized to the Debye-Hückel model as follows:

(8)

The appropriate boundary conditions for the dimensionless electrostatic potential Φ are

(9)

where ζ denotes the measurable electrical potential between the boundaries of the compact layer and the diffuse layer of the EDL at the shear plates.

When using the boundary conditions in Eq. (9), the closed-form solution of Eq. (8) is

(10)

3 Equations of fluid motion

To describe the fluid flow motion between two parallel flat plates driven by its upper stretching plate with the effects of the EDL, as displayed in Fig. 1, the following assumptions are imposed: the fluid is incompressible, the flow is laminar and steady, and the velocity of the flow away from the channel walls is independent of position along the x-direction. When the viscous force, the pressure force, and the electrical body force generated by the double layer electric field are all considered, the Navier-Stokes equations to describe this microchannel flow can be written as

(11)

(12)

(13)

The boundary conditions are

(14)

where u and v are the velocity components along the x-axis and y-axis, respectively, p denotes the pressure, ρ_f is the fluid density, ν_f is the kinematic viscosity, and E_x is the electric field strength.

Define the following similarity transformations:

(15)

where ψ(x, y) is the stream function,

(16)

Substituting Eq. (15) into Eqs. (11), (12), and (13), the continuity equation (11) is satisfied straightforwardly, and the rest of equations can be written, by emphasizing the pressure related terms, as follows:

(17)

(18)

Since the terms on the right-hand side of Eq. (18) are independent of x, which indicates

As a result, we obtain

(19)

where

is the Reynolds number, and K is the pressure constant to be determined.

To obtain the electric field strength E_x, we quantify the balance between the streaming current and the electrical conduction current^[7-8]. As it is known that in the absence of an applied electric field, when a dilute solution is forced through a channel under hydrostatic pressure, the streaming current is generated, which can be expressed as

(20)

where A_c is the cross-sectional area of the flow channel.

The streaming potential further produces a conduction current in the reverse direction which is

(21)

where λ₀ is the electrical conductivity of the fluid, E_s is the streaming potential, and L is a reference length of the channel.

At the steady state, the net electrical current is zero, i.e.,

(22)

Substituting Eq. (15) into Eqs. (20) and (21), we obtain

(23)

(24)

where W is the reference width of the channel and E_x=E_s/L. Solving Eq. (22), we obtain

(25)

Therefore, the last term on the left-hand side of Eq. (19) can be written as

(26)

where

is a dimensionless constant.

Using Eq. (26), Eq. (19) can then be reduced to

(27)

with the boundary conditions

(28)

Two physical parameters are important to quantify the flow phenomena in the microchannels and for the engineering design of microchannels, namely, the local skin frictional coefficients on both walls. We give the definitions of these two parameters before showing the detailed solution procedure of Eq. (27):

(29)

where

(30)

Substituting Eq. (30) into Eq. (29), we obtain

(31)

where ξ=x/H is the dimensionless length variable.

4 Solution method

It is worth noting that Eq. (27) contains two unknown constants K and which are difficult to be calculated accurately during the solution finding procedure by numerical or analytical approaches without particular treatments. Nevertheless, in this work, we employ and extend the HAM technique to obtain the exact solution of the nonlinear equation (27) subject to the boundary conditions in Eq. (28), in which the exact values of the unknown constants, K and , are calculated spontaneously and simultaneously as parts of the solution procedure without any particular tricks.

We first write the unknown constants K and as

(32)

and then express the functions F(η) and Φ(η) as

(33)

In the framework of the HAM, the kth order HAM deformation equations can be written as

(34)

(35)

with the boundary conditions

(36)

where ħ_f and ħ_ϕ are the HAM auxiliary parameters used for the convergence-control of the HAM analytical approximations. Also, R_{f, m}, R_{ϕ, m}, and χ_m are defined, respectively, by

(37)

(38)

and

(39)

The solutions of Eqs. (34) and (35) are in the forms of

(40)

(41)

where

(42)

(43)

and the constants C_{i, m} (i=0, 1, 2, 3, 4) are determined by the boundary conditions

(44)

The above HAM solution procedure can work accordingly after the initial approximations f₀(η) and ϕ₀(η) are chosen properly based on the boundary conditions (9) and (28), such as

(45)

During the HAM solution procedure, the unknown constant K is determined in the following way. Due to Eqs. (34) and (37), it is known that f_m'(η)=0 contains the unknown term ω_m-1, which can be obtained by using the boundary condition f_m'(1)=0 defined in Eq. (36). For example, for the zeroth order approximation for K when m=0 in Eq. (32), as f₁'(1)=0, we have

(46)

which leads to

(47)

Also, the constant integration term can be obtained from Eq. (45). In this way, the whole solution series in Eq. (32) for K and the constant integration term can be determined successively from m=1, 2, 3, ....

5 Results

To check the accuracy of our solution, we define the following error estimation formula:

(48)

Substituting mth order computational results into Eq. (48), the corresponding error can be obtained. For example, for prescribed values of Re=5 and G₁=10, by properly choosing ħ_f =-2 and ħ_ϕ=-1/2 (Paé technique is used for κ = 20), we find that the errors are perfectly small, as listed in Table 1. Also, we run a few validations between the analytic results obtained by the HAM and the results by the numerical technique, finite difference method. When applying the finite difference method to solve governing equation (27) together with the boundary conditions in Eq. (28), we first take one derivation of Eq. (27) with respect of η. Then, the central finite difference schemes are employed to represent the differential terms in the resulting governing equation, and a ghost node on each end of the computational zone is used to take care of the first order derivations in the boundary conditions. The integral term is represented by Simpson's rule. The nonlinearity of the governing equation is dealt with through iteration. As shown in Fig. 2, good agreement has been found between the analytical approximations by the HAM and the numerical results by the finite difference method.

As seen in Fig. 3, Φ(η) decreases monotonously as κ increases. Particularly, as κ becomes sufficiently large, Φ(η) diminishes to zero in the middle of the channel. This leads to the common assumptions of zeta potential^[22] in the application of capillary: (ⅰ) the capillary walls are flat, and (ⅱ) the double layer develops fully so that the potential in the middle of the capillary is zero. Traditionally, people took that κ=5 could be large enough to describe such electrokinetic phenomena in fluid flows since the potential in the mid plane is small compared with the potential near the wall. However, Burgreen and Nakache^[23] showed that the influence of the zeta potential still exists across the channel even when κ=10 in which the distance between the two plates is already considerably large. Also, it is known that the double layer field is significant near the channel walls. For instance, the double layer occupies only 8% of the channel cross-sectional area for κ=80.

In Fig. 4, the influences of κ on the pressure constant K are presented. From this figure, we can see that the pressure constant K decreases dramatically as κ increases when it is small. Once reaching the minimum value, K increases as κ evolves. The pressure constant K approaches to 1 as κ keeps continuously increasing. This indicates that the influence of the electrostatic potential on the pressure gradient is opposite to that caused by the upper stretching wall. κ has significant effects on K when the distance between two parallel plates is small. As the distance between the two walls consecutively increases, the electrostatic potential in the middle of the channel diminishes to zero. In this situation, the effect of the electrostatic potential on the pressure gradient becomes weaker and weaker, which results in that the effects on the pressure gradient due to the upper stretching wall of the channel become dominant.

It should be noted that, from Eq. (31), F"(0) and F"(1) are the core parts of the local skin friction coefficients C_fxl and C_fxu, respectively. Therefore, we consider the effects of κ on F"(0) and F"(1) that will straightforwardly reflect the effects of κ on C_fxl and C_fxu. As shown in Fig. 5, F"(0) decreases rapidly as κ increases when κ is not large. After passing the minimum value, F"(0) grows with κ increasing. Physically, the double layer electric field generates an electrical body force, which is opposite to the frictional resistance. As a result, its enhancement decreases the drag resistance, while as κ is large enough, the effect of electrostatic potential becomes negligible in most part of the fluid within the microchannel. It is shown in Fig. 6 that the variation of F"(1) with κ exhibits the totally reverse trend relative to that of F"(0). F"(1) rapidly enhances to its maximum value, then gradually reduces as κ continuously reduces. This means that the electrical body force exerts a positive influence on the skin friction coefficients for small κ, whose effect becomes weaker and weaker when the size of the channel increases. The streaming potential

is also an important quantity to affect the flow patterns in a microchannel. As shown in Fig. 7, E_s develops very quickly as κ evolves. When the electrostatic potential in the middle of the microchannel just diminishes to zero, the streaming potential reaches its peak value. After that, E_s reduces as κ increases. Physically, the reason is that the large distance of electrokinetic separation corresponds to a large volume transport so that more ions are carried to the end of the channel which leads to a higher charge accumulation.

It is a common sense that the Reynolds number Re plays an important role in the flow characteristics. As shown in Fig. 8, when κ is small, K decreases dramatically as Re enlarges at first, and then it enhances slightly as Re continuously evolves. However, when κ is large (i.e., κ=20), K decreases monotonously as Re increases. These trends show that Re plays a dominant role on K when it is small. When Re is large, the effect of the electrical body forces becomes important, especially when κ is not large. It is known from Eq. (31) that f"(0)/Re and f"(1)/Re are the main parts of the local skin friction coefficients C_fxl and C_fxu, respectively. As shown in Fig. 9, f"(0)/Re increases gradually as Re grows for both cases of κ=3 and κ=20. On the contrary, f"(1)/Re decreases consecutively as Re increases for both cases of κ=3 and κ=20, as displayed in Fig. 10. For both f"(0)/Re and f"(1)/Re, the effects of κ become significant when Re is sufficiently large. As seen in Fig. 11, E_s reduces gradually as Re enlarges, particularly, the reduction rate for κ=3 is faster than that for κ=20. Physically, when κ is appropriately small, the electrostatic potential affects the whole cross-section of the channel. In this case, κ plays an equally important role in the reduction of E_s as Re does. When κ is sufficiently large, its effect on the flow is negligible. Re becomes the dominant factor for the decrease of E_s.

The effects of G₁ on the various physical quantities of the flow have been shown in Table 2 and Table 3. It can be seen that K and f"(0) decrease as G₁ grows, but f"(1) and E_s rise as G₁ increases. Physically, G₁ is related to the electrical strength whose enlargement results in the increase of the electric potential. Since the pressure gradient caused by the electrostatic potential is opposite to that caused by the upper stretching wall, the increase of G₁ leads to the reduction of K. On the other hand, the electrical body force is along the same direction of the viscous force near the upper wall of the channel. As a result, f"(1) enhances as G₁ increases. Due to the conservation of mass flux, f"(0) has to diminish as G₁ grows. Obviously, the increase of the electrostatic potential is helpful to increase E_s. From these tables, we also notice that all these physical quantities vary almost linearly with G₁ from 0 to 200.

6 Conclusions

In this work, we proposed a mathematical model to describe the liquid flow in a microchannel driven by its upper stretching wall in the combination of electrokinetic effects. Comparing with the traditional models, the unphysical boundary condition of the electrostatic potential in the middle of the channel has been fully eliminated. Our model was solved using the similarity transformation. The accurate analytical approximations were given by the HAM, in which the pressure constant and the integral constant related to the electric potential were solved spontaneously by using the proper boundary conditions on the channel walls. It is expected that our model is consistent with the commonly accepted models in the field of fluid mechanics, and can offer the general and proper way to study the flow phenomena in microchannels.

Acknowledgements This work is supported in part by the Australian Research Council through a Discovery Early Career Researcher Award to Qiang SUN.