Coupled lattice Boltzmann method for simulating electrokinetic flows: A localized scheme for the Nernst–Plank model

Accepted Manuscript Coupled lattice Boltzmann method for simulating electrokinetic flows: a localized scheme for the Nernst–Plank model Hiroaki Yoshida, Tomoyuki Kinjo, Hitoshi Washizu PII: DOI: Reference:

S1007-5704(14)00105-1 http://dx.doi.org/10.1016/j.cnsns.2014.03.005 CNSNS 3125

To appear in:

Communications in Nonlinear Science and Numerical Simulation

Please cite this article as: Yoshida, H., Kinjo, T., Washizu, H., Coupled lattice Boltzmann method for simulating electrokinetic flows: a localized scheme for the Nernst–Plank model, Communications in Nonlinear Science and Numerical Simulation (2014), doi: http://dx.doi.org/10.1016/j.cnsns.2014.03.005

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Coupled lattice Boltzmann method for simulating electrokinetic flows: a localized scheme for the Nernst–Plank model Hiroaki Yoshidaa,b,∗, Tomoyuki Kinjoa,b , Hitoshi Washizua,b a

b

Toyota Central R&D Labs., Inc., Nagakute, Aichi 480-1192, Japan Elements Strategy Initiative for Catalysts and Batteries (ESICB), Kyoto University Katsura, Kyoto 615-8520, Japan

Abstract We present a coupled lattice Boltzmann method (LBM) to solve a set of model equations for electrokinetic flows in micro-/nano-channels. The model consists of the Poisson equation for the electrical potential, the Nernst– Planck equation for the ion concentration, and the Navier–Stokes equation for the flows of the electrolyte solution. In the proposed LBM, the electrochemical migration and the convection of the electrolyte solution contributing to the ion flux are incorporated into the collision operator, which maintains the locality of the algorithm inherent to the original LBM. Furthermore, the Neumann-type boundary condition at the solid/liquid interface is then correctly imposed. In order to validate the present LBM, we consider an electro-osmotic flow in a slit between two charged infinite parallel plates, and the results of LBM computation are compared to the analytical solutions. Good agreement is obtained in the parameter range considered herein, including the case in which the nonlinearity of the Poisson equation due to the large potential variation manifests itself. We also apply the method to a two-dimensional problem of a finite-length microchannel with an entry and an exit. The steady state, as well as the transient behavior, of the electroosmotic flow induced in the microchannel is investigated. It is shown that, although no external pressure difference is imposed, the presence of the entry and exit results in the occurrence of the local pressure gradient that causes ∗

Corresponding author. Email address: [email protected] (Hiroaki Yoshida)

Preprint submitted to Commun. Nonlinear Sci. Numer. Simulat.

March 8, 2014

a flow resistance reducing the magnitude of the electro-osmotic flow. Keywords: Nernst–Planck equation, Electro-osmotic flow, Electrical double layer, Lattice Boltzmann method 1. Introduction Recent remarkable developments in micro-/nano-fabrication techniques have prompted direct experimental observations of the electrokinetic phenomena important to micro-/nano-devices. The electrokinetic flow is a typical electrokinetic phenomenon, in which multiple transport processes, such as ionic diffusion, fluid flow, and electrostatic interaction, play important roles. In the thin layer adjacent to the interface between a charged solid wall and an electrolyte solution, referred to as the electrical double layer (EDL), the ion concentration is highly inhomogeneous [1–3], and the net charge within the EDL interacts with an externally applied electric field or pressure field, which causes various phenomena, such as electro-osmotic flow, streaming potential, and salt rejection from microchannels [4–12]. Although these phenomena have been known for over a century, active use of them in engineering applications has only recently attracted attention [13–19]. The most widely used model for analyzing electrokinetic flows consists of the Poisson–Boltzmann (PB) equation for the internal electrical potential and the Navier–Stokes equation for the flow field of electrolyte solutions. The ion concentration field is directly related to the electrical potential through the Boltzmann distribution, which is derived assuming thermodynamic equilibrium, in which the ionic distributions are not affected by the flows of the electrolyte solution. This model is pertinent under the condition that the external electric field is separate from the internal potential field and the system is in steady state. These conditions are satisfied in many applications, e.g., steady electro-osmotic flows through an infinitely long straight channel are correctly analyzed using this model. There are some important cases, however, which are not covered by the PB equation. For example, if the internal and external electric fields must be integrated into a single potential field because of strong nonlinearity of the electrical potential, it is not possible to assume a Boltzmann distribution with a unique reference (or bulk) potential. The PB equation should be replaced by the set of the Poisson equation and the Nernst–Planck equation for such cases, where the Nernst–Planck equation governs the transport processes of each ionic species 2

[20–27]. In the present paper, we propose a numerical framework to solve the model for electrokinetic flows based on the Nernst–Planck equation. More specifically, the model consists of the Poisson equation, the Nernst–Planck equation, and the Navier–Stokes equation. The method is also capable of analyzing transient behavior, which is not possible using existing methods that are based on the PB equation. The numerical scheme for solving each governing equation is based on the lattice Boltzmann method (LBM). The LBM [28, 29] was originally developed as an alternative numerical method for solving Navier–Stoke-type equations and has been extended to solve the convection–diffusion-type equations [30, 31]. Although there are a vast number of alternative schemes for these equations associated with finite element and finite-difference methods, the LBM is nevertheless attractive because it is easy to use for programming and is compatible to parallel and GPU computing. In addition, when we consider the electrokinetic phenomena in complex morphology, such as ion transport in fuel cells [32–34], the LBM is a promising tool in view of the success in flows through porous media. Several LBMs for analyzing electrokinetic flows have already been proposed [21, 22, 27, 35–39], and successfully applied, primarily to two-dimensional analysis. Most of these LBMs are based on the PB equation, whereas the method proposed by Wang and Kang [27] solves the Poisson, Nernst–Plank, and Navier–Stokes equations simultaneously and is capable of handling problems that are not covered by the PB approximation. However, the simplified Dirichlet boundary condition for the Nernst–Planck equation used in Ref. [27] still assumes the Boltzmann distribution near the boundary, which causes some problems when the external electric field is combined with the internal potential. Moreover, fully time-dependent problems are not solved accurately, because of the artificial flux across the boundary. In contrast, the method proposed by Capuani et al. [21] for fluid mixtures that covers the electrokinetic flows based on the Nernst–Plank equation as a special case [22], properly incorporates the original Neumann-type boundary condition. Although the locality of the original LBM is sacrificed because the scheme is based on the flux on the links between the lattice points, which is evaluated using the neighboring lattice points, the artificial flux across the boundary is avoided. Whereas in the method of Capuani et al. the LBM is supplemented with the discrete solution of the Poisson equation obtained using the finitedifference scheme rather than the LBM, in the method proposed herein, all the physical quantities are solved within the framework of the LBM. Partic3

ularly, the collision operator in the LBM for the Nernst–Planck equation is designed such that all the terms necessary for the ion flux are evaluated in a completely local manner, and communication with neighboring nodes is only through the streaming process. Since the method inherits the algorithmic simplicity of the original LBM, it is compatible with the massively parallel computing. The Neumann-type boundary condition at the solid/liquid interface is also imposed correctly, with employing the standard bounce-back procedure. The remainder of the present paper is organized as follows. In Section 2, the governing equations for the electrical potential, ionic transports, and flows of electrolyte solutions are stated. The numerical scheme based on the LBM is presented in Section 3, including the detailed numerical procedure. In Section 4, two specific problems are considered and investigated numerically using the present LBM. An electro-osmotic flow through a channel of infinite length is considered in Section 4.1, and the LBM results are compared with the analytical solution. In Section 4.2, a comparison with the previous method is also made. A two-dimensional microchannel of finite length is then considered in Section 4.3. The electro-osmotic flow induced by an external electrical field is investigated, and the ends of the finite-length microchannel were found to cause resistance due to the local pressure gradient. 2. Electrokinetic model 2.1. Governing equations In this section, we state the model equations describing the behavior of an electrolyte solution. The scale of our interest ranges from several tens of nanometers to several micrometers, which is sufficiently larger than atomic scale. Therefore, we assume that the continuum descriptions of the transport phenomena are still valid. The flow of the electrolyte solution is then governed by the Navier–Stokes equation: ∂uj = 0, ∂xj ∂ui ∂ui 1 ∂p ∂ 2 ui Fi =− +ν 2 + , + uj ∂t ∂xj ρ0 ∂xi ∂xj ρ0

