Extension of the Improved Bounce-back Scheme for Electrokinetic Flow in the Lattice Boltzmann Method

In this paper, an improved bounce-back boundary treatment for fluid systems in the lattice Boltzmann method [Yin, X.; Zhang J. extended to handle the electrokinetic flows with complex boundary shapes and conditions. Several numerical simulations are performed to validate the electric boundary treatment. Simulations are presented to demonstrate the accuracy and capability of this method in dealing with complex surface potential situations, and simulated results are compared with analytical predictions with excellent agreement. This method could be useful for electrokinetic simulations with complex boundaries, and can also be readily extended to other phenomena and processes.


Introduction
With growing interest in bio-Micro Electro Mechanical Systems (MEMS) and bio-Nano Electro Mechanical Systems (NEMS) applications and fuel cell technologies, electrokinetic flows have become one of the most important non-mechanical techniques in the application of microfluidics and OPEN ACCESS nanofluidics [1,2].Electro-osmotic flow (EOF) is a promising approach to drive the microfluidics under an external electric field, such as sample injection, chemical reaction, species separation and energy supply [3,4].Due to these important applications, EOF in microchannels has received interesting attention [5][6][7][8][9][10][11].
Electro-osmotic flow (EOF) is a basic electrokinetic phenomenon, where an electrical double layer (EDL) is formed due to the interaction between an electrolyte solution and a dielectric surface [12].From the macroscopic point of view, the EOFs are governed by the Navier-Stokes (NS) equations for fluid flow and the Poisson-Boltzmann equation for the electrical potential.Many studies have been carried out on electro-osmotic flow in microchannel.The lattice Boltzmann method (LBM) has been generally accepted as a useful simulation method for complex flows [13][14][15].The LBM approach is advantageous in dealing with complex boundary geometries [16,17] and could be potentially more efficient with advanced computational technologies.Because of its distinctive advantages over conventional numerical methods, the lattice Boltzmann method has introduced to simulate electro-osmotic flow in micro devices.Warren [18] made the first attempt to apply the LBM to solve the Navier-Stokes equations for the solution, while the conservation equation for each ionic species and the Poisson equation for the electrical potential were solved via the "moment propagation" method.He and Li [19] proposed a lattice Boltzmann scheme for analyzing the electrochemical processes in an electrolyte based on a locally electrically neutral assumption.With a multiple-component LBM model, this scheme has also been utilized to study the electrohydrodynamic drop deformation in an electric field [20].The electrokinetic flows in microchannels is simulated by the lattice Boltzmann method with one-dimensional linearized solution of the Poisson-Boltzmann equation [5,7,8].Melchionna and Succi [21] solved the nonlinear Poisson-Boltzmann equation by an efficient multi-grid technique and then predicted the flow behavior using a lattice Boltzmann scheme.The multi-grid technique has great efficiency to solve the nonlinear Poisson-Boltzmann equation; however it has rarely been extended for complex geometries [22].Recently, The LBM has been applied to study the mixing enhancement in heterogeneously charged microchannels [23][24][25][26][27][28] and the roughness and cavitation effects in electro-osmotic microfluidics [29,30].
As with other numerical methods, boundary conditions play crucial roles for the simulation validity and stability.However, unlike the tremendous efforts in developing accurate boundary treatments for LBM models for fluid flows and convection-diffusion systems [31][32][33][34][35][36], boundary methods for LBM models of electric field have not been addressed adequately.Typical electric field LBM simulations are performed in regular domains with flat boundaries aligned along the lattice grid lines.Several studies have considered rough surfaces [26,30]; however, the rough surfaces were actually modeled as flat, stair-like patches.Recently, Yin and Zhang [32] developed an improved bounce-back method for fluid flows, which is discussed in [37][38][39].
In this paper, an improved bounce-back method for fluid flows [32] is extended to simulate the electrokinetic flows with arbitrary boundary shapes.Numerical simulations demonstrate that the boundary treatments have accurately represented the spatial geometry as well as the surface potential.An example calculation is also performed to illustrate the application of our boundary treatment for electrokinetic studies.This study could be useful for LBM simulations of electric fields in systems with complex surface geometry and surface conditions.

Macroscopic Governing Equations for EOF
For incompressible EOF in microfluidic channel, the governing equations including the continuity equation and the momentum equation can be described as follows: where u is the velocity vector, ρ is the density of solution, P is the pressure, υ is the kinetic viscosity of the flow, F represents the external force and is given as: where e ρ is the net charge density, and E is the external electric field.
The drive force of the EOF is indicated by the body force term ( e ρ E ) in the momentum equation and is caused by the action of the induced electrical field on the net charge density in the EDL region.EDL theory [40] related the electrostatic potential and the ion distribution in the bulk solution can be well approximated by the Poisson equation as follows: where ψ is the electrical potential, ε and 0 ε are the dimensionless dielectric constant and permittivity of vacuum, respectively.For the flows over a non-conducting stationary surface, the ion distribution can be well approximated by the Boltzmann distribution: where i n is the ionic number concentration of the i -th species, , i n ∞ is the ion concentration in the bulk solution, i z is the valence of type-i ions, e is the charge of an electron, b k is the Boltzmann constant and T is the absolute temperature.For a symmetric electrolyte present study, the net charge density is given as follows: 2 sinh Combining Equations ( 4) and ( 6) yields the nonlinear Poisson-Boltzmann distribution for the EDL potential in the dilute electrolyte solution: For the surfaces with low surface electric potentials, the Debye-Hückel approximation ) can be applied and the Poisson-Boltzmann equation can be linearized to: where: and its reciprocal 1 κ − , the so-called Debye length, is usually used as a measure of the EDL thickness in Debye-Hückel theory.

