Numerical Study of Electro-Osmotic Fluid Flow and Vortex Formation

The phenomenon of electro-osmosis was studied by performing numerical simulations on the flow between parallel walls and at the nozzle microchannels. In this work, we propose a numerical approximation to perform simulations of vortex formation which occur after the passage of the fluid through an abrupt contraction at the microchannel. The motion of the charges in the solution is described by the Poisson–Nernst–Planck equations and used the generalized finite differences to solve the numerical problem. First, solutions for electro-osmotic flow were obtained for the Phan–Thien/Thanner model in a parallel walls channel. Later simulations for electro-osmotic flow were performed in a nozzle. The formation of vortices near the contraction within the nozzle was verified by taking into account a flow perturbation model.


Introduction
For more than 200 years, it has been known that certain liquid solutions can flow through a small channel by applying an electrical potential difference at both ends of the microchannel. This area of study, or technique of manipulation of the fluid, i.e., the transport of particles through application of an external electric field, is called electrokinetics. In this scope, there are several applications in research such as chemistry, engineering and biomedical [1,2]. Electrohydrodynamics is the connection between theories of electromagnetism and hydrodynamics. Thus, the dynamics of fluids, electrokinetics and electrochemistry form the basis for the study of electrohydrodynamics. A general review and further details on the physical principles of electrokinetic effects can be found in [3].
One of the first studies of electro-osmosis is the work of Reuss [4], in which it is studied the flow of a polar fluid by application of an external electric field which acts on ionic solution, driving the movement of charged particles near the walls. At 1879, Helmholtz [5] reported the electrical flow parameters for the electrokinetic transport and he named of the Electrical Double layer (EDL). Subsequently, Goy [6] and Chapman [7] contributed in the studies of the flow distribution of charges near the walls in a capillary tube. In 1918 von Smoluchowski [8] performed capillary studies, generalizing the theory proposed by Helmholtz. Another great contribution was given by Debye and Hückel in 1923 [9]. They determined the ionic concentration in the solution by linearizing the Boltzmann distribution to the energy. Subsequently several studies related to electrokinetic effects in Newtonian fluids were performed. Later [10] carried out studies in ultra-fine capillaries, observing the effect of the surface potential on the flow. One year later Rice and Whitehead published the work on the surface potential effects for cylindrical capillaries [11]. Theoretical studies and numerical simulations of electro-osmotic flows began to be published by Yang and Li [12] and Patankar and Hu [13]. In 2000, Bianchi et. al. [14] solved the problem of electro-osmotic flow in a T-junction using the finite element method. Analytical solutions for Newtonian flow in channels were demonstrated by Dutta and Beskok [15]. In 2002 Lin et al. [16] solved numerically the Nernst-Planck equation for the ionic transport together the Navier-Stokes equations for Newtonian fluids. The study of electro-osmotic flow for non-Newtonian fluids is more recent, due to the difficulty imposed depending on the constitutive model to be used. In 2008 Park and Lee [17] determined the electro-osmotic velocity for viscoelastic flows while Zhao et al. [18] carried out studies of electro-osmotic flows of fluids modeled by the power law. Tang et al. [19] also performed a numerical study of non-Newtonian fluids in electro-osmotic flows. Afonso et al. [20] demonstrated the analytical solution for viscoelastic flows in a parallel plates and circular section channel, taking into account a non-zero pressure gradient and the electrokinetic effect, where one of the constitutive models used was the simplified Phan-Thien/Thanner (sPTT). In 2015, Peng et al. [21] studied the effects of a concentration gradient within a channel, taking into account a mixture of different electrolytic properties. Recently, Song et al. [22] carried out studies of numerical instabilities in mixing a ferrous solution with water in a parallel walls channel. Niu et al. [23] studied a microfluidic pumping by physical and numerical experiments in which it was found that the 2D flow obeys power law decay. In other study Niu et al. [24] demonstrate that self-generated solvent flow can be used to generate long-range attractions on the colloidal scale. The dynamics of this system is governed by an effective conservative energy that for large separations depends on the inverse of the distance. In this way, associated investigation of pumping including electrokinetic phenomenon were performed by [25,26]. Recently Arcos et al. [27], studied the behavior of a sPTT fluid subject to time dependent zeta potentials. Pimenta and Alves [28], study the implementation of electrically driven flow models of viscoelastic fluids in the finite-volume framework of OpenFOAM, applied the induced-charge electro-osmosis around a conducting cylinder and charge transport across an ion-selective membrane. Pimenta and Alves [28] code is freely available as open-source code.
Microfluidics have found widespread application in science and technology, as in medicine. The flow in microfluidic devices is usually created either by applying a pressure gradient (e.g., using a syringe pump), or using electrokinetic effects. Although effective for laboratory experiments, using a syringe pump is not practical for development of portable equipment for lab-on-a-chip applications, such as biological sample analysis. In contrast, electrokinetic effects are particularly useful at micro-and nanoscales. As an example, electro-osmosis can be efficiently employed to induce fluid pumping using electric fields, therefore avoiding the need to integrate micro-pumps and mechanical valves for flow control, which increase the complexity and the cost of disposable microfluidic devices. Small scales typical of microfluidics, increase the relevance of surface forces and electrokinetics effects and, for non-Newtonian fluids, there is an enhancement of the role of fluid elasticity, beyond anything that can be attained at macroscale. Thus, given the importance of the topic, in this work computational experiments were performed for electro-osmotic flows between parallel plates and in a nozzle. The vortices are introduced into the flow by supposing the existence of a perturbation near the corners in the channel contraction. The effect of viscoelasticity was also verified on the nozzle with the imposed electrical perturbation.
In Section 2, we describe the governing equations of the problem to the channels. Section 3 we present the computational method used to perform the simulations. Section 4 reports the verification for the parallel plates channel. Moreover, we present the proposed formulation for the nozzle flow problem and computational results obtained. Conclusions are drawn in Section 5.

