# Analysis of the electric double layer effect on the microchannel flow stability

Sedat F. Tardu

► To cite this version:

Sedat F. Tardu. Analysis of the electric double layer effect on the microchannel flow stability. Microscale Thermophysical Engineering, 2004, 8 (4), pp.383-401. 10.1080/10893950490516938. hal-00262000

HAL Id: hal-00262000

https://hal.science/hal-00262000v1

Submitted on 18 Feb 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# ANALYSIS OF THE ELECTRIC DOUBLE LAYER EFFECT ON MICROCHANNEL FLOW STABILITY

**Sedat Tardu**

Laboratoire des Ecoulements: Géophysiques et Industriels, Grenoble Cedex,
LEGI BP 53 X France

It is shown that the microchannel flow under the electric double layer (EDL) effect is inviscidly unstable. A classical Orr–Sommerfeld analysis revealed that the critical Reynolds number of the primary (linear) instability may decrease by a decade, provided that the liquid contains a very small amount of ions and is associated with large enough zeta potential and low conductivity/ viscosity. When the EDL layer is thick, the inflexion point moves away from the viscous near wall region, and the inviscid instability mechanism becomes more dominant. The EDL modes are slower than the Poiseuille ones and the spectrum of the eigenvalues shows its strong destabilizing effect. There are some strong arguments suggesting that the nonlinear stability mechanism under the EDL effect is also more severe compared with the macroscale flow. The present investigation suggests that early transition in microchannel flows is plausible and can be checked through well-controlled experiments.

**Keywords** microfluidics, microchannels, EDL, linear stability, early transition

## INTRODUCTION

The “macro effects” should not be confused with real “micro effects” when explaining the deviations from the classical Navier–Stokes approach observed in microchannel flows. The roughness, for instance, is not a micro effect, since on the same scale it produces similar behavior in macro flows. Microscale effects can be neglected in conventional situations, but become increasingly important as the characteristic length of the flow decreases. Molecular, interfacial, and surface effects enter into this category.

Contrary to gas flows, the molecular effects in micro liquid flows are difficult to model. The slip length and the temperature jump can be related to the Knudsen number through the momentum and temperature accommodation coefficients and first- or second-order models in the rarefied regime of gas flows. Discontinuous boundary conditions in liquid flows and significant molecular effects hold only in extremely small devices having extremely high speeds [1]. Tardu [2] used the universal relation of slip length as a function of the shear rate given by Thompson and Troian [3] to estimate the critical width of microchannels to have slip liquid flows. He concluded that the hydraulic diameters of such devices have to be smaller than 0.52 µm. Even then, the slip velocity is only 0.04% of the centerline velocity at the transitional Reynolds number.
## NOMENCLATURE

<table><tr><td>a</td><td>channel half height</td></tr><tr><td>c</td><td>complex wave speed of the disturbance</td></tr><tr><td>-dP/dX</td><td>imposed streamwise pressure gradient</td></tr><tr><td>e</td><td>electron charge, 1.6021 × 10<sup>-19</sup> C</td></tr><tr><td>E<sub>x</sub></td><td>electrical field strength</td></tr><tr><td>E<sub>s</sub></td><td>streaming potential</td></tr><tr><td>EDL</td><td>electric double layer</td></tr></table>

$$G = \frac{(n_0 z \epsilon a)^2}{\lambda_0 \mu}$$
parameter appearing in the EDL velocity component
$$I_c$$
conduction current
$$I_s$$
streaming current
$$k_b$$
Boltzmann constant, 1.3805 × 10<sup>-23</sup> J/(mol K)

<table><tr><td>k</td><td>(2n<sub>0</sub>z<sup>2</sup>e<sup>2</sup>/εε<sub>0</sub>k<sub>b</sub>T)<sup>1/2</sup></td><td>inverse EDL thickness</td></tr><tr><td>l</td><td>length of the channel</td></tr><tr><td>L</td><td>Orr-Sommerfield operator</td></tr><tr><td>n<sub>0</sub></td><td>bulk ionic concentration</td></tr><tr><td>p</td><td>pressure disturbance</td></tr><tr><td>R</td><td>ratio of the EDL to the Poiseuille velocity components at the centerline of the channel</td></tr><tr><td>Re</td><td>Reynolds number based on the centerline velocity and half channel height</td></tr><tr><td>Re<sub>c</sub></td><td>critical Reynolds number</td></tr><tr><td>Re<sub>t</sub></td><td>transitional Reynolds number</td></tr><tr><td>T</td><td>absolute temperature</td></tr></table>

<table><tr><td>U, u</td><td>dimensional and nondimensional streamwise velocity component</td></tr><tr><td>u<sub>EDL</sub></td><td>nondimensional velocity induced by EDL</td></tr><tr><td>v</td><td>wall-normal velocity disturbance</td></tr><tr><td>X, x</td><td>dimensional and nondimensional streamwise coordinate</td></tr><tr><td>Y, y</td><td>dimensional and nondimensional wall normal coordinate</td></tr><tr><td>z</td><td>valence of positive or negative ions</td></tr></table>

### Greek symbols

<table><tr><td>α</td><td>wave number</td></tr><tr><td>α<sub>c</sub></td><td>critical wave number</td></tr><tr><td>ε</td><td>dielectric constant of the medium</td></tr><tr><td>ε<sub>0</sub></td><td>permittivity of vacuum, 8.854 × 10<sup>-12</sup> CV<sup>-1</sup> m<sup>-1</sup></td></tr><tr><td>ζ, ζ̄</td><td>dimensional and nondimensional zeta potential</td></tr><tr><td>κ = ka</td><td>Debye parameter</td></tr><tr><td>λ<sub>0</sub></td><td>electrical conductivity of the fluid</td></tr><tr><td>λ<sub>i</sub></td><td>eigenvalues</td></tr><tr><td>μ</td><td>dynamic viscosity of the fluid</td></tr><tr><td>ρ</td><td>charge density</td></tr><tr><td>ψ</td><td>electrostatic potential at any point in the electric double layer</td></tr><tr><td>Ψ</td><td>stream function of the wave-like disturbance used in the stability analysis</td></tr></table>

