An Exact Solution for Power-Law Fluids in a Slit Microchannel with Different Zeta Potentials under Electroosmotic Forces

Electroosmotic flow (EOF) is one of the most important techniques in a microfluidic system. Many microfluidic devices are made from a combination of different materials, and thus asymmetric electrochemical boundary conditions should be applied for the reasonable analysis of the EOF. In this study, the EOF of power-law fluids in a slit microchannel with different zeta potentials at the top and bottom walls are studied analytically. The flow is assumed to be steady, fully developed, and unidirectional with no applied pressure. The continuity equation, the Cauchy momentum equation, and the linearized Poisson-Boltzmann equation are solved for the velocity field. The exact solutions of the velocity distribution are obtained in terms of the Appell’s first hypergeometric functions. The velocity distributions are investigated and discussed as a function of the fluid behavior index, Debye length, and the difference in the zeta potential between the top and bottom.


Introduction
Recently, microfluidic device applications are increasing in the fields of chemical analysis, medical diagnostics, material synthesis, and others [1][2][3]. In the field of microfluidics, flow control in a microchannel is one of the most important issues. The problem with conventional pressure-driven flow is that as the channel size decreases, the hydraulic area becomes extremely small, resulting in a significant increase in the corresponding hydraulic resistance [4]. Electroosmotic flow (EOF) does not suffer from this problem because it is the motion of fluid that depends on the electric field across a microchannel [5][6][7].
Many efforts have been made to study electroosmotic flow (EOF) using Newtonian fluids. However, a few microfluidic devices are used more frequently for processing biological fluids such as blood, saliva, DNA, and polymer solutions, which cannot be treated as Newtonian fluids. To analyze the EOF of such fluids, an approach to non-Newtonian constitutive relations should be considered [8][9][10][11][12][13][14][15].
Among the various constitutive laws for non-Newtonian fluids, the power law model is the most popular because of its simplicity and suitability for analyzing a wide range of fluids. Thus, many researchers have conducted EOF studies using the power law model [16][17][18][19][20][21][22][23][24][25]. Zhao et al. analyzed the EOF of power-law fluids and obtained the approximate solution of the velocity field in a slit microchannel [16]. In addition, they studied a general Smoluchowski slip velocity over a surface [17] and provided an exact solution of the velocity distribution in a slit microchannel [18]. They also analyzed the EOF of power-law fluids in cylindrical [19] and rectangular [20] microchannels. Tang et al. conducted a numerical study of the EOF in microchannels of a power-law fluid using the lattice Boltzmann method [21]. Vasu and De analyzed a mathematical model of the EOF of power-law fluids in a rectangular microchannel at high zeta potential [22]. Babaie et al. and Hadigol et al. numerically analyzed the EOF of power-law fluids in a slit microchannel with pressure gradient [23,24]. Ng and Qi developed a simplified analytical model to describe the electrokinetic flow of a power-law fluid for varying wall potentials and channel heights in a slit channel [25].
Most previous studies have been performed on microchannels with the same zeta potential at the top and bottom walls. To our knowledge, however, many microfluidic devices are made with a combination of different materials, such as silicon dioxide (glass) as the base and polydimethylsiloxane as the top and side-walls. In these cases, asymmetric electrochemical boundary conditions should be applied for reasonable analysis of the EOF. Afonso et al. [26] and Choi et al. [27] analyzed the EOF of viscoelastic fluids in a microchannel with asymmetric zeta potential using the simplified Phan-Thien-Tannar model. Qi and Ng investigated the EOF of a power-law fluid through a slit channel where the walls were asymmetrically patterned with periodic variations in shape and zeta potential [28]. Hadigol et al. numerically investigated the characteristics of electroosmotic micromixing of power-law fluid in a slit microchannel with nonuniform zeta potential distributions along the walls of the channel [29]. Jiménez et al. investigated the start-up from rest of the EOF of Maxwell fluids in a rectangular microchannel with asymmetric high zeta potentials at the walls [30]. Peralta et al. conducted theoretical analysis of the start-up of oscillatory EOF in a parallel-plate microchannel under asymmetric zeta potentials [31]. Recently, Choi et al. presented the EOF in a rectangular microchannel using numerical analysis [32] and suggested an approximate solution for the EOF of power-law fluid with asymmetric zeta potential of a planar channel [33].
Obtaining the exact solution not only provides physical insight into the phenomena but can also serve as a benchmark for experimental, numerical, and asymptotic analyses. Zhao and Yang [18] have reported the exact solution for the EOF of a power-law fluid with a symmetric zeta potential. However, it remains a challenge to obtain the exact solution of the EOF in microchannels with asymmetric electrochemical boundary conditions.
In the present study, exact solutions for EOFs of power-law fluids in a slit microchannel with different zeta potentials at the top and bottom walls are presented. In addition, the key parameters affecting the velocity distribution of EOF, including the fluid behavior index, Debye length, and different zeta potentials at the top and bottom walls are analyzed. Figure 1 shows a two-dimensional EOF in a slit microchannel of height 2H. The top and bottom walls were charged with zeta potential ψ t and ψ b , respectively. An external electric field E 0 was applied to a power-law fluid with a constant density ρ and electric permittivity .  [23,24]. Ng and Qi developed a simplified analytical model to describe the electrokinetic flow of a power-law fluid for varying wall potentials and channel heights in a slit channel [25]. Most previous studies have been performed on microchannels with the same zeta potential at the top and bottom walls. To our knowledge, however, many microfluidic devices are made with a combination of different materials, such as silicon dioxide (glass) as the base and polydimethylsiloxane as the top and side-walls. In these cases, asymmetric electrochemical boundary conditions should be applied for reasonable analysis of the EOF. Afonso et al. [26] and Choi et al. [27] analyzed the EOF of viscoelastic fluids in a microchannel with asymmetric zeta potential using the simplified Phan-Thien-Tannar model. Qi and Ng investigated the EOF of a power-law fluid through a slit channel where the walls were asymmetrically patterned with periodic variations in shape and zeta potential [28]. Hadigol et al. numerically investigated the characteristics of electroosmotic micromixing of power-law fluid in a slit microchannel with nonuniform zeta potential distributions along the walls of the channel [29]. Jiménez et al. investigated the start-up from rest of the EOF of Maxwell fluids in a rectangular microchannel with asymmetric high zeta potentials at the walls [30]. Peralta et al. conducted theoretical analysis of the start-up of oscillatory EOF in a parallel-plate microchannel under asymmetric zeta potentials [31]. Recently, Choi et al. presented the EOF in a rectangular microchannel using numerical analysis [32] and suggested an approximate solution for the EOF of power-law fluid with asymmetric zeta potential of a planar channel [33].