(1) (2)

where t is the time and x is the spatial coordinate. (In the present paper, we either use boldface letters or assign indexes i, j, and k to designate the 4

vector element. We assume the summation convention for repeated indexes.) The functions u(t, x) and p(t, x) are the flow velocity and the pressure of the electrolyte solution, respectively, and F (t, x) is the body force acting on the unit volume. The constants ρ0 and ν are the density and the kinetic viscosity of the electrolyte solution. We next show the conservation equation for the ion species in the electrolyte solution. In the following, the subscript m designates the quantities of the mth species, and m = a and c correspond to the anion and the cation, respectively. With Cm and Jm denoting the concentration and the flux, respectively, the mass conservation of the mth species is given as ∂Cm ∂Jmj = 0, + ∂t ∂xj

Jmi = −

ezm Dm ∂φ ∂Cm − Dm + Cm ui , Cm kB T ∂xi ∂xi

(3)

where e is the unit charge, kB is Boltzmann’s constant, and T is the temperature. The constants zm and Dm in Eq. (3) are the valence and the diffusion coefficient of the mth species, respectively. Here, the first, second, and third terms in the right-hand side of the flux equation are the contributions of the electrochemical migration, the diffusion, and the convection of the electrolyte solution, respectively. This model of the ion flux is referred to as the Nernst–Planck model [4]. The electrical potential φ obeys the following Poisson equation: ∂2φ (4) εr ε0 2 = −ρe , ∂xj where ε0 is the permittivity of a vacuum and εr is the dimensionless dielectric constant. The net charge density ρe is related to the ion concentration via  F zm Cm , (5) ρe = m

where F is Faraday’s constant (F = eNA with NA being Avogadro’s number). In the present paper, we assume that the body force acting on the electrolyte solution is due only to the interaction of the net charge and the electric field: Fi = −ρe

∂φ . ∂xi

(6)

2.2. Boundary conditions on the solid/liquid interface The boundary condition for the flow velocity at the solid/liquid interface is the ordinary non-slip condition: ui = 0,

x ∈ ∂Ω, 5

(7)

where ∂Ω denotes the solid/liquid boundary. Although a slip at the solid/liquid boundaries manifests itself in nano-scale flows [40, 41], we adopt the simple non-slip condition, because the scale of the problems considered in the present paper is relatively large. Extension to the slip boundary condition is left for future research. The pertinent boundary condition for the ion concentration is the non-flux condition across the boundary: nj Jj = 0,

x ∈ ∂Ω,

(8)

where n is the unit normal vector pointing inward to the fluid region Ω. Substituting Eq. (3) into Eq. (8) yields −

ezm Dm ∂φ ∂Cm − Dm n j + Cm nj uj = 0. Cm nj kB T ∂xj ∂xj

(9)

This form of the Neumann-type boundary condition is rather complex to implement if the electrochemical migration term is regarded as the source term [27]. Instead of using Eq. (9), a simplified boundary condition in the following form is often used [20, 24, 25]:   ezm φ , x ∈ ∂Ω, (10) Cm = C0m exp − kB T where C0m is the reference (or bulk) concentration. This condition is derived by integrating Eq. (3), assuming u = 0 and a steady state. Equation (10) implies that the reference potential is fixed to φ = 0, where the concentration is C0m . Obviously, it is not appropriate to assume this one-to-one correspondence between the potential and the concentration if the internal and external potential fields are combined. Different values of the reference potential should be assigned at each point on the boundary. Furthermore, the time-dependent behavior is not captured correctly by Eq. (10), because the Dirichlet-type condition (10) allows ion flux across the boundary, which is clearly artificial on the non-reactive solid/liquid interface. Therefore, we use the original boundary condition (9) in the present paper. Finally, we show the boundary conditions for the electrical potential. In the present paper, we assume that either (i) the surface charge density σ or (ii) the zeta potential ζ is given at the solid/liquid interface. If the surface charge distribution is not known a priori, the Poisson equation should be solved over the domain including the dielectric solid bodies and/or solid 6

bodies with internal charge distribution. In the case of the specified surface charge density, the Neumann-type boundary condition originating from Gauss’ law is imposed: −ε0 εr nj

∂φ = σ, ∂xj

x ∈ ∂Ω.

(11)

At the boundary of the specified zeta potential, the Dirichlet-type boundary condition φ = ζ is imposed. 3. Lattice Boltzmann method In this section, we present the LBM for solving the model equations for the electrokinetic flows. The set of equations consists of (i) the Poisson equation (4) with Eq. (5), (ii) the Nernst–Planck equation (3), and (iii) the Navier–Stokes equation (2) with Eqs. (1) and (6). The lattice Boltzmann (LB) equations are defined for each of these equations. In Section 3.1, the LB equation commonly used for equations (i) through (iii) is described, and detailed definitions specific to each transport process are given in Sections 3.2 through 3.4. Finally, the computational procedure is presented in Section 3.5. 3.1. Lattice Boltzmann equation The LB equation governs the behavior of the distribution function fα (t, x), where α = 0, 1, 2, . . ., N . Physical quantities such as the electrical potential and the ion concentration are obtained as moments of the distribution function, as described below. Each of fα is transported over a regular spatial lattice with the assigned velocity. The direction of the velocity is defined in terms of the vector eα . The definition of the vector eα depends on the form of the partial differential equation to be solved, and the explicit expression will be given in the following subsections. In the present paper, we use Greek subscripts to indicate the quantities corresponding to the directions of the discrete velocities, as fα above. It is convenient to regard the quantities with the Greek index as vectors in RN . Thus, we introduce the following row and column vector notation: |q := (q0 , q1 , . . . , qN )T ,

q| := (q0 , q1 , . . . , qN ),

7

(12)

The LB equation then reads 1 |f (t + ∆t, x + eα ∆x) − |f (t, x) = (|f eq  − |f )(t, x) τ   ∆t |G(t, x) + |G(t + ∆t, x + eα ∆x) , + 2

(13)

where |f (t + ∆t, x + eα ∆x) is the column vector having components fα (t + ∆t, x + eα ∆x), and ∆t and ∆x are the time step and the grid interval, respectively. The function Gα corresponds to the forcing or source term. The first term in the right-hand side is the collision term, which defines the relaxation process during a time step, with τ and |f eq  being the relaxation coefficient and the equilibrium distribution function. The collision term used in the present paper is the single-relaxation-time (SRT) method, in which all the distribution functions relax with the common relaxation-time coefficient τ irrespective of α. Although SRT suffers from several drawbacks, which are eliminated by using the multiple-relaxation-time (MRT) method, we prefer the simplicity of the SRT in the present paper. The extended LBM, which introduces the MRT collision operator, is straightforwardly obtained applying the technique described in Refs. [31, 42, 43]. The definitions of Gα and |f eq  will be given in the following subsections. In the implementation of the LBM, the calculation of Eq. (13) is split into the collision process and the streaming process: • Collision process: 1 ∆t |G(t, x), (14) |fˆ(t, x) = |f (t, x) + (|f eq  − |f )(t, x) + τ 2 • Streaming process: ∆t |f (t + ∆t, x + eα ∆x) = |fˆ(t, x) + |G(t + ∆t, x + eα ∆x). 2

(15)

3.2. Scheme for the Poisson equation 3.2.1. LB equation We first present the definition of the velocity vector eα . Here, and in what follows, we restrict ourselves to the two-dimensional problem. In solving the Poisson equation, we use the five-velocity model defined as follows:   0 1 −1 0 0 [e0 , e1 , e2 , e3 , e4 ] = (16) 0 0 0 1 −1 8

The equilibrium distribution function in the collision term has the following form: |f eq  = |ωφ,

φ = 1|f .

(17)

The weight coefficient ωα satisfies the following conditions: 1|ω = 1,

e|ω = 0,

ei ej |ω = Λδij .

(18)

Here, the raw vectors, such as 1|, are, according to the rule in Eq. (12), ex | = (0, 1, −1, 0, 0).

1| = (1, 1, 1, 1, 1),

(19)

