Optimizing electroosmotic flow in an annulus from Debye Hückel approximation to Poisson–Boltzmann equation

Abstract

In this study, we consider steady and starting electroosmotic (EO) flow in an annulus channel with different zeta potentials (denoted by α, β, respectively) developed on the inner and outer channel walls. An analytical solution is first obtained under the linearized Debye–Hückel approximation (DHA), and then extended to the non-linear Poisson–Boltzmann equation (PBE) by the method of homotopy. The steady-state EO pumping rate for any given pair of (α, β) is optimized with respect to the ratio b of inner to the outer radii, and the corresponding temporal developments of the optimal EO flow are given detailed investigation of the effects of the electric double layer (EDL). The optimal EO pumping rates Q_M(α,β) are presented on the α–β plane for several electrokinetic widths (K) to illustrate their general trends versus the corresponding b (denoted by b_max), which serves as a useful guide for practical applications. Investigation is also given to the shifts of b_max and Q_M(α,β) with varying the parameter λ, which measures the nonlinearity of the PBE.

---

1 Introduction

In recent years, electro-osmosis has been widely used as a powerful method to transport fluids in micro-tubes.¹ The principle is to make use of the EDLs associated with the zeta potential in a micron- or nano-scale channel to transport, mix or separate electrolytes by externally applying a steady or time-varying voltage.^2–4 For example, Jia et al. presented a novel microfluidic device using induced-charge electroosmosis (ICEO) flow to focus and collect microparticles in the suspending medium.⁵ Das et al. presented a microfluidic device for trapping and detection of oil droplets in water by using an alternating current (AC) electroosmotic.⁶ Sasaki et al. fabricated a micromixer to rapidly mix the fluids by using AC electroosmotic flow.⁷ Jellema et al. developed a new microfluidic approach for charge-based particle separation using combined hydrodynamic and electrokinetic effects.⁸ Witek et al. evaluated the electrokinetic transport properties of the cells in microdevices fabricated in pristine and ultraviolet-modified poly(methyl methacrylate) (PMMA) and polycarbonate (PC).⁹

One out of the many important applications is the decontamination of porous media such as soils.¹⁰ It is a high-efficiency and quick method that decontaminate and repair saturated clay with low penetrability which has been polluted by heavy metal, and the heavy metal removal is over 90%.¹¹ Annamalai et al. effectively decontaminated and repaired the textile dye-contaminated field soil¹² and textile effluent-contaminated farming soil¹³ by using the electrokinetic remediation technology. On the theoretical side, Wu and Papadopoulos¹⁴ suggested that one model for porous media is the annular pore, and obtained an exact EO solution through the Debye–Hückel approximation. In 2000, the same solution was found independently by Tsao.¹⁵ The unsteady starting EO flow in an annulus was solved by Chang and Wang,¹⁶ and the oscillatory EO flow by Jian et al.¹⁷ Wu and Papadopoulos noted that under certain conditions altering the inner cylinder (core) radius may increase or decrease the EO flow. However, no details were supplied. On the other hand, the study of Chang and Wang¹⁶ is limited to annular tubes with the same zeta potentials on both inner and outer walls. In this study, we present a thorough investigation about how the core radius (with a different zeta potential) may influence the EO pumping rate.

Most solid surfaces are charged when in contact with a dissociating fluid. The surface charge varies greatly, and may assume different signs. Take proteins as an example, whose charge is of complicated nature. Many proteins such as the bacteriophage capsid have a transition from net-positive to net-negative charge depending on the solution pH, the salt concentration as well as the detailed distribution of amino acids.^18,19 Annular tubes represent an interesting example with two channel walls which may develop different zeta potentials. The purpose of this study is to provide semi-analytical solutions and investigate physical significance for both steady and transient EO flows in annular tubes from the linear D–H approximation to the nonlinear P–B equation.

2 Basic equations

Fig. 1a and b shows a schematic diagram and the cross section of the annular channel for EO flow considered in this study.

Fig. 1 (a) Schematic of the EO flow along an annular tube in the three-dimensional setting. (b) The cross section of the annulus channel. The radii of the inner and outer circles are bL (0 < b < 1) and L, and the zeta potentials ζ_i and ζ_o developing on the inner and outer walls are taken to be constant and uniform.

Assume that the electrolyte is a Newtonian fluid, and the fluid flow is solely caused by electro-osmosis. The fluid motion is governed by the Navier–Stokes (N–S) equation²⁰

where u′ is the fluid velocity, t′ the time, ρ the fluid density, μ the fluid viscosity, ρ_e the charge density, and E = (0,0,E(t′)) – externally applied longitudinal electric field. The flow of this study is parallel to the axis of the annular tube, thus eqn (1) with u′ = (0,0,u′) reduces towhich is subject to the no slip condition u′(bL) = u′(L) = 0.

The Poisson–Boltzmann equation for a symmetric electrolyte is^18,20