Governing Equations
We are assuming the fluid incompressible laminar and isothermal flow. Moreover, here the treatment of the governing equations will be given in dimensionless form: where u is the velocity field, t is time, p is the pressure, Re = ρUH/η 0 is the Reynolds number, U is the average velocity, H is the channel height, ρ the mass density and η 0 denotes the total shear viscosity η 0 = η s + η p . The rate of deformation tensor D = 1 2 ∇u + (∇u) T and T is the elastic stress. The dimensionless solvent viscosity coefficient is given by β = η s η 0 . The evolution in time of the polymeric stress tensor is related by where De = λU/H is the Deborah number and λ is the relaxation time of the fluid. Here we will be use a kernel-conformation tensor [29][30][31] and then determine the stress tensor. An alternative form to describe viscoelastic models is by using the conformation tensor, A. In general the equation for A can be written as ∂A ∂t where M(A) is define according to the viscoelastic model. By taking into account the decomposition of the velocity gradient ∇u T = Ω + B + NA −1 , it was possible reformulate the tensor conformation, where in Equation (6) both Ω and N are anti-symmetric tensors, B is symmetric and commutes with A. This decomposition enables the rewriting of the constitutive equation for tensor A One of interest subject is to solve Equation (4) or Equation (5) for high values of Deborah number De. Numerical methods are unstable for certain critical values of De. In order to overcome such failure, Fattal and Kupferman [29] proposed a reformulation of the differential constitutive equations into an equation for the matrix-logarithm of the conformation tensor. Extending the ideas proposed by Fattal and Kupferman [29,30], Afonso et al. [31] presented a generic kernel-conformation tensor transformation that allows apply different kernel functions to the matrix transformation, in which the evolution equation for k(A), can be expressed in its tensorial formulations as where B and M are symmetric tensors constructed by the orthogonalization of the diagonal tensors. Thus, the HiG-Flow system solves Equation (8) instead of Equation (4), for more details please see [32].
Tensor A is only used in the HiG-Flow system to perform kernel-conformation, which is intrinsic to the system when viscoelastic flow simulations are taking out. In fact, numerical simulations in HiG-Flow can be performed using numeric stabilizers, where the kernel conformation is on, allowing the user to establish a stabilizer in a friendly manner without having to make modifications to the code kernel, but only in the input file of simulation. Considering a Newtonian fluid flow, the tensor S is null and the velocity and pressure are only updated at each step of time.