In other words, they are the vectors for which the αth component is the quantity inside the brackets with the index α. If we choose, for instance, the weight coefficient to be  1/3, (α = 0) ωα = (20) 1/6, (α = 1, . . . , 4) then the coefficient of the tensor in Eq. (18) is Λ = 1/3. The vector Gα , which corresponds to the source term, is defined as |G = |ωρe .

(21)

If we set the value of the relaxation-time coefficient τ as τ=

1 ∆t εr ε0 , + 2 Λ∆x2

(22)

then φ obtained using the LBM converges to the solution of the following equation as ∆x → 0, keeping ∆t/∆x2 = const: ∂φ ∂2φ = ε r ε 0 2 + ρe . ∂t ∂xj We regard the solution in the limit of t → ∞ as the solution of Eq. (4).

9

(23)

3.2.2. Boundary condition Here, we explain how to implement the LBM on the boundary. If the boundary exists between x + eα ∆x ∈ Ω and x ∈ / Ω, it is not possible to carry out the streaming process (15). Therefore, we need to assign an appropriate value to fα . In the case of the Neumann-type boundary condition given by Eq. (11), we replace Eq. (15) by the following equation: fα (t + ∆t, x + eα ∆x) =fˆβ (t, x + eα ∆x) + ∆tσ/∆x ∆t + Gα (t + ∆t, x + eα ∆x). 2

(24)

Here, and in what follows, the index β indicates the direction opposite α, i.e., eα = −eβ . In the case of the Dirichlet-type boundary condition, i.e. the potential at the boundary is given as φ = ζ, we replace Eq. (15) by the following equation: fα (t + ∆t, x + eα ∆x) = − fˆβ (t, x + eα ∆x) + Λζ ∆t + Gα (t + ∆t, x + eα ∆x). 2

(25)

The boundary rules (24) and (25) properly reproduce the Neumann and Dirichlet boundary condition in the limit at which ∆x → 0, the proof of which follows from the results of the asymptotic analysis of boundary conditions in Ref. [31]. The periodic boundary condition with the potential bias, which is used in Section 4.1, is imposed based on the method described in Refs. [44, 45]. Specifically, the non-equilibrium part of the distribution function is translated periodically, while the equilibrium part is biased. The following equation describes the procedure in the x-direction, at the point where x is inside the domain but x ± ∆x is not: fα (t + ∆t, x ± ∆x ∓ Lx ) =fˆα (t, x) − ωα φ + ωα (φ ∓ ∆φ) ∆t + Gα (t + ∆t, x ± ∆x ∓ Lx ), 2

(26)

where the upper and lower signs apply to α = 1 and 2, respectively. Here, Lx is the length of the computational domain in the x-direction, and ∆φ is the potential bias.

10

3.2.3. Evaluation of derivative ∂φ/∂x Since the Navier–Stokes equation and the Nernst–Planck equation include the derivative of φ, we need to evaluate the value of ∂φ/∂x appropriately. The most commonly used method to evaluate the derivative in existing LBMs is the finite-difference approximation, in which the value of ∂φ/∂x is evaluated with the aid of the values of φ on the surrounding grid points [27, 46]. In the present study, however, we employ another method for evaluation of the derivative. The systematic asymptotic analysis detailed in Ref. [31] yields an expression of the distribution function expanded as a power series in the grid interval (∆x). Taking a product of the expansion with e| results in the formula: 1 ∂φ ei |f  + O(∆x2 ). =− (27) ∂xi τ Λ∆x Note the formula contains no finite-difference approximation, and thus the derivative is obtained in a completely local manner. As discussed in Ref. [31], the second-order accuracy is valid as long as ∆t/∆x2 is kept at constant. The numerical validation of the second-order accuracy has been made by Li et al., in the course of their study of the boundary condition for the thermal LBM [47]. The remarkable advantage of using this local evaluation (27) is that no special treatment is necessary near the boundary. Furthermore, compatibility with parallel computing is directly inherited from the original LBM. 3.3. Scheme for the Nernst–Planck equation 3.3.1. LB equation In solving the Nernst–Planck equation, we use the set of eα and ωα defined in the previous subsection. The flux Jm of the Nernst–Planck equation described in Eq. (3) consists of the electrochemical migration, diffusion, and convection terms. If we regard the coefficient of Cm in the migration term as a part of the convection velocity of the convection-diffusion equation, the equilibrium distribution function of the LBM for the general convection-diffusion equation described in Ref. [31] is transformed into:    ∆t ezm Dm ∂φ eq (28) |ej ω Cm , uj − |f  = |ω + ∆xΛ kB T ∂xj Cm = 1|f . (29)

11

The relation between the relaxation coefficient τ and the diffusion coefficient is 1 ∆t Dm . (30) τ= + 2 Λ∆x2 Since the electrochemical migration term is integrated into the equilibrium distribution function, there is no source term, and Gα = 0 in Eq. (13). Therefore, combined with the local evaluation of ∂φ/∂x in Eq. (27), the present scheme for the Nernst–Planck equation is also implemented locally, and the communication with the surrounding grid points is only through the streaming process (Eq. (15)). One of the advantages of this locality is that the original boundary condition for ion flux is easily implemented in a natural manner using the standard boundary rules, and the artificial flux in time-dependent problems is properly eliminated (Section 4.2). In parallel with the derivation of the convection-diffusion equation from the MRT-LBM detailed in Ref. [31], the numerical solution obtained using the present LB algorithm based on Eq. (28) is shown to converge to the solution of the Nernst–Planck equation with second-order accuracy with respect to the grid interval. 3.3.2. Boundary condition We consider the solid/liquid boundary at the midpoint between x + / Ω. Recall that, in the scheme for the Nernst–Planck eα ∆x ∈ Ω and x ∈ equation, all the contributions to the ion flux are integrated in the collision operator and thus there is no source term. Thanks to this feature, the nonflux boundary condition (9) of the Neumann type is easily imposed using the standard bounce-back procedure, which ensures no mass transport across the boundary: (31) fα (t + ∆t, x + eα ∆x) = fˆβ (t, x + eα ∆x). The Dirichlet boundary condition used in Section 4.3, i.e., the value-fixed boundary condition (Cm = C0 ), is imposed using the following equation: fα (t + ∆t, x + eα ∆x) = −fˆβ (t, x + eα ∆x) + ΛC0 .

(32)

The asymptotic analysis of the boundary conditions detailed in Ref. [31] is straightforwardly applied to prove that Eqs. (31) and (32) properly recover the Neumann and Dirichlet boundary conditions.

12

3.4. Scheme for the Navier–Stokes equation The set of the velocity vector eα used in solving the Navier–Stokes equation is different from that used in solving the Poisson and Nernst–Planck equations. Here, we use the nine-velocity model defined as   0 1 −1 0 0 1 −1 1 −1 [e0 , e1 , e2 , e3 , e4 , e5 , e6 , e7 , e8 ] = 0 0 0 1 −1 1 1 −1 −1 (33) The equilibrium distribution function is defined as   3 9 3 eq 2 2 |f  = ρ |ω + |ωej uj  + |ω(ej uj )  − |ωuj  , (34) Cu 2Cu2 2Cu2 ρ = 1|f , ρui = Cu ei |f . (35) Here, Cu = ∆x/∆t, and the weight coefficient ωα is defined as   4/9, (α = 0) 1/9, (α = 1, . . . , 4) ωα =  1/36. (α = 5, . . . , 8)

(36)

The body force is incorporated via

|G =

3Fj |ωej  , Cu2

(37)

where F is defined in Eq. (6). Note that the value of ∂φ/∂x, which is necessary in computing F using Eq. (6), is evaluated by means of Eq. (27). If we set the value of the relaxation-time coefficient as τ=

1 3∆t ν, + 2 ∆x2

(38)

then the flow velocity u and the pressure p ≡ Cu2 ρ/3 obtained by using the LBM converge to the solution of the Navier–Stokes equation in the limit at which ∆x → 0 [28, 29]. Strictly speaking, the numerical solution satisfies the Navier–Stokes equation with an artificial compressibility. However, since the variation of ρ (or p) is of higher order than ∆x and so the divergence of u is small, we apply this method to obtain the flow field. In order to impose the non-slip boundary condition (7), we use the simple bounce-back procedure in the same form as Eq. (31). 13