Mathematical Formulation
Obtaining the exact solution not only provides physical insight into the phenomena but can also serve as a benchmark for experimental, numerical, and asymptotic analyses. Zhao and Yang [18] have reported the exact solution for the EOF of a power-law fluid with a symmetric zeta potential. However, it remains a challenge to obtain the exact solution of the EOF in microchannels with asymmetric electrochemical boundary conditions.
In the present study, exact solutions for EOFs of power-law fluids in a slit microchannel with different zeta potentials at the top and bottom walls are presented. In addition, the key parameters affecting the velocity distribution of EOF, including the fluid behavior index, Debye length, and different zeta potentials at the top and bottom walls are analyzed. Figure 1 shows a two-dimensional EOF in a slit microchannel of height 2H. The top and bottom walls were charged with zeta potential ψt and ψb, respectively. An external electric field E0 was applied to a power-law fluid with a constant density ρ and electric permittivity ϵ. The velocity field in the microchannel is governed by the continuity and Cauchy momentum equations given as: The velocity field in the microchannel is governed by the continuity and Cauchy momentum equations given as:

Mathematical Formulation
where v = (u, v) is the velocity vector, p is the pressure, τ is the stress tensor, and F = (F x , F y ) is the body force. The stress tensor can be given by where µ is the effective viscosity. For a steady, fully developed, unidirectional flow with no applied pressure and negligible gravitational force, the body force acts only along the x-direction and the Cauchy momentum equation in Equation (2) can be simplified as: The effective viscosity of the power-law fluid can be expressed as: where m is the flow consistency index, and n is the flow behavior index. The body force along the x-direction is given by: The net charge density ρ e can be obtained by the Poisson equation, which takes the form of: With the assumption of Boltzmann distribution and small zeta potentials, the electrical potential profile in the electrical double layer (EDL) is governed by the linearized Poisson-Boltzmann equation expressed by: which is subject to the following boundary conditions: κ −1 is called the Debye length and is defined as κ −1 = (εk B T/2e 2 z 2 n ∞ ) 1/2 , where n ∞ and z are the bulk number concentration and the valence of ions, respectively, e is the fundamental charge, k B is the Boltzmann constant, and T is the absolute temperature. The solution for the electrical potential distribution is of the form: Then, the net charge density ρ e can be expressed as a function of the EDL potential ρ e (y) = −κ 2 εψ(y).
With all the aforementioned considerations, the Cauchy momentum equation in Equation (4) This equation is constrained by the following boundary conditions (no-slip conditions) Substituting Equation (10) into Equation (12) yields: Since most of the materials that make up the microchannels have negative zeta potential [34,35], the wall zeta potentials are assumed to be negative in the present study; thus, the flow occurs in Integrating Equation (14) with y leads to: where ψ m and R are the average zeta potential and the dimensionless zeta potential difference between the top and bottom walls, respectively, which are defined by: and C + and C − are integral constants. Both Equations (15a) and (15b) should be zero at y = y c . Therefore, Integrating Equation (15) with the corresponding boundary condition in Equation (13) leads to the velocity distribution: where By integrating Equation (19), the velocity distribution can be obtained as: where u s denotes the generalized Smoluchowski velocity for power-law fluids by employing the average zeta potential ψ m at the top and bottom walls on the basis of the work of Zhao et al. [16], which is expressed as: and n F 1 1 + 1 n ; 1 2 , 1 2 ; 2 + 1 n ; n+1 n F 1 1 + 1 n ; 1 2 , 1 2 ; 2 + 1 n ; The integral constant C can be obtained from the following equation: where J(x) are defined by: The details of the mathematical derivations are described in Appendix A. It is a challenge to obtain the explicit form for the integral constant C. Thus, in this study, a numerical method was used for evaluating C. F 1 (a; b 1 , b 2 ; c; x, y) in Equation (23) is the Appell's first hypergeometric function [36], which can be represented as a one-dimensional integral form [37]: where α, β, and γ are real values, and Γ(z) is a gamma function. It is evident from Equation (28) that, although the Appell's first hypergeometric function in Equation (23) has complex arguments, it always has a real value. Alternatively, Equation (21) can be expressed in the single form using Equation (15) as follows: where V(y) ≡ − 1 κ n n+1 1 √ C 2 +w 2 du dy [I(y) + C]F 1 1 + 1 n ; 1 2 , 1 2 ; 2 + 1 n ; Equation (21) is applicable to the EOF of power-law fluids with arbitrary zeta potentials at the top and bottom walls. If the top and bottom walls have the same zeta potential, then the velocity distribution is expressed as follows: where ; cos h 2 (κy) (32) which is identical to the result of Zhao and Yang [18] on the EOF of power-law fluid with a symmetrical zeta potential. The detailed derivations are described in Appendix B.

