Two-dimensional numerical simulation of microscale thermoelectrokinetic instability

Two-dimensional direct numerical simulation of thermoelectrokinetic instability in electrolyte bounded by semiselective surfaces was conducted for the first time. The results of numerical experiments show that the convective motion in the electrolyte takes place due to the stability loss of one-dimensional steady state solution. This motion manifested itself in twodimensional structures similar to electroconvective rolls. What fundamentally differs is formation of thin lines connecting the membranes, where the salt concentration reaches its maximum value. Electric current in the system, obviously, will firstly flow through these lines.


Introduction
Problems of electrokinetics have recently captured attention of scientists from different fields of knowledge because of rapid developments in micro-, nano-and biotechnology.It tends to be very promising for practical issues to use ion-exchange surfaces, permeable only for one type of ions, cation or anion.Among such surfaces are semiselective electric membranes, electrodes and systems of micro-and nanochannels.It is well known that electric current in electrolyte in microgap between semiselective membranes is proportional to voltage ¨V, if the values of ¨V are rather small.This is a manifestation of Ohm's law, so such a regime is usually called ohmic.With increase of ¨V transition to limited mode takes place (electric current weakly depends on voltage) [1].Further magnification of external electric field leads to reappear linear dependence of the electric current from the voltage.Transition from limited to overlimiting mode is caused by appearance of liquid movement owing to instability development, theoretically predicted in [2] and numerically confirmed in [3].However, it must be said, that numerical and experimental data used to match only qualitatively, but not quantitatively, e.g. the experimental critical value of potential drop was less than that obtained numerically.As it has been shown in [4], neglecting of thermic effects caused by Joule heating may be the reason of such an unpleasant mismatch.A new type of microscale instability, occurring then electric current at the layer is oppositely directed to gravity vector, was found out.The critical potential drop can be smaller than that one for electokinetic instability threshold, if special combination of new introduced parameters is considered.In this work we add to Nernst-Planck-Poisson-Stocks system of equation an energy equation with source term caused by Joule heating of the electrolyte owing to electric current in it.Also, an additional term, respected for bouncy force, occurs in the Stocks equation.Membranes are assumed to be partly insulted (Fig. 1).For two-dimensional direct numerical simulation of the problem, the special algorithm based on one suggested in [5] and developed in [3] was designed.