The interfacial effects are presumably most responsible for the deviations observed from the classical macro theory. One of the micro effects that may play an important role is the electric double layer (EDL). Most solid surfaces have an electric surface potential. The most common mechanism for the charging of surface layers in microfluidics is the depronotonation of surface groups including silica, glass, acrylic, and polyester [4]. The electrostatic charges present on the solid surface attract the counterions to establish an electrical field. The ions are immobile in the compact layer next to the wall and less than 1 nm thick. They can, however, be transported by an imposed flow field in the diffuse EDL layer. Consequently, the counterion concentration near the wall is larger than in the bulk of the fluid. The imposed pressure gradient accumulates the mobile ions downstream and sets up an electrical field whose potential is defined as the streaming or electrokinetic potential. The streaming potential and the net charge density induces a streamwise external force. In the steady state, the streaming current is in equilibrium with the conduction current in the opposite direction. This allows the determination of the streaming potential and the velocity profiles under the EDL effect. In macroscale flows, these effects are negligible because the thickness of the EDL is minor compared to the height of the channel. However, in micro flows the EDL plays a rather significant role.
Well-controlled recent experiments have clearly confirmed that the EDL can explain the behavior of the Poiseuille number in laminar regime, providing that the liquid contains a very small amount of ions. Kulinsky et al. [5], for example, reported an increase of 70% of the friction factor under the EDL effect in planar channels of 4 µm heights.

Early transition to turbulence in microchannel flows has been reported by several investigations in the past. Wang and Peng [6] indicated that the transition occurred at $Re < 300$ in microchannels 0.2–0.8 mm wide and 0.7 mm deep, which is significantly smaller than 750 corresponding to macroscale flows. (The Reynolds number throughout this article is based on the centerline velocity and half height of the channel.) A transitional number of about 500 has also been reported in the experiments conducted in 1-mm-wide trapezoidal microchannels with channel depths ranging between 79 and 325 µm [7–9]. The micro tube experiments of Mala and Li [10] indicate that there is an early transition from laminar to turbulent flows for $Re > 112–338$. Pointing at a plausible EDL effect, they have clearly stated that “There is an early transition from laminar to turbulent flow, and the range of $Re_t$ values varies somewhat, depending on the diameter and the material of the wall.” These effects may also be attributed to the roughness that may both influence the transition in the entry region, or directly affect the flow behavior through an implied roughness viscosity [11]. It may, however, also be presumed that direct interfacial effects play some role in early transition.

The main goal of this study is to directly investigate the EDL effect on the linear stability of planar channel flow and estimate indirectly the resulting transitional Reynolds numbers. It is asked here whether electrokinetic effects on liquid flows may cause the early transition or not and whether the linear stability EDL flow characteristics may explain the small transitional Reynolds numbers reported in some recent microchannel experiments.

## EQUATIONS GOVERNING 2-D CHANNEL FLOW UNDER THE EDL EFFECT

The ion distribution in a symmetric plane channel flow is given by the equilibrium Boltzmann equation:

$$n^{\pm} = n_0 \exp \left( \pm \frac{ze\psi}{k_b T} \right) \qquad (1)$$

where $n_0$ is the bulk ionic concentration, $\psi$ is the electrical potential, $z$ is the valence of positive or negative ions, $e$ is the electron charge, $k_b$ is the Boltzmann constant, and $T$ is the absolute temperature. The effects of finite ion size, the gradients of the dielectric strength and viscosity are neglected and the fluid properties (viscosity and permittivity) are constant.

The net charge density in a unit volume of fluid is:

$$\rho = (n^{+} - n^{-})ze = -2n_{0} \sinh \left( \frac{ze\psi}{k_{b}T} \right) \qquad (2)$$

The Poisson equation relating the electrical potential and the net charge density per unit volume reads as

$$\nabla^2 \psi = - \frac{\rho}{\varepsilon_0 \varepsilon} \qquad (3)$$
with $\epsilon$ and $\epsilon_0$ standing, respectively, for the dielectric constant of the medium and the permittivity of vacuum. Equation (3) combined with Eq. (2) may be written in nondimensional form as

$$ \nabla^2 \bar{\psi} = -\kappa^2 \bar{\rho}(y) = \kappa^2 \sinh(\bar{\psi}) \quad (4) $$

The scaling reference potential and charge density are, respectively, $E_S = \frac{k_b T}{ze}$ and $\rho_S = n_0 ze$. Equation (4) makes the nondimensional Debye parameter $\kappa = ka$ appear with:

$$ k = (2n_0 z^2 e^2 / \epsilon \epsilon_0 k_b T)^{1/2} \quad (5) $$

whose inverse is the EDL thickness, and $a$ is the channel half height.

Equation (4) may be linearized when the electrical potential is small compared with the thermal energy of ions, i.e., when $|ze\psi| \ll |k_b T|$, or $\bar{\psi} \ll 1$ to give:

$$ \nabla^2 \bar{\psi} = \kappa^2 \bar{\psi} \quad (6) $$

Equation (6) holds remarkably well for $\bar{\psi} < 2$. One has:

$$ \frac{d^2 \bar{\psi}}{dy^2} = \kappa^2 \bar{\psi} \quad (7) $$

for a two-dimensional channel, where $y = Y/a$ stands for the nondimensional wall normal distance.

The resulting charge density is $\bar{\rho}(y) = -\bar{\psi} = -\frac{1}{\kappa^2} \frac{d^2 \bar{\psi}}{dy^2}$. The boundary conditions for Eq. (7) are:

* At the wall $\psi = \bar{\psi}_0$. The electrical potential at the wall can be approximated by the zeta potential $\bar{\psi}_0 = \bar{\zeta}$. The general solution of the resulting EDL-velocity profile is not affected by this approximation. However, it is generally difficult to measure $\bar{\psi}_0$. It is easy to determine experimentally the zeta potential at the boundary between the compact and diffuse layers by determining the streaming potential and the pressure drop [12].

* At the centerline $\bar{\psi} = 0$, providing that the EDL diffuse layers developing on both sides of the channel do not overlap. The solution of Eq. (7) is:

$$ \bar{\psi} = \frac{\bar{\zeta}}{\sinh(\kappa y)} |\sinh(\kappa y)| \quad (8) $$

The charge density $\bar{\rho}(y)$ and the streaming electrical field strength caused by the ions mobilized under the imposed pressure gradient acts as an external force in the streamwise momentum equation of a fully developed 2-D channel flow:

$$ 0 = \mu \frac{d^2 U}{dY^2} - \frac{dP}{dX} + E_x \rho(Y) \quad (9) $$
The electrical field strength is related to the streaming potential by $E_x = \frac{E_s}{l}$, where $l$ is the length of the channel. The nondimensional form of Eq. (9) is:

$$0 = \frac{d^2 u}{dY^2} + 2 - \frac{2\Gamma \overline{E_s}}{\kappa^2} \frac{d^2 \overline{\psi}}{dY^2} \quad (10)$$

where $\Gamma = \frac{\zeta n_0 z e a^2}{l \mu U_s}$ and the scaling velocity is the centerline velocity of the Poiseuille component, i.e., $U_s = -\frac{a^2 dp/dx}{2\mu}$. Integrating Eq. (10) with the appropriate boundary conditions gives the velocity profile:

$$u = 1 - y^2 - 4 \frac{\Gamma \overline{E_s} \overline{\zeta}}{\kappa^2} \left\{ 1 - \left| \frac{\sinh \kappa y}{\sinh \kappa} \right| \right\} \quad (11)$$

The last equation is the nondimensional form of the velocity profile already obtained by Mala et al. [13]. The streaming current due to the flux of the net density charge is given by the following cross integral by definition:

$$I_s = \iint U \rho (Y) dY dZ \quad (12)$$

The streaming current produces the conduction current in the reverse direction:

$$I_c = E_x \lambda_0 S = \frac{E_s}{l} \lambda_0 S \quad (13)$$

where $\lambda_0$ is the electrical conductivity of the fluid, which is proportional to the ionic concentration. In the steady state there is no net current $I_s + I_c = 0$. This allows the closure of the equations and the determination of the streaming potential. The final version of the nondimensional velocity profile is:

$$u = 1 - y^2 - 4 \frac{I_1 - I_2}{\frac{\kappa^2 \sinh \kappa}{\overline{\zeta}^2 G} + 4 \left( I_3 - \frac{I_4}{\sinh \kappa} \right)} \left\{ 1 - \left| \frac{\sinh \kappa y}{\sinh \kappa} \right| \right\} \quad (14)$$

The parameter $G$ is given by $G = \frac{(n_0 z e a)^2}{\lambda_0 \mu}$. The quantities $I$ in Eq. (1) are:

$$I_1 = I_3 = \frac{\cosh \kappa - 1}{\kappa} \\ I_2 = \left( \frac{1}{\kappa} + \frac{2}{\kappa^3} \right) \cosh \kappa - \frac{2}{\kappa^2} \sinh \kappa - \frac{2}{\kappa^3} \quad (15)^1 \\ I_4 = \frac{\sinh \kappa \cosh \kappa}{2\kappa} - \frac{1}{2}$$

<footnote n="1">¹There is a small printing error in Mala et al. [13], p. 3082</footnote>
It is seen that the velocity profile can be decomposed into a macro Poiseuille component plus an EDL effect component. A closer inspection of Eq. (14) allows expressing the latter in another closed form as:

$$u_{\text{EDL}} = -4 \frac{I_1 - I_2}{\frac{\lambda_0 \mu \sinh \kappa}{\bar{\zeta}^2} + 4 \left( I_3 - \frac{I_4}{\sinh \kappa} \right)} \left\{ 1 - \left| \frac{\sinh \kappa y}{\sinh \kappa} \right| \right\} \quad (16)$$

The parameter $\bar{\lambda}_0\mu$ is the nondimensional electrical conductivity times the dynamic viscosity of the fluid. It can be seen that the scaling parameter for $\lambda_0\mu$ is the square of the scaling charge density times the scaling length, i.e., $(\rho_S/k)^2$ leading to:

$$\bar{\lambda}_0\mu = \frac{\lambda_0\mu}{(n_0\epsilon\epsilon_0k_bT/2)} \quad (17)$$

There are three parameters governing the EDL effect, namely

$$u_{\text{EDL}} = u_{\text{EDL}}(\kappa, \bar{\zeta}, \bar{\lambda}_0\mu) \quad (18)$$

Due to the grouping of the parameters this relation is equivalent to:

$$u_{\text{EDL}} = u_{\text{EDL}}(\kappa, G\bar{\zeta}^2) \equiv u_{\text{EDL}}\left(\kappa, \frac{\bar{\lambda}_0\mu}{\bar{\zeta}^2}\right) \quad (19)$$

as clearly suggested by Eqs. (14) and (16).

Gradients of the dielectric strength, viscosity, and conductivity should be incorporated in more realistic models. Yet these effects are presumably negligible in significantly dilute solutions, such as ultra-deionized, ultra-filtered water and pure organic liquids.

The main EDL effects may be summarized as follows:

* An increase of the friction constant and apparent viscosity.

* A decrease of the Nusselt number.

These effects are persistent, yet not significantly important, at least for large values of $\kappa$. The apparent viscosity increases, for example, by a factor of nearly 3 at $\kappa = 2$, but the EDL effect on the wall shear stress disappears quickly when $\kappa \ge 10$ [13]. The Nusselt number decreases by nearly 40% at $\kappa = 5$, and less than 5% at $\kappa = 50$ [2].

## INFLEXIONAL INVISCID INSTABILITY

One important aspect of the EDL velocity profile is the existence of an inflexional, or, in other words, maximum vorticity point. One has indeed:

$$\frac{\partial^2 u}{\partial y^2} = -2 + 4 \frac{\kappa^2 (I_1 - I_2)}{\sinh \kappa \left[ \frac{\bar{\lambda}_0 \mu \sinh \kappa}{\bar{\zeta}^2} + 4 \left( I_3 - \frac{I_4}{\sinh \kappa} \right) \right]} \sinh \kappa y \quad (20)$$
which can also be written as:

$$u'' = \frac{\partial^2 u}{\partial y^2} = -2 - \frac{R\kappa^2}{\sinh \kappa} \sinh \kappa y \quad (21)$$

where $R \le 0$ is the ratio of the EDL to the Poiseuille velocity components at the centerline of the channel, i.e., $R = \frac{u_{EDL}}{u_{Pois}}$ at $y = 0$. The maximum vorticity takes place at the inflexion point:

$$y_I + \pm \frac{1}{\kappa} \arcsin h \left( \frac{-2 \sinh \kappa}{R\kappa^2} \right) \approx \pm \left[ 1 + \frac{1}{\kappa} \ln \left( -\frac{2}{R\kappa^2} \right) \right] \quad (22)$$