3.5. Computational procedure The electrokinetic model considered herein is composed of the Poisson equation, which is independent of time, and the Navier–Stokes and Nernst– Planck equations, which are dependent on time. Since the numerical solution of φ is obtained as the long-time limit of the time-dependent solution, we introduce an artificial time axis besides t, denoted by t˜. The ion concentration (anion Ca and cation Cc ) and the flow velocity u of the electrolyte solution are obtained using the common axis t. However, since the timescales of each transport phenomenon are different, we need to assign different values to the time step ∆t. In the sequel, we denote the time steps for the anion, cation, and flow velocity by ∆ta , ∆tc , and ∆tu , respectively. Normally, ∆tu ∆ta , ∆tc because the kinetic viscosity is much larger than the diffusion coefficients of ion species. The iterative procedure is illustrated in Fig. 1 for the case of Nt ∆tu = ∆ta = ∆tc . The solution φ at each instance in t is obtained as the limit at which t˜ → ∞ for the t˜-dependent solution. In the actual implementation, the iteration with respect to t˜ is terminated when the difference between two successive values of φ reaches a certain tolerance, typically 10−7 [V], at all of the grid points in Ω. Although Fig. 1 shows the case of ∆ta = ∆tc , the procedure is readily extended to the case of ∆ta = ∆tc . In that case, the iteration with respect to the species with the smaller ∆t is repeated until the sum of ∆t reaches the larger time step. 4. Application to electro-osmotic flows in microchannels In this section, we apply the proposed LBM to specific problems. We first consider a simple one-dimensional problem of electro-osmotic flows between two parallel plates. The numerical results will be compared with the analytical solution in order to validate the present algorithm. We next consider a two-dimensional microchannel of finite length with an entry and an exit. It will be shown that the presence of the entry and exit results in the occurrence of the pressure gradient that causes as a flow resistance reducing the magnitude of the electro-osmotic flow. 4.1. Infinite-length channel We first consider the problem of the electro-osmotic flow in a channel of infinite length. More specifically, the slit microchannel formed between two parallel plates located at y = ±H is filled with a 1 : 1 electrolyte solution 14

(C0a = C0c ≡ C0 ). The surfaces of the channel are charged with a constant charge density σ, and thus the EDL is formed near the channel walls. If a potential bias is imposed in the x-direction, i.e., an external electric field is applied along the channel, Coulomb’s force acts on the electrolyte solution through the net charge within the EDL, and consequently an electro-osmotic flow occurs in the x-direction. In this subsection, the one-dimensional electroosmotic flow in steady state is investigated numerically using the proposed LBM algorithm. If the surface charge density is small and thus the potential variation is small, the problem is linearized, and an analytical solution can be derived. In the present subsection, we compare the results of the LBM computations with the analytical solution in order to validate the LBM algorithm. Furthermore, the nonlinear effect due to the large surface charge density is later discussed. Although the derivation of the analytical solution for similar problems has been described in many places (e.g., Refs. [6, 8]), it is worth outlining the derivation of the analytical solution in order to specify the analytical solution fitting the boundary condition considered herein. This also clarifies the assumptions of the analytical solution. Since the flow is assumed to be steady and one-dimensional in the ydirection, the Navier–Stokes equation combined with the Poisson equation reduces to −Px + ρ0 ν

d 2 ux d2 φ − E ε ε = 0, x r 0 dy 2 dy 2

(39)

where Px is the pressure gradient along the channel, and Ex is the applied electrical field strength. Integrating Eq. (39) leads to the following expression for the flow field:   Px 2 Ex εr ε0 ζ φ(y) 2 (y − H ) + −1 , (40) ux (y) = 2ρ0 ν ρ0 ν ζ where ζ is the potential at the wall (the zeta potential, ζ = φ(H)). Here, the non-slip boundary condition is also used to determine the integration constant. Equation (40) implies that the flow field is immediately obtained if the potential distribution is known. Under the assumptions of a steady-state and one-dimensional flow, the PB approximation is valid. Substituting the Boltzmann distribution Cm = C0 exp(−ezm φ/kB T ) into Eq. (5) transforms Eq. (4) into the nonlinear PB 15

equation:

2C0 F d2 φ = sinh 2 dy εr ε0



eφ kB T



.

(41)

Although there exist numerous methods by which to solve the nonlinear PB equation numerically (e.g., Refs. [35–38]), the PB equation cannot be solved analytically. In order to find an analytical solution, the Debye–H¨ uckel linear approximation, i.e., |eφ/kB T | 1, is used. Then, the right-hand side of the PB equation is approximated by κ2 φ, where κ is the inverse of the Debye length defined as  1/2 2C0 eF κ= . (42) εr ε0 kB T With the aid of the boundary condition of the surface charge density, i.e., dφ/dy = ±σ/εr ε0 at y = ±H, the linearized PB equation is analytically solved: σ cosh(κy) . (43) φ(y) = εr ε0 κ sinh(κH) Substituting Eq. (43) into Eq. (40) gives the following analytical expression for the flow field.   Px 2 Ex σ cosh(κy) − cosh(κH) 2 ux (y) = . (44) (y − H ) + 2ρ0 ν ρ0 νκ sinh(κH) The formula for the flow velocity averaged over the channel cross-section is obtained by integrating Eq. (44):  

1 H Px H 2 κH Ex σ av ux (y)dy = − 1− . (45) + ux = H 0 3ρ0 ν ρ0 νκ2 H tanh(κH) The following parameters are used in the numerical analysis. The parameters characterizing the electrolyte solution are ρ0 = 103 [kg/m3 ], ν = 0.889×10−6 [m2 /s], Da = Dc = 10−8 [m2 /s], and εr ε0 = 6.95×10−10 [C2 /Jm]. In addition, the temperature is T = 273 [K], and the channel height is fixed at H = 0.5 [µm]. The channel width of 2H = 1[µm], which is small compared with that in typical microfluidic devices, is chosen in the present paper to demonstrate the applicability to the cases in which the thickness of the EDL is important [10, 26, 48–51]. The strength of the electrical field along the channel is Ex = 104 [V/m], and we restrict ourselves to the case of no pressure gradient along the channel (Px = 0). The parameters necessary for the 16

LBM computations are ∆x = 0.01 [µm], and 100 × ∆tu = ∆ta = ∆tc = 0.001 [µs], which are chosen based on a careful examination of accuracy. The rectangular domain with length 2 [µm] is used in the LBM computation, and the periodic boundary condition is applied in the x-direction. The potential bias is imposed using the boundary rule based on the method described in Ref. [44]. In Fig. 2, the flow velocity averaged over the channel cross-section is shown as a function of the ion concentration C0 . The cases of σ = −2 × 10−5 , −10 × 10−5 , −20 × 10−5 , and −50 × 10−5 [C/m2 ] are shown, and the solid line in the figure indicates the analytical formula (45). The average flow velocity increases as the concentration decreases while keeping the surface charge density at constant. This is because the zeta potential ζ increases with decreasing concentration. In order to illustrate this, in Fig. 3, we show the zeta potential as a function of the concentration. The figure clearly shows that the zeta potential is large for low concentrations. We briefly mention the evaluation of the zeta potential. The zeta potential is defined only for the equilibrium state without the external electrical field. Under the assumptions made in deriving Eq. (45), since the external electrical field is separated form the PB equation (41), the analytical form of the zeta potential in is simply obtained from Eq. (43) with y = H. On the other hand, since the external electrical field is integrated into the Poisson equation in the LBM computation, it is not possible to evaluate the zeta potential directly from the potential field obtained numerically. Therefore, we estimate the zeta potential using the following relation, assuming a Boltzmann distribution for the ion concentration: 1/2   1 kB T  , ξ= ρe (y = H). (46) ln −ξ + ξ 2 + 1 ζ= e 2F C0 This relation is derived from Eq. (5) with Cm = C0 exp(−ezm φ/kB T ). In deriving the formula (45), the potential variation is assumed to be small, i.e., the PB equation is linearized assuming |eφ/kB T | 1. Therefore, the fact that the zeta potential is sufficiently small compared to the value of |ζ| kB T /e ∼ 23.5 = ζ∗ [mV] is an indication of the validity of the analytical formula. In Fig. 3, the zeta potential at C0 = 0.008, 0.016 [mol/m3 ] for the case of σ = −50 × 10−5 [C/m2 ] greatly exceeds the value of ζ∗ , and correspondingly the values of the average flow velocity exhibit a clear discrepancy from the analytical values because the Debye–H¨ uckel approximation is not valid. 17