Numerical Method
The numerical method adopted in this work requires the solution of the Navier-Stokes equations for fluid flow and the Poisson-Boltzmann equation for electric potential distribution.A lattice structure D2Q9 with complex boundary conditions is proposed to solve the governing equations using two LBM model, corresponding to the fluid flow and electric potential, respectively.It is necessary to introduce the LB evolution equations and the boundary treatments in this section.

Lattice Boltzmann Model for the NS Equations
The evolution equation corresponding to the NS equations with external force is given as: where ( , ) is the distribution function for the flow fields at location x and time t and the subscript α indicates the lattice direction.t Δ is the time step.The parameter f τ is the relaxation time for the density distribution function.F α is the forcing term corresponding to the applied external electric field.
In order to recover the correct NS equation, the density equilibrium distribution in this work can be typically expressed as: where c is the lattice speed defined as x t Δ Δ , x Δ the lattice grid size.For the D2Q9 lattice model, the lattice weight factors depends on the length of the corresponding lattice vector, and is given by The relaxation time for fluid flow is related with the fluid viscosity by: The forcing term caused by the interaction of the EDL field with the externally applied electrical field is incorporated into the discrete Boltzmann equation by following the method described: ( ) ( ) The Chapman-Enskog expansion can be used to recover the macroscopic NS equations, and Macroscopic quantities such as the density and fluid velocity can then be evaluated from the distribution functions as:

Lattice Boltzmann Model for Poisson-Boltzmann Equation
Here, to solve the Poisson-Boltzmann equation, we adopt an LBM algorithm proposed by Oulaid et al. [41] because of its good numerical accuracy.The Poisson-Boltzmann equation can be considered as a convection-diffusion equation at the steady state.A lattice distribution function g α is introduced, and its evolution is described by the following lattice Boltzmann equation: where ( , ) is the distribution function for electric potential.as a parameter for different schemes.Both the minimum and maximum values of θ (0 and 1) have been tested with diffusion and convection-diffusion systems, and no significant influence on the solution accuracy is found.In this work we use a forward scheme for the temporal derivative with 1 θ = for simplicity: The electric potential ψ can be calculated from the distribution functions by: and the equilibrium distribution eq g α of electrical potential evolution function is: Following the spirit for the fluid flow, the following differential equation through the Chapman-Enskog analysis can be derived: and the solution to the Poisson Equation ( 4) can be obtained at the steady state of the simulation when the partial differential term on the left-hand side approaches zero.The potential diffusivity in Equation ( 21) is defined as: ( )

Boundary Conditions
As with any other numerical methods, correct and accurate boundary treatments play a crucial role in LBM simulation.Many useful schemes for boundary condition have been developed for solving different physical problems.To model the fluid-solid interaction on the complex geometries, the mid-point bounce-back scheme [30] is used for the flow field.As shown in Figure 1  As we known, the evolution equations consist of two computational steps: ( , ) ( , ) with or g α .When implementing the boundary conditions with LBM, the difficulty is how to finish the collision and streaming steps at the boundaries.In Figure 1, after the collision step at the fluid node f x , the distribution function f α * leave f x , and is then assumed to be bounce-back at the midpoint where e e α α = − , m u is the midpoint velocity at the point m x to be determined.

For 1 2 Δ ≤
, the midpoint m x locates between b x and f x , and the midpoint velocity m u can be readily obtained with a linear interpolation: where b u is the imposed boundary velocity at the intersection point b x , and f u is the flow velocity calculated at the fluid node f x .For 1 2 Δ > , the midpoint m x is in the solid domain and therefore an extrapolation is needed to obtain velocity m u : where ff u is the velocity at the fluid point ff x .
Following the above velocity boundary treatment, here, we use the electric potential m ψ at the midpoint to calculate the electric potential distribution function at the bounce-back nodes: with the midpoint electric potential m ψ obtained via: where b ψ is the imposed boundary electric potential at the intersection point b x , f ψ is the electric potential calculated at the fluid node f x , and ff ψ is the electric potential at the fluid point ff x .