where ψ′ is the electric potential (associated with the EDL), ε is the electric permittivity of the fluid, z is the valence, e is the proton charge, n₀ is the bulk electrolyte concentration, k_b is the Boltzmann constant, and T is the absolute temperature. The boundary condition for ψ′ in eqn (3) is the zeta potential ζ given on the shear plane where the no-slip condition is imposed on the flow velocity u′ (see, e.g., ref. 21). For brevity, the interface at which we apply/prescribe the boundary conditions for both u′ and ψ′ is hereafter simply called a wall (either the outer or the inner one).

The flow is started by applying an external electric field suddenly at t′ = 0,

where E₀ is the strength of electric field.

In this study, we normalize the coordinates x′, y′ and radial distance r′ by setting x′ = Lx, y′ = Ly and r′ = Lr. Let ζ_M be the maximum zeta potential in the range of interest. The zeta potentials on the inner and outer walls are normalized to be α and β with α = ζ_i/ζ_M and β = ζ_o/ζ_M. Normalize all the variables in the N–S and P–B equations with appropriate reference scales: u = u′/(−E₀εζ/μ), t = t′/(ρL²/μ), ∇² = L²∇′², ψ = ψ′/ζ_M and H(t) = E(t′)/E₀. Eqn (2) and (3) can then be written in the non-dimensional form

N–S:

u_t = ∇²u + ∇²ψH(t) (4)

P–B:

where we have defined two non-dimensional parameters: the electrokinetic width
and

It is noted that λ measures the strength of the zeta potential relative to the thermal potential, representing the degree of nonlinearity of the PB equation. K is the ratio of the radius L and Debye length λ_D, therefore also called the non-dimensional electro-kinetic width. It indicates a relatively thin EDL if K is large and a thick one if K is small. For example, the capillary electrophoresis applications usually have K at the orders of about 10–100 with the Debye lengths (λ_D) in the range of about 1–10 nm in a microfluidic channel (with characteristic width, say, 100 nm). However, a smaller K (with overlapped EDLs) can be produced by using a smaller microchannel (i.e., the radius L on the scale of nanometers) or decreasing the bulk concentration (low n₀). In this study, we choose K > 3 in order that the result is not further complicated by the issues of overlapped EDLs. The range of λ values is of practical concern. Take, for example, T = 288 K (15 °C), k_bT = 0.0248 eV (k_b = 8.62 × 10⁻⁵ eV K⁻¹). The value of ζ = 25 mV (either positive or negative), can be taken as a useful reference that separates low-charged surfaces from highly-charged surfaces. In CE (capillary electrophoresis) applications, the wall zeta potential ζ may vary in a wide range of values (say, from 0 to ±75 mV, i.e., λ = 0–3z), which requires solution of the nonlinear Poisson–Boltzmann equation.

Let Q be the flow rate Q′ normalized by −E₀εζL²/μ, then we have

3 Debye–Hückel approximation

For low zeta potentials (λ → 0), eqn (5) is reduced to the linear Debye–Hückel approximation

∇²ψ = K²ψ (8)

subject to the boundary conditions ψ(b) = α and ψ(1) = β, whereas the flow velocity satisfies the no-slip condition u(b) = u(1) = 0.

3.1 Steady EO flow in an annulus

When the flow in the tube is steady (u_t = 0), we denote the steady-state flow velocity by u₀. Then we have from eqn (4)

∇²u₀ + ∇²ψ = 0 or ∇²(u₀ + ψ) = 0 (9)

Since both u₀ and ψ are radially symmetric, we can easily solve eqn (9)-Laplace equation for u₀ + ψ to obtain

so that the conditions u₀ + ψ = α for r = b and u₀ + ψ = β for r = 1 are satisfied. Moreover, the fact that eqn (8) is linear suggest that we represent the general solution of ψ in terms of the solutions of two fundamental problems

ψ = αψ₁ + βψ₂ (11)

where ψ₁ and ψ₂ satisfy eqn (8) and

(A1) ψ₁(b) = 1, ψ₁(1) =0 (12)

(A2) ψ₂(b) = 0, ψ₂(1) = 1 (13)

Note that ψ₁ is higher on the inner wall and ψ₂ is higher on the outer wall. Similarly we can decompose the velocity and the flow rate

u₀ = αu₁ + βu₂ (14)

Q₀ = αQ₁ + βQ₂ (15)

The solutions ψ₁ and ψ₂ are radially symmetric, and so we can recast eqn (8) in polar coordinates

This is the modified Bessel equation of order 0 (see, e.g., ref. 22, p. 71).

(A1) The solution ψ₁ to eqn (16) can be obtained in terms of the modified Bessel functions with the boundary condition (12)

where I₀ and K₀ are the modified Bessel functions of the first and second kind in the zeroth order, respectively. As a word of caution, K₀ (with the suffix 0) denotes the modified Bessel function, and should not be confused with the electro-kinetic with K.