In order to realize a more detailed comparison with the analytical formula, in Figs. 4 and 5, we show the profiles of the flow velocity and the ion concentration for the cases of σ = −2 × 10−5 and −50 × 10−5 [C/m2 ]. The profile for C0 = 0.002 [mol/m3 ] is not shown in Figs. 4(b) and 5(b). This is because the concentration is too low to compensate for the charge on the boundary, and no solution with a positive anion concentration (Ca > 0) is obtained. In the case of σ = −2 × 10−5 [C/m2 ], the analytical formula exhibits very good agreement because the assumption for the linearization is valid. On the other hand, the analytical formula overestimates the flow velocity for the cases of C0 = 0.008 and C0 = 0.032 [mol/m3 ], as discussed in the previous paragraph. Finally, we confirm that the present method is also valid for large zeta potentials, at which the Debye–H¨ uckel approximation does not apply. To this end, we consider the case of the electrolyte solution without an added salt. More precisely, only the cation is present in the electrolyte solution, which compensates the negative charge on the channel surface, and Ca = 0. In this case, the Nernst–Planck equation combined with the Poisson equation is solved analytically without the Debye–H¨ uckel approximation (e.g., Refs. [48, 52]): Cc (y) = φ(y) =

C∗ , 2 cos (κ∗ y/2)

2kB T ln cos(κ∗ y/2), e

(47) (48)

where κ∗ is defined by Eq. (42) with C0 replaced by C∗ . The value of C∗ indicates the concentration at y = 0, which is determined through the relation between the net charge in the electrolyte solution and that on the surface: (σe/εr ε0 kB T ) = −κ∗ tan(κ∗ H/2). The flow field without the pressure gradient is obtained by substituting the potential field given in Eq. (48) into Eq. (40):   2Ex εr ε0 kB T cos(κ∗ y/2) ux (y) = . (49) ln ρ0 νe cos(κ∗ H/2) Figure 6 shows the profiles of the flow velocity and the concentration of cation, obtained using the present LBM computation along with the analytical formulae (47) and (49). The absolute values of the zeta potential (Eq. (48) with y = H) are 47.5, 80.9, and 110 [mV], in the cases of σ = −20 × 10−5 , −50 × 10−5 , and −100 × 10−5 [C/m2 ], respectively, which are much larger 18

than ζ∗ . The excellent agreement exhibited in Fig. 6 thus confirms the validity of the present method for the parameter range in which the Debye–H¨ uckel approximation is invalid. 4.2. Channel with inhomogeneous zeta potential One powerful feature of the governing system based on the Nernst–Planck equation is that the time evolution of the flow and ion concentration is tracked. In this subsection, we show that the present model properly captures the transient behavior of the ion concentration. In order to demonstrate this point with a simple geometry, we consider an infinite channel with the nonuniformly distributed wall potential. More specifically, instead of assigning a surface charge density σ, we apply the potential distribution:   −ζ (0 < x < L/4), ζ (L/2 < x < 3L/4), φ(y = ±H) = (50)  0 (otherwise).

Here, we assume periodic boundary condition in the x-direction with the period L, and no external electric field. The inhomogeneous distribution of the zeta potential is chosen so that the effect of the counter electrode is incorporated in the closed domain with periodic system. The problem herein is preliminary to the investigation of AC electro-osmotic pumps with microelectrodes [53–58], in which the transient behavior during formation of the EDLs is important. In the LBM computation, we set L = 1 [µm] and C0 = 0.032 [mol/m3 ], and other parameters are the same as those in the previous subsection. Although the gap of the channel (H = 0.5 [µm]) considered herein is smaller than that in typical experimental setups for the AC electro-osmotic pumps, the important feature demonstrated using the problem, i.e., the correct imposition of the boundary conditions during formation of the EDL, is independent of the choice of the parameters. As mentioned in Introduction, the LBM for the Nernst–Plank equation proposed by Wang and Kang in Ref. [27] (WK method) is not sufficient when we consider the time-dependent behavior, because of the boundary treatment based on the Poisson–Boltzmann approximation (Eq. (19) of Ref. [27]). We therefore compare the simulation results with those obtained using the WK method, to show the difference clearly. Simulation starts with the initial ion concentration that is uniform at C0 . The potential ζ = 5 [mV] is turned on at t = 0. We plot in Fig. 7 the time 19

evolution of the ion concentration Cm at the point (x, y) = (L/8, H) on the channel surface. Whereas the present method does capture the migration process before the EDL is fully developed, the method in Ref. [27] skips the migration process, applying concentration obtained by the Boltzmann distribution. In order to see the transient behavior in greater details, in Fig. 8, we show the profiles of the potential φ and the net charge ρe , as well as the concentrations Cc and Ca , along the line at x = L/8 with y ∈ (0, H). The concentration profiles of cation and anion (panels (a) and (b)) obtained using the present method gradually develop to form the EDL. In contrast, the corresponding results of the WK method at short times show the very sharp gradients due to the surface concentration of the Boltzmann distribution. The behavior of the net charge (panel (d)) is almost identical to that of cation concentration. The relaxation process of the concentration profiles to the steady state of the WK method is faster than the that of the present method, which is also observed in the potential profiles (panel (c)). The contour plots in Fig. 9 show the difference in ρe (see Eq.(5)) at t = 3 × 10−7 [s] and 5 × 10−5 [s]. Although the two methods give almost the same results after a long-time has passed (panel (b)), the short-time behavior, which is important, e.g., in problems with AC voltage[54, 55, 58], is significantly different (panel (a)). The discrepancy in the time-dependent analysis observed above is caused by the approximated boundary condition for the Nernst–Planck equation described in Eq. (10), which allows the artificial flux across the boundary. In order to examine the flux, we show in Fig. 10 the current field defined as I = F (Jc − Ja ) (with Jm defined in Eq. (3)) at t = 3 × 10−7 [s]. The result of the present method shows the charge transport from the region near the boundary with positive zeta potential to that with negative zeta potential. On the other hand, the result of the WK method shows the unphysical current penetrating the boundary. Figure 11 compares the distribution of the current in the y-direction at various time instants. In panels (a) and (b), distributions along the lines at x = L/8 and x = 5L/8, respectively, with y ∈ (0, H) are shown. Before the turn-on of the zeta potential (t = 0), the flux is zero everywhere. Just after the turn-on of the zeta potential, at x = L/8, the positive current occurs toward the boundary, and the current in the opposite direction occurs at x = 5L/8. The generated current decays in the longtime regime and finally vanishes in the steady state. Clearly, the artificial currents (or ion flux) at y = H are observed at short times in Figs. 11(a) and (b) for the results of WK method. Particularly, at very shot times (≤ 10−8 20

[s]) the absolute value of the current at the boundary is extremely large. On the other hand, the present method properly imposes the no-flux condition (Eq. (8)) at all time instances. As explained in Section 3.3, the localized scheme for the Nernst–Planck model employing the equilibrium distribution function described in Eq. (28) realizes the no-flux boundary condition by means of the simple implementation of the standard bounce-back procedure. 4.3. Finite-length channel We next consider the two-dimensional problem of a finite-length channel with an entry and exit. The geometry of the problem is shown in Fig. 12. The length of the channel is L, and entry and exit regions of length LM are located next to the ends of the channel. The symmetric boundary condition in the y-direction means that the present problem is one unit of an array of channels with interval 2HM . The boundary conditions at x = 0 are φ = 0, p = p0 , and Cm = C0 , and those at x = 2LM + L are φ = ∆φ, p = p0 , and Cm = C0 . The boundary condition on the channel walls, i.e., LM < x < LM + L at y = ±H, is constant charge σ = 0 for the potential, no flux for the ion, and no slip for the flow. The boundary condition on the face, i.e., |y| > H at x = LM and at x = LM + L is the same as that on the channel walls, except that σ = 0. The values of the parameters characterizing the electrolyte solution are the same as those of Section 4.1, and the temperature is also fixed at T = 273 [K]. Since the problem is symmetric about y = 0, the upper half region is considered in the computation imposing the symmetric boundary condition, i.e., the derivative with respect to y is zero at y = 0. In Figs. 13 and 14, we show typical plots of the electrical potential φ, the net charge ρe , the pressure p−p0 , and the flow velocity vector u, in the steady state. The geometrical parameters are H = HM = 0.5 [µm], L = 4 [µm], and LM = 1 [µm]. The potential difference is ∆φ = −60 [mV] and the surface charge density is σ = −10−4 [C/m2 ]. The ion concentration is C0 = 0.032 [mol/m3 ] in Fig. 13, and C0 = 0.002 [mol/m3 ] in Fig. 14. In the case of a high concentration (Fig. 13), i.e., the Debye length is short (see Eq. (42)), the surface charge on the channel wall is shielded by a thin EDL formed adjacent to the channel wall. Therefore, the electrical potential and the net charge do not vary in the y-direction around the center of the channel (y = 0). On the other hand, in the case of a low concentration (Fig. 14), the Debye length is long and the EDL is much thicker. Therefore, the net charge spans over the entire width of the channel. Consequently, the shape of the flow velocity 21