The right-hand side of the last relationship is a good approximation to $y_I$ for $\kappa \ge 10$. Clearly, the EDL flow is inviscidly unstable. The inviscid instability is of Fjortoft type because after the inflexion point near the wall, $u'' > 0$ and the velocity is smaller than the velocity at the inflexion point $u_I$, resulting in $u''(u - u_I) < 0$ in this zone.

The present investigation is a parametrical one based on nondimensional quantities. Yet we will introduce a numerical case to fix the ideas. Thus, we choose a “reference EDL flow” with $\kappa = 41$, $G = 12,720$, and $\bar{\zeta} = 2.1254$. This case corresponds to the flow of an infinitely diluted KCl aqueous solution ($n_0 = 3.764 \times 10^{19} \text{ m}^{-3}$) through a microchannel of height 100 µm subject to a zeta potential of 50 mV. The Debye length is 1.2 µm, and the bulk conductivity is $\lambda_0 = 8 \times 10^{-9} \text{ S/m}$. The parameters $G$ and $\bar{\zeta}$ and are identical to those of Mala et al. [10] wherein detailed electro-viscous effects can be found. More severe cases, leading to large EDL effects have been treated in the literature [14].

Figure 1 shows the distribution of the ratio of the EDL to Poiseuille velocity components at the centerline, together with the location of the inflexion point, versus $\kappa$. It appears that $R$ increases when $\kappa$ decreases, showing that the decelerating EDL effect becomes important with either a small separation distance between the plates or a large

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Distribution of $R$, the ratio of the EDL to Poiseuille velocity components (open circles and ordinate at left) and of $y_I$, the inflexion point position (bold circles and ordinate at right) vs. the nondimensional Debye-Hückel parameter $\kappa$. The zeta potential is $\bar{\zeta} = 2.1254$. The parameters $\kappa$ and $G$ change as the half width $a$ of the channel is varied (see the text). One has $G = 12,720$ at $\kappa = 41$ (reference case).</figcaption></figure>
EDL thickness. At the same time the inflexion point moves away from the wall toward a region wherein viscous effects are less significant. Thus, it may be expected that the microflow becomes increasingly susceptible to inviscid instability with decreasing $κ$. These arguments are logical, but a definite picture of the EDL linear stability mechanism can only be given by a detailed Orr–Sommerfeld analysis.

## STABILITY ANALYSIS

The linear hydrodynamic stability under the EDL effect is studied through classical methods. The Orr–Sommerfeld equation is solved by a Galerkin-like procedure [15]. The normal mode solutions of the disturbance equation are:

$$ (u, v, p) = R(\hat{v}, \hat{v}, \hat{p}) \exp(i\alpha x) \quad (23) $$

where R is the real part, $α$ is the dimensionless wave number of the disturbance, and $x$ is the streamwise coordinate. Introducing the stream function:

$$ \Psi(x, y, t) = \phi(y, t) \exp(i\alpha x) \quad (24) $$

the Orr–Sommerfeld equation takes the form:

$$ \frac{\partial}{\partial t} L \phi = \frac{1}{\text{Re}} L^2 \phi - i\alpha \left( u L \phi - \frac{d^2 u}{dy^2} \phi \right) \quad (25) $$

with the boundary conditions $\phi = \frac{\partial \phi}{\partial y} = 0$ at $y = \pm 1$. The operator L is $L = \frac{\partial^2}{\partial y^2} - \alpha^2$. The stream function is expanded in a Chebyshev polynomial series:

$$ \phi(y, t) = \sum_{n=1}^{N} a_n(t) T_{2n-2}(y) \quad (26) $$

where $T_m(y) = \cos(m \cos^{-1} y)$ denotes the Chebyshev polynomials of the first kind. We took $N = 256$ through this study. Equation (21) takes the form:

$$ \mathbf{Q} \frac{d\mathbf{a}}{dt} = (\mathbf{P} - i\alpha\mathbf{J}) \cdot \mathbf{a} \quad (27) $$

and the matrices are determined by using the $τ$-method described by Orszag [16]. The last equation can be written as:

$$ \frac{d\mathbf{b}}{dt} = \mathbf{D} \cdot \mathbf{b} \quad (28) $$

by introducing $\mathbf{B}$ the matrix that diagonalizes $\mathbf{Q}^{-1} \cdot (\mathbf{P} - i\alpha\mathbf{J})$, $\mathbf{b} = \mathbf{B}^{-1} \cdot \mathbf{a}$, and the diagonal matrix $\mathbf{D} = \lambda_i \delta_{ij} = \mathbf{B}^{-1} \cdot \mathbf{Q}^{-1} \cdot (\mathbf{P} - i\alpha\mathbf{J}) \cdot \mathbf{B}$. The eigenvalues are denoted by $\lambda_i$ and the flow is unstable when. The disturbances with symmetric streamfunctions are considered only. The method had similar results when compared to Grosch and Salwen [17], who used different sets of expansion functions.
Effect of the Debye–Hückel Parameter

The EDL flow depends on the groups of the parameters $κ$ and $G\bar{\zeta}^2$ (Eq. (19)). When the effect of the nondimensional Debye–Hückel length is first analyzed, the pa- rameter $G\bar{\zeta}^2$ is kept constant for a given $κ$. The sensitivity of the hydrodynamic stability characteristics to the zeta potential and the parameter $G$ will be analyzed in the next sec- tion. Recall that $G\bar{\zeta}^2 = 57,460$ and $κ = 41$ in the reference case introduced previously. For subsequent computations, the zeta potential, the ionic concentration, and the physical characteristics of the fluid (the permittivity, conductivity and viscosity) are kept constant and the half channel height $a$ is varied alone. Therefore, $κ$ and $G$ decrease proportionally with $a$ and $a^2$, respectively. Thus, $G\bar{\zeta}^2$ is 8740 and 2185 at, for example, $κ = 16$ and $κ = 8$, respectively.