Eqn (10) and (12) give

because this case corresponds to (α,β) = (1,0).

From eqn (7), (17) and (18), we may derive the flow rate,

Since Q₁ is zero at b = 0 and b = 1. There must be a maximum flow rate as the inner radius increases from zero. We see from eqn (19) that given K, Q₁ is a function of b. Fig. 2a shows the EO flow rate Q₁ versus b at different values of K. Note the logarithmic rise in flow as inner radius b increases from zero. Two velocity distributions (cf. eqn (18)) are presented in Fig. 2b and c. Fig. 2a shows that for all K values, the flow rate Q₁ rapidly rises to the maximum as b increases from zero, then the flow rate decreases to zero with continuing increase of b. And the optimum inner radius b for maximizing Q₁ obviously becomes larger at larger K. Fig. 2b shows that the maximal velocity occurs near the inner wall (i.e. the wall which has the higher zeta potential). Fig. 2b and c show that for larger K, the entire velocity profile is higher, and the location of the maximum velocity is closer to the inner wall. Note that in Fig. 2b and c, these plots are drawn for ru₁ instead of u₁ for easily seeing the contribution of the velocity to the EO pumping rate.

Fig. 2 (a) Q₁ versus b (from top: K = 100, 50, 20, 10, 5, 3). The maximum for each K curve is denoted by a red point. The red spot shifts to larger b with increasing K. (b) Velocity ru₁ versus radius r corresponding to optimum inner radius b = 0.397 for K = 100, and b = 0.03, 0.6. (c) Velocity ru₁ versus radius r corresponding to optimum inner radius b = 0.2988 for K = 10, and b = 0.02, 0.6.

(A2) Next we study the case when the outer wall has unity zeta potential and the inner wall zeta potential is zero. The solution for ψ₂ to eqn (16) with the boundary condition (13) is given by

From eqn (10) and (13), we have

because this case corresponds to (α,β) = (0,1).

From eqn (7), (20) and (21), we can obtain the flow rate,

Fig. 3a shows the variations of the flow rate Q₂ with respect to b for different K. The introduction of even a small core causes a logarithmic decrease of the flow rate. Fig. 3b shows the corresponding velocity distribution that the larger K, the larger velocity and closer the location of the maximum velocity to the outer wall (i.e. b = 1, the wall which has the higher zeta potential).

Fig. 3 (a) Q₂ versus b (from top: K = 100, 50, 20, 10, 5). (b) The velocity distribution in the tube for b = 0.5 (from top: K = 100, 50, 20, 10, 5, 3).

The general case. After solving the two fundamental problems (A1) and (A2), we can readily obtain the general solutions from eqn (11), (14) and (15) for differentially charged walls. Fig. 4a and b shows the flow rate Q₀ with varying b for different K with α < 1 and β = 1. Recall that the subscript “0” denotes steady state for the general case.

Fig. 4 (a) Q₀ versus b for α = 0.1, β = 1 (from top: K = 100, 50, 20, 10, 5, 3). (b) Q₀ versus b for α = 0.5, β = 1 (from top: K = 100, 50, 20, 10, 5, 3).

Fig. 5a–d shows the flow rate Q₀ with varying b for different K, α = 1 and β ≤ 1. From the results, we find that if the zeta potential on the outer wall is equal or larger than that on the inner wall, the maximum flow always occurs when b = 0, or without the inner cylinder. If the zeta potential on the inner wall is larger than that on the outer wall, adding an interior core may increase the flow rate, but only if the value of K is large enough. It is noted that if both walls have the same zeta potential (α = β = 1), a core only hinders the flow. If α = 1, β < 0.3, then adding a core is beneficial for the values of K ≥ 3.

Fig. 5 (a) Q₀ versus b for α = 1, β = 0.1 (from top: K = 100, 50, 20, 10, 5, 3). (b) Q₀ versus b for α = 1, β = 0.3 (from top: K = 100, 50, 20, 10, 5, 3). (c) Q₀ versus b for α = 1, β = 0.5 (from top: K = 100, 50, 20, 10, 5, 3). (d) Q₀ versus b for α = 1, β = 1 (from top: K = 100, 50, 20, 10, 5, 3).

Fig. 6a–c presents selected maps for the maximal EO flow rate Q_M(α,β)-specified by the level curves, and its corresponding optimal radius b_max for K = 100, 20, 5. The level of Q_M increases toward the right upper corner as both α and β are increased. It is noted that in the upper left of b = 0 of these maps, adding a core would effectively increase the flow rate. Conversely, in the bottom right of b = 0, a core would absolutely hinder the EO flow. If α < β, we always have the maximal EO flow rate when there is no inner core (b = 0). However, if α > β, we may optimize the EO pumping rate with respect to b. The most interesting finding is that for large K (=100) the b_max that optimizes the EO pumping shifts in the manner of an expansion fan of straight lines from b_max = 0 at α/β = 1 to 0.4 at α/β = ∞ (β = 0). The pattern of straight lines in b_max of these maps is reasonable because Q_M(α,β) should depend on the ratio of α and β, not on their individual magnitudes. For a smaller K, the fan starts at a larger α/β (=1.22 for K = 5), and terminates at a smaller b_max (=0.22 for K = 5) at α/β = ∞ (β = 0). As we increase the nonlinearity λ, b_max moves to a larger value; this trend is more significant for smaller K (=5), and is very little for large K (=100).