Phan-Thien/Tanner Model
Here we are interested in use the PTT model to solve the constitutive equation and then determine the velocity field. For this model, the right-hand side of Equation (4) can be written as: The dimensionless parameter ε is related to the steady-state elongational viscosity in extensional flows and ξ is a parameter related with the molecular slip. If ξ is null, the model reduces to the simplified PTT (sPTT). On the other hand if ξ is not null, there will be a non-zero second normal-stress difference in shear, leading to secondary flows in ducts with non-circular cross-sections, which is superimposed on the streamwise flow [33]. In fact the right-hand side for the conformation tensor Equation (5) is given by In this way the equations of motion can be solved for the PTT model fluid flow. Therefore, the electrical contribution must be taken into account for fully description of flow. In the next section we will show how the electrical source term was obtained.

Electro-Osmotic Force
Electro-osmotic fluid flows were first studied in the flow between a parallel wall microchannel and secondly in a contraction microchannel.
The schematic representation of electro-osmotic flow is show in Figure 1. The applied potential along the axis of the channel provides the driving force necessary to occur the electro-osmotic flow. The PTT constitutive model [34] was used for studies involving electro-osmotic flow of non-Newtonian fluids. Numerical simulations are performed in two dimensions, taking into account L, w >> H, where H is the channel height, L is channel length and w the channel width. Moreover, due to the symmetry of the problem, we analyzed only half of the channel, i.e., 0 ≤ y ≤ H.
For the problems subjected to the electro-osmotic forces, there exists a source term in Equation (2), F = ρ e E, where E is the electric field. The electric field appears due to two contributions; one is the applied potential φ and the other due to the induced potential ψ which changes in the transversal direction to the channel walls. Thus, E = ∇φ + ∇ψ. The formation of the Debye layer occurs due to the spontaneous movement of the charged species near the channel wall, causing a charge redistribution in the fluid that originates the electrical double layer [35]. Therefore, the equations to be solved for these potentials are given by where ρ e = δ(n + − n − ) is the charge density, δ = n 0 ezH 2 / ζ 0 with n 0 the reference concentration, e is the elementary charge, z the charge valence, the dielectric constant and ζ 0 the potential on the wall. Moreover, the density charge is related by n + and n − which are the positive and negative ionic concentration, respectively. In this work we solve numerically the Nersnt-Planck equation for the ionic concentration: where Pe = UH/D is the Peclet number that depends on ionic diffusion D. The potential to thermal energy ratio is given by α = ezζ 0 /k B T with k B the Boltzmann's constant and T is the absolute temperature. The set of Equations (11)-(13) are solved to obtain the electrical force.

Computational Procedures
The HiG-Flow system [32,36] was used to obtain numerical solutions reported in this work. Computational domain to the simulation is obtained through HiG-Tree, which generates a hierarchical mesh (for further information, refer to references [32,36]). For bi-dimensional case, this mesh is a generalized quad-tree [37]. Hierarchical meshes impose difficulties in the numerical scheme based on Cartesian approximations, and requires the use of spatial interpolations at unknown points of the stencil. The interpolations of the properties in the center of the faces and in the center of the cells are made by the technique of moving least squares, which uses a given set of points where the property is known to estimate a unknown value in a neighbor point. Differential equations are discretized by the method of finite differences. Solvers using the PETSc library (Portable, Extensible Toolkit for Scientific Computation) [38] are used to solve linear systems [32,36]. The machine used to perform numerical simulations has a Core i7 2.0 GHz CPU, 16 Gb memory. The mesh with refinements along the channel were obtained using HiG-Tree. The total length of the mesh is 20H. Better results were obtained with a maximum ∆x/∆y = 4. Near the walls the minimum size is ∆y min /H = 7.8125 × 10 −4 and ∆x min /H = 3.1250 × 10 −3 . Figures 2 and 3 shows the mesh used to simulate the electro-osmotic flow in parallel plates and in a nozzle, respectively.   In order to illustrate the velocity profile convergence we were used three different grids. The main difference is the grid refinement, increasing the refinement near the walls to improve convergence. Effect of mesh refinement is show in Figure 5. Please note that for the mesh with ∆y min /H = 2.5000 × 10 −2 represented by circles, the nearest point the wall is clearly outside the expected curve. Increasing refinement, ∆y min /H = 6.2500 × 10 −3 represented by squares, all points remain within the analytical curve, but there are spaces with no points due the sharp effect near the wall and for ∆y min /H = 7.8125 × 10 −4 represented by triangles, these spaces were filled and then the simulated curve agrees with analytical curve. In the next section we report the results.