Formulation of the problem
In dimensionless form the problem is described by the Nernst-Planck-Poisson-Stokes (NPPS) system Here, ܿ േ are the molar concentrations of cations and anions; ‫܃‬ is the fluid velocity vector; ߔ is the electrical potential in the liquid phase; ܲ is the pressure.The Joule heating effect is described by energy equation The boundary conditions are as follows: The above problem is described by seven dimensionless control parameters: ߥ is the Debye numberǡ ߢ is a coupling coefficient between the hydrodynamics and the electrostatic (it is essential that the coupling coefficient depends only upon the physical properties of the electrolyte), ߂ܸ is the potential drop between the membranes, ܴܽ is the Rayleigh number, ‫݅ܤ‬ is the Biot number (characterizing the system's thermal insulation with respect to environment), ‫݁ܮ‬ is Lewis number, ‫‬ is the membrane interface concentration in the boundary conditions.The characteristic quantities and value of dimension parameters can be found in [4].We set in our calculations ț=0.2, Bi=0.01, p =5, Ȟ=0.1,Le=0.013unless otherwise indicated.

Theoretical flashback
It is worthwhile to remind here some basic results from [4]: the system is stabilized by the thermal effects at ܴܽ Ͳ and destabilized at ܴܽ ൏ Ͳ. Changing of the sign of ܴܽ is nothing more than the rotation of the system at ͳͺͲι relative to the direction of gravity.In the problem, there are two The first one is the Zaltzman-Rubinstein electrokinetic instability, characterized by the first term in the right-hand side of (7).The second is connected with the Joule heating, and is characterized by the second term in the right-hand side of (7) there ݇ is a wave number of superimposed small perturbations.
For the case when the heating is absent (ܴܽ ൌ Ͳ) the instability is caused by the Coulomb force in the space charge region and the corresponding formation of a slip velocity at its border.The nonuniformity of the slip velocity leads to the electrokinetic instability [2] The competing mechanism is associated with the Joule heating, mainly in the electroneutral region.The second mechanism is connected with heating and the thermic expanding of the liquid, but it is totally different from the Rayleigh-Benard convection and the instability is caused by an induced nonuniformity of the conductivity in the electrolyte.

Direct numerical simulation
In this paper we present the results of two-dimensional direct numerical simulation of the problem.The designed algorithm is based on one suggested by Nikitin in [5] and adopted for electrokinetic problems in [3].A finite-difference method with second-order accuracy is applied for the spatial discretization.A uniform grid is used in the homogeneous tangential ‫ݔ‬ direction; the grid is stretched in the normal ‫ݕ‬ direction via a ‫݄݊ܽݐ‬ stretching function in order to properly resolve the thin double layers attached to the membrane surfaces.For solving stiff problems implicit methods require the inversion of rather large matrices and thus are extremely costly, while explicit methods of time advancement require a very small time step and, hence, are prohibitively ineffective.A semi-implicit method is found to be a reasonable compromise: only a part of the right-hand side of the system is treated implicitly.A special semi-implicit ͵ ଵ ଷ -step Runge-Kutta scheme is used for the eventual integration in time.The details of the numerical scheme will be presented elsewhere.
For the natural "room disturbances", the infinite spatial domain is changed to a large enough finite domain that has lengths ‫ݔ݈‬ ൌ ݈ in both spatial dimensions, with the corresponding wave number ݇ ൌ ʹߨȀ݈.The condition that the solution be bounded as ‫ݔ‬ ՜ λ is changed to periodic boundary conditions.The length of the domain l has to be taken large enough to make the solution independent of the domain size.

Simulation results
Two-dimensional modeling confirms the results of investigation conducted in [4]: two competing mechanisms of instability are determined by the parameters ߢ and Ra.The relation between these parameters determines which of the instability mechanisms will be decisive for the destabilization of the system.The case Ra = 0 separates the destabilizing and stabilizing effects of the Joule heating.In contrast to the Rayleigh-Benard convection, the second mechanism stabilizes the system with heating at the bottom (ܴܽ Ͳ) and destabilizes it with heating at the top (ܴܽ ൏ Ͳ).
For ܴܽ ൏ Ͳ, with decreasing or increasing ȁܴܽȁǡ the heat effects prevail over the electrokinetic effects and a drastic change of instability modes occurs: the critical voltage ܸ ‫כ‬ decreases dramatically.Moreover, the short-wave instability changes to a long-wave instability.
The convective motion in the electrolyte takes place due to the stability loss of one-dimensional steady state solution.Electroconvective vortices observed in [3] differ from those obtained during our calculations: they are short-wave and their centres are situated on the border of electric double layer, and our vortices are long-wave and situated in the centre of the channel.(Fig. 2).If the heat-insulation of the system is good enough, overlimiting mode might occur at essentially smaller values of potential drop than in isothermal model.Two-dimensional modelling shows that for supercritical values of potential drop two-dimensional structures take place.These structures (Fig. 3 (b)), are similar to electrokinetic "spikes" (Fig. 3 (a)), but have a long-wave behaviour.The interesting feature of such structures is formulation of thin lines with high ion concentration between two membranes.Because of higher ion concentration these lines have a higher conductivity, and, hence, the electric current in the system goes mostly through these lines.

Conclusion
A new kind of instability caused by Joule heating near charge-selective surfaces and its influence on the electrokinetic instability are investigated theoretically by two-dimensional direct numerical simulation.The numerical experiments confirm the results of one-dimensional analysis [4] and allow clarifying the physical mechanism of thermoelectrokinetic instability.This work was supported, in part, by the Russian Foundation for Basic Research (Projects No. 14-08-00789-a, No. 15-58-45123-IND_a and 16-08-00643-a)

Figure 1 .
Figure 1.Schematics of the system.