Fig. 6 (a) Optimal EO pumping rate Q_M(α,β) on the α–β plane for K = 100. Level curves of constant Q_M and b_max are also specified in the plot. (b) Optimal EO flow rate Q_M(α,β) on the α–β plane for K = 20. Level curves of constant Q_M and b_max are also specified in the plot. (c) Optimal EO pumping rate Q_M(α,β) on the α–β plane for K = 5. Level curves of constant Q_M and b_max are also specified in the plot.

3.2 Unsteady starting EO flow in an annulus

Following the approach of Chang and Wang,¹⁶ we may decompose the velocity for the transient state,

u = u₀ − ũ (23)

where u₀ is the steady-state velocity and ũ is the transient flow velocity deficit, decaying to zero as time evolves.

Substituting eqn (23) and (9) in eqn (4), we have

ũ_t = ∇²ũ, t > 0 (24)

Let

Eqn (24) gives

The solution is

f_n(r) = Y₀(γ_n)J₀(γ_nr) − J₀(γ_n)Y₀(γ_nr) (27)

where J₀ and Y₀ are the Bessel functions of the first and second kind in the zeroth order and γ_n are the eigenvalues, satisfying

Y₀(γ_n)J₀(γ_nb) − J₀(γ_n)Y₀(γ_nb) = 0 (28)

From the asymptotic expansion of Bessel function,²² we have

Table 1 shows that the approximation (29) compares well with those computed with great accuracy, in particular, at larger b or greater n.

Table 1 Accurate values of γ_n for various radius ratios b of the annular channel. The approximate values from eqn (29) are also shown in the parentheses

n\γn	b = 0.1	b = 0.25	b = 0.5
1	3.315(3.491)	4.098(4.189)	6.246(6.283)
2	6.858(6.981)	8.324(8.378)	12.55(12.57)
3	10.38(10.47)	12.53(12.57)	18.84(18.85)
4	13.89(13.96)	16.73(16.76)	25.12(25.13)
5	17.39(17.45)	20.92(20.94)	31.41(31.42)
13	45.35(45.38)	54.44(54.45)	81.68(81.68)
14	48.84(48.87)	58.65(58.64)	87.96(87.97)
15	52.34(52.36)	62.82(62.83)	94.25(94.25)

---

The orthogonality of the Bessel functions f_n gives

with

c_n = {[Y₀(γ_n)J₁(γ_n) − J₀(γ_n)Y₁(γ_n)]² − b²[Y₀(γ_n)J₁(γ_nb) − J₀(γ_n)Y₁(γ_nb)]²/2} (31)

Recall the two fundamental problems for steady flow where we have set

(A1) u₀ = u₁ for (α,β) = (1,0)

(A2) u₀ = u₂ for (α,β) = (0,1)

The relationships: eqn (11), (14) and (15) are also valid for transient solution. Let

ũ = αũ₁ + βũ₂ (32)

(A1) For ψ = ψ₁, eqn (18) and the D–H approximation (8) give

Note that ũ₁ is the velocity deficit for ψ = ψ₁; then we recall eqn (25) and write

At t = 0, we have

Substituting (36) in eqn (34), and using eqn (30), we obtain

(A2) For ψ = ψ₂, eqn (21) and the D–H approximation (8) give

Now that ũ₂ is the velocity deficit for ψ = ψ₂; then we recall eqn (25) and write

At t = 0, we have

Substituting (40) in eqn (38), and using eqn (30), we obtain

After solving (A1) and (A2) for the transient problems, we may obtain the general solution for the velocity deficit from eqn (32)

and the EO flow rate deficit from eqn (33)where , are given by eqn (37) and (41), respectively, f_n are given by eqn (27), and γ_n are given by eqn (29).

Then we have the transient velocity and EO flow rate for the general pair of zeta potentials (α,β)

u = u₀ − ũ (44)

Q = Q₀ − Q̃ (45)

Fig. 7a and b shows some typical velocity distributions u at various times while different zeta potentials develop on the two walls for K = 10, b = 0.5. It is seen that the entire velocity increases as time evolves, and finally reaches the steady state. The electric double layer (EDL), a thin layer next to the solid walls, serves as the sole pumping unit of the EO flow when the external voltage is applied. Hence the EO pumping first appears near the walls at initial times, then the pumping gradually extends to the center of the gap by diffusion of momentum. The location of the maximal velocity is close to the wall that has the higher zeta potential at all times, but shifts away from that wall as time evolves.