The neutral curves deduced from the hydrodynamic stability analysis are summa- rized in Figure 2. The friction factor increases by only 10% at $κ = 41$, but the critical Reynolds number decreases by a factor nearly equal to 2: the critical wave and Reynolds numbers of the micro flow are, respectively, $α_c = 1.10$ and $Re_c = 3190$, to be compared with $α_c = 1.02$ and $Re_c = 5772$ of the conventional Poiseuille flow. The destabiliza- tion under EDL disappears quickly as $κ$ increases. Thus, at $κ = 164$ (a minichannel of 400 µm height) the critical Reynolds number is $Re_c = 5730$ (triangle in Fig. 2), which only slightly differs from $Re_c = 5772$. The micro EDL effects for the micro/minichannels of height larger than typically 100 µm are not significant. Note, however, that, in a purely

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Neutral curves of the EDL flow compared with the Poiseuille flow. The open circles correspond to Poiseuille flow with $κ = ∞$ and/or $\bar{\zeta} = 0$. The results of Grosh and Salwen [17] are shown by small bold squares. Bold circles correspond to $κ = 41$, $G = 12,720$, and $\bar{\zeta} = 2.1254$. The rest of the results are obtained by changing the microchannel height and keeping constant the rest of the parameters. The triangle is obtained for $κ = 164$ and $G = 203,520$. The neutral curves for $κ = 8$, $G = 484$ (bold squares), $κ = 16$, $G = 1937$ (open squares), and $κ = 41$, $G = 12,720$ (bold circles) are compared with the macroscale flow (open circles).</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Distribution of the critical wave (open symbols, ordinate at left) and Reynolds numbers (closed symbols, ordinate at right) deduced from Figure 2 vs. $\kappa$.</figcaption></figure>

mathematical sense, the macro Poiseuille flow is reached only when $\kappa \to \infty$ and/or the zeta potential is $\bar{\zeta} \to 0$.

There is a considerable destabilization of the flow at smaller Debye–Hückel parameters such as $\kappa = 16$, and $\kappa = 8$. In order to fix the ideas, these cases would correspond to channels of heights 40 µm and 20 µm, respectively, subject to the reference ionic concentration, zeta potential viscosity, and conductivity. The critical Reynolds number decreases to $Re_c = 817.5$ for $\kappa = 16$. The decrease is more than a decade at $\kappa = 8$ with $Re_c$ as small as $Re_c = 328.75$. The critical wave number increases at the same time and reaches $\alpha_c = 1.7$ at $\kappa = 8$.

The distribution of $Re_c$ and $\alpha_c$ as a function of $\kappa$ is shown in Figure 3. The critical Reynolds number changes linearly versus $\log(\kappa)$ in a large range of the Debye–Hückel parameter before reaching a plateau region near $\kappa = 4$ ($Re_c = 279$, $\alpha_c = 1.9$). The critical wave number, on the other hand, increases linearly with $\log(\kappa)$ when $\kappa < 41$. The analysis is deliberately stopped at $\kappa = 4$, since there may be some adverse effects due to the overlap of the EDL diffuse layers developing on both sides of the channel for $\kappa < 4$.

There are two main features of the neutral curves shown in Figure 2. First, the band of unstable wave numbers of the EDL micro flow is significantly larger compared with the Poiseuille macro flow due to the inviscid inflexional instability (unstable for $Re \to \infty$ for a given $\alpha$). The smaller the parameter $\kappa$, the broader the neutral curve is. Second, the upper branch of the neutral curves diverges slowly as the Reynolds number increases, while the lower branch of the micro and macro flows have similar behaviors. This can be seen with a careful look at Figure 2 for $\kappa = 16$ and $\kappa = 8$. The increase of the critical wave number corresponding to the upper branch of the neutral curves is seen more clearly in Figure 4. The neutral curve in the macro Poiseuille flow is closed for large Reynolds numbers, whereby the upper branch joins the lower one when $Re$ is large enough. This is simply because the Poiseuille flow is inviscidly stable. Under the presence
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Increase of the wave number along the upper branch of the neutral curves for three different values of $\kappa$ deduced from Figure 2.</figcaption></figure>

of the EDL, the upper and lower branches diverge so that the neutral curves are open. This means that the inviscid instability overcomes the shear layer instability at $Re \to \infty$ as expected. When $\kappa$ is low, the inflexion point leaves the viscosity-dominated near wall region, so that the inviscid instability mechanism can become more operative. The so-called divergence is slow. The wave number begins to increase at $Re = 200 \times 10^3$ for $\kappa = 41$, and earlier when the EDL layer is thicker, i.e., at approximately $Re = 10 \times 10^3$ for $\kappa = 16$, and $\kappa = 8$.

## Spectrum of Eigenvalues

Figure 5 shows the spectrum of eigenvalues as a function of $Re$ for a single fixed wave number $\alpha = 1$ and $\kappa = 16.4$. Only the first three eigenvalues corresponding to the complex wave speed $c = -\lambda/i\alpha$ are shown for the sake of clarity. The modes are ordered in such a way that ascending values correspond to descending values of $Im\{c\}$. Remember that the flow is stable (unstable) for $R\{\lambda\} < 0$ ($> 0$), or equivalently when $Im\{c\} < 0$ ($> 0$). The $Re$ was varied from 300 to $20 \times 10^3$. The curves corresponding to Poiseuille and EDL flows are denoted respectively by $P - i$ and $E - i$, where $i$ is the mode number.

Figure 5b shows that only the first mode has a different $Im\{c\}$ under the EDL effect. The imaginary parts of higher modes collapse entirely with those of the macroscale flow. Note that the imaginary part of the first mode in micro flow denoted by $Im\{c_{E-1}\}$ is larger than that of the macroscale Poiseuille flow $Im\{c_{P-1}\}$, showing the strong destabilizing effect of the EDL.

The EDL modes are slower than the Poiseuille ones, in the sense that $R\{c_{E-i}\}$ is smaller than $R\{c_{P-i}\}$ for all modes (Fig. 5a). At moderate Reynolds number, they may
<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Real (a) and imaginary parts (b) of the first three eigenvalues vs. Re, $\alpha = 1$, and $\kappa = 16$. The eigenvalue $i$ in EDL and Poiseuille flow is denoted by $E - i$ and $P - i$, respectively.</figcaption></figure>

be classified as slow modes since the real part of the complex wave speed is smaller than the average velocity of the channel, i.e., $R\{c_{E-i}\} < \frac{2}{3}$ for roughly $Re < 5 \times 10^3$. Modes 2 and 3 of macroscale flow are clearly fast modes, since their phase speed is close to the centerline velocity.