field is similar to a Poiseuille flow because of the driving force over the entire channel width, whereas, in Fig. 13, the flow velocity is plug-like due to the driving force localized near the channel wall. Moreover, in Fig. 14, the twodimensionality near the entry and exit of the channel, which is observed in the plots of the electrical potential and the net charge, is stronger than that observed in Fig. 13. In both Figs. 13 and 14, the pressure variation occurs in the x-direction, which counteracts the electro-osmotic flow, although no external pressure difference is imposed (p = p0 at x = 0 and x = 2LM + L). The pressure gradient in the x-direction will be further investigated later in this subsection. We next investigate the transient behavior before the steady state is established. The initial condition for the transient analysis is the equilibrium state in which ∆φ = 0. More precisely, we first carry out LBM computation while setting ∆φ = 0 until the electrical potential field converges. Then, using the obtained equilibrium state as the initial condition at t = 0, the LBM computation with ∆φ = 0 is performed in order to track the transient behavior of t > 0. Figure 15 shows the time evolution of the flow velocity averaged over the cross-section at the middle of the channel (x = LM + L/2). The results obtained for several ion concentrations are shown. The steady state is established approximately 10−6 [s] after the electrical filed (or ∆φ) is turned on. This coincides with the time scale for momentum dissipation over the channel width, i.e., (2H)2 /ν ∼ 10−6 [s]. As is clear for the case of C0 = 0.128 [mol/m3 ], an oscillation takes place, which results from the sudden turn-on of the driving force localized near the channel wall. Although not clear from the figure, in the case of C0 = 0.002 [mol/m3 ], a very slow oscillation still exists at t = 10−4 [s]. In this subsection, however, we regard the solution at t = 10−4 [s] as the steady-state solution because the amplitude is negligible (less than 1% for the average flow velocity). Figures 16 and 17 show the time evolution of the flow velocity profile at the entry (x = LM ), at the middle (x = LM + L/2), and at the exit (x = LM + L). At the exit of the channel for C0 = 0.032 [mol/m3 ] (Fig. 16 (c)), the fluid near the channel wall (y = H) is pushed out by the force acting in the EDL. Instead, drawing of the fluid, i.e., the negative flow velocity, occurs around the center of the channel (y = 0). On the contrary, at the entry of the channel (Fig. 16 (a)), the force drawing the fluid into the channel is relatively strong, and thus a positive flow is generated along the entire channel width. In the case of a low concentration of C0 = 0.002 [mol/m3 ] (Fig. 17), since the EDL is not localized near the channel wall, the driving force acts in a wider region, including 22

the channel center. Therefore, the negative flow velocity does not occur even at the exit of the channel (Fig. 17(c)). We further investigate the flow velocity of the steady state. Figure 18 shows the flow velocity averaged over the cross-section at the middle of the channel (x = LM + L/2) as a function of the ion concentration. The results for different values of HM are shown. In addition, for comparison, we plot the prediction using the formula for the infinite-length channel (45) setting Px = 0 and ∆φ Ex = . (51) L + 2LM H/(HM + H) Equation (51) is derived assuming that the channel is sufficiently long compared to the width (and thus the electric field is one dimensional), and is valid for the case in which the effect of the entry and exit is negligible. Clearly, the analytical formula (45) for the infinite-length channel overestimates the average flow velocity. It is also observed the dependence on HM is suppressed in the LBM results. These differences result from two factors. One is that the estimation of the electrical field using Eq. (51) is not sufficiently accurate because of the effect of the entry and exit. The other is that, as mentioned above, the pressure gradient created in the two-dimensional problem causes a driving force in the direction opposite the electro-osmotic flow. In order to elucidate these factors, in Fig. 19, we show the electrical potential and the pressure distribution along the centerline (y = 0) for the case of HM = H = 0.5 [µm]. In Fig. 19(a), the electrical potential distribution obtained using Eq. (51) (and the corresponding equations for the entry and exit regions) is also indicated by the dashed line. Although Eq. (51) gives a slightly steep gradient, i.e., the electrical field is overestimated, the concentration dependence of the degree of the overestimation is not explained by this discrepancy, because the potential gradients for C0 = 0.008, 0.032, and 0.128 [mol/m2 ] are almost identical. In contrast, Fig. 19(b) reveals a clear dependence of the pressure gradient on the concentration. The pressure gradient is larger for lower concentrations. Therefore, we conclude that the differences in Fig. 18 are due primarily to the exclusion of the pressure gradient in the analytical formula (45). Finally, Fig. 20 shows the flow velocity profiles at the middle of the channel (x = LM + L/2) in the steady state for the case of HM = H = 0.5 [µm]. The dashed line indicates the profiles expressed by Eq. (44) corresponding to the analytical formula shown in Fig. 18, whereas the solid line indicates the 23

profiles expressed by the same equation (44) but with the values of Ex and Px computed by taking the derivative of the potential and pressure distributions of Fig. 19. Equation (44) is capable of capturing the results of numerical analysis if the values of the potential and pressure gradient are appropriate. In the case of C0 = 0.002 [mol/m3 ], however, the analytical prediction still shows a discrepancy. This is because, as discussed in Section 4.1, the variation of the electrical potential (or the magnitude of the zeta potential) is so large that the nonlinear effect is not negligible, and the assumption regarding the derivation of the analytical formulae (Eqs. (44) and (45)) is invalid. Therefore, in order to precisely investigate the electro-osmotic flow with a large variation in electrical potential, the numerical analysis of the nonlinear model equation described in Section 2 is necessary. 5. Summary In the present paper, a framework of numerical analysis for the electrokinetic model describing micro-/nano-flows has been proposed. The Nernst– Planck equation is adopted in order to incorporate the electrochemical migration and the convection of the electrolyte solution, in addition to the diffusion of the ion species. The governing equations are solved using the LBM. In the scheme for solving the Nernst–Planck equation described in Section 3.3, all the contributions to the ion flux are included in the collision process that is implemented locally. In addition, since the gradient of the electrical potential, which appears in the migration term of the Nernst–Planck equation and in the body-force term of the Navier–Stokes equation, is evaluated completely locally at each grid point, the coupling algorithm never violates the locality of the collision process of the LBM. Therefore, simplicity in implementing the boundary condition and compatibility with parallel computing are directly inherited from the original LBM. The present LBM was applied to two specific problems, namely, the onedimensional problem of the electro-osmotic flow in the infinite-length channel and the two-dimensional problem consisting of a finite-length channel with an entry and exit. The results of the LBM computation were compared to the analytical solution derived assuming the linearized Poisson–Boltzmann equation, as well as to the analytical solution of the Nernst–Planck equation without linearization, and the method was validated for the cases of small and large zeta potentials. A comparison with the existing method was also carried out, to show that the artificial flux across the boundary was properly 24

eliminated by using the present method. For the case of the finite-length channel, the presence of the entry and exit caused a local pressure gradient along the channel that reduced the magnitude of the electro-osmotic flow, even if no global pressure gradient was imposed. Since the set of model equations is general, i.e., no approximation of linearization, and the LBM is capable of handling complicated geometries, the numerical method in the present paper has potential applications for wide range of electrokinetic flows. Therefore, the next stage of our study will be engineering applications that include more complicated geometries, such as the electro-osmotic flow through a channel of varying cross-section and the electro-osmotic flow and streaming potential in porous media. Acknowledgements The present work was partially supported by MEXT program “Elements Strategy Initiative to Form Core Research Center” (since 2012). (MEXT stands for Ministry of Education Culture, Sports, Science and Technology, Japan.) References [1] J. Israelachvili, Intermolecular and surface forces 3rd Edition, Academic press, 2011. [2] B. Honig, A. Nicholls, Classical electrostatics in biology and chemistry, Science 268 (1995) 1144–1149. [3] H. Washizu, K. Kikuchi, Electric polarizability of DNA in aqueous salt solution, J. Phys. Chem. B 110 (2006) 2855–2861. [4] J. Newman, Electrochemical Systems, Prentice-Hall, Englewood Cliffs, NJ, 1991. [5] H. Li, R. J. Gale, Hydraulic and electroosmotic flow through silica capillaries, Langmuir 9 (1993) 1150–1155. [6] D. Li, Electrokinetics in Microfluidics, Volume 2, Interface Science and Technology, Elsevier, 2004.