Boundary Conditions on the Walls
On the channel walls, Neumann conditions apply for p and φ. The potential ψ assumes the potential reference value ζ 0 . The concentrations at the walls are given by the Boltzmann distribution Equation (17) and the non-slip condition is imposed for velocity.

Inflow and Outflow
Dirichlet conditions are applied for pressure, so the channel inlet and outlet have exactly the same pressure, resulting in a pure electro-osmotic flow, i.e., p in equal to p out . Similar conditions were applied for the electrostatic potential φ, i.e., φ in = φ out agreeing with ∇φ = 0. Concentration assumes the n 0 reference value at inlet, and the outlet Neumann conditionn · ∇n ± = 0 is imposed to obtain a fully developed flow. It was taking Neumann condition for the velocity at the inlet and outlet; thus, the flow is developed under the action of the electric body force.

PTT Model
The PTT model was used to solve the pure electro-osmotic fluid flow through parallel plates. We will first show the results potential and ionic concentration which are coupled properties of the flow, according to set of Poisson-Nernst-Planck equations given by Equations (11)- (13). Figures 6 and 7 show the behavior of ψ and n + , n − respectively. The distribution of charged particles causes the potential variation near the wall, and the ionic concentration near the wall is related by changes on potential ψ in that neighborhood. For κ = 10 potential changes are seen at greater distances to the wall if compared to the curve show for κ = 20. As the Debye layer becomes thinner, κ = 50, the potential changes quickly near the channel walls so the charges are approaching to the walls as shown for κ = 300. Therefor the velocity depends on y and κ it will also vary according to the relaxation time of the fluid. In this way, the velocity profile is affected by the Deborah number De κ = λκu sh [20,39], as shown in Figure 8. For this problem we set De κ = 0.5 and De κ = 2.0. In addition, the parameters ξ = ε = 0.01 were fixed. The Debye parameter assumed two different values, κ = 10 and κ = 100. For low Deborah number De κ = 0.5, the velocity profile indicates behavior similar to that of a Newtonian fluid. Increasing the number of Deborah is apparent the changes on velocity profile respecting to the Newtonian profile due to the low viscosities near the wall, in fact influenced by the appearance of the shear stress of the PTT model. The normal and shear stress are shown in Figures 9 and 10, in which we can observe that the numerical results are similar to the analytical solution.

sPTT Model
In this model the parameter ξ = 0 at the expression Equation (10). Moreover, instead of considering the pure electro-osmotic flow, we will impose a non-zero pressure gradient on the channel. In this sense, by taking into account the steady-state flow, the solutions for the velocity and the viscoelastic tensor for the sPTT model depend on the parameter Γ defined as: If the pressure gradient is greater than zero, dp/dx > 0, the fluid is pushed in the opposite direction of flow as shown in Figure 11 for Γ = 1.0, Γ = 2.0 and Γ = 2.778, represented by downward triangles, diamonds and left-pointer triangles respectively. On the other hand, if dp/dx < 0, the fluid is pushed in the same direction of flow, resulting in the profiles shown in Figure 11 for Γ = −0.5 e Γ = −1.778 represented by circles and squares respectively. The Figure 12 shows the velocity profile of the electro-osmotic flow for sPTT fluid. For Γ = −1.0 represented by right-pointer triangles and Γ = 2.77 by left-pointer triangles. As with the full PTT fluid, the velocity profile in the sPTT model is influenced by the Deborah number, which in these simulations is De κ = 2.5. In fact there exists a dependence on De 2 κ . More details, please see [20]. This effect is understand for the curves when Γ = 0 showed in the Figure 12 represented by solid lines for De κ = 0 and De κ = 2.5, where it seen changes relative to the Newtonian profile. Normal and shear stress are shifted relative to the corresponding curves for pure electro-osmotic fluid flow which Γ = 0 represented by solid lines as shown in Figure 13. For the normal component T xx , the curves indicate the decrease of the stress near the channel wall when the pressure gradient is positive, Γ = 2.77, represented by upwards triangles. On the other hand, if the pressure gradient is negative, Γ = −1.0, so T xx increases near the wall, as can be seen the curve represented by downwards triangles. Similar behavior is observed for shear stress T xy .