Fig. 7 (a) u versus r (b = 0.5) at different times with α = 1, β = 0.5 and K = 10 (from bottom t = 0.001, 0.01, 0.03, 0.05, 0.1, 10). The maximum for each t curve is denoted by a red point. (b) u versus r (b = 0.5) at different times with α = 0, β = 1 and K = 10 (from bottom t = 0.001, 0.01, 0.03, 0.05, 0.1, 10).

Fig. 8a and b shows some typical velocity distributions u at various times while we have the same zeta potential on the two walls for K = 20 and 3 with fixed b = 0.5. A large K indicates a relatively thin EDL, therefore for large K the EO flow only appears near the walls, and the electrolyte near the channel center would not be influenced at initial times. From Fig. 8a and b, we can see that for large K there are two maxima of the velocity near the walls at small times after the flow is started, but eventually there is only one maximal velocity as both maxima shifts towards the channel center as time evolves. On the other hand, for small K there is only one maximal velocity soon after the flow is started.

Fig. 8 (a) u versus r (b = 0.5) at different times with K = 20 and α = β = 1 (from bottom, t = 0.001, 0.01, 0.03, 0.05, 0.1, 10). (b) u versus r (b = 0.5) at different times with K = 3 and α = β = 1 (from bottom, t = 0.001, 0.01, 0.03, 0.05, 0.1, 10).

Fig. 9a and b shows the temporal evolutions of the flow rate Q for different K and b with α = 0.5 and β = 1. As shown in these plots, the flow rate reaches the steady state at the same time for all K values. This is consistent with eqn (42) which shows that the transient time is dominated by 1/γ₁² and is independent of K, but almost linear in (1 − b) (cf. eqn (29)), namely, the smaller b, the greater the transition time. In addition, the EO flow rate for larger K increases more rapidly than that for smaller K at initial stages of flow development though the transient-time scales are the same.

Fig. 9 (a) The transient EO flow rate versus time for the annulus channel with α = 0.5, β = 1 and b = 0.5. (b) The transient EO flow rate versus time for the annulus channel with α = 0.5, β = 1 and b = 0.75.

4 Extension to the Poisson–Boltzmann equation

For higher zeta potentials (λ ⟪ 0), the DHA is no longer valid, the analysis must be performed based on the nonlinear P–B equation.

It is difficult to solve nonlinear equations, but once we have solved the linear DHA, there are several approaches to extend the results to the nonlinear P–B equation. Here we employ the method of homotopy. The method of homotopy is powerful, accurate and suitable for many highly nonlinear problems.²³ Recently, homotopy methods have been widely applied to the fields of science, economics²⁴ and engineering.²⁵ This method proves to be very efficient in leading to convergent results compared to some traditional numerical methods.²⁶ In this study, we combine Newton's and homotopy methods to solve the nonlinear Poisson–Boltzman equation.

4.1 Method of homotopy

The boundary conditions are

ψ(b) = α, ψ(1) = β (46)

Define

L(ψ) = ∇²ψ − K²ψ (48)

Introduce the homotopy function

H(ψ,s) = sF(ψ) + (1 − s)L(ψ) (49)

The “end conditions” are given by eqn (47) and (48)

H(ψ,0) = L(ψ) = 0 (50)

H(ψ,1) = F(ψ) = 0 (51)

The general idea of homotopy is to start with the known solution for eqn (50) (from DHA), and progressively approach to the final solution for the P–B eqn (51) by solving the homotopy eqn H(ψ,s) = 0 with gradually increasing the homotopy parameter s to ensure convergence of the solution at each increment of s.

Let 0 = s₀ < s₁ < … < s_M = 1, M is an arbitrary integer. In this study, it is sufficient to set M = 5, and take s₁, s₂, s₃ and s₄ are 0.2, 0.4, 0.6 0.8, respectively.

(A) For s = s₀, there is H(ψ₀,s₀) = L(ψ₀) = 0, with ψ₀(b) = α, ψ₀(1) = β.

And the initial solution is

(B) For s = s₁, there is

H(ψ₁,s₁) = s₁F(ψ₁) + (1 − s₁)L(ψ₁) = 0 with ψ₁(b) = α, ψ₁(1) = β (53)

where

The nonlinear equation H(ψ₁,s₁) = 0 is then discretized by the central finite difference. Let h = (1 − b)/N, i = 1 ∼ (N − 1), eqn (54) becomes

Or, we can simply put it in the compact form

H = Aψ₁ + G = 0 (55)

where ψ₁ = [ψ₁₁, ψ₁₂, …, ψ_1(N−1)]^T, and A is the coefficient matrix of ψ₁. Then let ψ₁(0) = ψ₀, where ψ₀ is the initial solution. Apply Newton's method, we obtain the iterated solution till convergence

ψ₁(j + 1) = ψ₁(j) − J⁻¹H[ψ₁(j)] (56)