The eigenvalue spectrum is further analyzed for a more destabilizing case at $\kappa = 8$. Figure 6 shows the resulting distribution of the complex wave speed real and imaginary
<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Real (a) and imaginary parts (b) of the first three eigenvalues vs. Re, $\alpha = 1$, and $\kappa = 8$. The eigenvalue $i$ in EDL and Poiseuille flow is denoted by $E - i$ and $P - i$, respectively.</figcaption></figure>
parts. It is seen in Figure 6a that all EDL modes are now slow modes in the range of Reynolds numbers investigated here. The differences between $R\{c_{E-i}\}$ and $R\{c_{P-i}\}$ are larger compared to the case $\kappa = 16$. The parameter $\kappa$ does not appreciably affect the imaginary part of the spectrum as can be seen in Figure 6b, except a slight increase with decreasing $\kappa$ in the first mode.

Stability Sensitivity to the Zeta Potential and the Parameter G

We have previously shown that the EDL velocity component depends on the group of parameter $G\bar{\zeta}^2$ for a given Debye–Hückel parameter $\kappa$ (Eq. (14)). The parameter $G$ can be expressed as a function of the nondimensional conductivity and viscosity product $G = \frac{\kappa^2}{\overline{\lambda_0\mu}}$, showing that $G$ varies like the inverse of $\overline{\lambda_0\mu} = \frac{\lambda_0\mu}{(n_0\varepsilon\varepsilon_0k_bT/2)}$ for fixed $\kappa$. Thus, $G$ increases with decreasing conductivity and/or viscosity.

To each couple $(\kappa, G\bar{\zeta}^2)$ corresponds a different set of the critical wave and Reynolds numbers $(\alpha_c, Re_c)$. It is difficult to illustrate the interdependency of $Re_c$ via $G\bar{\zeta}^2$, since $\alpha_c$ changes continuously as $G\bar{\zeta}^2$ varies. Therefore, we fixed the wave number at its macroscale flow critical value $\alpha_c = 1.02$ to give a general tendency and extract clear information.

Figure 7 shows the distribution of $Re_c$ versus $G\bar{\zeta}^2$ for four different values of $\kappa$. Consider first the results corresponding to the smallest value $\kappa = 8$. The effect of the inviscid instability is so important in this case that the critical Reynolds number changes rapidly even for insignificantly small values of the zeta potential. For instance, $Re_c$ varies by 25% with $Re_c = 4300$, already at $G\bar{\zeta}^2 = 10$. That would correspond to a zeta potential of only 0.66 mV for the reference case. Increasing the zeta potential by 10 from this value decreases $Re_c$ by nearly 3. There is a range of $G\bar{\zeta}^2$ in which the critical Reynolds number decreases sharply before reaching a plateau region wherein $Re_c$ is close to the critical number deduced from the neutral curves. This range varies

<figure><img src="image_10.png" /><figcaption>Figure 7. Distribution of the critical Reynolds number vs. $G\bar{\zeta}^2$ for $\alpha = 1.02$.</figcaption></figure>
<table><caption>**Table 1.** Sensitivity of the critical Reynolds number to $G\bar{\zeta}^2$ variations</caption><thead><tr><th>κ</th><th>Range of Gζ̄² within Re<sub>c</sub> decreases sharply<br>Lower value–upper value</th></tr></thead><tbody><tr><td>41</td><td>10<sup>3</sup>–10<sup>5</sup></td></tr><tr><td>20</td><td>10<sup>2</sup>–10<sup>4</sup></td></tr><tr><td>16.4</td><td>10<sup>2</sup>–10<sup>4</sup></td></tr><tr><td>8.2</td><td>10–10<sup>3</sup></td></tr></tbody></table>

<table><caption>**Table 2.** Experimental conditions related to the results presented in Figure 8</caption><thead><tr><th>Authors</th><th>Configuration</th><th>Hydraulic radius (μm)</th></tr></thead><tbody><tr><td>Mala and Li [10]</td><td>Microtubes</td><td>102.5–25</td></tr><tr><td>Qu et al. [11]</td><td>Trapezoidal silicon microchannels</td><td>71–23.5</td></tr><tr><td>Gao et al. [27]</td><td>2-D microchannels</td><td>100–53</td></tr></tbody></table>

with $κ$ (Table 1) and it is approximately two decades large. The lower and upper limits of the $G\bar{\zeta}^2$ range increase with $κ$. It is clear from both Figure 7 and Table 1 that $Re_c$ is more sensitive to $G\bar{\zeta}^2$ when the Debye–Hückel parameter is small enough, i.e., for $κ ≤ 40$. The first and main condition to expect significant EDL effects on hydrodynamic linear stability is to deal with liquids of low ionic concentration. The second condition requires large enough $G\bar{\zeta}^2$ value. This is fulfilled by low conductivity, which is generally associated with wall/liquids interactions that lead to high zeta potential $\bar{\zeta}$. The second condition is less critical than the first. To show this, consider a median value of $κ = 20$. Figure 7 shows that $Re_c$ decreases by a factor of 6 at $G\bar{\zeta}^2 = 10^3$, corresponding, for example, to $G = 250$ and $\bar{\zeta} = 2$; that is, a liquid that is 50 times more conductive than the reference case.

DISCUSSION AND CONCLUSION

The electro-viscous effects modify considerably the linear stability characteristics of microchannel flows. Two main conditions are required for a significant alteration of the critical Reynolds number:

1. The nondimensional Debye–Hückel parameter should be smaller than $κ ≤ 40$. There is no effect on stability for liquids with high ionic concentration and/or minichannels.

2. Large enough zeta potential and liquids with low conductivity/viscosity leading typically to $G\bar{\zeta}^2 ≥ 5000$ for $\frac{1}{k} ≈ 1$ μm are necessary. This condition is less critical. With a Debye length $\frac{1}{k} = 2$ μm, which is just two times larger than
the reference case, the critical Reynolds number decreases by nearly 80% in a
32 µm width channel at $G\bar{\xi}^2$ as small as $G\bar{\xi}^2 = 100$. The key question: Can
this effect be experimentally verified and validate early transition reported in
some recent investigations?

To give preliminary answers to these questions and discuss related phenomena, we
first have to note that the hydrodynamic linear stability and transition mechanisms are
different in plane Poiseuille flows. What we are interested in from an applications point
of view is indeed the transitional Reynolds number $Re_t$ related but not strictly equal
to $Re_c$.

According to Orszag and Patera [18], the instability process of wall bounded
macroscale shear flows involves three steps:

1. Primary (linear) instability of the basic shear flow

2. Nonlinear saturation of the primary instability and formation of a secondary flow