Fluid Flow in a Nozzle
Channels with contractions and expansions are generally used when one is interested in mixing different fluids in a microchannel. Here the idea is trying to mimic the behavior of the fluid flow for this type of geometry, specifically regarding to the electrical effects near the corners. Physical description of these effects related to the movement of the charges near the corners is very complex, but we can try in some way to find a coherent approach linking the phenomenon to numerical simulation. As experimentally found by [40], near the corners there are fluctuations in zeta potential on the walls, causing a change on the velocity in that neighborhood. Here, initially the proposal is to put on an external perturbation on the potential, which we will consider to be due to punctual charges located in the nozzle near the corners, as show in Figure 14. Let us consider that the charges are very distant from each other, so we will be not taken into account the isolated interaction between them. In addition, we consider that the interaction of the perturbation charges with the part before the contraction of the channel is negligible, assuming that the concentration of ions in this place remains stable, i.e., the ionic depletion starts when the fluid arrives inside the contraction of the channel, in the region with intense gradient in Figure 14. It is known from the electromagnetism theory that the potential ϕ(x, y) due to an elementary point charge is given by where P(x 0 , y 0 ) is the place of the point charge is located, K is the dielectric constant and e the elementary charge. In this way, the potential ψ p on the walls near the corners is the sum of reference zeta potential with perturbation potential, ψ p = ζ 0 + ϕ, and one can be write: where we defined ω = KeH∇φ/ζ 2 0 as a constant which depends on external applied field ∇φ, i.e., increasing ∇φ, increases the electro-osmotic velocity and the perturbation effect on the standard potential ψ is consequently accentuated. On the other hand, increasing the distance from the wall, decreases the potential due the perturbation resulting ψ p → ζ 0 from Equation (21). Additionally, the distance between parallel plates must be greater than the distance of action of the imposed perturbation, i.e., (x − x 0 ) 2 + (y − y 0 ) 2 H. By making ψ p = ψ * p ζ 0 , x = x * H e ω = ω * H, one can write the expression for dimensionless potential: Figure 14. Illustration of imposed perturbation near the corners.
According to [41], the potential on the corners have an inverse square root of ionic concentration dependence, ψ p ∼ n −1/2 . In this way, the ionic concentration depends on inverse square of potential, lead to write the following relation: where n bc is the value of ionic concentration on the channel wall without perturbation effect. When ψ p approaches the zeta potential, perturbation effects are negligible then the ionic concentration is given by n p = n bc . On the other hand, if ψ p increases due the perturbation effect, the ionic concentration decreases representing depletion of charges near the corners. Using the expression for ψ p one can writing Equation (23) in dimensionless form: For this proposal, the Equations (22) and (24) Figures 15 and 16 show the ionic concentrations n + and n − respectively, in (a) for zero perturbation ω * = 0 and (b) for strong corner effect imposed by ω * = 6.0 × 10 −3 , where accentuated charge changes was observed after contraction in the channel. These changes are directly linked to the ψ changes and in fact occurs due to the coupling of Poisson-Nernst-Planck equations. The perturbation effect increases as ω * is increased, as show in Figure 17. Therefore, the perturbation make changes on the pressure and velocity near the corners as show in Figures 18 and 19 . Increasing the perturbation parameter ω * , the flow is concentrate in the center of the channel as show in Figure 20, reducing the cross-section area of the fluid flow, contributing in the case of fluid mixing to be more efficiently process.  Figure 21 shows the curves of potential ψ as a function of the applied perturbation. For ω * = 0, there is no perturbation and the potential curves represented by squares are similar to that curve obtained for the parallel plate channel. By imposing ω * = 1.0 × 10 −3 , curves represented by circles, noticed an increase in the absolute value of ψ in the center of the contraction and after the expansion, ψ approaches zero again indicating the decrease of perturbation effect out of the contraction. Increasing the perturbation parameter, ω * = 3.0 × 10 −3 , curves represented by upwards triangles, can be seen an more accentuated increment in the potential within the contraction and after the expansion of the channel. For ω * = 6.0 × 10 −3 , curves represented by downward triangles, perturbation effect changes the potential to have approximately the same value in the plateau within and after the contraction. The behavior of the ionic concentration curves is shown in Figure 22. For ω * = 0, curves represented by squares, the behavior is similar to that obtained for parallel plates channel. If ω * = 1.0 × 10 −3 , curves represented by circles, we observe a decrease in concentration n + and an increase of n − within the contraction and this effect is smaller after the channel expansion. When ω * = 3.0 × 10 −3 the concentration n + decreases while n − increases and the curves are inverted relative to the reference concentration n 0 , as shown in Figure 23. This effect is even more pronounced for ω * = 6.0 × 10 −3 . In this way, the coupling between the potential and ionic concentrations is related in Figures 21-23, indicating that the ionic concentration variation is affected by the potential increment due to perturbation. It is natural to expect that the perturbation imposed on the nozzle will cause numerical instabilities in the ionic concentration and potential ψ, since we are forcing a new distribution of charges for the problem, but if the perturbation is small enough, the numerical method is stable. The velocity profiles are affected by both the perturbation and the variation of the pressure within the contraction in the channel. These curves were obtained at the center of the contraction and Figure 24 shows the Newtonian and viscoelastic profiles. The curves represented by squares correspond to ω * = 0, i.e., velocity profile variations occur due to the variation of the pressure in the contraction, and tends to form a crest before to the velocity curve decay and then decreases until zero on the wall, where full squares correspond to the velocity for the Newtonian fluid. This effect is most accentuated for the viscoelastic fluid, as can be seen in the curve represented by empty squares. Increasing the perturbation effect, ω * = 1.0 × 10 −3 , the crests are suppressed as shown the curves represented by full and empty circles corresponding to the Newtonian and viscoelastic fluid, respectively. If ω * = 3.0 × 10 −3 the profile forms a crest in the center of the channel and increases the distance from the center, the velocity tends to zero quickly, as shown the curves with upwards triangles. This effect is similar to that applying a pressure gradient on the channel ends, as seen in Section 4.2. Finally, if ω * = 6.0 × 10 −3 , it is observed the existence of negative velocities near the wall, represented by the downwards triangles. Behavior of the polymeric tensor, Figures 25 and 26, corroborates with the results obtained in Section 4.2 in order to show the similar effect to of applying a negative pressure gradient when the fluid enters in the contraction, represented by the fully points and a positive pressure gradient after the expansion of the channel represented by the empty points.

Conclusions
Computational experiments were performed to obtain the solution of electro-osmotic flows. In order to validate the HiG-Flow, first the results for the pure electro-osmotic fluid flow of viscoelastic fluids (using the full PTT) through parallel plates were obtained and compared with satisfactory accuracy against analytical results. Later, the sPTT model was studied and it was seen the influence of an applied gradient pressure on the flow for a parallel plates channel. Furthermore, it was observed the vortex formation in a nozzle by imposed perturbation near the corners. Results were obtained and the contraction in the channel naturally imposes a pressure gradient on that neighborhood and then it was found that the imposed perturbation near the corners causes a similar effect. It is notable that the imposition of the perturbation causes instabilities on the solution, as can be seen in the concentration and potential curves. It is expected because one is imposing conditions unknown to the flow, but for weak perturbation the numerical solution keep stable.