where and

(C) Repeat the process B for s = s₂, s = s₃… until s_M = 1 to obtain the nonlinear solution ψ = ψ_M for the P–B equation.

Fig. 10a–d shows the iterated solutions for the EDL potentials of s₀ ∼ s₅ with K = 5, b = 0.5 for several pairs of (α,β). As shown in these figures, the linear D–H approximation remains a good approximation for small λ. Yet, the discrepancy between DHA and PBE is getting large as λ is further increased. In addition, from these figures we can also see clearly that a larger λ corresponds to a stronger zeta potential in the entire channel, especially at the center, and the EDL is more tightly restricted in a thin layer next to the channel walls.

Fig. 10 (a) The sequence of homotopy solutions for the PBE (λ = 2) with α = 0.8, β = 0.3, K = 5, b = 0.5 (from top: s₀ ∼ s₅). (b) The sequence of homotopy solutions for the PBE (λ = 10) with α = 0.8, β = 0.3, K = 5, b = 0.5 (from top: s₀ ∼ s₅). (c) The sequence of homotopy solutions for the PBE (λ = 2) with α = 1, β = 1, K = 5, b = 0.5 (from top: s₀–s₅). (d) The sequence of homotopy solutions for the PBE (λ = 10) with α = 1, β = 1, K = 5, b = 0.5 (from top: s₀–s₅).

4.2 Steady flow

From eqn (4) we still have ∇²(u₀ + ψ) = 0 for the steady state of P–B equation, and the general solution for eqn (9) with u₀ + ψ = α for r = b and u₀ + ψ = β for r = 1 is given by

As an illustration, let us fix b = 0.5. Fig. 11a–d shows the velocity distributions of the P–B equation with various λ for K = 10 or 5. The result for the λ = 0 case is from the D–H approximation.

Fig. 11 (a) Steady velocity distributions associated with the EDL potential of the PBE (K = 10) with α = 1, β = 1 (from bottom: λ = 0, 0.1, 1, 2, 3, 5). (b) Steady velocity distributions associated with the EDL potential of the PBE (K = 5) with α = 1, β = 1 (from bottom: λ = 0, 0.1, 1, 2, 3, 5). (c) Steady velocity distributions associated with the EDL potential of the PBE (K = 10) with α = 0, β = 1 (from bottom: λ = 0, 0.1, 1, 2, 3, 5). (d) Steady velocity distributions associated with the EDL potential of the PBE (K = 5) with α = 0, β = 1 (from bottom: λ = 0, 0.1, 1, 2, 3, 5).

As expected, the maximal velocity occurs at the channel center for the symmetry condition (α = 1, β = 1). The location of the maximal velocity is closer to the wall that has the higher zeta potential for the unsymmetric condition (α = 0, β = 1). The EO pumping rate increases with increasing λ, and the increase is more efficient at smaller K.

4.3 Unsteady starting flow

Again we write u = u₀ − ũ, at t = 0. Eqn (57) gives

Using eqn (30), we find

Then we have the velocity and flow rate in unsteady state from eqn (25) and (7). Fig. 12 shows typical velocity distributions at various times while the two channel walls develop the same zeta potentials (α = β = 1) on two walls for λ = 5 with K = 10 and b = 0.5. As shown in Fig. 12, there are two maxima of the velocity near the walls at small times after flow is started, but there is only one maximal velocity near the center at the steady state.

Fig. 12 The transient velocity distributions at various times associated with PBE (λ = 5) with α = β = 1, K = 10 and b = 0.5 (from bottom, t = 0.001, 0.01, 0.03, 0.05, 10).

Fig. 13 shows typical velocity distributions at various times while there are different zeta potentials on the two channel walls for λ = 2 with K = 10 and b = 0.5. It is seen that the velocity increases gradually as time evolves, and finally reaches the steady state, and the location of the maximal velocity is at all times close to the wall which has the higher zeta potential.

Fig. 13 The transient velocity distributions at various times associated with PBE (λ = 2) for α = 0.8, β = 0.3 with K = 10 and b = 0.5 (from bottom, t = 0.001, 0.01, 0.03, 0.05, 10).

Fig. 14a and b shows variations of the flow rate Q₀ with increasing t for different K and λ, and α = β = 1, b = 0.5. From these curves, we see that the transient time scale is independent of λ though the flow rate is larger for larger λ. In addition, for a smaller K the flow rate increases rapidly as λ increases, but the influence of λ will be diminishing for larger K.

Fig. 14 (a) K = 5, α = β = 1, b = 0.5 (from bottom, λ = 0.1, 1, 2, 3, 5, 8, 10). (b) K = 10, α = β = 1, b = 0.5 (from bottom, λ = 0.1, 1, 2, 3, 5, 8, 10).