Results and Discussions
The key parameters that affect velocity distribution are the fluid behavior index n, electrokinetic parameter κH, and dimensionless zeta potential difference R = (ψ t − ψ b )/(ψ t + ψ b ) between ψ t and ψ b . In this section, the effects of these parameters on velocity distribution are investigated. Figure 2 shows the dimensionless velocity (u/u s ) distributions from Equation (21) for different values of fluid behavior index n at a fixed κH of 15. Figure 2a represents the velocity distributions with same zeta potentials (R = 0) at the bottom and top, while Figure 2b indicates those of asymmetric zeta potentials with R of 0.2 (ψ t /ψ b = 1.5). In both cases of symmetric and asymmetric zeta potentials, as the fluid behavior index n decreases, the velocity gradient near the wall increases, and the plug-like characteristics of velocity distribution are enhanced. This is because the fluid with smaller fluid behavior index is less viscous, and the velocity can easily change from zero at the wall to the Smoluchowski velocity at the core region. which is identical to the result of Zhao and Yang [18] on the EOF of power-law fluid with a symmetrical zeta potential. The detailed derivations are described in Appendix B.

Results and Discussions
The key parameters that affect velocity distribution are the fluid behavior index n, electrokinetic parameter κH, and dimensionless zeta potential difference R = (ψt-ψb)/(ψt+ψb) between ψt and ψb. In this section, the effects of these parameters on velocity distribution are investigated. Figure 2 shows the dimensionless velocity (u/us) distributions from Equation (21) for different values of fluid behavior index n at a fixed κH of 15. Figure 2a represents the velocity distributions with same zeta potentials (R = 0) at the bottom and top, while Figure 2b indicates those of asymmetric zeta potentials with R of 0.2 (ψt/ψb = 1.5). In both cases of symmetric and asymmetric zeta potentials, as the fluid behavior index n decreases, the velocity gradient near the wall increases, and the plug-like characteristics of velocity distribution are enhanced. This is because the fluid with smaller fluid behavior index is less viscous, and the velocity can easily change from zero at the wall to the Smoluchowski velocity at the core region.    Figure  3a shows the velocity distributions of the shear thinning fluid (n = 0.8) and Figure 3b shows those of   Figure 3a shows the velocity distributions of the shear thinning fluid (n = 0.8) and Figure 3b shows those of the shear thickening fluid (n = 1.2). In both cases, as κH increases, the velocity distribution changes from parabolic type to plug-like type. The increase of κH means a decrease in Debye length. In other words, the EDL thickness, on which the electrostatic body force is applied, decreases and the velocity distribution changes to a plug-like type.
Micromachines 2018, 9, x FOR PEER REVIEW 7 of 11 the shear thickening fluid (n = 1.2). In both cases, as κH increases, the velocity distribution changes from parabolic type to plug-like type. The increase of κH means a decrease in Debye length. In other words, the EDL thickness, on which the electrostatic body force is applied, decreases and the velocity distribution changes to a plug-like type.
(a)  The velocity distributions near the top and bottom walls develop from zero (on the wall) to close to the generalized Smoluchowski velocity determined by the corresponding zeta potentials; in the core region, these two velocity distributions near the walls are almost linearly connected. Therefore, as the difference in zeta potential between the top and bottom walls increases, the velocity gradient in the core region increases. The velocity gradient in the core region decreases and increases the viscosity of the shear thinning fluid and shear thickening fluid, respectively. As a result, as the dimensionless zeta potential difference R increases, the velocity at the center (y/H = 0) increases for shear thinning fluids and decreases for shear thickening fluids.  The velocity distributions near the top and bottom walls develop from zero (on the wall) to close to the generalized Smoluchowski velocity determined by the corresponding zeta potentials; in the core region, these two velocity distributions near the walls are almost linearly connected. Therefore, as the difference in zeta potential between the top and bottom walls increases, the velocity gradient in the core region increases. The velocity gradient in the core region decreases and increases the viscosity of the shear thinning fluid and shear thickening fluid, respectively. As a result, as the dimensionless zeta potential difference R increases, the velocity at the center (y/H = 0) increases for shear thinning fluids and decreases for shear thickening fluids.

Conclusions
In this study, the exact solutions are proposed for fully developed two-dimensional steady unidirectional EOFs of power-law fluids with different zeta potentials at the top and bottom walls. The exact solutions are expressed in terms of Appell's first hypergeometric functions. The effects of parameters such as the fluid behavior index n, electrokinetic parameter κH, and zeta potential ψt and ψb on the velocity distribution are investigated.   (A2) where ψ(y), I(y) and w are defined in Equations (10), (20), and (25), respectively. The term |ψ(y)| ψ(y) means the sign of the electrical potential. Since the wall zeta potentials are assumed to be negative (ψ t , ψ b < 0), ψ(y) always has a negative value; therefore, the term |ψ(y)| ψ(y) is −1. With all the aforementioned considerations, the velocity distribution is obtained as Equation (21).
At point y = y c , the two equations according to the interval in Equation (21) should have the same value. Since I(y c ) + C = 0 by Equation (18), V + (y c ) = V − (y c ) = 0. Therefore, the integral constant C can be obtained from the following equation which is simplified as Equation (27).