25

[7] J. C. T. Eijkel, A. Berg, Nanofluidics: what is it and what can we expect from it?, Microfluid. Nanofluid. 1 (2005) 249–267. [8] J. H. Masliyah, S. Bhattacharjee, Electrokinetic and Colloid Transport Phenomena, John-Wiley & Sons, Inc., Hoboken, New Jersey, 2006. [9] R. B. Schoch, J. Han, P. Renaud, Transport phenomena in nanofluidics, Rev. Mod. Phys. 80 (2008) 839. [10] L. Bocquet, E. Charlaix, Nanofluidics, from bulk to interfaces, Chem. Soc. Rev. 39 (2010) 1073–1095. [11] H. Daiguji, Nanofluidics: High mobility in tight spaces., Nature Nanotechnology 5 (2010) 831. [12] L. Wang, J. Fan, Nanofluids research: Key issues, Nanoscale Res. Lett. 5 (2010) 1241–1252. [13] F. H. J. van der Heyden, D. J. Bonthuis, D. Stein, C. Meyer, C. Dekker, Power generation by pressure-driven transport of ions in nanofluidic channels, Nano Lett. 7 (2007) 1022–1025. [14] C. Davidson, X. Xuan, Electrokinetic energy conversion in slip nanochannels, J. Power Sources 179 (2008) 297–300. [15] A. Achilli, T. Y. Cath, A. E. Childress, Power generation with pressure retarded osmosis: An experimental and theoretical investigation, J. Membr. Sci. 343 (2009) 42–52. [16] E. Kjeang, N. Djilali, D. Sinton, Microfluidic fuel cells: A review, J. Power Sources 186 (2009) 353–369. [17] R. Chein, H. Chen, C. Liao, Analysis of electro-kinetic pumping efficiency through finite-length nano-scale surface-charged capillaries, J. Electroanal. Chem. 630 (2009) 1–9. [18] J. M. Perry, K. Zhou, Z. D. Harms, S. C. Jacobson, Ion transport in nanofluidic funnels, ACS Nanoano 4 (2010) 3897–3902. [19] A. Siria, P. Poncharal, A.-L. Biance, R. Fulcrand, X. Blase, S. T. Purcell, L. Bocquet, Giant osmotic energy conversion measured in a single transmembrane boron nitride nanotube, Nature 494 (2013) 455–458. 26

[20] R. J. Yang, L. M. Fu, C. C. Hwang, Electroosmotic entry flow in a microchannel, J. Colloid Interf. Sci. 244 (2001) 173–179. [21] F. Capuani, I. Pagonabarraga, D. Frenkel, Discrete solution of the electrokinetic equations, J. Chem. Phys. 121 (2004) 973–986. [22] I. Pagonabarraga, F. Capuani, D. Frenkel, Mesoscopic lattice modeling of electrokinetic phenomena, Comput. Phys. Commun. 169 (2005) 192– 196. [23] A. Mansouri, C. Scheuerman, S. Bhattacharjee, D. Y. Kwok, L. W. Kostiuk, Transient streaming potential in a finite length microchannel, J. Colloid Interf. Sci. 292 (2005) 567–580. [24] E. Y. K. Ng, S. T. Tan, Study of EDL effect on 3-D developing flow in microchannel with Poisson–Boltzmann and Nernst–Planck models, J. Numer. Meth. Eng. 71 (2007) 818–836. [25] H. M. Park, J. S. Lee, T. W. Kim, Comparison of the Nernst–Planck model and the Poisson–Boltzmann model for electroosmotic flows in microchannels, J. Colloid Interf. Sci. 315 (2007) 731–739. [26] Y. S. Choi, S. J. Kim, Electrokinetic flow-induced currents in silica nanofluidic channels, J. Collid Interf. Sci. 333 (2009) 672–678. [27] M. Wang, Q. Kang, Modeling electrokinetic flows in microchannels using coupled lattice Boltzmann methods, J. Comput. Phys. 229 (2010) 728– 744. [28] S. Chen, G. D. Doolen, Lattice Boltzmann method for fluid flows, Annu. Rev. Fluid Mech. 30 (1998) 329–364. [29] S. Succi, The lattice Boltzmann equation for fluid dynamics and beyond, Oxford Univ. Press, New York, 2001. [30] I. Ginzburg, Equilibrium-type and link-type lattice Boltzmann models for generic advection and anisotropic-dispersion equation, Adv. Water Resour. 28 (2005) 1171–1195. [31] H. Yoshida, M. Nagaoka, Multiple-relaxation-time lattice Boltzmann model for the convection and anisotropic diffusion equation, J. Comput. Phys. 229 (2010) 7774–7795. 27

[32] P. Asinari, M. C. Quaglia, M. R. von Spakovsky, B. V. Kasula, Direct numerical calculation of the kinematic tortuosity of reactive mixture flow in the anode layer of solid oxide fuel cells by the lattice Boltzmann method, J. Power Sources 170 (2007) 359–375. [33] Y. Suzue, N. Shikazono, N. Kasagi, Micro modeling of solid oxide fuel cell anode based on stochastic reconstruction, J. Power Sources 184 (2008) 52–59. [34] K. Suga, S. Takenaka, T. Ito, M. Kaneda, T. Kinjo, S. Hyodo, Evaluation of a lattice Boltzmann method in a complex nanoflow, Phys. Rev. E 82 (2010) 016701. [35] Z. Guo, T. S. Zhao, Y. Shi, A lattice Boltzmann algorithm for electroosmotic flows in microfluidic devices, J. Chem. Phys. 122 (2005) 144907. [36] Z. Chai, B. Shi, Simulation of electro-osmotic flow in microchannel with lattice Boltzmann method, Phys. Lett. A 364 (2007) 183–188. [37] M. Wang, J. Wang, S. Chen, Roughness and cavitations effects on electro-osmotic flows in rough microchannels using the lattice PoissonBoltzmann methods, J. Comput. Phys. 226 (2007) 836–851. [38] J. Wang, M. Wang, Z. Li, Lattice evolution solution for the nonlinear Poisson-Boltzmann equation in confined domains, Commun. Nonlinear Sci. Numer. Simulat. 13 (2008) 575–583. [39] J. Zhang, Lattice Boltzmann method for microfluidics: models and applications, Microfluid. Nanofluid. 10 (2011) 1–28. [40] H. Washizu, T. Ohmori, Molecular dynamics simulations of elastohydrodynamic lubrication oil film, Lubr. Sci. 22 (2010) 323–340. [41] A. E. Kobryn, A. Kovalenko, Slip boundary conditions in nanofluidics from the molecular theory of solvation, Molecular Simulat. 37 (2011) 733–737. [42] P. Lallemand, L.-S. Luo, Theory of the lattice Boltzmann method: Dispersion, dissipation, isotropy, Galilean invariance, and stability, Phys. Rev. E 61 (2000) 6546–6562.

28