Fig. 15a and b shows variations of the flow rate Q₀ versus b for different λ with K = 5, α = 1 and β = 0.1. It is seen that the flow rate Q rapidly rises to the maximum as b increases from zero, then decreasing to zero with further increasing b. From Fig. 15a and b, we can also see that for this small K = 5, the flow rate Q could be markedly increased by increasing λ, but the influence of the nonlinearity is diminishing for large K. In addition, it is indicated in Fig. 15a that increasing the nonlinearity λ shifts the optimum inner radius toward a larger value.

Fig. 15 (a) Q versus b for K = 5 and α = 1, β = 0.1 (from top: λ = 10, 5, 3, 2, 1, 0.1). The maximum for each λ curve is denoted by a red point. (b) Q versus b for K = 50 and α = 1, β = 0.1 (from top: λ = 10, 5, 3, 2, 1, 0.1).

Table 2 shows the optimum sizes of the core for various λ and K. Apparently, increasing λ would increase the optimum radius b_max and the corresponding flow rate Q_M. For a small K, the nonlinearity markedly shifts the optimum radius b_max and increases the flow rate Q_M, yet this effect is small for larger K. In addition, at smaller λ (≤2), the nonlinearity has little effect on the optimum radius b_max and the flow rate Q_M, and conversely for large λ that effect is very much significant.

Table 2 Optimal core radii b_max for α = 1, β = 0.1. The corresponding EO pumping rate Q_M are listed in the parentheses. The results for the λ = 0 case are obtained from the Debye–Hückel approximation

λ\K	5	10	20	100
0	0.1779	0.2575	0.309	0.352
	(0.585)	(0.8286)	(0.9677)	(1.0821)
0.1	0.178	0.2575	0.309	0.352
	(0.5855)	(0.8287)	(0.9678)	(1.0821)
1	0.1804	0.2588	0.3097	0.3522
	(0.5883)	(0.8308)	(0.969)	(1.0823)
2	0.1877	0.2628	0.3117	0.3526
	(0.5981)	(0.8374)	(0.9728)	(1.0831)
3	0.1991	0.2685	0.3148	0.3532
	(0.614)	(0.8476)	(0.9785)	(1.0842)
5	0.2269	0.2825	0.3219	0.3546
	(0.6587)	(0.8741)	(0.9927)	(1.0867)
10	0.2842	0.312	0.3368	0.3572
	(0.7743)	(0.9345)	(1.0228)	(1.0915)

---

5 Concluding remarks

In this study, we presented semi-analytical solutions of steady and unsteady starting EO flow for an annulus channel. The two channel walls may develop different zeta potentials (α at the inner r = b and β at the outer r = 1). The problem for this simple geometry is actually complicated by five physical parameters, the pair of zeta potentials (α, β), b – the ratio of the inner to outer radii, K – the electrokinetic width as well as λ which measures the nonlinearity of the P–B equation. It is certainly worthy of having some analytical solutions for such a problem with many parameters.

Our closed-form solutions for steady EO flow under DHA are more compact and convenient than previously reported by other authors. In addition, we obtained the analytical series solution for the transient starting EO flow, and numerically extended the results to the P–B equation by the method of homotopy. The EO flow in an annulus channel has some distinguished features, which could be summarized as follows.

(1) If the zeta potential on the outer wall is equal or larger than that of the inner wall, introduction of an interior core would always hinder the flow. If the zeta potential on the inner wall is larger, an interior core may increase the flow rate, but only if the value of K is large enough or β < 0.3.

(2) For large K, the EDL is a thin layer next to the solid walls so that the maximal velocity occurs near the wall right after the flow is started, but its location eventually shifts towards the center of the channel. For the steady EO flow with α = β, the maximal velocity has a location very close to the channel center, but the location is closer to the wall that has the higher zeta potential if α ≠ β.

(3) A small K indicates a relatively thick EDL that the charge distributes more widely across the channel. For the case of α = β, the location of the maximal velocity is at the channel center soon after the flow is started. The unsteady velocity profile has a more uniform distribution compared to the case of a larger K.

(4) The EO flow reaches the steady state at the same transient time for all K values with a given channel geometry. But the smaller b, the greater the transition time. This is consistent with the analysis that the transient time is dominated by 1/γ₁² ∼ (1 − b)². In addition, the EO flow rate for larger K increases faster than that for smaller K at the initial stage of flow development though the transient times are the same.

(5) A larger λ means a stronger zeta potential with EDL tightly restricted in a thin layer next to the channel wall. There are two maxima of the velocity near the walls right after the flow is started, and then the location of the maximal velocity is shifting towards the channel center. For α ≠ β, the maximal velocity occurs at a location closer to the wall of higher zeta potential, and the lager λ, the closer the location is to the wall. The transient-time scale is independent of K and λ though the EO flow rate is larger at a larger K or λ. In addition, for a smaller K (=5, 10) the EO pumping rate increases significantly with increasing λ, but the effect of λ is diminishing for larger K (=20, 100).