3. Secondary instability, which is the linear stability of the secondary flow

The first step gives $Re_c = 5772$ for macroscale plane Poiseuille flows, but transitional Reynolds numbers as low as $Re_t = 400$ have been reported in experiments. This is because the second step results in a Landau-type mechanism leading to subcritical instability or the metastability. The instability occurs with finite amplitude when all infinitesimal disturbances are stable at a significantly lower critical value compared to $Re_c$ [19]. The quasi-steady nonlinear analysis conducted by Orszag and Patera [18] suggests a critical Reynolds number for the secondary instability of $R \approx 400$. Three-dimensional disturbances decay for $R < 1000$; however, 1000 cannot be taken as a precise number because the three-dimensional instability requires that a threshold two-dimensional amplitude be achieved, which depends on the perturbation level. The experiments indicate that the ratio of the transitional Reynolds number to the critical Reynolds number of the primary instability is approximately $r = \frac{Re_t}{Re_c} \approx \frac{750}{5772} \approx \frac{1}{8}$ [20] but this is not an exact universal value either.

Some arguments suggest that the EDL can reinforce the subcritical nature of macro
Poiseuille flow through the nonlinear saturation of the primary instability (step 2). Stuart
[21] has shown that there are two main mechanisms generating the subcritical equilib-
rium state. The first mechanism prevents the subcritical instability through the distortion
of the mean motion. This effect is proportional to $Re_c \alpha_c^2$, which is significantly small
in EDL flow. The second mechanism reinforces the metastability by the wall normal
distortion of the fundamental and it is inversely proportional to $Re_c$. EDL directly en-
hances the subcritical state through the second mechanism and indirectly by retarding
the first mechanism. One of the simple explanations of the secondary instability (step 3)
is the creation of local inflexional streamwise profiles with strong shear susceptible to
inviscid instability. The basic flow in EDL is already inflexional. Interactions of the first
and second steps may lead to significantly smaller transitional Reynolds numbers in the
presence of EDL. The reasonable conclusion of these arguments is that the EDL flow is
subcritical and that presumably $r = \frac{Re_t}{Re_c} \le \frac{1}{8}$.

It is somewhat difficult to reach high Reynolds numbers in microchannels because
high pressure can cause channels to break. We show in Figure 8 the estimated transi-
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Estimated values of the transitional numbers under the EDL effect compared with some experiments and the range of Reynolds numbers that can be reached in experiments in microchannels. Experimental configurations are resumed in Table 2.</figcaption></figure>

transitional Reynolds as a function of $\kappa$, together with some microchannel and microtube data reported in the literature. Recall that Hagen-Poiseuille flow is linearly stable, and some aspects of the transition in pipe flows can be explained only by secondary instabilities. We supposed that $r = \frac{Re_t}{Re_c} = \frac{1}{8}$ as in macroscale flow. We need a Debye length to express data in terms of $Re - \kappa$ distribution. We took $\frac{1}{k} = 1~\mu\text{m}$ to be consistent. The open circles in Figure 8 correspond to the reference case. A second case, shown by bold circles, is also introduced. It consists of a liquid that is 11.5 times less conductive than the reference case leading to $G\bar{\zeta}^2 = 5000$ at $\kappa = 41$. It is clearly seen in Figure 8 that the experimental Reynolds numbers are larger than $Re_t$, concluding that experimental verification of the present study is plausible.

Early transition in microchannel flows has been reported in several experimental studies as we already indicated in the introduction. There are two data points related to transition in Figure 8. The first is from Mala and Li [10], who indicated $Re_t \approx 300$ in microtubes with hydraulic diameters $D_h = 130~\mu\text{m}$. The second is from trapezoidal silicon microchannel experiments of Qu et al. [11, Fig. 8b, p. 362] and is shown by a bold square in Figure 8. The results obtained by these authors at $D_h = 142~\mu\text{m}$ suggest a transitional number of 400. There is a satisfactory agreement between the estimated transitional numbers and the experiments in Figure 8, but that may be coincidental. Controlled experiments in a way similar to those reported by Ren et al. [22] have to be conducted. That may be achieved, for instance, by keeping the same channel with the same roughness distribution and by changing the ionic concentration of the liquid. Somewhat unfortunately, the maximum Debye length in Ren et al.'s [22] experiments is four times smaller than the reference case taken here. Ren et al. [23] and Sze et al. [24]
have measured the zeta potential ($\zeta$ = 245 mV), Debye length ($\frac{1}{k}$ = 0.305 $\mu$m), and the bulk conductivity ($\lambda_0$ = $1.053 \times 10^{-4}$ S/m) in microchannels of height varying from 14.1 $\mu$m to 40.5 $\mu$m. Note that the conductivity is $1.3 \times 10^3$ times larger than the reference case. The present estimation is only qualitatively valid in this case because the zeta potential is too high for the Poisson equation linearization. Furthermore, the important contribution of the corners to the EDL field [25, 26] and to the local flow instability has to be considered to model Ren et al.'s [22] experiments. With these restrictions in mind, one has $\kappa$ = 23.6 and $G\bar{\zeta}^2$ = 500 for the 14.1 $\mu$m height channel reported by these authors. Up to 20% higher flow resistance is found and the transitional Reynolds number decreases by some 63%, with $Re_t$ ≈ 270 in this particular situation. There is a significant effect on the transition, but the resulting $Re_t$ is presumably too high to be experimentally reached in such a small device. In Wang and Peng [6] the value of $\kappa$ is presumably larger than 120 and the low transitional Reynolds numbers of approximately 300 they report cannot be explained by the EDL effect. Strong effect of roughness in the inlet region is suspected in these experiments.

Early transition in microchannels can have interesting applications for enhancement of mixing together with heat and mass transfer processes. A simple estimation shows that the ratio of the Nusselt numbers in EDL and macroscale flow is roughly inversely proportional to the ratio of the transitional Reynolds numbers. A decrease of $Re_t$ by a factor $f$ therefore implies an enhancement of the heat transfer by the same $f$ in some giving working conditions.

Fully developed turbulent EDL microchannel flow reveals several questions. A transitional Reynolds number of $Re_t$ = 817.5 ($\kappa$ = 16, $G$ = 1937) would result in a half height channel in wall units (i.e., scaled by the shear velocity and viscosity) of $a^+ \approx 30$. Recalling that the logarithmic equilibrium layer begins at 30 wall units in macroscale conventional turbulent channel flows, the structure of “micro turbulence” in such situations is open to question. Our team tries to clarify these points by using direct numerical simulations.