[43] Z. Guo, C. Zheng, Analysis of lattice Boltzmann equation for microscale gas flows: Relaxation times, boundary conditions and the Knudsen layer, Int. J. Comput. Fluid Dyn. 22 (2008) 465–473. [44] Q. Zou, X. He, On pressure and velocity boundary conditions for the lattice Boltzmann BGK model, Phys. Fluids 9 (1997) 1591–1598. [45] A. J. Wagner, I. Pagonabarraga, Lees-Edwards boundary conditions for lattice Boltzmann, J. Stat. Phys. 107 (2002) 521–537. [46] T. Inamuro, T. Ogata, S. Tajima, N. Konishi, A lattice Boltzmann method for incompressible two-phase flow with large density differences, J. Comput. Phys. 198 (2004) 628–644. [47] L. Li, R. Mei, J. F. Klausner, Boundary conditions for thermal lattice Boltzmann equation method, J. Comput. Phys. 237 (2013) 366–395. [48] H. Daiguji, P. Yang, A. J. Szeri, A. Majumdar, Electrochemomechanical energy conversion in nanofluidic channels, Nano Lett. 4 (2004) 2315– 2321. [49] D. Stein, M. Kruithof, C. Dekker, Surface-charge-governed ion transport in nanofluidic channels, Phys. Rev. Lett. 93 (2004) 35901. [50] Q. Pu, J. Yun, H. Temkin, S. Liu, Ion-enrichment and ion-depletion effect of nanochannel structures, Nano Lett. 4 (2004) 1099–1103. [51] S. J. Kim, Y. C. Wang, J. H. Lee, H. Jang, J. Han, Concentration polarization and nonlinear electrokinetic flow near a nanofluidic channel, Phys. Rev. Lett. 99 (2007) 044501. [52] B. Rotenberg, I. Pagonabarraga, D. Frenkel, Dispersion of charged tracers in charged porous media, Europhys. Lett. 83 (2008) 34004. [53] A. Ramos, H. Morgan, N. G. Green, A. Castellanos, AC electrokinetics: a review of forces in microelectrode structures, J. Phys. D: Appl. Phys. 31 (1998) 2338–2353. [54] A. Ajdari, Pumping liquids using asymmetric electrode arrays, Phy. Rev. E 61 (2000) 45–48.

29

[55] A. Ajdari, Electrokinetic ‘ratchet’ pumps for microfluidics, Appl. Phys. A 75 (2002) 271–274. [56] N. A. Mortensen, L. H. Olesen, L. Belmon, H. Bruus, Electrohydrodynamics of binary electrolytes driven by modulated surface potentials, Phys. Rev. E 71 (2005) 056306. [57] M. Lian, J. Wu, Ultrafast micropumping by biased alternating current electrokinetics, Appl. Phys. Lett. 94 (2009) 064101–064101. ˇ ˇ [58] J. Hrdliˇcka, P. Cervenka, M. Pˇribyl, D. Snita, Zig-zag arrangement of four electrodes for ac electro-osmotic micropumps, Phys. Rev. E 84 (1) (2011) 016307.

30

Figure captions

Figure 1: Flow chart of the LBM for the electrokinetic flows.

31

Figure 2: Flow velocity ux averaged over the channel cross-section as a function of the ion concentration C0 for various values of the surface charge density σ on the channel wall. The symbol indicates the results of the LBM computation, and the solid line indicates the analytical formula (45).

32

Figure 3: Zeta potential ζ as a function of the ion concentration for various values of the surface charge density σ on the channel wall. The symbol indicates the zeta potential estimated from the LBM computation results using Eq. (46), and the solid line indicates the analytical formula (43) with y = H. The value of ζ∗ is indicated by the dashed line.

Figure 4: Flow velocity profiles ux for various values of the ion concentration C0 in the cases of (a) σ = −2 × 10−5 [C/m2 ] and (b) σ = −50 × 10−5 [C/m2 ]. The symbol indicates the results of the LBM computation, and the solid line indicates the analytical formula (44).

33

.5

Figure 5: Profiles of the ion concentration Cm normalized by C0 for various values of C0 in the cases of (a) σ = −2 × 10−5 [C/m2 ] and (b) σ = −50 × 10−5 [C/m2 ]. The symbols indicate the results of the LBM computation ( : cation and : anion), and the lines indicate the profiles of the Boltzmann distribution with Eq. (43) (the solid line: cation, and the dashed line: anion).

◦

Figure 6: Profiles of (a) flow velocity ux and (b) ion concentration Cc normalized by C∗ for various values of the surface charge density σ. The symbol indicates the results of the LBM computation, and the solid line indicates the analytical formula: Eq. (49) in panel (a) and Eq. (47) in panel (b).

34

Figure 7: Time evolution of the ion concentration Cm at the surface (x, y) = (L/8, H) for C0 = 0.032 [mol/m3 ]. The solid line indicates the result obtained with the present method, and the dashed line indicates the result obtained with the method proposed by Wang and Kang [27].

35

Figure 8: Time evolution of the profiles of (a) cation concentration Cc , (b) anion concentration Ca , (c) potential φ, and (d) net charge ρe , along the line at x = L/8 with y ∈ (0, H). The solid line indicates the results obtained with the present method, and the dashed line indicates the results obtained with the method proposed by Wang and Kang [27].

36

Figure 9: Contour plots of the net charge ρe in the case of ζ = 5 [mV] and C0 = 0.032 [mol/m3 ], at (a) t = 3 × 10−7 [s] and (b) 5 × 10−5 [s]. The solid line indicates the results obtained with the present method, and the dashed line indicates the results obtained with the method proposed by Wang and Kang [27].

37

Figure 10: Vector plots of the current density I in the case of ζ = 5 [mV] and C0 = 0.032 [mol/m3 ], at t = 3 × 10−7 [s]. (a) the result obtained with the present method, and (b) the result obtained with the method proposed by Wang and Kang [27].

38

Figure 11: Time evolution of the current in the y-direction Iy . The profiles along the line (a) at x = L/8 and (b) at x = 5L/8, with y ∈ (0, H). The solid line indicates the results obtained with the present method, and the dashed line indicates the results obtained with the method proposed by Wang and Kang [27].

Figure 12: Geometry of the channel of finite length.

39

Figure 13: Two-dimensional plots of (a) the electrical potential φ, (b) the net charge ρe , (c) the pressure p − p0 , and (d) the vector field u, for the case of C0 = 0.032 [mol/m3 ].

40

Figure 14: Two-dimensional plots of (a) the electrical potential φ, (b) the net charge ρe , (c) the pressure p − p0 , and (d) the vector field u, for the case of C0 = 0.002 [mol/m3 ].

41

Figure 15: Time evolution of the flow velocity ux averaged over the channel cross-section at x = LM + L/2 for various values of ion concentration C0 . The sampled values are indicated by the symbol connected by the solid line.

Figure 16: Flow velocity profiles ux (a) at the entry (x = LM ), (b) at the middle (x = LM + L/2), and (c) at the exit (x = LM + L) of the channel, for the case of C0 = 0.032 [mol/m3 ]. In each panel, the profiles at t = 10−9 , 10−8 , 10−7 , and 10−6 [s] are indicated by the solid line, and the profiles in the steady state are indicated by the dashed line.

42

Figure 17: Flow velocity profiles ux (a) at the entry (x = LM ), (b) at the middle (x = LM + L/2), and (c) at the exit (x = LM + L) of the channel, for the case of C0 = 0.002 [mol/m3 ]. In each panel, the profiles at t = 10−9 , 10−8 , 10−7 , and 10−6 [s] are indicated by the solid line, and the profiles in the steady state are indicated by the dashed line.

Figure 18: Flow velocity ux averaged over the channel cross-section as a function of the ion concentration C0 for the cases of HM = 0.5H, H, and 1.5H. The symbols connected by the dotted lines indicate the results of the LBM computation, and the other lines indicate the analytical formula (45) with Px = 0 and Eq. (51).

43

Figure 19: Distribution of (a) the electrical potential φ and (b) the pressure p − p0 along the center line (y = 0) for various values of ion concentration C0 . The electrical potential distribution obtained using Eq. (51) (and the corresponding equations for the entry and exit regions) is also indicated in panel (a) by the dashed line.

44

Figure 20: Flow velocity profiles ux at the middle of the channel (x = LM + L/2) for (a) C0 = 0.002 and 0.008 [mol/m3 ], and (b) C0 = 0.032 and 0.128 [mol/m3 ]. The symbol indicates the results of the LBM computation. The dashed line indicates Eq. (44) with Px = 0 and Eq. (51), and the solid line indicates the same equation (51) with the values of Ex and Px computed by taking the derivative of the potential and pressure distributions of Fig. 19.

45

The highlights of the submitted manuscripts are the following: • A numerical scheme based on the lattice Boltzmann method for electrokinetic flows. • The scheme includes a simplified algorithm for solving the Nernst―Planck equation. • The boundary condition for ion flux is correctly implemented. • The scheme is validated comparing with the one-dimensional analytical solution. • A two-dimensional electro-osmotic flow is analyzed as an application of the method.