Electric Potential with Flat Surface
First, we apply the improved bounce-back scheme to calculate the potential profile between two parallel plates, both of constant potential, immersed in an electrolyte solution.Near the charged surfaces, ions in the electrolyte solution will be redistributed and the electric diffuse layer will be established.The ion charge density can be related to the local potential via the Boltzmann distribution.For surfaces with low surface potentials, the Debye-Huckel approximation can be applied and the Poisson-Boltzmann equation can be solved by Equation (8).The solution of this linearized Poisson-Boltzmann equation between two identical plates of a separation H is: where x is the transverse location across the gap with 0 x = at the centerline.We use 0 1   ; green symbol =0.2 Δ ).
As for typical LBM boundary models, the numerical accuracy is studied.We choose different channel heigh 16, 32, 64, 128 H = and calculate the error between the LBM results and theoretical solutions.
The errors are plotted in Figure 4, and linear fitting are conducted in the logarithmic graph.The fitting slope is usually considered as the accuracy order of a numerical model.As show in Figure 5, the accuracy order is 1.97, about 2, indicating a second-order for this system by the improved bounce-back treatment.

Electric Potential with Complex Geometry
Next, to examine the performance of our method for more complex boundary shapes and conditions, we consider the electric field between two coaxial circular surfaces with inner and outer radii in R and out R , respectively.With no net charge, the general solution is given as: ( ) ln The corresponding exact solution for this case is: In the simulation, we set 40 surfaces are put at the center of the domain.Simulation result for this system is plotted in Figure 3.We also find that the results agree well with analytical solution.

Application in Electro-Osmotic Flows
In this section, a charged spherical particle immersed in an electrolyte solution is considered with a side length of 2 μm.A 100 100 100 × × uniform grid is used and the particle center ( ) x y z locates at the center of the cubic domain.Detailed numerical results on electro-osmotic flow and the effects of variation of ionic concentration, the sphere radius, external electric field and electric potential on velocity profile are presented.The numerical results are also compared with analytical solutions.In the simulation, the Poisson-Boltzmann equation is solved to obtain a steady solution firstly.And then the Navier-Stokes equations with the external force term is solved.We select a symmetric solution with : 1:1 z z = (for example, KCl, NaCl, etc.) and assume the solution has similar physical properties.The parameters are the ionic molar concentration 0.01 , where A N is Avogadro's number, the dielectric constant of the solution set to be 1.0.Periodic boundary conditions are applied in all the three directions, and hence the simulated system actually represents a cubic array of spheres uniformly distribution in space.
The algorithm and boundary treatment described in previous sections have been used to simulate the electric flow with curved boundary.For the purpose of validation, the solution of the Poisson-Boltzmann equation around a spherical particle with thin EDL layers is given as: where r is the distance to the spherical center.We use 0 10 mV ψ = and 500 V/ m x E = in our simulation.The particle has a radius of R = 0.6 μm.
Figure 6a shows the electric potential as a function of the distance to the particle center.The red circles are from our LBM calculation and the black curve is the theory solution according to Equation (32).Good agreement can be observed between them.The potential distribution at c y y = is presented in Figure 6b.The distribution appears circularly symmetric and isotropic, and this is confirmed by the fact that all the simulated ψ ~r data points fall approximately on a single curve in Figure 6a.
When an external electric field is applied, an electric force F will be generated in the electrolyte solution near the surface due to non-zero charge in that region, and the electrostatic force can thus induce fluid flows along the electric field direction.This phenomenon is called the electro-osmosis.

Conclusions
In this paper, we have extended the improved bounce-back boundary treatment for LBM flow simulations to electric field simulations.Several simulations have also been performed to examine our boundary methods in term of ability to deal with complex boundary situations.An example simulation of electro-osmotic flow with a charge sphere particle immersed in an electrolyte solution has also been presented.Comparisons with theoretical predictions show excellent agreement for all simulations, and our method therefore could be useful for future electrokinetic simulations with complex boundary geometries.Furthermore, the boundary treatment in this work can be applied to LBM simulations for other processes and phenomena that can be described by similar differential equations.
g τ is the relaxation time for the electric potential distribution function.G α is related to the net charge term in Equation (4) by: , the link between the fluid node f x and the solid node s x intersect the physical boundary at b x .The fraction of the intersected link in the solid domain region is

Figure 1 .
Figure 1.Schematic illustrations for the mid-point bounce-back scheme.
channel heigh16,32, 64, 128 H = in our simulation.The calculated potential profile is plotted in Figure2.The symbols are results from LBM simulations, and the black lines are theoretical solutions predicted from Equation(29).Excellent agreement can be seen with different channel height between the simulation results and theory.

Figure 2 .
Figure 2. Electric potential distributions from our LBM simulation (symbols) and the analytical solution (black lines) with different height between two identically charged plates in an electrolyte solution.

Figure 4 .
Figure 4.The error for electric potential between two identically charged plates in an electrolyte solution.

Figure 5 .
Figure 5. Electric potential distributions between two coaxial circular surfaces.

and 2 C
can be determined by the boundary conditions on the surfaces.For the Dirichlet boundary conditions, the electric potential on both surfaces are

, and the electric potential with 0 ψ
as a constant.The external electric field is only applied in x -direction, i.e.,

Figure 6 .
Figure 6.The electric potential distribution (a) and (b) and electro-osmotic flow (c) and (d) around the spherical particle in the y c y = plane.