(6) If α < β, we always have the maximal EO flow rate when there is no inner core (b = 0). However, if α > β, we may optimize the EO pumping rate with respect to b. The several maps on the α–β plane (Fig. 6a–c) present the important interesting results, and provides a useful guide for practical design. The main finding is that for large K (=100) the b_max that optimizes the EO pumping shifts in the manner of an expansion fan of straight lines from b_max = 0 at α/β = 1 to 0.4 at α/β = ∞ (β = 0). For a smaller K, the fan starts at a larger α/β (=1.33 for K = 5), and terminates at a smaller b_max (=0.21 for K = 5) at α/β = ∞ (β = 0). As we increase the nonlinearity λ, b_max moves to a larger value; this trend is more significant for smaller K (=5), and is very little for large K (=100).

As a final remark, the results of this study are mostly presented in dimensionless quantities. Yet they are readily transformed to their dimensional quantities with referring to the reference scales. For example, the reference time scale is ρL²/μ, which is about 10⁻⁸ s = 10 ns for L = 100 nm if we use the density and viscosity of water. This indicates that transport of electrolyte by EO pumping is very efficient for the duct on the size of nanoscales. Although the present results for the Debye–Hückel approximation are extended to the nonlinear Poisson–Boltzmann equation, in other realistic applications, one may need more accurate models to account for the physics and chemistry such as finite size effects, overlapped EDLs, etc.

Acknowledgements

This work was completed while the corresponding author (C. C. Chang) was visiting Guangxi University. The authors would like to thank the financial support by the Guangxi University and the Ministry of Science and Technology (Taiwan) under Project No. MOST 102-2221-E-002 -169 -MY3.

References

1. D. P. J. Barz and P. Ehrhard, Lab Chip, 2005, 5, 949–958.
2. P. Watts and C. Wiles, Chem. Commun., 2007, 443–467,  10.1039/b609428g.
3. L. OuYang, C. Wang, F. Du, T. Zheng and H. Liang, RSC Adv., 2014, 4, 1093–1101.
4. Z. Yuan, A. L. Garcia, G. P. Lopez and D. N. Petsev, Electrophoresis, 2007, 28, 595–610.
5. Y. Jia, Y. Ren and H. Jiang, RSC Adv., 2015, 5, 66602–66610.
6. D. Das, Z. Yan, N. V. Menon, Y. Kang, V. Chan and C. Yang, RSC Adv., 2015, 5, 70197–70203.
7. N. Sasaki, T. Kitamori and H. B. Kim, Lab Chip, 2006, 6, 550–554.
8. L. C. Jellema, T. Mey, S. Koster and E. Verpoorte, Lab Chip, 2009, 9, 1914–1925.
9. M. A. Witek, S. Wei, B. Vaidya, A. A. Adams, L. Zhu, W. Stryjewski, R. L. McCarley and S. A. Soper, Lab Chip, 2004, 4, 464–472.
10. R. Iyer, Part. Sci. Technol., 2001, 19, 219–228.
11. M. A. Hashim, S. Mukhopadhyay, J. N. Sahu and B. Sengupta, J. Environ. Manage., 2011, 92, 2355–2388.
12. S. Annamalai, M. Santhanam, S. Sudanthiramoorthy, K. Pandian and M. Pazos, RSC Adv., 2016, 6, 3552–3560.
13. S. Annamalai, S. Selvaraj, H. Selvaraj, M. Santhanam and M. Pazos, RSC Adv., 2015, 5, 81052–81058.
14. R. C. Wu and K. D. Papadopoulos, Colloids Surf., A, 2000, 161, 469–476.
15. H. K. Tsao, J. Colloid Interface Sci., 2000, 225, 247–250.
16. C. C. Chang and C. Y. Wang, Electrophoresis, 2008, 29, 2970–2979.
17. Y. Jian, L. Yang and Q. Liu, Phys. Fluids, 2010, 22, 207–216.
18. H. J. Butt, K. Graf and M. Kappl, Physics and Chemistry of Interfaces, Wiley, New York, 2003.
19. R. J. Nap, A. L. Božič, I. Szleifer and R. Podgornik, Biophys. J., 2014, 107, 1970–1979.
20. H. Bruus, Theoretical Microfluidics, Oxford University Press, 2008.
21. J. G. Santiago, Anal. Chem., 2001, 73, 2353–2365.
22. M. Abramowitz and I. A. Stegun, Handbook of mathematical functions with formulas, graphs and mathematical tables, U.S. Govt. Print. Off., 1965.
23. C. M. Bender and S. A. Orszag, Advanced mathematical methods for scientists and engineers, McGraw-Hill, New York, 1978.
24. K. L. Judd, Numerical Method in Economics, MIT Press, 1998.
25. F. Soltanian, M. Dehghan and S. M. Karbassi, Int. J. Comput. Math., 2010, 87, 1950–1974.
26. T. M. Wu, Appl. Math. Comput., 2005, 168, 1169–1174.

---