## REFERENCES

1. M. Gad-el-Hak, The Fluid Mechanics of Microdevices, *Journal of Fluids Engineering*, vol. 121, pp. 5-33, 1999.

2. S. Tardu, Transferts Thermiques Dans Les Microcanaux, *Traité EGEM*, Tome 6, Ch. 6, Micro-fluidique, Hermès, p. 6, 2004.

3. P. A. Thompson and S. M. Troian, A General Boundary Condition for Liquid Flow at Solid Surfaces, *Physical Review Letters*, vol. 63, pp. 766-769, 1997.

4. K. V. Sharp, R. J. Adrian, J. G. Santiago, and J. I. Molho, Liquid Flows in Microchannels, In M. Gad-el-Hak (ed.), *The MEMS Handbook*, CRC Press, Boca Raton, FL, 2002.

5. L. Kulinsky, Y. Wang, and M. Ferrari, Electroviscous Effects in Microchannels, *Proceedings of SPIE, The International Society of Optical Engineering*, vol. 3606, pp. 158-168, 1999.

6. B. X. Wang and X. F. Peng, Experimental Investigation on Liquid Forced Convection Heat Transfer through Microchannels, *International Journal of Heat and Mass Transfer*, vol. 37, pp. 73-82, 1994.

7. M. M. Rahman and F. Gui, Experimental Measurements of Fluid Flow and Heat Transfer in Microchannel Cooling Passages in a Chip Substrate, *ASME EEP*, vol. 4, pp. 685-692, 1993.

8. M. M. Rahman and F. Gui, Design, Fabrication, and Testing Microchannel Heat Sinks for Aircraft Avionics Cooling, *Proc. 28th Intersociety Energy Conversion Engineering Conference*, vol. 1, pp. 1-6, 1993.
9. F. Gui and R. P. Scaringe, Enhanced Heat Transfer in the Entrance Region of Microchannels, *Proceedings of the 30th Intersociety Energy Conversion Engineering Conference*, vol. 2, pp. 284–294, 1995.

10. G. M. Mala and D. Li, Flow Characteristics of Water in Microtubes, *International Journal of Heat and Fluid Flow*, vol. 20, pp. 142–148, 1999.

11. W. Qu, G. M. Mala, and D. Li, Pressure-Driven Water Flows in Trapezoidal Silicon Microchannels, *International Journal of Heat and Mass Transfer*, vol. 43, pp. 353–364, 2000.

12. J. Lykelma and M. Minor, On Surface Conduction and Its Role in Electrokinetics, *Colloids and Surfaces A*, vol. 140, pp. 33–41, 1998.

13. G. M. Mala, D. Li, and J. D. Dale, Heat Transfer and Fluid Flow in Microchannels, *International Journal of Heat and Mass Transfer*, vol. 40, pp. 3079–3088, 1997.

14. C. Yang, D. Li, and J.-H. Masliyah, Modeling Forced Liquid Convection in Rectangular Microchannels with Electrokinetic Effects, *International Journal of Heat and Mass Transfer*, vol. 41, pp. 4229–4249, 1998.

15. C. H. Von Kerczek, The Instability of Oscillatory Poiseuille Flow, *Journal of Fluid Mechanics*, vol. 116, pp. 91–114, 1982.

16. S. A. Orszag, Accurate Solution of Orr-Somerfeld Stability Equation, *Journal of Fluid Mechanics*, vol. 50, pp. 689–703, 1971.

17. C. E. Grosch and H. Salwen, The Stability of Steady and Time-Dependent Plane Poiseuille Flow, *Journal of Fluid Mechanics*, vol. 34, pp. 177–205, 1968.

18. S. A. Orszag and A. T. Patera, Secondary Instability of Wall Bounded Shear Flows, *Journal of Fluid Mechanics*, vol. 128, pp. 347–385, 1983.

19. P. G. Drazin and W. H. Reid, *Hydrodynamic Stability*, Cambridge University Press, Cambridge pp. 370–464, 1981.

20. H. Schlichting, Boundary Layer Theory, 7th edition, New York: McGraw–Hill, 1979.

21. J. T. Stuart, On the Non-Linear Mechanics of Wave Disturbances in Stable and Unstable Parallel Flows. Part 1. The Basic Behaviour in Plane Poiseuille Flow, *Journal of Fluid Mechanics*, vol. 9, pp. 353–370, 1960.

22. L. Ren, D. Li, and W. Qu, Electro-Viscous Effects on Liquid Flow in Microchannels, *Journal of Colloid and Interface Science*, vol. 233, pp. 12–22, 2001.

23. L. Ren, W. Qu, and D. Li, Interfacial Electrokinetic Effects on Liquid Flow in Microchannels, *Interational Journal of Heat and Mass Transfer*, vol. 44, pp. 3125–3134, 2001.

24. A. Sze, D. Erickson, L. Ren, and D. Li, Zeta-Potential Measurement Using Smoluchowski Equation and the Slope of the Current-Time Relationship in Electroosmotic Flow, *Journal of Colloid and Interface Science*, vol. 261, pp. 402–410, 2003.

25. C. Yang and D. Li, Electrokinetic Effects on Pressure-Driven Liquid Flows in Rectangular Microchannels, *Journal of Colloid and Interface Science*, vol. 194, pp. 95–107, 1997.

26. C. Yang and D. Li, Analysis of Electrokinetic Effects on the Liquid Flow in Rectangular Microchannels, *Colloids and Surfaces A*, vol. 143, pp. 339–353, 1998.

27. P. Gao, S. Le Person, and M. Favre Marinet, Scale Effects on Hydrodynamics and Heat Transfer in Two-Dimensional Mini and Microchannels, *International Journal of Thermal Science*, vol. 41, p. 10, 2002.

28. N. Itoh, Spatial Growth of Finite Wave Disturbances in Parallel and Nearly Parallel Flows. Part 1. The Theoretical Analysis and the Numerical Results for Plane Poiseuille Flow, *Transactions of the Japan Society for Aeronautical and Space Science*, vol. 17, pp. 160–174, 1974.

29. X. F. Peng, G. P. Peterson, and B. X. Wang, Heat Transfer Characteristics of Water Flowing through Microchannels, *Experimental Heat Transfer*, vol. 7, pp. 265–283, 1994.