Titre : Électrophorèse en non-équilibre d'une microparticule selective en ions Résumé :

Les travaux de cette these sont consacrés à l'étude théorique et numérique du mouvement d'une particule sphérique sélective de cations dans un électrolyte sous l'action d'un champ électrique externe d'intensités diverses. La région externe de l'électrolyte n'est pas limitée, par conséquent, l'influence des parois n'est pas prise en compte.

Les principaux résultats suivants ont été obtenus dans les travaux:

1. La solution asymptotique et numérique du problème de l'électrophorèse dans un champ électrique d'intensité faible.

2. Dans le cas limite d'une faible intensité de champ électrique, la dépendance de la vitesse électrophorétique de l'intensité du champ électrique est établie, ce qui est une généralisation de la formule de Helmholtz-Smoluchowski pour une particule diélectrique (condition d'équilibre) à une particule sélective d'ions (condition non-équilibre).

3. Une solution analytique du problème est obtenue dans le cas limite d'un champ électrique de haut intensité dans chacune des couches limites minces. Aussi bien la compréhension de la dépendance de la vitesse des particules au champ électrique que le phénomène d'instabilité électrocinétique dans la région de charge sont d'une importance pratique. Cet effet peut être utilisé pour créer de nouveaux dispositifs microfluidiques pour mélanger des liquides à la micro-échelle. Les données de simulation numérique peuvent être utiles pour évaluer l'efficacité de l'utilisation de l'instabilité électrocinétique dans la création et l'utilisation de tels dispositifs. The work is devoted to the theoretical and numerical study of the motion of a spherical cation-selective particle in an electrolyte under the action of an external electric field of various strengths. The outer electrolyte region is not limited, so the problem is completely free from the possible influence of the walls.

The following main results were obtained in the work: 1. Asymptotic and numerical solution of the electrophoresis problem in a weak electric field.

2. In the limiting case of the weak electric field, the dependence of the electrophoretic velocity on the electric field strength is derived, which is a generalization of the Helmholtz-Smoluchowski formula for a dielectric particle (an equilibrium process) to the ion-selective particle (non-equilibrium process).

3. An analytical expansion of the problem is obtained for the limiting case of a strong electric field in each of the thin boundary layers. 4. The emergence of electrokinetic instability in the region of the incoming cations flux was found. A scenario of the flow transition from regular to chaotic is obtained.

Both the understanding of the dependence of the particle velocity on the electric field strength and the discovered phenomenon of electrokinetic instability of the space charge region are of practical importance. The electrokinetic instability can be used to create new microfluidic devices for mixing liquids on a micro-scale. Numerical simulations data can be useful for evaluating the efficiency of using electrokinetic instability in the creation of such devices.

Résumé

Les phénomènes électrocinétiques sont couramment utilisés dans le domaine de la production de systèmes de petite taille (de l'ordre du micro ou nanomètre) et des perspectives technologiques de création de "laboratoires sur chips" microfluidiques. La plupart de ces dispositifs ont un fluide de travail conduction qui doit être mis en mouvement. 2. Trouver une solution analytique au problème posé dans le cas limite d'une faible intensité de champ électrique ; obtenir la relation analytique de la vitesse de la microparticule (basée sur l'état de l'équilibre des forces agissant sur la particule) en fonction de l'intensité du champ électrique appliqué.

3. Trouver une solution analytique dans le cas d'un champ électrique de grande amplitude avec une formation de couches limites imbriquées les unes après les autres et dans le volume externe de l'électrolyte ; obtenir une expression analytique de la vitesse de glissement électroosmotique à la limite de la couche de diffusion et de la vitesse électrophorétique. 4. Développer un algorithme numérique pour résoudre le problème complet non-linéaire dans une formulation axisymétrique en coordonnées sphériques et effectuer des simulations numériques du problème pour une large gamme de paramètres. 5. Etudier l'effet de la charge de surface d'une particule sur la solution du problème ; établir les limites d'applicabilité des approches analytiques pour les champs électriques de faible et forte amplitude. [START_REF] Wang | Million-fold preconcentration of proteins and peptides by nanofluidic filter[END_REF]. Identifier les principales transitions de bifurcation lorsque l'on modifie les paramètres du problème. Le problème est formulé pour le cas axisymétrique dans un système de coordonnées sphériques (r, θ), se déplaçant avec la vitesse de la microparticule avec l'origine au centre de la particule. Le système d'équations est transformé en une forme sans dimension dans laquelle apparaît deux paramètres : le nombre de Debye ν et le coefficient de couplage des propriétés électrostatiques et hydrodynamiques de l'électrolyte, κ.

Les conditions limites suivantes sont définies à la surface de la microparticule : le potentiel est égal à une constante qui, sans perte de généralisation, peut être égal à zéro ; les composantes de la vitesse du fluide satisfont aux conditions d'imperméabilité et d'adhérence ; la concentration en ions positifs est fixée (égale à une constante empirique p), le flux d'ions négatifs à la surface de la particule est nul.

Loin de la particule, à l'infini, les conditions aux limites suivantes sont définies : le champ électrique externe -E ∞ est appliqué, les composantes de la vitesse à l'infini tendent vers la vitesse de la particule U ∞ , la concentration en ions tend vers celle de l'équilibre.

Les conditions initiales suivantes sont considérées : t = 0 : c + = c -= 1, où c + et c -sont les concentrations en ions chargés positivement et négativement, respectivement.

Ce problème est résolu à l'aide de méthodes analytiques et numériques. Il reste néanmoins non linéaire à cause des équations de transport ionique. Le système reste fortement associé à un petit paramètre apparaissant dans l'équation de Poisson -le nombre de Debye ν. Physiquement, cela signifie la présence de fines couches ioniques près de la surface de la microparticule. La solution numérique du problème repose sur la méthode des différences finies.

La simulation numérique directe (par différences finies) du problème est complétée par une analyse analytique par le biais de la stabilité linéaire d'un écoulement stationnaire.

Cette solution analytique est basée sur une approche asymptotique et est conçue pour deux cas limites -un champ électrique de faible amplitude (E ∞ → 0) et un champ électrique de très forte amplitude (E ∞ → ∞). La solution analytique est également basée sur le fait que l'épaisseur de la double couche électrique est faible par rapport au rayon de la microparticule, ce qui signifie que ν → 0.

Pour un champ électrique de faible amplitude, la méthode asymptotique raccordée a été utilisée avec le petit paramètre E ∞ → 0 et les solutions à l'ordre zéro et à l'ordre 1 ont été obtenues. Comme ν << 1, le domaine a été subdivisé en deux : interne et externe. Près de la surface de la particule, la décomposition interne est utilisée, le développement d'une nouvelle variable est introduit. Quant à la solution externe, elle est prise électriquement neutre, ρ = c + -c -= 0 ou c + = c -. Cette analyse permet de trouver ainsi les champs de concentration en ions et la distribution du potentiel électrique. Pour déterminer la vitesse électrophorétique U ∞ , en l'absence d'inertie, on écrit l'équilibre entre les forces de contraintes visqueuses et celles de contraintes électrostatiques de Maxwell. Les calculs prédisent une réponse linéaire de la vitesse des microparticules au champ électrique appliqué (de faible intensité) et est analogue à la formule de Helmholtz-Smoluchowski pour les particules diélectriques.

La prochaine étape est consacrée à la solution semi-analytique du problème de l'électrophorèse d'une particule sélective en ions dans un champ électrique de forte intensité, i.e.

E ∞ → ∞. La principale différence par rapport au cas précédemment étudié de faible intensité réside dans la modification de la structure de la région chargée électriquement près de la surface de la particule. Dans le cas d'un champ électrique de faible intensité, la base du nuage ionique chargé autour de la particule est une double couche électrique (DCE) et le flux ionique à travers la particule peut être négligé. Dans le cas où E ∞ → ∞, la situation change radicalement : maintenant, la charge dans la DCE peut être négligée par rapport à la charge déposée par les ions. Une structure complexe de couches imbriquées est formée. Une double couche électrique mince d'épaisseur O(ν), ν → 0, est adjacente directement à la surface de la microparticule. Au-dessus de la DCE, il y a la zone dite de charge d'espace, qui a une épaisseur de l'ordre de ν 2 3 log(ν), beaucoup plus épaisse que la double couche électrique avec, de plus, une charge beaucoup plus grande (ρ > 0) dans cette couche. Cette couche est hors équilibre et est formée par le flux d'ions entrants. Dans cette zone de charge d'espace, la concentration en sel est très faible, K = c + + c -≈ 0 (c'est pourquoi, dans la théorie des membranes, on parle de zone de solution appauvrie).

Cette zone de charge d'espace est suivie d'une mince couche de diffusion électriquement neutre (ρ ≈ 0), la concentration en sel K varie de K ≈ 0 (à la frontière avec la zone de charge d'espace) à la valeur d'équilibre K = 2. L'analyse asymptotique est basée sur la non prise en compte de la double couche électrique. Ces 3 couches dépendent de l'angle et disparaissent en un point correspondant à l'angle critique θ = θ 0 . L'analyse de stabilité linéaire a été réalisée pour chaque couche séparément. Pour la région de charge d'espace, il a été possible d'obtenir une solution totalement analytique et une formule permettant d'estimer l'épaisseur de cette couche. La solution dans la région de la couche de diffusion a conduit à un système non linéaire complexe, qui a été résolu numériquement. Il n'a pas été possible d'obtenir une expression analytique pour l'épaisseur δ de cette couche, mais il s'est avéré qu'elle était de l'ordre de O(E -1 ∞ ). Dans la région externe de ces 3 couches (dans l'électrolyte), on suppose la neutralité électronique, c'est-à-dire ρ = 0, mais aussi l'équilibre de la solution saline, c'est-à-dire K = 2.

Chapter 1

Introduction to Electrokinetics

Electrokinetics utilize the electric fields to cause electrostatic forces acting on charged or polarizable fluids and suspended particles. These electrostatic forces induce the motions of fluids and particles.

Electrokinetic phenomena are currently widely used. This is primarily due to the latest developments in the production of micro-scale devices and the technological capabilities of creating microfluidic "labs on a chip". Most of these devices have a working fluid that needs to be set in motion.

Methods that are used in practice to bring fluid into motion on a conventional scale are often based on fluid instability caused by inertia. However, at the microscale, this phenomenon is suppressed by viscous forces. Thus, the liquid should not be set in motion due to turbulence, but due to the forces of molecular diffusion. On a very small scale, molecular diffusion is fast enough, however, in microfluidic devices with a typical size in the range from 10 to 100 µm, it will take about 100 seconds to mix liquids with a diffusion coefficient of the order of 10 -10 m 2 /s [START_REF] Squires | Induced-charge electro-osmosis[END_REF]. Therefore, there is a need for new methods of pumping, mixing, manipulating, and separating liquids and particles at the microscale.

The development of microfabrication technology in the last decades makes it possible to scale down the electric fields to micro and even nanoscale. At these scales, the effects of electrokinetics become dominant, which makes it one of the most promising methods for manipulating particles. One of the remarkable distinction of electrokinetics is that it does not require the creation of mechanical moving parts, which is extremely difficult on a small scale. Since the propulsion of both liquid and particles occurs due to the electric field, the devices based on electrokinetic effects are much reliable than mechanical devices.

In fact, the electrokinetic phenomena provide one of the most popular and effective nonmechanical ways to bring fluid into motion at the micro and nanoscale. The main idea is as follows: a liquid in which there are ions of a dissolved substance (electrolyte) adjoins a charged solid surface. A diffusion cloud of oppositely charged ions gets collected near this surface. An external electric field produces a force acting on the charged diffusion layer, which causes the fluid to flow relative to the particle or solid surface. The electrokinetic flow near a fixed surface is called the electroosmotic flow, and the electrokinetic movement of particles suspended in a fluid is called electrophoresis. Electroosmosis and electrophoresis have a wide range of applications in the field of analytical and colloidal chemistry, mixing the liquids, separating particles moving in a liquid by size and other properties.

In many microfluidic devices, transport is caused by using electric current through the ion-selective material placed in aqueous electrolytes. This type of material is a principal part in many applications such as water desalination [START_REF] Biesheuvel | Nonlinear dynamics of capacitive charging and desalination by porous electrodes[END_REF][START_REF] Kwak | Microscale electrodialysis: Concentration profiling and vortex visualization[END_REF]4], separation of synthetic particles and biological cells [5], preconcentration of samples [START_REF] Wang | Million-fold preconcentration of proteins and peptides by nanofluidic filter[END_REF], and fuel cells [START_REF] Pourcelly | Ionic Interactions in Natural and Synthetic Macromolecules. Applications of charged membranes in separation, fuel cells and emerging processes[END_REF]. The most commonly used in practice with selective property is permselective membranes [START_REF] Xu | Ion exchange membranes: State of their development and perspective[END_REF].

They can be classified into cation-selective (cation-exchange) and anion-selective (anionexchange) membranes. Since negative-charged groups are fixed inside the cation exchange membrane, anions are rejected by the negative charge and cannot permeate through the cation exchange membrane (see Fig. After 1798, when A. Volta created the first stable source of electric current, many experiments with an electric field followed. One of them was the experiments conducted in Moscow by F. Reuss, who published his discovery in 1809 [START_REF] Reuss | Sur un novel effet de l'électricité galvanique[END_REF]. The two types of experiments were conducted. In the first one, the clay plug was placed in a U-shaped tube and under applied voltage, the water level in one part of the tube rises [START_REF] Wall | The history of electrokinetic phenomena[END_REF]. In the second experiment, the silica sand particles migrated in the electric field. So the phenomena of electroosmosis and electrophoresis were discovered, respectively.

After the Reuss experiments, there were plenty of experiments related to electroosmosis. The first quantitative experiments were carried out by G. Wiedemann, whose first work [START_REF] Wiedemann | Ueber die bewegung von flüssigkeiten im kreise der geschlossenen galvanischen säule[END_REF] on electroosmosis was published in 1852. He studied electroosmosis in tubes and found out, for example, that the ratio of the mass of water transferred to the electric current does not depend on the applied voltage and the tube radius. This is in full accordance with the Smoluchowski equation obtained 60 years later [START_REF] Smoluchowski | Contribution a la théorie de l'endosmose électrique et de quelques phénoménes corrélatifs[END_REF]. Wiedemann's research was an important step towards a quantitative understanding of electroosmosis.

The phenomenon inverse to electroosmosis was discovered experimentally and reported by G. Quincke [START_REF] Quincke | Ueber die fortführung materieller theilchen durch strömende elektricität[END_REF] in the work published in 1859. In this experiment, the water was pumped through the tube and the potential difference between its ends was observed.

This phenomenon was called the streaming potential (or the Quincke effect). Quincke has found that the potential difference varies for different systems, but it has the same sign. It follows from the fact that many natural materials on the phase interface acquire a negative charge in an aqueous solution. Another fact discovered by Quincke is that the potential difference at water pumped through the tube does not depend on the crosssection and thickness of the plug and, furthermore, streaming potential depends linearly on the applied pressure. As a result, Quincke postulated the existence of a charge outside the charged solid surface. The sign of this charge is opposite to the sign of the surface charge, i.e. it screens the surface charge from the liquid side. The presence of this charge was a prerequisite for a qualitative explanation of both electroosmosis and the streaming potential. The introduction of this charge was crucial for colloidal science and essentially meant the discovery of the electric double layer (EDL). On the basis of Quincke's EDL idea, H. Helmholtz in 1879 derived a quantitative theory of electrokinetic phenomena.

The study of electrokinetic phenomena was continued by E. Dorn, who publish results in 1880 [START_REF] Dorn | Ueber die fortführung der electricität durch strömendes wasser in röhren und verwandte erscheinungen[END_REF]. The experimental setup was composed of a vertical glass tube with electrodes at two vertical levels, filled with water and grains of sand that settled in the tube. The potential difference between the two electrodes was measured by Dorn during the process of sand grains sedimentation. This potential difference is called the sedimentation potential or the Dorn effect.

By 1880, four classical electrokinetic phenomena were discovered, namely electrophoresis, electroosmosis, streaming potential (Quincke effect), and sedimentation potential (Dorn effect).

From the Reuss experiments, it was clear that the clay particles were charged, but the origin of this charge was unclear until 1886 when it was shown by J. Hoff that the osmotic pressure of a solution of a monovalent salt was actually twice the expected value for an undissociated salt. Thus, it became clear that dissociation of salts into ions occurs in aqueous solutions.

The theory of EDL continued to be developed and in 1910 Gouy has proposed [START_REF] Gouy | Sur la constitution de la charge électrique à la surface d'un électrolyte[END_REF] the structure of EDL which included a diffuse layer of charge near solid surface. In 1913, Chapman independently formulated an equivalent theory [START_REF] Chapman | LI. a contribution to the theory of electrocapillarity[END_REF]. Although Gouy's theory of EDL has been mentioned by Smoluchowski in the footnote, it has not been used in the derivation of the Helmholtz-Smoluchowski formula for electroosmotic velocity since it does not depend on the exact structure of EDL. Otherwise, Helmholtz and Smoluchovsky could not have received a sufficiently correct relation between velocity and electric field strength. The theory of Helmholtz was published in 1879 and improved by Smoluchowski in 1903 [START_REF] Smoluchowski | Contribution a la théorie de l'endosmose électrique et de quelques phénoménes corrélatifs[END_REF]. Thus, the Helmholtz-Smoluchowski theory precedes the Gouy-Chapman theory of EDL.

The Gouy-Chapman model of EDL introduces an important parameter, which characterizes the thickness of EDL and is usually called Debye length. It gives the thickness estimation of the ionic shell surrounding each sufficiently large particle in the electrolyte and this thickness depends on the ions concentration and the temperature,

λ D = εRT z 2 F 2 c ∞ (2.1)
where ε is the dielectric permittivity; R is the gas constant; T is the absolute temperature;

z is the valence of ions; F is Faraday's constant and c ∞ is the concentration of ions in the bulk. Due to a lack of knowledge about the EDL structure, the particles were considered to be electrically neutral (without surface charge), in spite of the fact they can migrate in the electric field. It has been proven later [START_REF] Billitzer | Elektrische doppelschicht und absolutes potential. kontaktelektrische studien i[END_REF][START_REF] Billitzer | Über die elektrizitätserregung durch die bewegung fester körper in flüssigkeiten. kontaktelektrische studien ii[END_REF] that dielectric particles cannot be electrically neutral.

EDL

--- - - - - - - - - -

Modern theoretical and experimental works

Electrophoresis is a fundamental phenomenon not only for electrokinetics but also for theoretical physics. Significant efforts have been made to understand this phenomenon, starting with the works of Smoluchowski and Helmholtz [START_REF] Smoluchowski | Contribution a la théorie de l'endosmose électrique et de quelques phénoménes corrélatifs[END_REF][START_REF] Helmholtz | Studien über electrische grenzschichten[END_REF]. The ultimate goal of these works is to find the dependence of the particle velocity on the applied external electric field. Such knowledge is necessary, in particular, for the design of new methods of characterizing, separating, and transporting particles [START_REF] Bazant | Induced-charge electrokinetic phenomena[END_REF].

The electrophoretic velocity strongly depends on both the physical properties of the particle itself and the electrolyte in which it is located. Distinctive results were obtained for different types of particles studied in the literature: dielectric, metal, and ion-selective (see Fig. There is also a relatively new type of particle for colloidal science -Janus particles with non-uniform distributed surface or bulk properties in an electric field. Currently, they are widely used in electronic displays, medicine and diagnostics, chemical catalysis.

The classical case of dielectric particles has been studied for over a hundred years. The mathematical description of the electrophoresis of particles with a thin electric double layer (EDL) formed near its surface is based on the classical theory of EDL, developed by Helmholtz and Smoluchowski [START_REF] Smoluchowski | Contribution a la théorie de l'endosmose électrique et de quelques phénoménes corrélatifs[END_REF][START_REF] Helmholtz | Studien über electrische grenzschichten[END_REF],

U ∞ = εζ µ E ∞ , (2.2)
that is, the electrophoretic velocity U ∞ is linearly proportional to the strength of the applied electric field E ∞ and the zeta potential ζ. Moreover, it is independent of the particle radius a. The velocity also depends on the properties of the electrolyte, namely, the dielectric permittivity ε and the dynamic viscosity of the liquid µ.

The relation (2.2) describes the so-called "linear" electrophoresis or electrophoresis of the first kind since the velocity is linearly proportional to E ∞ . This formula is valid for the weak electric field. With a sufficiently large electric field strength, the validity of this formula is violated due to the non-linearity of the original equations.

There are many experimental and theoretical works devoted to the study of the motion of dielectric particles suspended in the electrolyte under the influence of an external electric field. The effects of surface conductivity, the thickness of EDL and surface charge were studied by Van Der Put et al. [START_REF] Van Der Put | Electrical conductivity of dilute and concentrated aqueous dispersions of monodisperse polystyrene particles. influence of surface conductance and double-layer polarization[END_REF] and Lyklema [START_REF] Lyklema | Fundamentals of Interface and Colloid Science[END_REF]. Wiersema et al. [START_REF] Wiersema | Calculation of the electrophoretic mobility of a spherical colloid particle[END_REF] obtained a numerical solution for a spherical particle. Assuming the symmetry of the problem, O'Brien and White [START_REF] O'brien | Electrophoretic mobility of a spherical colloidal particle[END_REF] numerically calculated the dependence of the electrophoretic mobility of the particle on the zeta potential on its surface.

The special mention among the modern works deserves the theoretical studies conducted by Yariv's group. Hamed and Yariv [START_REF] Yariv | An asymptotic derivation of the thin-debye-layer limit for electrokinetic phenomena[END_REF] and Yariv [START_REF] Hamed | Induced-charge electrokinetic flows about polarizable nano-particles: the thick-debye-layer limit[END_REF] used the expansion in a small parameter (electric field strength) to analyze the stationary electrokinetic flow near a perfectly polarized spherical nanoparticle. The small value of this parameter means a significant Debye length in comparison with the particle radius.

Yariv and Davis [START_REF] Yariv | Electro-osmotic flows over highly polarizable dielectric surfaces[END_REF] used the macroscale model proposed by Yossifon et al. [START_REF] Yossifon | Symmetry breaking in induced-charge electroosmosis over polarizable spheroids[END_REF] to analyze the electrokinetic flow around the dielectric surface for the case of thin EDL. The authors developed and analyzed electrical effects where polarization changes the zetapotential distribution. They have found that the non-linearity has a slowing down effect, i.e. the increase of the external electric field makes the microparticles to move slower than predicted by the Helmholtz-Smoluchowski theory.

According to Smoluchowski's theory, the particle velocity linearly depends on the value of the zeta potential (Eq. (2.2)). But this dependence is valid only for small and moderate values of the zeta potential and becomes invalid for high values [START_REF] Schnitzer | Dielectric-solid polarization at strong fields: Breakdown of smoluchowski's electrophoresis formula[END_REF][START_REF] Schnitzer | Macroscale description of electrokinetic flows at large zeta potentials: Nonlinear surface conduction[END_REF]. Nevertheless, Schnitzer and Yariv [START_REF] Schnitzer | Dielectric-solid polarization at strong fields: Breakdown of smoluchowski's electrophoresis formula[END_REF] have analytically proven that this relationship is also valid for moderate zeta potentials. In their next work, Schnitzer and Yariv [START_REF] Schnitzer | Macroscale description of electrokinetic flows at large zeta potentials: Nonlinear surface conduction[END_REF] presented a general analysis of the electrokinetic transport near a strongly charged dielectric particle under the assumption of a thin EDL, which is not limited to the case of the weak electric field.

In a strongly non-linear mode, another formula must be used for electrophoretic mobility (not the Smoluchowski formula). Schnitzer et al. [START_REF] Schnitzer | Weakly nonlinear electrophoresis of a highly charged colloidal particle[END_REF] generalized the results of [START_REF] Schnitzer | Macroscale description of electrokinetic flows at large zeta potentials: Nonlinear surface conduction[END_REF] to the case of large zeta potential values, when the surface conductivity becomes significant.

In all problems associated with dielectric particles, the charge on the particle is as-sociated only with EDL. In this case, the effect of the electric field is neglected. Such processes are called processes of the "first kind" by Dukhin (Rubinstein and Shtilman [START_REF] Rubinstein | Voltage against current curves of cation exchange membranes[END_REF] call these processes "underlimiting"). Since the influence of electric current can be neglected, these processes are in electrochemical equilibrium.

Electrokinetic phenomena for metal particles have also been studied for a long time both theoretically and experimentally [START_REF] Gamayunov | Pair interaction of particles in electric field. 1. features of hydrodynamic interaction of polarized particles[END_REF][START_REF] Gamayunov | Study of flows induced in the vicinity of conducting particles by an external electric field[END_REF][START_REF] Murtsovkin | Nonlinear flows near polarized disperse particles[END_REF][START_REF] Murtsovkin | Steady flows in the neighborhood of a drop of mercury with the application of a variable external electric field[END_REF][START_REF] Barinova | Electroosmosis at spherical and cylindrical metal surfaces[END_REF]. The behavior of metal particles strongly depends on the external field. In weak fields, particles behave like dielectric ones. However, in strong fields, these particles behave very differently from the dielectric particles. In the high-intensity electric field, a qualitatively new effect arises: an electric charge in a liquid near a particle is created by an electric current passing through the particle, and not by the surface charge. Such processes in electrochemistry are called non-equilibrium, and Dukhin called this electrokinetic phenomenon of the "second kind".

In this case, according to Dukhin, a quadratic dependence of the particle velocity on E ∞ appears, as well as the dependence on the particle radius a. The latter is very important for practical applications, as it make possible to sort particles by size. There is currently no mathematical models describing the movement of metal particles.

The third type of particles is ion-selective particles. The electrokinetic flow near an ion-selective particle is associated with the non-equilibrium behavior of the space charge.

The surface of an ion-selective particle allows one type of ion to penetrate through it and hence to create an electric current. At the weak electric field strength, the ion flux through the surface is negligible and therefore a quasi-equilibrium EDL is formed. The electroosmosis associated with the sliding of the electrolyte, which is the result of the tangential component of the electric field on the charge of the quasi-equilibrium EDL, is called the electroosmosis of the "first kind" [START_REF] Lyklema | Fundamentals of Interface and Colloid Science[END_REF], which is similar to the classical electroosmosis for dielectric particles (Fig. 2.2a). With an increase in the electric field strength, the ion flux through the surface of the particle becomes significant and can no longer be neglected.

The structure of the EDL at the interface is transformed from quasi-equilibrium to nonequilibrium. A distinctive feature of this new structure is that an expanded space charge region (SCR) is added to the charge of the quasi-equilibrium EDL. This non-equilibrium electroosmosis associated with the space charge region is called electroosmosis of the "second kind" and is described in the works of Dukhin and Mishchuk [START_REF] Dukhin | Electrokinetic phenomena of the second kind and their applications[END_REF][START_REF] Mishchuk | Electroosmosis of the second kind[END_REF]. It refers to a similar action of a tangential electric field on the charge inside SCR. In accordance with the Newton's third law, electroosmosis ultimately leads to the movement of a microparticle and electrophoresis of the first or second kind, respectively.

Experimental studies of electrophoresis and electroosmosis of an ion-selective particle were performed by Mishchuk and Takhistov [START_REF] Mishchuk | Electroosmosis of the second kind[END_REF], Barinova et al. [START_REF] Barinova | Electroosmosis at spherical and cylindrical metal surfaces[END_REF], Mishchuk and Barinova [START_REF] Mishchuk | Theoretical and experimental study of nonlinear electrophoresis[END_REF], as well as Barany [START_REF] Barany | Electrophoresis in strong electric fields[END_REF]. The main result of these studies is that at the transition from electrophoresis of the first kind to electrophoresis of the second kind, the linear dependence of the electrophoretic velocity on the electric field is transformed to a quadratic one

U ∞ ∼ εa µ E 2 ∞ . (2.3)
In addition, in accordance with the work of Dukhin and Mishchuk, the dependence of the electrophoretic velocity on the particle radius a also arises. In contrast to the case of the dielectric particle, non-linear effects increase the velocity of electrophoresis. The formula (2.3) was derived from dimensional analysis.

The mathematical formulation of the problem near an ion-selective flat surface was carried out by Rubinstein and Shtilman [START_REF] Rubinstein | Voltage against current curves of cation exchange membranes[END_REF]. This statement was successfully used in the works of Rubinstein and Zaltzman [START_REF] Rubinstein | Electro-osmotically induced convection at a permselective membrane[END_REF], Zaltzman and Rubinstein [START_REF] Zaltzman | Electro-osmotic slip and electroconvective instability[END_REF], Demekhin et al. [START_REF] Demekhin | Linear and nonlinear evolution and diffusion layer selection in electrokinetic instability[END_REF][START_REF] Demekhin | Direct numerical simulation of electrokinetic instability and transition to chaotic motion[END_REF][START_REF] Demekhin | Three-dimensional coherent structures of electrokinetic instability[END_REF], Pham et al. [START_REF] Pham | Direct numerical simulation of electroconvective instability and hysteretic current-voltage response of a permselective membrane[END_REF], Druzgalski et al. [START_REF] Druzgalski | Direct numerical simulation of electroconvective instability and hydrodynamic chaos near an ion-selective surface[END_REF] for a theoretical study of electrokinetic phenomena near a flat ion-selective surface. Unlike plane geometry, the case of curved microparticle interface is one of the conditions for the non-linear electrokinetic phenomena emergence [START_REF] Dukhin | Electrokinetic phenomena of the second kind and their applications[END_REF]. The experimental works [START_REF] Mishchuk | Electroosmosis of the second kind[END_REF][START_REF] Mishchuk | Theoretical and experimental study of nonlinear electrophoresis[END_REF][START_REF] Barany | Electrophoresis in strong electric fields[END_REF] are devoted to electrokinetic phenomena that arise near the curved ion-selective surface of a spherical particle. These works also give estimated ratios (for example, for the electrophoretic velocity) for some limiting cases. However, some theoretical aspects of the motion of such particles remain unclear.

Mathematical formulation

The principles of conservation of momentum, mass, and charge are used to define the state of a problem. The conservation laws are applied in an assumption of a continuum fluid.

The electrolyte is assumed to be a dielectric Newtonian fluid in which ions of a certain substance (salt, acid or alkali) are dissociated. An example of such an electrolyte is distilled water with sodium chloride (N aCl) dissolved in it. The N aCl easily dissociates into positively charged sodium ions N a + (cations) and negatively charged chlorine ions

Cl -(anions).

Mass transport in electrolyte: diffusion, convection and electromigration

Suppose that some substance is dissolved in a volume V of the liquid. The closed surface of this volume denotes as S. The concentration c(x, y, z, t) of this substance depends on space and time. Since the flow of substance through the surface in general is non-zero, we have to write the following relation,

S j • n ds + ∂ ∂t V c(x, y, z, t) dv = 0 (2.4)
According to Ostrogradsky-Gauss theorem,

S j • n ds = V ∇ • j dv (2.5) It means, V ∂c ∂t + ∇ • j dv = 0 (2.6)
This equality must be true for any volume V , so the integrand function must be zero,

∂c ∂t + ∇ • j = 0 (2.7)
The ionic flux j can be represented by three mechanisms of mass transfer: diffusion, convection and electro-migration. All of them must be taken into account.

Diffusion

The transport of mass or diffusion takes place in an electrolyte with two or more species if there exists a concentration gradient. Mass diffusion in the electrolyte is a consequence of the Brownian motion of ions.

Consider an electrolyte solution consisted of i ion species. The two basic concentration units -mass concentration and molar concentration can be defined, respectively, by,

ρ i = m i V = mass of species i volume of solution (kg • m -3 ) (2.8) c i = n i V = number of moles of species i volume of solution (mol • m -3 ) (2.9)
In the present work we will use only molar concentration c i .

The diffusion flux of ions arises due to the inhomogeneous distribution of ion con-centrations in the volume of the electrolyte. At constant temperature and pressure, the diffusive flux is related to the concentration through Fick's first law as,

j D = -D i ∇c i (2.10)
where D i is the diffusion coefficient of i-th species with units m 2 • s -1 . In general for i > 2 the system will have several concentration gradients and the diffusive flux of each species can be affected by other concentration gradients. The exception is the case of a dilute solution (or weak solution) for which flux for each species is unaffected by the presence of the other. So (2.10) is valid only for weak solutions.

Convection

The next mechanism of the mass transfer is convection. It is the mass transfer due to the bulk motion of a fluid. The velocity of ions undergoing mass transfer has both diffusive and convective components. The difference between these two is the following:

the convection is a mass transfer due to the average velocity of all ions, and diffusion is due to the instant velocity of a particular ion, compared to the average velocity of the fluid as a whole.

By means of an appropriate choice of the reference velocity, it is possible to define the convective flux. For mass transfer of dilute species, the solvent (water for aqueous solutions) prevails in the momentum of the system and it is natural to take the velocity of the liquid as a reference velocity. Thus, the convective flow is due to ion transport by fluid motion and is proportional to the fluid velocity,

j C = U • c i (2.11)
where U is the velocity of the liquid.

Electromigration

For the electro-migration according to Coulomb's law the electric field E is the negative gradient of the electrostatic potential E = -∇Φ. The force exerting on a particle is the magnitude of the particle charge multiplied by the charge number z i and the electric field E. The force per one mole should be written as -z i N A e ∇Φ, where N A is the number of ions in one mole of substance; e is elementary charge for one ion; F = N A • e is the Faraday's constant.

Electro-migration flux is proportional to the force acting on the particle multiplied by the concentration of ions c i . The proportionality factor ν i is a transport property called the molar mobility of the i th ionic species (ν i = D i /R T ).

j E = -ν i z i F c i ∇Φ (2.12)
Combining all three transport mechanisms,

j = -D i ∇c i -ν i z i F c i ∇Φ + U c i (2.13) ∂c i ∂t + ∇ • j = ∂c i ∂t + ∇ • (-D i ∇c i -ν i z i F c i ∇Φ + U c i ) = 0 (2.14) or ∂c i ∂t + U • ∇c i = D i ∇ 2 c i + z i D i F R T ∇ • (c i ∇Φ) (2.15)
Equation (2.15) is called the Nernst-Planck equation or the ion transport equation.

In the general case, this equation should be supplemented by the source term, which is responsible for the ion flux due to chemical reactions with the generation of new ions.

In this work, we restrict ourselves to the following assumption: chemical reactions with the generation of new ions occur only on the electrodes, and since the case of infinitely distant electrodes is considered, the source term in the volume of the electrolyte is zero.

Flow velocity and electrostatic fields

It is assumed that the electrolyte is symmetric, binary and monovalent, i.e. it contains only two types of ions (with concentrations c + and c -) and the charge numbers (valence) of ions are equal in modulus (z + = -z -= 1), the diffusion coefficients of cations and anions are equal (D + = D -= D). The liquid is regarded as incompressible and has a dynamic viscosity µ and electrical permittivity ε. Then the ion transport in the electrolyte is described by the following equations,

∂c ± ∂t + U • ∇c ± = D ± F RT ∇ • (c ± ∇Φ) + ∇ 2 c ± . (2.16)
Since the problem has an electric field, the system must include the electrostatic part, which is described by the Poisson equation,

ε∇ 2 Φ = -F (c + -c -).
(2.17)

where F (c + -c -) -is the electric charge density.

For the velocity field consider the Navier-Stokes equation describing the motion of a viscous fluid in the general form,

ρ ∂ U ∂t + ( U • ∇) U -∇P + µ∇ 2 U + f = 0 (2.18) ∇ • U = 0, (2.19)
where f is the volume force and P is the pressure.

Let us simplify this equation in accordance with the conditions of the problem. First of all we need to estimate the Reynolds number. In our case, the characteristic length is the microparticle radius. From the experimental work [START_REF] Mishchuk | Electroosmosis of the second kind[END_REF], devoted to electroosmosis of the second kind around ion-selective particles, it follows that the electroosmotic velocity near sufficiently large particle with a radius of 210µm in an electric field strength 1000

V/m reaches 3.5 •10 -4 m/s. Then the Reynolds number is estimated as Re ≈ 0.08, which is a negligible value, and therefore the inertial terms in the Navier-Stokes equation can be neglected and the equation can be taken in the Stokes approximation [START_REF] Schlichting | Boundary Layer Theory[END_REF]. The volume force is the electric force acting on the liquid, which is equal to F (c + -c -)∇Φ [START_REF] Probstein | Physicochemical Hydrodynamics[END_REF]. Thus, the equation for hydrodynamics, taking into account the volumetric electric force, takes the form,

-∇P + µ∇ 2 U = F (c + -c -)∇Φ, (2.20)
∇ • U = 0 (2.21)
The system of equations has to be converted in dimensionless form. To do that the following characteristic quantities were used, a: characteristic length (particle radius); Typical values of the physical quantities and constants of the problem are given in Table 2.1.

When the system becomes dimensionless, two additional parameters appear: ν is the Debye number, which is the ratio of the Debye length λ D to the particle radius a (ν 1 is a small parameter of the problem),

ν = λ D a , λ D = εΦ 0 F c ∞ 1 2 , Parameter Symbol Value Faraday constant F 96485.33 A • s/mole Universal gas constant R 8.314 m 2 • kg/s 2 • K • mole Absolute temperature T 298.15 K Thermal potential Φ 0 2.569 ×10 -2 m 2 • kg/s 3 • A Particle radius a 5 ×10 -6 m Ion diffusion coefficient D 1.994 ×10 -9 m 2 /s Dynamic viscosity µ 8.93 ×10 -4 kg/m • s Absolute permittivity ε 7.08 ×10 -10 s 4 • A 2 / m 3 • kg Concentration c ∞ 0.1 mole/m 3 Debye length λ D 4.34 •10 -8 m Electric field strength E ∞ 10 2 -10 4 V/m Table 2.
1 -Typical values of the physical quantities and constants of the problem. The ion diffusion coefficient and concentration is given for sodium chloride (NaCl).

and κ -parameter connecting the hydrodynamic and electrostatic parts of the problem.

κ = εΦ 2 0 µD .
The parameter κ characterizes the physical properties of the electrolyte solution and does not change for a fixed electrolyte.

Finally, the system of dimensionless equations have the following form,

∂c ± ∂t + U • ∇c ± = ±∇ • (c ± ∇Φ) + ∇ 2 c ± , (2.22)
ν 2 ∇ 2 Φ = c --c + , (2.23)
-∇P + ∇ 2 U = (c + -c -) κ ν 2 ∇Φ = -κ∇ 2 Φ∇Φ, (2.24)
∇ • U = 0. (2.25)
Hereinafter, all values will be given in dimensionless form, unless otherwise indicated.

Since it is assumed that the microparticle has a spherical shape, it is convenient to carry out our calculations in a spherical coordinate system, the origin of which is at the particle center, i.e. the motion is considered in the coordinate system associated with the particle. The problem statement for spherical coordinates (r, θ) in the axisymmetric case (no dependence on the azimuthal angle) is written as follows,

∂c + ∂t + U 1 r ∂c + ∂θ + V ∂c + ∂r = 1 r 2 sin θ ∂ ∂θ sin θc + ∂Φ ∂θ + 1 r 2 ∂ ∂r r 2 c + ∂Φ ∂r + + 1 r 2 sin θ ∂ ∂θ sin θ ∂c + ∂θ + 1 r 2 ∂ ∂r r 2 ∂c + ∂r , (2.26
)

∂c - ∂t + U 1 r ∂c - ∂θ + V ∂c - ∂r = - 1 r 2 sin θ ∂ ∂θ sin θc -∂Φ ∂θ + 1 r 2 ∂ ∂r r 2 c -∂Φ ∂r + + 1 r 2 sin θ ∂ ∂θ sin θ ∂c - ∂θ + 1 r 2 ∂ ∂r r 2 ∂c - ∂r , (2.27)
ν 2 r 2 ∂ ∂r r 2 ∂φ ∂r + 1 sin θ ∂ ∂θ sin θ ∂φ ∂θ = c --c + ;
(2.28)

- 1 r ∂P ∂θ + ∂ 2 U ∂r 2 + 2 r ∂U ∂r + 1 r 2 ∂ 2 U ∂θ 2 + cot θ r 2 ∂U ∂θ - u r 2 sin 2 θ + 2 r 2 ∂V ∂θ = = (c + -c -) 1 r κ ν 2 ∂Φ ∂θ ; (2.29)
- ∂P ∂r + ∂ 2 V ∂r 2 + 2 r ∂V ∂r + 1 r 2 ∂ 2 V ∂θ 2 - 2V r 2 + cot θ r 2 ∂V ∂θ - 2U r 2 cot θ - 2 r 2 ∂U ∂θ = = (c + -c -) 1 r κ ν 2 ∂Φ ∂r .
(2.30)

∂ ∂θ (sin θ r U ) + ∂ ∂r sin θ r 2 V = 0. (2.31)
To completely close the mathematical formulation, the system of equations (2.22)-(2.25) must be supplemented with boundary conditions.

Boundary conditions 2.4.1 Boundary conditions at solid/liquid interface

For definiteness, we assume that the particle is cation exchange, that is, it does not pass anions. On the surface of the particle with r = 1, the boundary conditions are:

r = 1 : c + = p; c -∂Φ ∂r - ∂c - ∂r = 0; Φ = 0; U = 0. (2.32)
The first BC in Eq. (2.32) was introduced in the paper by Rubinstein and Shtilman [START_REF] Rubinstein | Voltage against current curves of cation exchange membranes[END_REF] (see also Rubinstein and Zaltzman [START_REF] Zaltzman | Electro-osmotic slip and electroconvective instability[END_REF], Zaltzman and Rubinstein [START_REF] Rubinstein | Electro-osmotically induced convection at a permselective membrane[END_REF], and Chapter 6 in the book by Probstein [START_REF] Probstein | Physicochemical Hydrodynamics[END_REF]). In order to better understand this BC, let us consider the structure of the membrane. The cation-exchange membrane is an organic polymer consisting of a matrix and pores. In the matrix, the anions (c a ) are fixed and immobile [START_REF] Zabolotsky | Ion transport in membranes[END_REF], which create a fixed membrane charge p (c a = p). When the membrane is placed in the electrolyte without any electrical field, the pores are filled with this electrolyte and ions of the opposite sign (c + ) get accumulated. Moreover, their number is practically equal to the charge of the membrane (c + = p), that is, the membrane, as a whole, is screened from the inside. If the membrane charge is large enough (p 1), it is more difficult for external forces to change the amount of cations inside the membrane and c + = p is assumed inside the membrane and at its surface with the electrolyte (see Eq.

(2.32)). Studies of planar membranes [START_REF] Rubinstein | Electro-osmotically induced convection at a permselective membrane[END_REF][START_REF] Zaltzman | Electro-osmotic slip and electroconvective instability[END_REF][START_REF] Demekhin | Linear and nonlinear evolution and diffusion layer selection in electrokinetic instability[END_REF][START_REF] Demekhin | Direct numerical simulation of electrokinetic instability and transition to chaotic motion[END_REF][START_REF] Demekhin | Three-dimensional coherent structures of electrokinetic instability[END_REF][START_REF] Pham | Direct numerical simulation of electroconvective instability and hysteretic current-voltage response of a permselective membrane[END_REF][START_REF] Druzgalski | Direct numerical simulation of electroconvective instability and hydrodynamic chaos near an ion-selective surface[END_REF] show that the solution is practically independent of the value of p for p 3 (see also [START_REF] Frants | Nonequilibrium electrophoresis of an ion-selective microgranule for weak and moderate external electric fields[END_REF]).

The second boundary condition in (2.32) means the absence of anions through the particle surface. The third condition defines a constant potential on the surface, which can be set to zero without loss of generality. The last condition in (2.32) is the no-slip condition a solid surface, i.e. flow velocity is equal to zero.

Note that since positively charged ions can pass through the surface of a particle, their flow creates an electric current through the surface, r = 1 :

j = c + ∂Φ ∂r + ∂c + ∂r .
(2.33)

We suppose the membrane to be ideal, i.e. it is able to withstand the accumulation of ions inside. Therefore the current density is zero over the surface, π 0 j sin θdθ = 0.

(2.34)

Boundary conditions away from the interface

Far from the particle, at r → ∞, the electrolyte solution is electrically neutral, and the salt concentration is in equilibrium. The electric field tends to the external value E ∞ , the fluid velocity tends to the electrophoretic velocity U ∞ (see 2.3),

c ± → 1; ∂Φ ∂r → -E ∞ cos θ; U → -U ∞ sin θ; V → U ∞ cos θ. (2.35) x y E ∞ U θ r = 1 V ∞ U θ = 180 o θ = 0 o Figure 2.
3 -Schematic of the flow near the granule under the external electric field, E ∞ and the fluid flow velocity U ∞ at infinity which are directed along the xcoordinate. The electrophoretic velocity of the particle has an opposite direction. In the spherical polar system, x = r cos θ and y = r sin θ, U = U θ is the tangential velocity at the particle surface and V = U r is the normal velocity at the particle surface. At the co-moving reference frame, the far-field velocity condition is U → -U ∞ sin θ and V → U ∞ cos θ.

The system is closed by imposing conditions on the concentration at the initial time:

t = 0 : c + = c -= 1 (2.36)
Chapter 3

Numerical methodology

The solution of the system is sought on a non-uniform staggered grid using the finite difference method. The scalar values like electrostatic potential Φ and the ion concentrations c ± are taken at the centers of the grid cells, and the vector values -the velocity components U and V , in the middle of their lateral faces. Derivatives are calculated using second-order two-point difference schemes with an offset. The numerical method is based on the works of Nikitin [START_REF] Nikitin | Third-order-accurate semi-implicit runge-kutta scheme for incompressible navier-stokes equations[END_REF][START_REF] Nikitin | Finite-difference method for incompressible navier-stokes equations in arbitrary orthogonal curvilinear coordinates[END_REF][START_REF] Nikitin | On the rate of spatial predictability in near-wall turbulence[END_REF][START_REF] Nikitin | Statistical characteristics of wall turbulence[END_REF][START_REF] Nikitin | A spectral finite-difference method of calculating turbulent flows of an incompressible fluid in pipes and channels[END_REF].

Discretization

Let us give an example of a program code that defines a grid. Consider a grid of Im cells along the variable r = 1 ÷ R max and Jm cells along θ = 0 ÷ π. Cell nodes are numbered from i = 0 to Im and from j = 0 to Jm, respectively. Further we will call them "nodes" for simplicity. The centers of the cells correspond to "center nodes", which in the code are numbered from 1 to Im and Jm, respectively.

The stretching of the mesh is set, respectively, by the functions ξ and ζ so that the nodes are

r i = rn(i) = ξ i Im and θ j = tn(j) = ζ j Jm . Accordingly, ξ(0) = 1, ξ(1) = R max , ζ(0) = 0, ζ(1) = π. For center nodes, the relations take the form r i-1/2 = rm(i) = ξ i-1/2 Im and θ j-1/2 = tm(j) = ζ j-1/2 Jm .
The indices i and j actually define a uniform grid in the new variables x and z: r = ξ(x) and θ = ζ(z). Numerical differentiation of scalar functions is carried out as follows,

∂c ∂r = ∂c ∂x dξ dx = (c(i + 1, j) -c(i, j)) /Im / dξ dx + O 1 Im 2 ,
where the indices i and j correspond to center nodes, and the value of the derivative itself is determined in the node with the index i. Array indexing is hereinafter specified by parentheses. The value Im/ dξ dx taken in the same center node i is precomputed and stored in the array element rn1(i) (The number 1 hereinafter means that the array stores the distance between two adjacent points since the grid is irregular). To differentiate by θ, a similar array tn1 is created. For velocity functions, the difference scheme changes slightly,

∂V ∂r = (V(i, j) -V(i -1, j)) /rm1(i) + O 1 Im 2 ,
where the i index for the V array now corresponds to the node, and for the rm1 arrayto the center node. For θ, the tm1 array is created in the same way.

It should be noted that linear interpolation is used to shift the functions themselves to another grid. Thus, the concentration values c(i, j) are defined as (c(i + 1, j) + c(i, j)) /2, moreover, this is the second order approximation.

For example, the first derivative of the function defined in the nodes:

f (i) = 1 rm1(i) f(i) -f(i -1) h r
The first derivative of the function defined in center nodes (i -1/2):

g (i) = 1 rm1(i) g(i + 1) -g(i) h r

Numerical method for Poisson equation

Poisson equation in spherical coordinates has the form,

∂ ∂r r 2 ∂Φ ∂r + 1 sin θ ∂ ∂θ sin θ ∂Φ ∂θ = r 2 ν 2 (c --c + ) (3.1)
In this equation first term depends only on r and the second -only on θ, therefore we can write,

a i Φ i-1,j + b i Φ i,j + c i Φ i+1,j + d j Φ i,j-1 + e j Φ i,j + f j Φ i,j+1 = g i,j (3.2)
where Φ i,j = Φ(r i , θ j ) a i , b i , c i -coefficients depending on r, and d j , e j , f j -coefficients depending on θ. The details of the calculation of the coefficients are presented in Appendix A.

Independence of the coefficients a(i), b(i), c(i) from the angle θ is violated only in the last layer of nodes at i = Im, because the boundary condition Φ = -E ∞ r cos θ is dependent on θ. In order to eliminate this feature we make the following substitution,

Φ = Φ cos θ , ⇒ ∂ Φ ∂r Rmax = -E ∞ r (3.3)
and we solve the equation with respect to Φ. In both parts of equation, we multiply by 1 cos θ and we get,

1 cos θ ∂ ∂r r 2 ∂( Φ cos θ) ∂r + 1 cos θ sin θ ∂ ∂θ sin θ ∂( Φ cos θ) ∂θ = r 2 ν 2 cos θ (c --c + ) (3.4) ∂ ∂r r 2 ∂ Φ ∂r + 1 cos θ sin θ ∂ ∂θ sin θ ∂( Φ cos θ) ∂θ = r 2 ν 2 cos θ (c --c + ) (3.5)
with the boundary conditions,

Φ r=1 = 0, ∂ Φ ∂r Rmax = -E ∞ r (3.6)
This equation is solved using spectral decomposition (in terms of eigenvectors and eigenvalues). Let us introduce the notation for the operators regarding variables r and θ,

T r = ∂ ∂r r 2 ∂ ∂r (3.7) T θ = 1 sin θ ∂ ∂θ sin θ ∂ ∂θ (3.8)
then we can write the Poisson equation in the form,

(T r + T θ )[ Φ] = R (3.9)
where R is the right part of Poisson equation.

We need to simplify one of these operators to the diagonal form. After that, we will have a system of equations with the second operator. For this purpose, we use the spectral decomposition of the matrix of the operator T θ ,

T θ = X Φ Λ Φ X -1 Φ (3.10)
where Λ -diagonal matrix, the elements of which are the eigenvalues of T θ . They need to be real. The column vectors of X are the eigenvectors of T θ .

T r [ Φ] + T θ [ Φ] = R (3.11)
using relation (3.10), we get,

T r [ Φ] + X Φ Λ Φ X -1 Φ [ Φ] = R (3.12)
Multiply this equation by X -1 Φ from the right side and we get,

T r [X -1 Φ Φ] + Λ Φ [X -1 Φ Φ] = X -1 Φ R (3.13)
which leads to,

(T r + Λ Φ )[X -1 Φ Φ] = X -1 Φ R (3.14)
This equation is solved with respect to X -1 Φ and then this solution is multiplied by X and cos θ to find the final solution for Φ. The details of the calculation of matrix for operator (T r + Λ Φ ) see in Appendix A.

Numerical method for Stokes equation

For the Stokes equation,

-∇P + ∇ 2 U = κ ν 2 (c + -c -)∇Φ
we will use the following property,

curl curl U = ∇(∇ • U ) -∇ 2 U
Let us express the second term from the Stokes equation,

∇ 2 U = ∇(∇ • U ) -curl curl U
and take into account that from the continuity equation

∇ • U = 0. Denote ω = curl U .
Then,

-∇P -curl ω = κ ν 2 (c + -c -)∇Φ Consider the general case curl of U = (u 1 , u 2 , u 3 ), ∇ • U and ∇P in arbitrary curvi- linear orthogonal coordinates, ω = curl U = = h 1 h ∂(h 3 u 3 ) ∂x 2 - ∂(h 2 u 2 ) ∂x 3 ; h 2 h ∂(h 1 u 1 ) ∂x 3 - ∂(h 3 u 3 ) ∂x 1 ; h 3 h ∂(h 2 u 2 ) ∂x 1 - ∂(h 1 u 1 ) ∂x 2 ∇ • U = 1 h ∂ ∂x 1 hu 1 h 1 + ∂ ∂x 2 hu 2 h 2 + ∂ ∂x 3 hu 3 h 3 ∇P = 1 h 1 ∂P ∂x 1 ; 1 h 2 ∂P ∂x 2 ; 1 h 3 ∂P ∂x 3 ; where h = h 1 • h 2 • h 3 ; h 1 , h 2 , h 3
are the Lame coefficients. For the spherical coordinates the Lame coefficients are:

h 1 = h r = 1, h 2 = h θ = r, h 3 = h ϕ = r sin θ, and h = h r h θ h ϕ =
r 2 sin θ. For our case of axisymmetrical flow, azimuthal velocity component u 3 = 0, and there is no dependence on angle ϕ. Denote U = u 1 and V = v 1 , then,

ω = 0; 0; r sin θ r 2 sin θ ∂(rV ) ∂r - ∂U ∂θ = 0; 0; 1 r ∂(rV ) ∂r - ∂U ∂θ ∇ • U = 1 r 2 sin θ ∂ ∂r r 2 sin θ • U + ∂ ∂θ (r sin θ • V ) ∇P = ∂P ∂r ; 1 r ∂P ∂θ ; 0
Since ω 1 = ω 2 = 0 we can rewrite curl ω in the following way,

curl ω = h r h ∂(h ϕ ω 3 ) ∂θ ; - h θ h ∂(h ϕ ω 3 ) ∂r ; 0
or substituting the Lame coefficients,

curl ω = 1 r sin θ ∂(sin θω 3 ) ∂θ ; - 1 r ∂(rω 3 ) ∂r ; 0
Let us omit the index 3 for ω 3 and substitute curl, divergence and gradient into Stokes and continuity equations. Then one can get the following system,

- ∂P ∂r - 1 r sin θ ∂(sin θ ω) ∂θ = κ ν 2 (c + -c -) ∂Φ ∂r (3.15) - 1 r ∂P ∂θ + 1 r ∂(rω) ∂r = κ ν 2 (c + -c -) 1 r ∂Φ ∂θ (3.16) ω = 1 r ∂(rv) ∂r - ∂u ∂θ 1 r 2 ∂(r 2 u) ∂r + 1 r sin θ ∂ ∂θ (sin θ • v) = 0 (3.17)
Denote the right parts of Stokes equation as,

F r = κ ν 2 (c + -c -) ∂Φ ∂r ; F θ = κ ν 2 (c + -c -) 1 r ∂Φ ∂θ
∂ 2 P ∂r∂θ + 1 r 2 ∂ ∂θ 1 sin θ ∂(sin θω) ∂θ = - 1 r ∂F r ∂θ (3.18) - 1 r ∂ 2 P ∂r∂θ + 1 r ∂ ∂r ∂(rω) ∂r = 1 r ∂(rF θ ) ∂r (3.19)
- 1 r sin θ ∂ 4 Ψ ∂r 4 - 1 r ∂ 2 ∂r 2 1 r 2 ∂ ∂θ 1 sin θ ∂Ψ ∂θ + 1 r 3 ∂ ∂θ 1 sin θ ∂ ∂θ ∂ 2 Ψ ∂r 2 - - 1 r 5
∂ ∂θ

1 sin θ ∂ ∂θ sin θ ∂ ∂θ 1 sin θ ∂Ψ ∂θ = Ω (3.21)
where

Ω = 1 r ∂(rF θ ) ∂r - 1 r ∂F r ∂θ
. For this equation it is possible to separate variables r and θ by multiplying Eq. (3.21) by r 5 sin θ, which leads to,

r 4 ∂ 4 Ψ ∂r 4 + sin θ ∂ ∂θ 1 sin θ ∂ ∂θ r 4 ∂ 2 ∂r 2 1 r 2 + r 2 ∂ 2 ∂r 2 Ψ+ + sin θ ∂ ∂θ 1 sin θ ∂ ∂θ sin θ ∂ ∂θ 1 sin θ ∂ ∂θ Ψ = -Ω r 5 sin θ
For the numerical solution we use the same method as for Poisson equation -the spectral decomposition. We end up with the same operator T θ regarding the angle as for Poisson equation,

T θ = sin θ ∂ ∂θ 1 sin θ ∂Ψ ∂θ then, r 4 ∂ 4 Ψ ∂r 4 + r 4 ∂ 2 ∂r 2 1 r 2 + r 2 ∂ 2 ∂r 2 T Ψ + T 2 Ψ = R (3.22)
where T θ = XΛX -1 , Λ is a diagonal matrix. By multiplying Eq. (3.22) by X -1 leads to,

r 4 ∂ 4 ∂r 4 (X -1 Ψ) + r 4 ∂ 2 ∂r 2 1 r 2 + r 2 ∂ 2 ∂r 2 Λ(X -1 Ψ) + Λ 2 (X -1 Ψ) = X -1 R (3.23)
The discretization of Eq. (3.23) is detailed in Appendix B.

Time advancement (Nernst-Plank equation)

The spatial discretization with sufficiently fine resolution leads to stiff problems and requires implicit methods for time advancement. Semi-implicit methods, where only a part of the Nernst-Plank operator is treated implicitly, are most efficient.

To solve the Nernst-Planck equation, a semi-implicit Runge-Kutta scheme of the third order of accuracy [START_REF] Nikitin | Third-order-accurate semi-implicit runge-kutta scheme for incompressible navier-stokes equations[END_REF][START_REF] Nikitin | Finite-difference method for incompressible navier-stokes equations in arbitrary orthogonal curvilinear coordinates[END_REF] was compiled. An implicit operator is a parameter to this method. After each step, its error is estimated and the step size is adjusted.

dw dt = F (t, w)
where w(t) is an unknown vector function, t is an independent variable (time), and F is a non-linear operator. A semi-implicit method to advance from w n at time t n to w n+1 at time t n+1 = t n + τ is based on the following explicit third-order accurate Runge-Kutta method,

w * -w n τ = 2 3 F n (3.24) w * * -w n τ = 1 3 F n + 1 3 F * (3.25) w n+1 -w n τ = 1 4 F n + 3 4 F * * (3.26)
where F n ≡ F (t n ; w n ), F * ≡ F (t n + 2τ /3; w * ), and F * * ≡ F (t n + 2τ /3; w * * ).

This particular method is selected from a 2-parameter family of the third-order accurate Runge-Kutta methods according to some requirements. First, the third step of the method (3.26) does not include F * . This is essential for the following construction of the semi-implicit scheme. Second, the method has the coincident abscissas (F * and F * * are evaluated for the same time moment t n + 2τ /3).

The semi-implicit scheme is constructed by perturbing (3.24)-(3.26) taking into account the following relations,

w * = w(t n + 2τ /3) + O(τ 2 ) w * * = w(t n + 2τ /3) + O(τ 3 ) w n+1 = w(t n+1 ) + O(τ 4 )
We assume, that the implicit terms are linear and denote the corresponding operator L (implicit operator, some approximation to the Jacobian ∂F/∂w) and γ -positive number.

Consider the following implicit scheme,

w * -w n τ = 2 3 F n + γ L(w * -w n ) (3.27) w * * -w n τ = 1 3 F n + 1 3 F * + γ L(w * * -w * ) (3.28) w n+1 -w n τ = 1 4 F n + 3 4 F * * + γ L(w n+1 -wn+1 ) (3.29)
where wn+1 is some O(τ 2 ) approximation to w(t n+1 ) (i.e. wn+1 = w(t n+1 ) + O(τ 3 )).

It is easy to see, that the right-hand side of (3. Thus, the local error of the scheme w n+1 -w(t n+1 ) is still O(τ 4 ), so that the scheme retains the third order of accuracy. Note that the absence of F * in the right-hand side of (3.29) is a necessary condition for this.

The specific form of the implicit operator does not affect the order of accuracy of the scheme but determines its stability properties. In the trivial case of L ≡ 0, the scheme turns into a simple third-order accurate explicit Runge-Kutta scheme. Generally speaking, the stability is the higher, the closer L is to the Jacobian of the Nernst-Plank

operator. An important criterion for choosing L is how effectively may be solved the corresponding set of Eq. (3.27)-(3.29).

The Nernst-Plank equation in the spherical coordinate system has the following form,

∂c ± ∂t + U 1 r ∂c ± ∂θ + V ∂c ± ∂r = 1 r 2 1 sin θ ∂ ∂θ sin θ ∂c ± ∂θ ± c ± ∂Φ ∂θ + + 1 r 2 ∂ ∂r r 2 ∂c ± ∂r ± r 2 c ± ∂Φ ∂r
We suppose that F (w) can be split in two operators F (w) = G(w) + H(w), where operator H(w) creates instability,

w 1 -w 0 ∆t = G(w 0 ) + H(w 1 )
w 1 -∆tH(w 1 ) = w 0 + ∆tG(w 0 )

and the system w 1 -∆tH(w 1 ) can be easily solved. The simplest case is when H(w) = A w with a band matrix A. If H(w) = A w then (I -∆tA)w 1 = w 0 + ∆tG(w 0 ).

The main problem here is to choose the right operator, which will be treated implicitly because it has a strong influence both on the stability and on the time step.

The proposed Runge-Kutta scheme has 4 steps. At each step will be solved the linear system associated with the implicit operator and will calculate the values of Φ, U , V .

We don't describe here in detail the formation of the operator of the Nernst-Planck equation but indicate which part of this operator is treated implicitly. The following terms of Nernst-Planck equation is treated implicitly,

U 1 r ∂c ± ∂θ ; V ∂c ± ∂r ; 1 r 2 1 sin θ ∂ ∂θ ± sin θc ± ∂Φ ∂θ ; 1 r 2 ∂ ∂r r 2 ∂c ± ∂r ± r 2 c ± ∂Φ ∂r
The matrix of this operator is block-tridiagonal; therefore, the solution of the corresponding system is quite simple. Nevertheless, this part of the algorithm takes up to 80% of the whole calculation time.

Chapter 4

Electrophoresis of an ion-selective particle in a weak electric field It was experimentally observed that at low electric field strength E ∞ , the ion-selective particle behaves qualitatively similar to dielectric impermeable particle and its electrophoretic velocity is proportional to the external electric field strength (U ∞ ∼ E ∞ ). In this thesis, the analysis of the ion-selective particle behavior begins with asymptotic analysis for the case of the weak electric field when E ∞ << 1. The expansion of unknown functions in a series with regards to the small parameter E ∞ will be regular.

f = f 0 + E ∞ f 1 + E 2 ∞ f 2 + ...
At the same time, there is another small parameter in the problem -the Debye number ν. Moreover, the analysis at ν << 1 leads to a singular decomposition since ν is in front of the highest derivative in the Poisson equation (2.23). Physically, this is associated with the formation of a thin EDL. The singular decomposition implies the separation of the solution region into two: the inner (corresponding to a thin layer in which the functions change dramatically) and the outer one (corresponding to the bulk of electrolyte), followed by the matching of the inner and outer solutions [START_REF] Van Dyke | Perturbation Methods in Fluid Mechanics[END_REF].

Let us consider in more detail the mechanism of the formation of a thin boundary layer arising near a ion-selective particle. Since the cation selective particle is a membrane with the deficit of cations in its structure, it has a negative charge in total (due to deficit of cations and excess of anions). When this particle is placed in an electrolyte, which contains positive and negative ions, the positively charged ions are attracted to the particle surface in order to compensate for the net negative charge inside the particle. These ions gather around the particle and a well-ordered and immovable layer covers the surface of the particle. This layer is called the stationary layer or the Stern layer. Cations in this layer neutralize some part of the particle negative charge. However as ions are often surrounded by water molecules (in the case of N aCl/water solution); these are rather big and cannot neutralize the surface charge completely. The remaining charge attracts additional ions from the bulk of the electrolyte so that a second layer forms around the particle. This layer is located farther from the surface of the particle. The Coulomb force becomes weaker with distance from the particle and therefore, the second layer, called the diffuse or Gouy-Chapman layer, is much less ordered and mobile. These two layers together form the Electric Double Layer (EDL). At the interface between the liquid and solid, EDL is always present, even in the absence of an external applied electric field (Fig. 4.1). This structure of EDL is in equilibrium since it is not related to the passing of electric current through the EDL. When an external electric field is applied, the structure of EDL is changed. In a weak electric field, EDL does not undergo significant qualitative changes, since the current through the particle surface is negligible.

The existence of two small parameters requires to define the ratio between them. In our case, the relation ν/E ∞ << 1 naturally arises, which allows us, at the beginning, to consider the regular expansion regarding E ∞ and to analyze separately the case ν << 1 for each approximation.

The purpose of the asymptotic study is to describe the stationary regime. Therefore, in this section, we consider the case ∂ ∂t = 0. We perform an asymptotic expansion for ν → 0 (small Debye numbers) for a weak electric field, E ∞ → 0.

The zeroth-order steady-state approximation with respect to E ∞

In the absence of an external electric field (E ∞ = 0), there is no particle motion and the electric current through the particle surface is zero. Therefore, the decomposition of the main functions can be represented as,

c ± = c ± 0 + E ∞ c ± 1 + O(E 2 ∞ ), Φ = Φ 0 + E ∞ Φ 1 + O(E 2 ∞ ), (4.1)
U = E ∞ U 1 + O(E 2 ∞ ), V = E ∞ V 1 + O(E 2 ∞ ), j = E ∞ j 1 + O(E 2 ∞ ). (4.2)
The particle velocity U ∞ is also an unknown parameter and, according to experimental observations, it can be assumed that U ∞ = O(E ∞ ) for E ∞ << 1. This velocity will be found at the end of the present analysis from the balance of forces acting on the particle.

Using the expansions (4.1) -(4.2) after substituting them in (2.22) -(2.25), it is possible to solve the system one by one starting from the zeroth-order approximation and then moving on to terms of a higher-order with respect to E ∞ . The main result of this section is to look for the particle velocity U ∞ , so for this purpose, it will be enough to consider the zero and first-order approximations.

E ∞ u 1 • ∇(c ± 0 + E ∞ c ± 1 + O(E 2 ∞ )) = (4.3) = ±∇ (c ± 0 + E ∞ c ± 1 + O(E 2 ∞ ))∇(Φ 0 + E ∞ Φ 1 + O(E 2 ∞ )) + ∇ 2 (c ± 0 + E ∞ c ± 1 + O(E 2 ∞ )) ν 2 ∇ 2 (Φ 0 + E ∞ Φ 1 + O(E 2 ∞ )) = c - 0 + E ∞ c - 1 + O(E 2 ∞ ) -c + 0 -E ∞ c + 1 -O(E 2 ∞ ) (4.4)
For the zeroth-order approximation, only terms of the zeroth-order with respect to E ∞ are taking into account.

The equations (4.3) -(4.4) with the corresponding boundary conditions in the zerothorder steady-state approximation take the form,

± ∇ • (c ± 0 ∇Φ 0 ) + ∇ 2 c ± 0 = 0 (4.5) ν 2 ∇ 2 Φ 0 = c - 0 -c + 0 (4.6) r = 1 : c + 0 = p, c - 0 ∂Φ 0 ∂r - ∂c - 0 ∂r = 0, Φ 0 = 0, (4.7) r → ∞ : c + 0 → 1, c - 0 → 1, ∂Φ 0 ∂r → 0. (4.8)
The system of equations above do not contain the applied electric field E ∞ , therefore ∂Φ 0 /∂θ = 0. After applying the Laplace operators and divergence in spherical coordinates and in steady state,

1 r 2 ∂ ∂r ±r 2 c ± 0 ∂Φ 0 ∂r + 1 r 2 ∂ ∂r r 2 ∂c ± 0 ∂r = 0 (4.9) ν 2 1 r 2 ∂ ∂r r 2 ∂Φ 0 ∂r = c - 0 -c + 0 (4.10)
with the same boundary conditions (4.7) and (4.8).

The solution of this system will be carried out based on a singular asymptotic approach for ν << 1. And, as it was mentioned earlier, it will require the separation of regions into inner and outer regions. General ideas of the method are presented in [START_REF] Van Dyke | Perturbation Methods in Fluid Mechanics[END_REF].

Inner expansion (ν → 0)

In order to find the inner solution, it is necessary to stretch the thin EDL region near the granule surface r = 1 by the change of variable, z = r -1 ν . Then,

∂ ∂z = 1 ν ∂ ∂r (4.11)
The system for the inner expansion, where r ≈ 1, is as follows,

∂ ∂z c + 0 ∂Φ 0 ∂z + ∂c + 0 ∂z = 0, (4.12)
∂ ∂z -c - 0 ∂Φ 0 ∂z + ∂c - 0 ∂z = 0, (4.13)
∂ 2 Φ 0 ∂z 2 = c - 0 -c + 0 , (4.14)
z = 0 : Φ 0 = 0, c + 0 = p, -c - 0 ∂Φ 0 ∂z + ∂c - 0 ∂z = 0. (4.15)
Here and in what follows, the upper bar denotes the inner solution (solution inside the EDL). Thus for this inner solution, only the boundary condition at z = 0 can be used.

In the zeroth-order approximation with respect to ν there is no influence of the spherical geometry (no dependence on the angle θ) and the system (4.12) -(4.15) can be solved or the solution for the flat geometry [START_REF] Zaltzman | Electro-osmotic slip and electroconvective instability[END_REF] can be used:

Φ 0 = ζ + 2 ln e -ζ/2 + 1 + (e -ζ/2 -1)e -z √ K ∞ 0 e -ζ/2 + 1 -(e -ζ/2 -1)e -z √ K ∞ 0 , (4.16
)

∂Φ 0 ∂z = 4 √ 2 e -z √ K ∞ 0 (e -ζ -1) e -2z √ K ∞ 0 (e -ζ/2 -1) 2 -(e -ζ/2 + 1) 2
, (4.17)

c + 0 = K ∞ 0 2 e -ζ/2 + 1 -(e -ζ/2 -1)e -z √ K ∞ 0 e -ζ/2 + 1 + (e -ζ/2 -1)e -z √ K ∞ 0 2 , (4.18)
c - 0 = K ∞ 0 2 e -ζ/2 + 1 + (e -ζ/2 -1)e -z √ K ∞ 0 e -ζ/2 + 1 -(e -ζ/2 -1)e -z √ K ∞ 0 2 , (4.19)
the K ∞ 0 and ζ should be taken at the outer edge at z → ∞:

K ∞ 0 = c + 0 + c - 0 , ζ = Φ 0 .
The unknown constants must be defined by matching inner and outer solutions. The idea of matching is that the inner and outer solutions should agree on an intermediate Therefore the matching condition can be written as,

K 0 | r=1 = K 0 z→∞ = K ∞ 0 (4.21) Φ 0 | r=1 = Φ 0 z→∞ = ζ (4.22) ρ 0 | r=1 = ρ 0 | z→∞ = 0 (4.23)

Outer expansion (ν ≡ 0)

Since ν → 0, the first approach is to make the approximation ν ≡ 0, and hence find the

solution of the problem, 1 r 2 ∂ ∂r r 2 c + 0 ∂Φ 0 ∂r + r 2 ∂c + 0 ∂r = 0 (4.24) 1 r 2 ∂ ∂r -r 2 c - 0 ∂Φ 0 ∂r + r 2 ∂c - 0 ∂r = 0 (4.25) c - 0 -c + 0 = 0 (4.26)
Adding and subtracting Eq. (4.24) and (4.25), we obtain,

1 r 2 ∂ ∂r r 2 (c + 0 -c - 0 ) ∂Φ 0 ∂r + r 2 ∂(c + 0 + c - 0 ) ∂r = 0 (4.27) 1 r 2 ∂ ∂r r 2 (c + 0 + c - 0 ) ∂Φ 0 ∂r + r 2 ∂(c + 0 -c - 0 ) ∂r = 0 (4.28)
The change of variables has been made for the convenience of solving the problem

c + 0 -c - 0 = ρ 0 and c + 0 + c - 0 = K 0 , then, 1 r 2 ∂ ∂r r 2 ρ 0 ∂Φ 0 ∂r + r 2 ∂K 0 ∂r = 0 (4.29) 1 r 2 ∂ ∂r r 2 K 0 ∂Φ 0 ∂r + r 2 ∂ρ 0 ∂r = 0 (4.30)
From the Eq. (4.26) it follows that ρ 0 ≡ 0. Physically that means absence of charge in the bulk of electrolyte, i.e. electroneutrality,

∂ ∂r r 2 ∂K 0 ∂r = 0 (4.31) ∂ ∂r r 2 K 0 ∂Φ 0 ∂r = 0 (4.32)
The system is complemented by the following boundary conditions:

r → ∞ : K 0 → 2, ∂Φ 0 ∂r → 0 (4.33)
From the equation (4.31) function

K 0 = c 1 r + c 2 .
From the solution for inner region and matching conditions (4.21) -(4.23) one can get the following zeroth-order solution:

ρ 0 ≡ 0, K 0 ≡ 2, Φ 0 ≡ ζ, ζ = ln p (4.34)

The first-order approximation with respect to E ∞

The first-order approximation with respect to E ∞ is described by the following system,

u 1 • ∇c ± 0 = ±∇ • (c ± 0 ∇Φ 1 + c ± 1 ∇Φ 0 ) + ∇ 2 c ± 1 , (4.35)
ν 2 ∇ 2 Φ 1 = c - 1 -c + 1 , (4.36)
-∇P 1 + ∇ 2 u 1 = -κ(∇ 2 Φ 0 ∇Φ 1 + ∇ 2 Φ 1 ∇Φ 0 ), ∇ • u 1 = 0 (4.37)
with the boundary conditions:

r = 1 : c + 1 = 0, c - 1 ∂Φ 0 ∂r + c - 0 ∂Φ 1 ∂r - ∂c - 1 ∂r = 0, Φ 1 = 0, (4.38
)

u 1 = 0, (4.39) r → ∞ : c + 1 → 0, c - 1 → 0, ∂Φ 1 ∂r → -cos θ, (4.40
)

U 1 → - U ∞ E ∞ sin θ, V 1 → U ∞ E ∞ cos θ. (4.41)
Since the solution of the system depends on E ∞ , the charge distribution inside the EDL will be inhomogeneous with respect to the angle θ (Fig. 4.1). Consequently, an electroosmotic flow is formed near the particle u 1 = 0. Note that the system (4.35) -(4.41) is a linear system of equations taking into account the previously obtained functions of zeroth-order c + 0 , c - 0 , Φ 0 .

Inner expansion (ν → 0)

The system (4.35) -(4.41) with respect to the variable z = r-1 ν in the inner expansion can also be simplified by cancelling terms of a higher-order in ν. As a result, the system takes the form,

∂ ∂z c + 0 ∂Φ 1 ∂z + c + 1 ∂Φ 0 ∂z + ∂c + 1 ∂z = 0, (4.42)
∂ ∂z c - 0 ∂Φ 1 ∂z + c - 1 ∂Φ 0 ∂z - ∂c - 1 ∂z = 0. (4.43)
After integrating and using the condition (4.38) and defining

j 1 = c + 0 ∂Φ 1 ∂r + c + 1 ∂Φ 0 ∂r + ∂c + 1 ∂r , c + 0 ∂Φ 1 ∂z + c + 1 ∂Φ 0 ∂z + ∂c + 1 ∂z = νj 1 → 0 (4.44) c - 0 ∂Φ 1 ∂z + c - 1 ∂Φ 0 ∂z - ∂c - 1 ∂z = 0 (4.45) ∂ 2 Φ 1 ∂z 2 = c - 1 -c + 1 (4.46) z = 0 : c + 1 = 0, c - 1 ∂Φ 0 ∂z + c - 0 ∂Φ 1 ∂z - ∂c - 1 ∂z = 0 (4.47) Φ 1 = 0. (4.48)
From the outer electrically neutral solution, we can write,

z → ∞ : c + 1 -c - 1 → 0. (4.49)
Note that unknown variables in the system depend on the angle θ, but this dependence is parametric.

For the uniqueness of the solution of the system (4. By introducing new variables,

E 0 = ∂Φ 0 ∂z , E 1 = ∂Φ 1 ∂z , ρ 1 = c + 1 -c - 1 , K 1 = c + 1 + c - 1 . (4.50)
Adding and subtracting the equations (4.42) and (4.43), the following system can be obtained,

K 0 E 1 + K 1 E 0 + ∂ρ 1 ∂z = 0, (4.51)
ρ 0 E 1 + ρ 1 E 0 + ∂K 1 ∂z = 0, (4.52
)

∂E 1 ∂z = -ρ 1 . (4.53)
The equation (4.52) can be represented in the form (using also Eq. (4.14)):

∂K 1 ∂z = ∂E 0 ∂z E 1 + ∂E 1 ∂z E 0 = ∂E 1 E 0 ∂z , (4.54)
after a single integration, one can get,

K 1 = E 1 E 0 + A, (4.55)
where A is the integration constant. For z → ∞ E 0 → 0 and therefore,

A = K 1 z→∞ ≡ K 1 ∞ . ( 4
, ∂ 2 E 1 ∂z 2 -(E 0 2 + K 0 )E 1 = K 1 ∞ E 0 . (4.57)
The boundary conditions (4.47) -(4.49) on the particle can be transformed into the form,

z = 0 : ∂E 1 ∂z -E 0 E 1 = K ∞ 1 , z → ∞ : ∂E 1 ∂z → 0 (4.58)
After introducing the new unknown

F (z) = E 1 /K ∞ 1 , the boundary-value problem (4.57) -(4.58) turns into, d 2 F dz 2 -(E 2 0 + K 0 )F = E 0 , (4.59)
z = 0 : dF dz -E 0 F = 1, z → ∞ : dF dz → 0. (4.60)
Integrating E 1 = K ∞ 1 F from zero to infinity and using the boundary conditions (4.58), one can obtain,

∞ 0 ∂Φ 1 ∂z dz = Φ 1 z→∞ = K ∞ 1 ∞ 0 F dz, Φ 1 z=0 = 0
The solution of the system (4.59) -(4.60) leads to,

∞ 0 F dz = - 1 2 , (4.61)
and the matching conditions impose the following equality on the boundary between the inner and outer decompositions,

z → ∞ : K 1 = -2Φ 1 . (4.62)
The solution for the inner expansion is not given explicitly since it is obtained numerically by integrating the equations (4.59) -(4.60).

Outer expansion

The system (4.35) -(4.41) with respect to K 1 = c + 1 + c - 1 together with the matching condition for the decompositions (4.62) can be written as,

u 1 • ∇K 0 = ∇ 2 K 1 , ∇ • (K 0 ∇Φ 1 + K 1 ∇Φ 0 ) = 0, (4.63) r = 1 : K 1 ∂Φ 0 ∂r + K 0 ∂Φ 1 ∂r - ∂K 1 ∂r = 0, (4.64) 2Φ 1 = -K 1 , (4.65) r → ∞ : K 1 → 0, ∂Φ 1 ∂r → -cos θ. (4.66)
Taking into account that from zeroth-order solution K 0 ≡ 2 and Φ 0 ≡ ζ,

∇ 2 K 1 = 0, ∇ 2 Φ 1 = 0, (4.67) r = 1 : 2 ∂Φ 1 ∂r = ∂K 1 ∂r , 2Φ 1 = -K 1 , (4.68) r → ∞ : K 1 → 0, ∂Φ 1 ∂r → -cos θ. (4.69)
The solution to the system (4.67) -(4.69) can be easily found,

Φ 1 = -r + 1 4r 2 cos θ, K 1 = 3 2r 2 cos θ. (4.70)
and the following relations are true,

r = 1 : Φ 1 = - 3 4 cos θ, ∂Φ 1 ∂θ = 3 4 sin θ, K 1 = 3 2 cos θ. (4.71)
The electric current j(θ) on the particle's surface at r = 1,

j = j 1 E ∞ = E ∞ 2 ∂K 1 ∂r + K 0 ∂Φ 1 ∂r + K 1 ∂Φ 0 ∂r r=1 + O(E 2 ∞ ) = = -3 E ∞ cos θ + O(E 2 ∞ ). (4.72)
In a first approximation, the current is symmetrical regarding the y axis; j disappears at the poles of the particle at θ = ± π 2 . Thus, we can write out the first order approximation for the functions K and Φ according to the expansions (4.1):

K = K 0 + K 1 E ∞ = 2 + 3 2r 2 E ∞ cos θ (4.73) Φ = Φ 0 + Φ 1 E ∞ = ζ + -r + 1 4r 2 E ∞ cos θ (4.74)
We do not write out the approximation for ρ explicitly, but note that the expansion inside the EDL (zeroth-order) can be obtained from the equations (4.18) -(4. [START_REF] Helmholtz | Studien über electrische grenzschichten[END_REF], and in the outer region ρ = 0. These approximations (4.73) -(4.74) will be used in the future for comparison with the results of numerical simulations.

Electroosmotic flow and slip velocity

From the solution in the previous paragraph, it was shown that the electrostatic problem can be solved independently of the hydrodynamic part.

The Stokes equation (4.37) with the boundary conditions (4.40) and (4.41) can now be solved. For inner decomposition, the system (4.37) has the form,

-ν 2 ∂P 1 ∂θ + ∂ 2 U 1 ∂z 2 = -κ ∂ 2 Φ 0 ∂z 2 ∂Φ 1 ∂θ , ∂P 1 ∂z = κ ν 2 ∂ ∂z ∂Φ 1 ∂z ∂Φ 0 ∂z , V 1 = 0, (4.75)
with boundary conditions,

z = 0 : U 1 = 0, z → ∞ : ∂U 1 ∂z = 0. (4.76)
After eliminating the pressure P 1 from the system, one can obtain the following equation,

∂ 2 U 1 ∂ 2 z = -κ ∂ 2 Φ 0 ∂z 2 ∂Φ 1 ∂θ + κ ∂ 2 Φ 1 ∂θ∂z ∂Φ 0 ∂z , (4.77)
After integrating this expression and using the second boundary condition (4.76), one can obtain,

∂U 1 ∂z = -κ ∂Φ 0 ∂z ∂Φ 1 ∂θ -2κ ∞ z ∂ 2 Φ 1 ∂θ∂s ∂Φ 0 ∂s ds.
The second integration with respect to z will lead to,

∞ 0 ∂U 1 ∂z dz = -κ   Φ 0 ∂Φ 1 ∂θ ∞ 0 - ∞ 0 Φ 0 ∂ 2 Φ 1 ∂θ∂z dz   -2κ ∂ ∂θ ∞ 0 ∞ z E 1 E 0 ds dz.
Taking into account that U 1 z=0 = 0, Φ 0 z=0 = 0 and Φ 0 z→∞ → ζ, one can find the following expression for the slipping velocity U s ,

U s = U 1 z→∞ = -κ ζ ∂Φ 1 ∂θ z→∞ + κ ∞ 0 Φ 0 ∂ 2 Φ 1 ∂θ∂s ds -2κ ∂ ∂θ ∞ 0 ∞ z E 1 E 0 dsdz = = -κ ζ ∂Φ 1 ∂θ z→∞ + κ ∂ ∂θ   ∞ 0 Φ 0 E 1 ds -2 ∞ 0 ∞ z E 1 E 0 dsdz   = = -2κ   ζ 2 + ∞ 0 Φ 0 F ds -2 ∞ 0 ∞ z F E 0 dsdz   ∂Φ 1 ∂θ z→∞ .
From the equation (4.71),

∂Φ 1 ∂θ z→∞ = ∂Φ 1 ∂θ r=1 = 3 4 sin θ
denoting the integrals,

S 1 = ∞ 0 Φ 0 F ds, S 2 = -2 ∞ 0 ∞ z F E 0 dsdz,
the slip velocity is,

U s = -κ ζ δ(p) ∂Φ ∂θ r=1 (4.78)
In our case, the equation turns into,

U s = - 3 4 ζκ sin θ δ(p). (4.79)
Here The outer solution for the velocity field is described by the system,

δ(p) = 1 + 2(S 1 + S 2 )/ζ
-∇P 1 + ∇ 2 u 1 = 0, ∇ • u 1 = 0 (4.81)
with boundary conditions

r = 1 : V 1 = 0, U 1 = U s = - 3 2 κ δ sin θ, (4.82) r → ∞ : U 1 = - U ∞ E ∞ sin θ, V 1 = U ∞ E ∞ cos θ. (4.83)
After introducing the stream function Ψ in spherical coordinates,

U = - 1 sin θ 1 r ∂Ψ ∂r , V = 1 sin θ 1 r 2 ∂Ψ ∂θ (4.84)
and corresponding transformations, the system (4.81) -(4.83) turns into the following,

∇ 4 Ψ 1 = 0, (4.85) r = 1 : Ψ 1 = 0, ∂Ψ 1 ∂r = 3 4 κ ζ δ sin 2 θ, (4.86) r → ∞ : ∂Ψ 1 ∂θ = r 2 U ∞ E ∞ sin θ cos θ, ∂Ψ 1 ∂r = r U ∞ E ∞ sin 2 θ. (4.87)
The general solution of the equation (4.85) can be written as,

Ψ 1 = ∞ n=1 A n r n+3 + B n r n+1 + C n r 2-n + D n r -n Q n (cos θ), (4.88)
where Q n are the Gegenbauer polynomials (see [START_REF] Happel | Low Reynolds number hydrodynamics[END_REF][START_REF] Leal | Laminar Flow and Convective Transport Processes[END_REF]), the coefficients A n , B n , C n and D n can be found from the boundary conditions (4.86) -(4.87). The solution of the system (4.85) -( 4.87) has finally the following form,

Ψ 1 = -r 2 + 3 2 r - 1 2r U ∞ E ∞ -κ 3 4 ζδ r - 1 r Q 1 (cos θ), (4.89)
where Q 1 is the first Gegenbauer polynomial,

Q 1 (cos θ) = 1 2
(cos 2 θ -1).

Electrophoretic velocity

To find the electrophoretic velocity U ∞ , the general formulation (2.26) -(2.35) must be supplemented by imposing the condition of zero force acting on the particle [START_REF] Happel | Low Reynolds number hydrodynamics[END_REF][START_REF] Leal | Laminar Flow and Convective Transport Processes[END_REF]. The balance of viscous stress forces F V and Maxwell stresses F M on the particle surface Σ has the form,

F V + F M = Σ (-P n + τ nn n + τ nl l + τ M nn n + τ M nl l)dσ = 0 (4.90)
The projection of the normal unit vector n on the x axis is n x = cos θ, and the projection of the tangent unit vector l is l x = -sin θ, see fig. 4.2. Hence,

τ nn n = τ rr n x + τ rθ l x ; τ M nn n = τ M rr n x + τ M rθ l x (4.91) n x = cos θ, l x = -sin θ, dσ = 2πr 2 sin θdθ, dl = rdθ,
τ rr = 2 ∂V ∂r , τ rθ = 1 r ∂V ∂θ + ∂U ∂r - U r , (4.92)
τ M rr = κ 2 ∂Φ ∂r 2 - 1 r 2 ∂Φ ∂θ 2 , τ M rθ = κ 1 r ∂Φ ∂r ∂Φ ∂θ . ( 4
F V = π π 0 -sin θ ∂ 3 Ψ ∂r 3 -sin θ ∂ 3 Ψ ∂r∂θ 2 + 2 sin θ ∂ 2 Ψ ∂r 2 + 5 cos θ ∂ 2 Ψ ∂r∂θ -4 sin θ ∂Ψ ∂r -8 cos θ ∂Ψ ∂θ dθ, (4.94) F M = -πκ π 0 cos θ sin θ ∂Φ ∂r 2 + sin 2 θ ∂Φ ∂θ ∂ 2 Φ ∂r 2 + sin 2 θ ∂Φ ∂θ ∂ 2 Φ ∂θ 2 dθ. (4.95)
The balance of viscous and electric forces should be performed on any closed curve around the particle. But it is more convenient to take it at the edge of the EDL, which corresponds to the outer solution. Then,

F V + F M = π 3κ ζδ -6 U ∞ E ∞ = 0. (4.96)
Electrophoretic velocity is found from this expression, namely,

U ∞ = κ ζδ 2 E ∞ , (4.97)
Taking into account the expression (4.80), one can get,

U ∞ = 1 2 κζ(1 -0.11ζ)E ∞ (4.98)
that in dimensional form has the form,

U ∞ = εζ 2µ 1 -0.11 ζF RT E ∞ . (4.99)
One can emphasize that this formula differs from the well-known Helmholtz-Smoluchowski formula for dielectric particles,

U ∞ = κζE ∞ (4.100)
which in dimensional form has the expression,

U ∞ = εζE ∞ µ . (4.101)

Numerical results

The numerical solution of the time-dependent system (2.22) -(2.25) is presented in this section. It is assumed that the electrolyte is an aqueous solution of sodium chloride (NaCl). For this solution the parameter κ is fixed and κ = 0.26. The solution to the problem depends on the Debye number, but the main results are presented for ν = 0.0087 and ν = 0.002, which corresponds to a particle size radius of 5 µm and 21 µm correspondingly. The problem substantially depends on the external electric field E ∞ and the results are divided into three cases: for a weak, moderate and strong electric field [START_REF] Frants | Nonequilibrium electrophoresis of an ion-selective microgranule for weak and moderate external electric fields[END_REF][START_REF] Ganchenko | The movement of an ionexchange microparticle in a weak external electric field[END_REF][START_REF] Ganchenko | Extreme nonequilibrium electrophoresis of an ion-selective microgranule[END_REF][START_REF] Ganchenko | Instabilities, bifurcations, and transition to chaos in electrophoresis of charge-selective microparticle[END_REF]. In this section, the results for the weak and moderate electric field are presented. The numerical results for the weak electrical field is also compared with the analytical solution and experiments.

A numerical simulation of the complete system of equations was carried out. The results presented for κ = 0.26 and ν = 0.0086, which corresponds to a solution of N aCl with an equilibrium concentration of c∞ = 0.1 mol/m 3 and the particle radius of ã = 5 µm. where the cross-sections for the angles θ = 0 • and θ = 180 • are given, differs. This difference in the behaviour of the functions can be explained as follows. In the first-order

approximation, ρ = ρ 0 + ρ 1 E ∞ + O(E 2 ∞ )
and ρ 1 has the same order of magnitude as E ∞ , which is small (E ∞ = 0.05). Therefore, ρ ρ 0 , which can be seen in Fig. 4

.3 (b).

At the same time for electrical conductivity, the first-order approximation has the form

K = K 0 + K 1 E ∞ + O(E 2
∞ ) and is no longer a small value and therefore K = K 0 , as can be seen from Fig. 4.4.

With an increase in the electric field strength E ∞ , the nonequilibrium effects become more significant and the results of numerical analysis begin to deviate from the analytical one derived for small values of E ∞ . In Figs. In the region of the incoming ion flow, at 90 • < θ < 180 • , the charge density ρ, and the EDL thickness are increasing, while in the region of the outgoing ion flow, at 0 • < θ < 90 • , both the charge density and the EDL thickness are reduced. In this case, the charge in the region 90 • < θ < 180 • remains positive, while in the region 0 • < θ < 90 • it changes sign and becomes negative. This fact is associated with the influence of the ions flow into and out of the surface, that is, it is a typical non-equilibrium effect.

With the increase in the electric field strength, the symmetric electrically neutral clouds of the enriched and depleted electrolyte lose their symmetry (see Fig. 4.6(a)). The zone of depleted electrolyte also called the desalination zone, due to salt concentration here is much lower than in the bulk or it even tends to zero. In addition to the loss of symmetry, these clouds expand, and the zone of the depleted solution expands faster. In general, however, the law of conservation of mass is still valid, (equilibrium value) to a very small value, K ≈ 0 (Fig. 4.6(b)). With increasing E ∞ , the region with K ≈ 0 expands, as in the case of flat membranes [START_REF] Rubinstein | Electro-osmotically induced convection at a permselective membrane[END_REF][START_REF] Zaltzman | Electro-osmotic slip and electroconvective instability[END_REF][START_REF] Demekhin | Direct numerical simulation of electrokinetic instability and transition to chaotic motion[END_REF], and the electrical resistance of this zone increases. In the outer region of the electrolyte, K slowly increases with distance from the surface of the particle and reaches the equilibrium value of K = 2. Inside the EDL of the enriched electrolyte region, the total ion concentration slowly decreases and also reaches the equilibrium value of K = 2 in the outer region of the electrolyte.

Recall that for small electric field strength E ∞ , the charge inside the EDL is positive and does not depend on the angle θ (see Fig. In addition, for small values of the angle θ, the electroosmotic velocity changes direction for relatively small values of E ∞ = 0.05 and 3.

For small values of E ∞ , the electrokinetic effects are in a quasi-equilibrium state and are completely controlled by the concentration of positive ions on the particle surface (parameter p).

The problem has a parameter p representing the concentration of ions inside the pores of the membrane near the surface. It is only known that p > 1 [START_REF] Rubinstein | Voltage against current curves of cation exchange membranes[END_REF]. The value of p is difficult to measure experimentally. Therefore, the fact of independence (or weak dependence) on p for sufficiently large values of E ∞ is very important. controlled by the flow of ions (and not by concentration of positive ions p on the surface, as it was for small values of E ∞ ), and non-equilibrium effects become predominant. from some point [START_REF] Frants | Nonequilibrium electrophoresis of an ion-selective microgranule for weak and moderate external electric fields[END_REF]. Thus, it can be assumed that for moderate electric field strength, the particle velocity U ∞ does not depend on this parameter for sufficiently higher values of p.

Chapter 5

Electrophoresis of an ion-selective particle in strong electric field

As was noted earlier, the Smoluchowski theory predicts a linear increase in the particle velocity with an increase in the ζ potential. However, this dependence was deduced under the assumption of small values of the zeta potential and a weak electric field and becomes invalid for a sufficiently large value of the zeta potential or electric field [START_REF] Schnitzer | Dielectric-solid polarization at strong fields: Breakdown of smoluchowski's electrophoresis formula[END_REF][START_REF] Schnitzer | Macroscale description of electrokinetic flows at large zeta potentials: Nonlinear surface conduction[END_REF].

In a high-intensity electric field, the electric current passing through the particle becomes significant, and the EDL changes its structure qualitatively, acquiring a nonequilibrium state (Fig. 5.1).

Figure 5.1 -The structure of a non-equilibrium EDL near a cation-selective particle in a strong electric field. Fig. 5.1 illustrates the behavior of an ion-selective particle when a quite strong electric field is applied. First, consider the area to the left of the particle. Cations moving in a direction of the electric field (from left to right), when they reach the left side of the particle's surface, they can enter inside it and exit from the backside, thus creating a positive electric current through the particle. At the same time, the anions move in the opposite direction. Thus, a low concentration of both cations and anions is formed in this area, i.e., local desalination of the electrolyte occurs.

Cations do not accumulate inside the particle despite the force of attraction to the negatively charged functional groups of the particle's material since the strength of the applied external electric field significantly exceeds the forces of attraction between ions of different signs. Thus, the electro-migration mechanism is dominant. Cations move out from the particle from the right side and continue to move in the direction from the particle, i.e. the concentration of cations is increased here. Anions, in turn, move towards the particle's surface, and since they cannot pass through it, some of them "stick" to the particle's surface, despite the repulsive forces with the negatively charged particle. The remaining anions bypass the particle and continue to move from right to left. Thus, in the region to the right of the particle, there is an increased concentration of both anions and cations (enriched electrolyte solution).

For a flat ion-exchange membrane, the regions of the incoming and outgoing ion flux are separated by an electrolyte layer (see Fig. 1 in [START_REF] Zaltzman | Electro-osmotic slip and electroconvective instability[END_REF]). An interesting feature of the problem with the ion-selective particle is that the topology of the problem changes and both regions are located on the same boundary, and the flows of incoming and outgoing ions are separated by points 1 and 1 , in which the electric current disappears, j(θ 0 ) = 0 (Fig. 5.2). These points separate the EDL with a positive charge for the region of the incoming stream and the EDL with a negative charge for the outgoing stream of ions.

The space charge zone exists only for the region of the incoming flow and disappears when approaching the point θ = θ 0 . The electroosmotic velocity at the boundary of the space charge region creates a thin diffusion layer. At θ = θ 0 , the diffusion boundary layer is expected to detach (see Levich [START_REF] Levich | Physicochemical hydrodynamics[END_REF]).

Let us give the estimation of the thickness of thin inner layers at ν → 0 and E ∞ → ∞, without considering the relationship between the order of ν and E ∞ . The thickness of the EDL is estimated as O(ν), the estimate of the thickness of the SCR y m , taken from [START_REF] Rubinstein | Electro-osmotically induced convection at a permselective membrane[END_REF],

is y m = O(ν 2/3 E 1/ 3
∞ ). It will be shown later that the thickness of the diffusion layer is estimated as δ = O(E -1 ∞ ). Note that all three estimations of the thicknesses of the layers are much smaller than the particle radius. The thickness of the EDL, in turn, is much smaller than the thickness of the SCR y m (this is true if E ∞ ≤ O(ν -2 )) and y m /δ is of the order of O(ν 2/3 E 4/3 ∞ ), where δ is the thickness of the diffusion layer (see Fig. 5.

2).

In what follows, it will be shown that both layers (SCR and diffusion layer) have the same order of magnitude, and the last ratio is O(1). For convenience, it can be introduced a new variable y = r -1 which is reckoned from the interface between the liquid and the particle surface. In Fig. 5.2, the region under consideration is 0 < y < ∞ and 0 < θ < 180 • which in turn is divided into several subregions (layers). Thin layers of I, II and III are described by the inner solution, while the outer solution is associated with region IV. The region I is EDL and it forms on the interface between solid particle and electrolyte. The region II is a space charge region with a thickness of y m , it follows just after EDL. The next region III following up SCR is a diffusion layer with a thickness of δ, which is the last in the structure of boundary layers nested in each other. After this structure is located the electrically neutral outer region IV of the electrolyte with dimensionless equilibrium electrical conductivity K = 2.

Note that in the region of outgoing ion flow (see later Fig. 5.6 (a)) forms a jet with high electrical conductivity where the advection of salt ions becomes comparable with diffusion. A similar form of the far-field singularity was found for the case of a strongly charged colloidal particle [START_REF] Schnitzer | Weakly nonlinear electrophoresis of a highly charged colloidal particle[END_REF].

The earliest analysis of the Nernst-Planck-Poisson system was carried out by Felici [START_REF] Felici | Recent advances in the analysis of d.c. ionized electric fields -part i[END_REF][START_REF] Felici | Recent advances in the analysis of d.c. ionized electric field -part ii[END_REF]. Asymptotic analysis for boundary value problems of the Nernst-Planck-Poisson equation system with regards the small parameter ν was carried out by Urtenov [START_REF] Ch | Boundary value problems for systems of nernst-planck-poisson equations (asymptotic expansions and related issues)[END_REF].

5.1 Analytical solution in space charge region

Electrostatic part

Now the case of high electric field strength, E ∞ → ∞ is under consideration. An electric charge near the particle surface is created due to the "sticking" of cations at a nonequilibrium ion flow through the particle. The first layer is formed near the charged surface (EDL), and a space charge region is formed behind it. The net charge in SCR is much greater than the net charge of EDL, so it is possible to neglect the EDL during the asymptotic analysis.

Since the thickness of the SCR is small compared to the radius of the particle, one can take r ≈ 1. The solution can be considered as time-independent, i.e. ∂/∂t = 0, and the change in the solution with respect to the radial direction occurs much more strongly than in the angular direction, therefore, the derivatives with respect to θ can be neglected, i.e.

∂/∂r ∂/∂θ.

All these considerations could be obtained exactly by stretching the variables in a thin layer, as was done in the previous chapter. However, here, due to the complex structure of the layers, it is more convenient not to introduce a multitude of stretched variables, but to operate in the usual system (r, θ) only by compiling the equations applicable in the corresponding layers.

Thus, from the Nernst-Planck and Poisson equations one can obtain the following system, 1

r 2 ∂ ∂r r 2 ∂c ± ∂r ± r 2 c ± ∂Φ ∂r = 0 ν 2 r 2 ∂ ∂r r 2 ∂Φ ∂r = c --c +
Since it was assumed that the influence of spherical geometry is insignificant within the space charge region (where r 1), then one can obtain,

∂ ∂r ∂c ± ∂r ± c ± ∂Φ ∂r = 0 ν 2 ∂ ∂r ∂Φ ∂r = c --c + Denote ∂Φ ∂r = E, then, ∂ ∂r c + E + ∂c + ∂r = 0 ∂ ∂r c -E - ∂c - ∂r = 0 ν 2 ∂E ∂r = c --c + (5.1)
From the boundary condition (2.32),

c -E - ∂c - ∂r = 0 (5.2)
And from condition (2.33),

c + E + ∂c + ∂r = j (5.3)
Adding and subtracting equations (5.2) and ( 5.3), one can get,

(c + + c -)E + ∂ ∂r (c + -c -) = j (5.4) (c + -c -)E + ∂ ∂r (c + + c -) = j (5.5)
and make the following substitution for the equations (5.2) and (5.3),

c + -c -= -ν 2 ∂E ∂r (5.6)
Introducing a new variable y = r -1, then ∂ ∂r = ∂ ∂y and using Eq. (5.1), Eq. ( 5.4) can be written as,

(c + + c -)E -ν 2 ∂ 2 E ∂y 2 = j (5.7)
and Eq. (5.7) leads to,

-ν 2 ∂E ∂y E + d dy (c + + c -) = j or j = d dy - 1 2 ν 2 E 2 + c + + c - (5.8)
From the equation (5.8) by means of integration one can obtain,

- 1 2 ν 2 E 2 + c + + c -= j(y -y m ) (5.9)
where y m is an unknown integration constant. By substituting (5.9) into (5.7),

j(y -y m ) + 1 2 ν 2 E 2 E -ν 2 ∂ 2 E ∂y 2 = j ν 2 ∂ 2 E ∂y 2 -j(y -y m ) + 1 2 ν 2 E 2 E + j = 0 or ν 2 ∂ 2 E ∂y 2 + j(y m -y) - 1 2 ν 2 E 2 E + j = 0 (5.10)
Since ν → 0 and d dy = O(1), the first term in (5.10) has a higher order of smallness.

At the same time, the term 1 2 ν 2 E 3 cannot be neglected, because E > O(ln 1 ν ) [START_REF] Zaltzman | Electro-osmotic slip and electroconvective instability[END_REF], which means that E → ∞. It follows from (5.8) that j = O(1) and j can be neglected with respect to the second term in the equation (5.10). As a result, one can obtain the cubic algebraic equation in E,

j(y m -y) - 1 2 ν 2 E 2 E = 0, (5.11)
which was previously obtained and analyzed for the case of flat membranes by Urtenov and Babeshko [START_REF] Ch | Boundary value problems for systems of nernst-planck-poisson equations (asymptotic expansions and related issues)[END_REF].

The first trivial solution E 1 = 0 makes sense only if y > y m . When solving the equation,

- 1 2 ν 2 E 2 + j(y m -y) = 0
we get two more roots: E 2 and E 3

E 2 = - 1 ν 2j(y m -y), E 3 = 1 ν 2j(y m -y) E 2 < 0 has no physical meaning.    E = 1 ν 2j(y m -y), 0 < y < y m E = 0, y ≥ y m (5.
12)

It follows from the system above that the integration constant y m is the boundary of the SCR. From the solution in SCR, i.e. at 0 < y < y m ,

ν dΦ dy = √ 2
= ν(Φ(y) -Φ(0) =0 ) = - 2 √ 2 3j (j(y m -y)) 3/2 y 0 = - 2 √ 2 3j (j(y m -y)) 3/2 + 2 √ 2 3j (jy m ) 3/2 ∆Φ = 2 √ 2 3jν (jy m ) 3/2 - 2 √ 2 3jν
(j(y m -y)) 3/2 (5.13)

At the SCR boundary (y = y m ), the equation above turns into,

∆Φ| y=ym = 2 √ 2 3ν j 1/2 y 3/2 m = ∆Φ m (5.14)
and we can express the thickness of SCR,

y 3 m = 9 8 ν 2 ∆Φ 2 m j (5.15)
Here ∆Φ m = Φ(y m ) -Φ(0) is the potential drop in the SCR.

It is possible to find ρ and K (or c + and c -). From (5.1) ρ = -ν 2 ∂E ∂y , and from (5.9) K = j(y -y m ) + 1 2 ν 2 E 2 . Thus, from Eq. (5.12),

ρ = -ν 2 ∂ ∂y 1 ν 2j(y m -y) = -ν √ 2 ∂ ∂y (jy m -jy) 1/2 = = -ν √ 2 1 2 (-j)(jy m -jy) -1/2 = ν √ j √ 2 (y m -y) -1/2 K = j(y -y m ) + 1 2 ν 2 1 ν 2 (2j(y m -y)) = 0 So finally,      c + -c -= ν √ j √ 2 (y m -y) -1/2 c + + c -= 0
From the system above we can express the concentrations,

c + = ν √ j 2 √ 2 (y m -y) -1/2 and c -= - ν √ j 2 √ 2 (y m -y) -1/2
Finally, we present the solution in the space charge region,

∆Φ = 2 √ 2 3jν (jy m ) 3/2 - 2 √ 2 3jν (j(y m -y)) 3/2 ρ = ν √ j √ 2 (y m -y) -1/2 and K ≈ 0

Hydrodynamic part

Let us consider the Stokes equation in spherical coordinates in accordance with the following assumptions: ∂ ∂θ ∂ ∂r and V U (see Fig. 2.3). Then,

- 1 r ∂P ∂θ + ∂ 2 U ∂r 2 = ρ 1 r κ ν 2 ∂Φ ∂θ - ∂P ∂r = ρ 1 r κ ν 2

∂Φ ∂r

Using the following relations:

X = π -θ, ∂ ∂X = -∂ ∂θ , y = r -1,, ∂ ∂r = ∂ ∂y , ρ = -ν
- ∂ ∂X κ 2 E 2 + c = ∂ 2 U ∂y 2 -κ ∂E ∂y ∂Φ ∂X - κ 2
∂E 2 ∂X = ∂ 2 U ∂y 2 -κ ∂E ∂y ∂Φ ∂X 1 κ ∂ 2 U ∂y 2 = ∂E ∂y ∂Φ ∂X - 1 2 ∂E 2 ∂X
To find the velocity, the above equations is integrated twice,

1 κ ym(X) y ∂ 2 U ∂y 2 dy = ym(X) y ∂E ∂y ∂Φ ∂X dy - 1 2 ym(X) y ∂E 2 ∂X dy (5.17)
After substituting expressions for Φ and E,

Φ = 1 ν 2 √ 2 3j (jy m ) 3/2 -(j(y m -y)) 3/2 , E = √ 2 ν j(y m -y)
When calculating the solution for electrostatics, j and y m are constants because ∂ ∂r >> ∂ ∂θ . These quantities are slowly varying functions of the angle θ. In the analysis of electrostatics, this did not affect the solution, while for hydrodynamics this becomes important since otherwise U ≡ 0.

The equation (5.17) will turn into,

1 κ ym(X) y ∂ 2 U ∂y 2 dy = 4 3ν 2 ym(X) y ∂ ∂y ( j(y m -y)) ∂ ∂X (jy m ) 3/2 -(j(y m -y)) 3/2 dy- - 1 ν 2 ym(X) y ∂ ∂X (j(y m -y))dy (5.18)
Let us consider each integral separately,

1 κ ym(X) y ∂ 2 U ∂y 2 dy = 1 κ     ∂U ∂y (y m ) =0 - ∂U ∂y (y)     = - 1 κ ∂U ∂y - 1 κ ym(X) 0 ∂U ∂y dy = - 1 κ (U (y m ) -U (0) =0 ) = - 1 κ U m
After calculating the right-hand side of Eq. ( 5.18), we obtain a formula similar to that obtained for a flat membrane 0 < y < y m (X),

U m = κ 8 ∆Φ 2 1 j ∂j ∂X + κ 2 ∂ ∂X ∆Φ 2 (5.19)
In contrast to the case of a flat membrane, for the spherical particle, it is also necessary to take into account the component V of the velocity, which can be obtained from the continuity equation, ∂ ∂y

(V sin X) - ∂ ∂X (U sin X) = 0 ym(X) 0 ∂ ∂y (V sin X)dy - ym(X) 0 ∂ ∂X (U sin X)dy = 0 V (y m (X)) sin X = ∂ ∂X   ym(X) 0 U (y) sin Xdy   - -U (y m (x)) sin X ∂ ∂X y m (X)
We use the boundary conditions U (0) = V (0) = 0 and V (y m ) = V m , U (y m ) = U m , then,

V m sin X = ∂ ∂X (U (y m (X)) sin X y m ) -U m sin X ∂y m ∂X or V m + U m ∂y m ∂X = 1 sin X ∂ ∂X (U m sin X y m ) (5.20)
substitute U m from Eq. (5.19), ∆Φ and y m from Eq. (5.14) - (5.15),

U m = κ 8 ∆Φ 2 1 j ∂j ∂X + κ 2 ∂ ∂X ∆Φ 2 , y 3 m = 9 8 ν 2 ∆Φ 2 j ∆Φ 2 = 8 9 1 ν 2 y 3 m j U m = κ 8 8 9 1 ν 2 y 3 m j 1 j ∂j ∂X + κ 2 ∂ ∂X 8 9 1 ν 2 y 3 m j = κ ν 2

Self-similar solution in diffusion layer

A thin electrically neutral diffusion layer is located behind the space charge zone for y > y m and has a characteristic thickness δ(X). In this layer, occurs the convection of the concentration of salt K along the particle surface with an electroosmotic velocity and at the same time K diffuses in this layer from K ≈ 0 in the depleted region to the equilibrium value K = 2 in the outer region of the electrolyte. Adding the equations (2.26) - (2.27) and assuming that the thickness of the diffusion layer is small compared to the radius of the particle, one can obtain the convection-diffusion equation in the approximation of the boundary layer, which implies the relation for V ,

∂K ∂t + U ∂K ∂X + V ∂K ∂y = ∂ 2 K ∂y 2 . ( 5
V sin X = V m sin X + η ∂ ∂X U m sin X , (5.28)
where η ≡ y -y m , and the velocity component V is formed by the electroosmotic part of V m (in fact, this is the intake of liquid from the space charge region) and by the change in the longitudinal component U m along the particle surface.

The equation (5.23) after substituting (5.28) turns into a linear partial differential

equation for K, ∂K ∂t sin X -U m sin X ∂K ∂X + η ∂ ∂X U m sin X ∂K ∂η + V m sin X ∂K ∂η = ∂ 2 K ∂η 2 sin X, (5.29)
where U m and V m are determined by the relations (5.21) and (5.22), with the boundary conditions,

η = 0 : K = 0, ∂K ∂η = j; η/δ → ∞ : K = 2.
(5.30)

The necessity for a solution (5.29) -( 5

j δ = 4 √ π exp (-b 2 /δ 2 ) 1 -erf (-b/δ) (5.34)
Let us stretch the unknowns so that they have O(1) in a new form,

u m = E 2 ∞ U m , v m = E 2 ∞ V m , ψ = E 2 ∞ Ψ, u ∞ = E 2 ∞ U ∞ , ∆φ = E ∞ F, φ = E ∞ Φ, j = E ∞ J, ∂/∂t = E 2 ∞ ∂/∂T, (5.35) b = ν 2/3 E 1/3 ∞ B, y m = ν 2/3 E 1/3 ∞ Y m , δ = E -1 ∞ ∆.
The system (5.32) -(5.34) after substituting the expressions (5.21) and (5.22) for u m and v m turns to the following system,

∂∆ ∂T sin X -κ ∂ ∂X F 2 ∆ 8 J ∂J ∂X + ∆ 2 ∂F 2 ∂X = 2 ∆ sin X, (5.36) ∂B ∂T sin X -κ ∂ ∂X F 2 8J ∂J ∂X + 1 2 ∂F 2 ∂X B sin X = = κ 3 8/3 16 ∂ ∂X 7 20
F 8/3 J 4/3 ∂J ∂X + J 5 ∂ ∂X F 8/3 J 4/3 sin X (5.37) - 3 2/3 2 sin X ∂ ∂T F 2/3 J 1/3 , J ∆ = 4 √ π exp (-χ 2 B 2 /∆ 2 ) 1 -erf (-χB/∆) , (5.38)
with respect to four unknown functions J, F , ∆ and B. The boundary conditions for the unknown are the conditions of symmetry at the pole,

∂∆ ∂X = 0, ∂F ∂X = 0, ∂J ∂X = 0, ∂B ∂X = 0 at X = 0. (5.39)
At X → X 0 , which corresponds to θ → θ 0 , the electric current J → 0 by definition. After a brief analysis of the equation (5.38), it can be shown that ∆ → ∞ for X → X 0 , because the right side of the equation (5.38) is always positive. Thus, the left side of the equation cannot tend to zero for X → X 0 , therefore ∆ should tend to infinity. In turn, this means the separation of the diffusion boundary layer. The separation of the boundary layer of this kind was originally predicted by Levich [START_REF] Levich | Physicochemical hydrodynamics[END_REF] and is shown in Fig. 4.5 (a) and 5.6 (a)

obtained by direct numerical simulation. Thus, theoretical justification now exists for the problem of separation of the diffusion boundary layer.

Two small parameters of the problem -ν and A -1 ∞ -can be combined into one pa-

rameter χ = O(1), χ = ν 2/3 E 4/3 ∞ = O(y m ) O(δ) .
Note that χ is the ratio of the characteristic length of the space charge region,

ν 2/3 E 1/3 ∞
and the characteristic length of the diffusion layer, E -1 ∞ (see stretching parameters in the equations (5.35)). This ratio χ is estimated in table 1 with λD = 10 nm, ã = 25 µm -500 nm and Table 5.1 -Typical parameter values.

E ∞ = 10-10 3 . ã E ∞ ν χ 25

Solution in the bulk of electrolyte

The derived system (5.36) -(5.38) of the three equations has four unknown variables: J, F , ∆ and B. The fourth equation will be deduced from the solution in the outer region 1 < r < ∞. In this region, K = 2 and ρ = 0, so the Poisson equation,

ν 2 ∇ 2 Φ = c --c +
turns into the Laplace equation (since c + = c -) with two boundary conditions,

∇ 2 ϕ = 0 (5.40)
with two boundary conditions. Here ϕ = Φ/E ∞ . The first boundary condition for r = 1 should be taken from the solution for the inner region (inside the space charge region).

The electric potential ϕ in the inner region varies from ϕ = 0 on the particle surface to ϕ = F at the boundary of the diffusion layer. This implies the boundary condition,

r = 1 : ϕ = F (5.41)
The far field boundary condition, where P k (cos θ) are Legendre's polynomials.

r → ∞ : ∂Φ ∂r → -E ∞ cos θ turns into, r → ∞ : ∂(ϕE ∞ ) ∂r → -E ∞ cos θ ⇒ ∂ϕ ∂r → -cos θ
The boundary condition follows from,

c -∂ϕ ∂r - ∂c - ∂r = 0 c + ∂ϕ ∂r + ∂c + ∂r = j
Add the last two equations,

K ∂ϕ ∂r + ∂ρ ∂r = j (5.44)
The relation for ion flux (5.44) is need to be used to get the dependency between J and F . In the electroneutral region of the electrolyte K = 2 and ρ = 0, then,

r = 1 : J = 2 ∂ϕ ∂r (5.45)
The dependency between J and F can be obtained if Eq. ( 5.43) substitute into Eq. (5.45),

J = -2 cos θ - π 0 (F (s) + cos s) ∞ k=0 (2k + 1)(k + 1)P k (cos s)P k (cos θ) sin s ds. (5.46) r = 1 : J = 2 ∂ϕ ∂r .
(5.47)

The system of equations (5.36) -(5.38) with boundary conditions (5.39) is closed with respect to the unknown functions J, F , ∆ and B.

Electrophoretic velocity

Consider the hydrodynamic part at ρ = 0 and K = 2. Now one can solve the flow problem, consisting of homogeneous Stokes equations taken for the equation of the stream function Ψ (see [START_REF] Leal | Laminar Flow and Convective Transport Processes[END_REF]),

D 2 (D 2 Ψ) = 0, D 2 = ∂ 2 ∂r 2 + sin θ r 2 ∂ ∂θ 1 sin θ ∂ ∂θ , (5.48)
together with the slip condition for r = 1, which is derived from, r = 1 :

U = V = 0 U at y m is equal to, U m = - 1 r sin θ ∂Ψ ∂r ⇒ ∂Ψ ∂r = -U m sin θ, r ≈ 1 in SCR r = 1 : Ψ = 0 r = 1 : Ψ = 0, ∂Ψ ∂r = -U m sin θ, (5.49)
and the condition for the velocity at infinity, which is obtained from,

r → ∞ : u → -U ∞ sin θ; v → U ∞ cos θ. u = - 1 r sin θ ∂Ψ ∂r = -U ∞ sin θ ⇒ ∂Ψ ∂r = U ∞ r sin 2 θ v = 1 r 2 sin θ ∂Ψ ∂θ = U ∞ cos θ ⇒ ∂Ψ ∂θ = r 2 U ∞ cos θ sin θ r → ∞ : ∂Ψ ∂θ = r 2 U ∞ sin θ cos θ, ∂Ψ ∂r = rU ∞ sin 2 θ. (5.50)
Since the problem (5.48) -(5.50) is linear, it is convenient to present the solution as a superposition of two solutions, Ψ = Ψ ∞ + Ψ 1 , where Ψ ∞ is a solution with zero slip velocity, U m = 0, and Ψ 1 is a solution with zero velocity for r → ∞, U ∞ = 0. The first solution has a simple form,

Ψ ∞ = -U ∞ r 2 - 3 2 r + 1 2r sin 2 θ 2 . ( 5.51)
The second solution can be represented as,

Ψ 1 = -(r 2 -1) π 0 U m (s) ∞ k=1 k(k + 1)(k + 2) 4r k Q k (cos s)Q k (cos θ)ds, (5.52)
where Q k (cos θ) are the Gegenbauer's polynomials. The distribution of U m is described by the equation (5.21).

The value of U ∞ is not determined a priori. It is determined from the condition of the absence of forces per particle (or balance of forces), taking into account the internal decomposition. Despite its simplicity, this flow problem is technically complex, and one can deform the integration surface to the outer side of the space charge region, where ρ = 0. For a particle that is not affected by volume forces, it is convenient to apply Leel's theorem instead of this restriction (see [START_REF] Leal | Laminar Flow and Convective Transport Processes[END_REF], formula 4-180): to obtain the condition for the absence of forces acting on the particle, the stream function Ψ = Ψ ∞ + Ψ 1 must be orthogonal to the first Gegenbauer polynomial Q 1 = cos θ. In our case, the consequent of this theorem finally lead to the following dependence,

U ∞ = - 1 2 π 0 U m (θ) sin 2 θdθ, (5.53)
which for fixed κ is a function of only the variable χ.

Numerical results for moderate and strong electric fields

The electrophoresis velocity with increasing E ∞ begins to deviate very much from the linear behavior shown in the previous chapter.

Let us first present the results for the salt concentration K (or electrical conductivity). As shown by the theory of Rubinstein-Zaltzman [START_REF] Rubinstein | Electro-osmotically induced convection at a permselective membrane[END_REF] for flat membranes, when E ∞ increases, the effect of EDL becomes insignificant, and the role of electric current (ion flux through the particle surface) is the main one in the formation of a special electric layerthe space charge region in Rubinstein's terminology [START_REF] Rubinstein | Voltage against current curves of cation exchange membranes[END_REF].

On the left side of the microparticle in the region of the incoming ion flux, an electrically neutral desalination zone with K ≈ 0 is formed. In this region, K varies from K ≈ 0 near the particle surface to K = 2 at the outer boundary of the diffusion layer (denoted by the dashed curve in Fig. 5.4(a)). The diffusion layer increases downstream and at θ ≈ θ 0 its separation occurs, as predicted by the analytical model at E ∞ → ∞ (such kind of separation was first predicted by Levich [START_REF] Levich | Physicochemical hydrodynamics[END_REF]).

On the right side of the particle in the region of the outgoing ion flux, near θ = 0 • a region of high salt concentration K is formed. Note that if the distance y = r -1 from the particle surface increases, then the salt concentration tends to the equilibrium value K = 2 in the bulk. Moreover, on the right side of the particle, this tendency to the equilibrium value is much sluggish than on the left side (see Fig. 5.4(b)). The decrease in salt concentration in the depleted region is compensated by the increase of salt concentration in the enriched electrolyte region, so the expression (4.102) is still valid. The salt concentration in the electrically neutral region to the right of the particle increases rapidly with increasing E ∞ . Moreover, the shape of this region is also strongly deformed from the cloud-like structure at E ∞ ≤ 5. As E ∞ increases to a gradually elongated structure resembling to a jet (E ∞ ∼ 20), in which the salt is pushed out of the particle at high velocity.

The next important function is the charge density ρ(r, θ) = c + -c -. It is the density of the charge under the action of the electric field that brings the liquid into hydrodynamic motion. At moderate values of the electric field strength, E ∞ = 5 the positive charge region is formed to the left of the particle in the region of the incoming ion flow θ 0 < 180 • .

The angle θ 0 corresponds to the limit point when the SCR disappears with an electric current j = 0 (see Eq. The formation of the SCR for limiting currents is a well-known phenomenon, which was theoretically predicted by Rubinstein and Shtilman [START_REF] Rubinstein | Voltage against current curves of cation exchange membranes[END_REF] (see also [START_REF] Rubinstein | Electro-osmotically induced convection at a permselective membrane[END_REF] and [START_REF] Zaltzman | Electro-osmotic slip and electroconvective instability[END_REF]). Note that the total charge in the SCR is much higher than in the EDL. The SCR coincides with the region of the depleted electrolyte solution, where K ≈ 0. Far from the SCR ρ → 0 and the electrolyte solution becomes electrically neutral. To the right of the particle surface, in the region 0 • < θ < θ 0 , where cations go out of the particle, the EDL with a negative charge density is formed. There is no SCR in this region (see Fig. 5.7(b)). Let us proceed to consider the flow pattern near the surface of an ion-selective particle.

The stream function for an axisymmetric flow is defined as, When the applied electric field is strong enough, this creates a more sophisticated electrokinetic flow mechanism called "electroosmosis of the second kind" [START_REF] Dukhin | Electrokinetic phenomena of the second kind and their applications[END_REF]. Dukhin in his work [START_REF] Dukhin | Electrokinetic phenomena of the second kind and their applications[END_REF] suggested that the well-known Smoluchowski slip formula can be applied to electroosmosis of the second kind under extreme conditions if the zeta potential is taken into account as the voltage drop inside the SCR. Using these heuristic arguments, Dukhin obtained [START_REF] Dukhin | Electrokinetic phenomena of the second kind and their applications[END_REF] an electrophoretic velocity, which is proportional to the square of the applied electric field. This analysis has been improved in [START_REF] Mishchuk | Electroosmosis of the second kind[END_REF]. In another recent work, asymptotic methods were used to obtain approximate current-voltage characteristics of a one-dimensional model of the problem. Using these characteristics to describe the transfer process in the cation-selective particle, the authors obtained a more complex model for solving the problem of electrophoresis. In their recent work, Ben et al. [START_REF] Ben | Nonlinear electrokinetics and "superfast" electrophoresis[END_REF] have found that

u = - 1 sin θ 1 r ∂ψ ∂r , v = 1 sin θ 1 r 2 ∂ψ ∂θ . a) E f -2 -1 0 1 2 -2 -1.5
U ∞ ∼ E 2/3
∞ , in the high Peclet number limit. Both of these predictions of the velocity proportionality degree with the applied electric field strength are in some qualitative agreement with the existing experimental data for different modes [START_REF] Dukhin | Electrokinetic phenomena of the second kind and their applications[END_REF][START_REF] Barany | Superfast electrophoresis of conducting dispersed particles[END_REF]. There is another work by Rubinstein and Zaltzman [START_REF] Rubinstein | Electro-osmotically induced convection at a permselective membrane[END_REF], in which they show that

U ∞ ∼ E 1/3 ∞ .

Electrokinetic instability and the transition to chaotic regime

When the external electric field becomes greater than a certain critical value

E ∞ > E *
∞ (E * ∞ ≈ 27 for ν = 0.0087), the numerical solution of the non-stationary problem effectively time dependent as t → ∞. Small perturbations superimposed on the initial fields, simulating external noise, grow and the solution is changed to some non-stationary solution: the electrokinetic instability manifests itself. In flat membranes, this instability was discovered by Rubinstein and Zaltzman [START_REF] Rubinstein | Voltage against current curves of cation exchange membranes[END_REF][START_REF] Rubinstein | Electro-osmotically induced convection at a permselective membrane[END_REF]. For a cation-exchange membrane, this instability manifests itself in the region of incoming cation flux, where ions enter the particle. The opposite side of the particle surface is stable.

Three factors make the case of the spherical particle more complicated than the case of a flat membrane: 1) the stationary solution is not one-dimensional; 2) convection of salt and liquid takes place; 3) both surfaces corresponding to the incoming and outgoing ion fluxes are parts of the same surface and the place of their separation θ = θ 0 is a singular point (see the analytical solution).

Instability is possible only in the region of the SCR and the adjacent diffusion layer (see the results of linear stability analysis in [START_REF] Ganchenko | Extreme nonequilibrium electrophoresis of an ion-selective microgranule[END_REF]). The numerical solution shows that at low supercriticality E ∞ > E * ∞ in the vicinity of the angle θ = 180 • the flow loses stability and small sinusoidal disturbances propagate in the direction of smaller angles and disappear at θ ≈ θ 0 . The flow to the right of the particle for these values of parameters remains stable [START_REF] Ganchenko | Instabilities, bifurcations, and transition to chaos in electrophoresis of charge-selective microparticle[END_REF]. the propagating perturbations of the charge density ρ take the form of spikes, as it was first described in the papers [START_REF] Demekhin | Linear and nonlinear evolution and diffusion layer selection in electrokinetic instability[END_REF][START_REF] Demekhin | Direct numerical simulation of electrokinetic instability and transition to chaotic motion[END_REF][START_REF] Demekhin | Three-dimensional coherent structures of electrokinetic instability[END_REF][START_REF] Chang | Nanoscale electrokinetics and microvortices: How microhydrodynamics affects nanofluidic ion flux[END_REF] devoted to the study of flat membranes. The boundary between the SCR and the diffusion region is quite sharp. The salt concentration K of traveling waves forms "cloud-like" structures. Fig. 5.12 compares the small neighbourhood of region II of the SCR shown in Fig. 5.11(a)

with the ρ distribution calculated numerically for a flat membrane (see [START_REF] Demekhin | Direct numerical simulation of electrokinetic instability and transition to chaotic motion[END_REF][START_REF] Chang | Nanoscale electrokinetics and microvortices: How microhydrodynamics affects nanofluidic ion flux[END_REF]).

A striking similarity of space charge density distributions is seen: inside the spikes, the charge density ρ is lowered (the spikes are "empty" of charge). Spike-like structures are connected by thin areas with a large space charge. The work of Shelistov et al. [START_REF] Shelistov | Self-similar solution to the problem of electrokinetic instability in semipermeable membranes[END_REF] showed that the opening angle of the spike does not depend on the parameters of the problem and is about θ s ≈ 120 • . This estimation is applicable not only for flat membranes but also for the more complex case of a pleated membrane (see [START_REF] Chang | Competition between dukhin's and rubinstein's electrokinetic modes[END_REF], Fig. 4). In this case, the spike-like coherent structures are not stationary: they propagate in the direction of smaller angles, and, therefore, they must be distorted by convection. Nevertheless, the opening angle of Without loss of generality, the initial time can be considered equal to zero, t = 0.

At t = 0, the space charge (see Fig. 5.13(a)) with a large amplitude is located near the interface at r = 1 and moves in the direction of increasing r, i.e. away from the surface of the particle. At t = 0.005, the charge propagates in space, it is located far from the interface and its amplitude decreases significantly. At t = 0.03 the charge with a small amplitude comes back to the particle surface, and at t = 0.07 it returns to the same position and restores its profile and value. This time is, in fact, a period, since the process is repeated periodically.

At t = 0, the salt concentration (see Fig. the concentration shock wave unfolds and eventually returns to the same position and restores its profile and value. This process is also repeated periodically.

If the particle radius changes from 25 µm to 250 µm, then for an aqueous solution The dependence of the electric current j created by the cation flow through the particle surface and normalized to its maximum j max on the angle θ is shown in Fig. 5.15 for various q q 0 q 0 q 0 q 0 t = t 0 since there is a period-doubling on the velocity graph. Note that, due to the high viscosity, the oscillations of velocity U ∞ , both periodic and chaotic, have quite small amplitude, no more than 5% -6% of the average velocity. Therefore, their detection in experiments is difficult.

Another important characteristic of the problem is θ 0 -the diffusion layer separation angle. Fig. 5.17 compares the angle θ 0 obtained using the semi-analytical approach with similar results of numerical simulation.

For the parameter χ > 1.28 the stationary semi-analytical solution is in good agreement with the numerical results. As noted above, the oscillations of unknowns, both periodic and chaotic, have a small amplitude, which is no more than 6% of the average value. In Fig. 5.17, the averaged values of θ 0 = 1/(t 2 -t 1 )

t 2
t 1 θ 0 dt are shown for a sufficiently long averaging interval. The discrepancy between the numerical and analytical approaches for χ < 1.28 can be explained by the fact that the analytical results go far beyond the limits of numerical modeling and relate to the regime of overlimiting currents (for a fixed ν and for small values of χ, the value E ∞ will be small, which contradicts the hypothesis E ∞ → ∞ used in the analytical approach). The dashed line stands for the DNS results and the solid line stands for the analytical solution.

of the electrophoretic velocity U ∞ normalized by the parameter κ vs parameter χ. For the semi-analytical approach, the dependencies are shown for κ = 0.1 and 1 (curves 1 and 2, respectively). This range of the parameter κ was chosen because it covers the most commonly used in practice in electrolyte solutions. Fig. 5.21 also shows the dependence obtained from the numerical solution (curve 3 for κ = 0.26). There is good agreement between the numerical and analytical approaches, except for the range of small values of χ. But, as noted above, this is since the assumption of the analytical approach E ∞ → ∞ is not fulfilled (the interval of small values of χ, according to Dukhin's terminology, corresponds to electrophoresis of the first kind).

Based on physical considerations, Dukhin derived the well-known formula for electrophoresis of the second kind [START_REF] Dukhin | Electrokinetic phenomena of the second kind and their applications[END_REF], which in dimensional form is the following:

u ∞ = 2 εaE 2 ∞ µ .
(5.54)

In our dimensionless formulation, this formula (5.54) becomes the following form: It should be noted that in experimental works a fairly wide range of electrolyte concentrations from 0.01 to 1 mol/m 3 is taken. However, at such low concentrations the electrolyte can still be considered highly diluted and the changes in dynamic viscosity μ and absolute electrical permittivity ε, compared to the water, can be neglected.

U ∞ κ = 2. ( 5
The most important quantities -the electrophoretic velocity U ∞ and the particle mobility U ∞ /E ∞ are shown in Fig. 6.2.

For E ∞ < 1, the analytical dependence gives a good approximation with the experimental data. The significant discrepancy between the numerical and analytical curves can The solid line corresponds to the numerical solution, the dashed line to the analytical solution (4.98). Experimental data from [START_REF] Mishchuk | Theoretical and experimental study of nonlinear electrophoresis[END_REF] are labeled with round markers; triangular markers stand for data from [START_REF] Barany | Electrophoresis in strong electric fields[END_REF].

be easily explained: it is obvious that the assumption E ∞ → 0, proceeding from which the analytical solution was obtained, is not satisfied.

In addition to the comparison of the particle velocity itself, we also compared the electroosmotic slip velocity, which is shown in Fig. in the outgoing ion flux region, it becomes negative. As a result, the velocity profile of electroosmotic slip changes dramatically: the maximum velocity shifts towards large angles θ and is reached at θ ≈ 135 • . Moreover, for small angles the velocity changes its direction at E ∞ = 3 and 5.

Strong electric field

The results of direct numerical simulation and semi-analytical analysis for the strong electric field were compared with experimental data. Experimental data from the works of Barani et al. [START_REF] Baran | Superfast electrophoresis of ion-exchanger particles[END_REF], Mishchuk and Takhistov [START_REF] Mishchuk | Electroosmosis of the second kind[END_REF], Mishchuk and Dukhin [START_REF] Mishchuk | Electrokinetic phenomena of the second kind[END_REF], Mishchuk and Barinova [START_REF] Mishchuk | Theoretical and experimental study of nonlinear electrophoresis[END_REF] were taken and analyzed. In these works, the radius of the cationexchange particle, ã, varies in the range from 0.5 to 600 µm; the concentration of the NaCl solution is 10 -4 mol/l. The external electric field strength ranges from 1 V/m to 100 kV/m. In addition, the ion-exchange particles in the experiments were fabricated from various types of ion-selective materials.

The dimensional slip velocity, ũm , along the particle surface is shown in Fig. 6.4. The velocity field from the numerical simulation is taken at the point ỹm . The experimental points are taken from the work of Mishchuk and Takhistov [START_REF] Mishchuk | Electroosmosis of the second kind[END_REF]. In the incoming ion flux region, the maximum velocity is reached at an angle θ ≈ 60 • . The theoretical maximum velocity at this point, as well as the profile of the dependence, are in good agreement with the experimental results. In the outgoing ion flux region, the electroosmotic velocity is zero and this velocity cannot be measured experimentally. Here, the angle is counted backwards to match the experiments.

The dependence of the electrophoretic velocity on the electric field is of great practical interest, therefore, all data on particle radii, electric field strengths, and works from which these values were taken were grouped in Table 6.1.

∞ ∼ const 1 E ∞ , u ∞ ∼ const 2 E 4/3 ∞ and u ∞ ∼ const 2 E 2 ∞ ,
∞ = 2 κ E 2
∞ , which corresponds to the line 3 in Fig. 6.5. The experimental data deviation from the Dukhin's prediction is seen. The dependence

u ∞ ∼ const 2 E 4/3
∞ (line 2) from our work fits much better the experimental data.

For a better understanding, let us present the normalized experimental and theoretical electrophoretic velocity

U ∞ = u ∞ /E 2 ∞ as a function of the universal variable χ = ν 2/3 E 4/3
∞ . This dependence is shown in Fig. 6.6.

As in the previous graph, the vertical dotted lines represent the conventional boundaries of low I, medium II, and high III electric field strength regions. The numerical solution does not work at a sufficiently high electric field strength, so in Fig. 6.6 the solid line breaks off at χ ≈ 8. Meanwhile, the semi-analytical solution complements the numerical one and can be extended to any arbitrarily large value of parameter χ. The discrepancy between the analytical and numerical solutions for small values of χ is almost two times, but as χ grows, this difference rapidly decreases to almost zero. For χ > 8, only the analytical method works, which is in very good agreement with the experimental The markers for the experimental points are given in Table 6.1; the solid line stands for the results of our DNS at ν = 0.002, and the dashed line stands for the semi-analytic results. Inset A: behavior of u ∞ for small E ∞ , the solid line stands the DNS, and the dashed line stands for the analytic solution. Inset B: behavior of u ∞ for large E ∞ .

data.

As shown earlier, the entire semi-analytical solution and, in particular, the electrophoretic velocity, strongly depends on the phenomenon of diffusion-convection in the diffusion layer. As a result, for a fixed value of the parameter κ, the normalized electrophoretic velocity U ∞ is a function of the χ parameter,

U ∞ = f (χ). ( 6.1)
The semi-analytical approach for large values of χ show that for χ → ∞ the following expression holds:

U ∞ ∼ 1 √ χ . ( 6.2)
Unfortunately, it was not possible to derive a simple analytical expression relating the ion-selective particle velocity to the strength of the applied electric field. In Dukhin's approach, the influence of the diffusion layer is neglected and there is no dependence on χ (see equation (5.55)). This may be the main reason for the discrepancy between Dukhin's prediction (5.55) and experiments. In the experiments shown in Fig. 6.6, one can see that U ∞ is not constant at χ → ∞, but decreases as χ increases, which is in full compliance with the semi-analytical approach. The expression (6.2) can also be written in the following form,

U ∞ ∼ E -2/3 ∞ or u ∞ ∼ E 4/3 ∞ ,
It is in good agreement with experiments (see inset B of Fig. 6.6).

Thus, it can be concluded that the electrophoretic velocity is proportional to E 4/3 ∞ , and not to E 2 ∞ , as it is assumed in Dukhin's works.

Conclusion

In the present work, the electrostatics and hydrodynamics near an ion-selective microparticle in an electrolyte solution are studied analytically for the limiting cases of a weak and strong electric field, and numerically in a wide range of parameters. Mathematically, this problem is described by a complex system of partial differential equations - The list of main new results presented in this work are as follows:

1. Analytical solutions of the problem applicable for weak electric fields (the first kind phenomena) have been found.

2. An analytical formula for the electrophoretic velocity for the first kind phenomena is derived. This formula is a generalization of the well-known Helmholtz-Smoluchowski formula for the case of ion-selective microparticles.

3. The semi-analytical solution of the electrophoresis problem in a strong electric field (the second kind phenomena) has been found.

4.

A correction to the semi-empirical dependence of the electrophoretic velocity vs electric field strength (Dukhin's quadratic dependence) was obtained for the limit E ∞ → ∞.

5. Exact numerical solutions of the problem applicable to a wide range of parameters have been found.

A previously unknown result of a special type of electrohydrodynamic instability

on the surface of a microparticle was analysed.

7. Successful comparisons of analytical models with the numerical solutions and experimental data have been carried out.

Prospects for further development of the topic

The next stage in the study of ion-selective particle electrophoresis can be a theoretical study of the effect of increased ion concentration in the region of outgoing ion flux. This effect was discovered during this research and appears only at a sufficiently high electric field strength (above a critical one). However, its impact on the problem was not taken into account. This effect is interesting since it can find direct application in the field of medical diagnostics, namely, it can be used for the preconcentration of samples and subsequent detection of substances in trace amounts. Such substances, in particular, can be markers of various diseases, which are actually charged particles and at an early stage of the development of the disease, their concentration is excessively low. Thus, this effect is an directly related to this potential application.

It is also linked to the desalination problem as at the rear pole of the sphere, the concentration is very high in opposite to the front pole of the particle. ψ(i + 1, j) -ψ(i, j) rm1(i + 1) -ψ(i, j) -ψ(i -1, j) rm1(i) = = α(i) ψ(i -1, j) + β(i) ψ(i, j) + γ(i)ψ(i + 1, j) α(i) = 1 rn1(i) rm1(i)

List of Figures

; γ = 1 rn1(i) rm1(i + 1)

; β = -α(i) -γ(i)

∂ 4 Ψ ∂r 4 i = ∂ 2 ∂r 2
∂ 2 Ψ ∂r 2 = α(i) ∂ 2 ∂r 2 (ψ(i -1, j)) + β(i) ∂ 2 ∂r 2 (ψ(i, j))+ +γ(i) ∂ 2 ∂r 2 (ψ(i + 1, j)) = = α(i)[α(i -1) ψ(i -2, j) + β(i -1) ψ(i -1, j) + γ(i -1) ψ(i, j)]+ +β(i)[α(i) ψ(i -1, j) + β(i) ψ(i, j) + γ(i) ψ(i + 1, j)]+ +γ(i)[α(i + 1) ψ(i, j) + β(i + 1) ψ(i + 1, j) + γ(i + 1) ψ(i + 2, j)] = = α(i) α(i -1) • ψ(i -2, j) + (α(i) β(i -1) + β(i) α(i)) • ψ(i -1, j)+ +(α(i) γ(i -1) + β 2 (i) + α(i + 1) γ(i)) • ψ(i, j)+ +(β(i) γ(i) + β(i + 1) γ(i)) • ψ(i + 1, j) + γ(i) γ(i + 1) • ψ(i + 2, j)

So the four-order term r 4 ∂ 4 ∂r 4 can be expressed in the following way, a 5 (i) ψ(i -2, j) + b 5 (i) ψ(i -1, j) + c 5 (i) ψ(i, j) + d 5 (i) ψ(i + 1, j)+ +e 5 (i) ψ(i + 2, j)

where coefficients defined as follows, a 5 (i) = rn 4 (i) • α(i) α(i -1) b 5 (i) = rn 4 (i) • α(i) (β(i -1) + β(i))

c 5 (i) = rn 4 (i) • (α(i) γ(i -1) + β 2 (i) + α(i + 1) γ(i))

d 5 (i) = rn 4 (i) • γ(i) (β(i) + β(i + 1))
e 5 (i) = rn 4 (i) • γ(i) γ(i + 1)

Now consider the second-order term r 4 ∂ 2 ∂r

1 r 2 + r 2 ∂ 2 ∂ Ψ,
rn 4 (i) α(i) ψ(i -1, j) rn 2 (i -1) + β(i) ψ(i, j) rn 2 (i) + γ(i) ψ(i + 1, j) rn 2 (i + 1) + +rn 2 (i) [α(i) ψ(i -1, j) + β(i) ψ(i, j) + γ(i) ψ(i + 1, j)] = = rn 2 (i) α(i) 1 + rn 2 (i) rn 2 (i -1) ψ(i + 1, j) + 2rn 2 (i) β(i)ψ(i, j)+ +rn 2 (i) γ(i) 1 + rn 2 (i) rn 2 (i + 1) ψ(i + 1, j) = = b 3 (i) ψ(i -1, j) + c 3 (i) ψ(i, j) + d 3 (i) ψ(i + 1, j)

where, b 3 (i) = rn 2 (i) α(i) 1 + rn 2 (i) rn 2 (i -1) c 3 (i) = 2 rn 2 (i) β(i)

d 3 (i) = rn 2 (i) γ(i) 1 + rn 2 (i) rn 2 (i + 1)
The boundary conditions at the outer edge of calculation domain were found from the

  4. L'émergence d'une instabilité électrocinétique dans la région du flux entrant de cations a été trouvée. Un scénario de transition de flux d'un régime régulier à chaotique est obtenu.

  Mots clés : Microparticule ionique sélective, électrophorèse, électroosmose, électrolyte, équations de Nernst-Planck-Poisson-Stokes, vitesse électrophorétique, électrophorèse du second type, électrophorèse non linéaire, gradient de conductivité électrique, instabilité électrohydrodynamique Institut de Mécanique et d' Ingénierie de Bordeaux (I2M, CNRS UMR 5295) ENSCBP bat A, 16 avenue Pey-Berland, 33607 PESSAC Cedex Title : Non-equilibrium electrophoresis of ion-selective microparticle Abstract :

Keywords:

  Ion-selective microparticle, electrophoresis, electroosmosis, electrolyte, Nernst-Planck-Poisson-Stokes equations, electrophoretic velocity, electrophoresis of the second kind, non-linear electrophoresis, electric conductivity gradient, electrohydrodynamic instability Institut de Mécanique et d' Ingénierie de Bordeaux (I2M, CNRS UMR 5295) ENSCBP bat A, 16 avenue Pey-Berland, 33607 PESSAC Cedex

  Les méthodes utilisées pour mettre en oeuvre un mouvement fluide à une échelle conventionnelle sont souvent basées sur les forces d'inertie. A l'échelle microscopique, les forces visqueuses dominent et le mouvement du fluide a plus de mal à être initié. Les phénomènes électrocinétiques constituent l'un des moyens non mécaniques les plus répandus et les plus efficaces pour faire bouger un fluide à l'échelle micro et nanométrique. L'idée principale est la suivante : un liquide dans lequel se trouvent des ions d'une substance dissoute (électrolyte) est adjacent à une surface solide chargée. Un nuage de diffusion d'ions de charge opposée s'accumule près de cette surface. Un champ électrique externe produit alors une force agissant sur la couche de diffusion chargée, ce qui provoque un écoulement de fluide par rapport à la surface chargée. L'écoulement électrocinétique près d'une surface fixe (paroi) est appelé écoulement électroosmotique, et le mouvement électrocinétique des particules chargées en suspension dans un fluide est appelé électrophorèse. L'électroosmose et l'électrophorèse ont de nombreuses applications en chimie analytique, dans les problèmes de mélange de liquides, de modification des propriétés thermoconductrices d'une substance, de séparation des particules se déplaçant dans un liquide par les taille et par d'autres propriétés ou dans le domaine biomédical pour n'en citer que quelques-unes. L'objectif dans cette thèse est d'étudier théoriquement et numériquement le mouvement d'une microparticule sphérique sélective aux ions dans un champ électrique constant dans une large gamme d'intensités. Pour atteindre cet objectif, il était nécessaire de : 1. Formuler théoriquement le problème du mouvement d'une microparticule sélective d'ions (perméable à un type d'ions) dans une solution d'électrolyte binaire (ions positifs et négatifs) lorsqu'un champ électrique externe est appliqué.

7 .

 7 Tester les relations théoriques obtenues : les comparer avec les résultats obtenus par simulation numérique directe du problème et faire également une comparaison avec les données expérimentales disponibles. Au début, une étude bibilographique de l'état de l'art actuel du problème est faite du point de vue de la théorie et des expériences. Ensuite, la formulation mathématique du problème est discutée. Le modèle géométrique consiste en une microparticule sphérique, échangeuse de cations à travers sa surface, de rayon a, située dans un réservoir infiniment étendu rempli d'un fluide visqueux de type Newtonien incompressible, considéré comme un électrolyte hautement dilué. Dans ce cas, un électrolyte est un liquide diélectrique dans lequel des ions d'un certain sel sont dissous (par exemple, NaCl ou KCl). Le problème est étudié sur la base des hypothèses suivantes : un électrolyte binaire (ne contenant que 2 types d'ions) et symétrique (les nombres de charges ou les valences d'ions sont identiques en module, z + = -z -= 1), les coefficients de diffusion des ions positifs et négatifs sont les mêmes (D + = D -= D), la particule est idéalement sélective (le flux d'anions à travers la particule est totalement absent). Dans ce système, un champ électrique d'intensité constante E ∞ agit, ayant des composantes normales et tangentielles par rapport à la surface de la particule. La composante normale provoque une augmentation de la charge d'espace sur la surface et la composante tangentielle provoque le mouvement du fluide le long de la surface de la particule. Lorsqu'un champ électrique est appliqué, la particule commence à se déplacer et acquiert la vitesse électrophorétique U ∞ . Cette vitesse peut être obtenue à partir de l'équilibre des forces agissant sur la particule, à savoir les forces visqueuses et électriques. Il convient de noter ici qu'un courant électrique traversant la particule est provoqué par le flux d'ions, en l'occurrence de cations. La description mathématique du problème est basée sur un modèle continu. Le modèle mathématique du problème correspond à un système non linéaire d'équations aux dérivées partielles. L'équation de Poisson décrit la distribution du potentiel électrique. Le mou-vement d'un fluide est décrit par l'équation de Stokes (très faibles vitesses). Le transfert d'ions est décrit par l'équation de Nernst-Planck, qui inclut les trois mécanismes de transfert d'ions possibles : diffusion, convection et électromigration.

3 ∞

 3 La vitesse électrophorétique est trouvée sur la base de l'équilibre des contraintes visqueuses et des contraintes électrostatiques de Maxwell agissant sur la particule comme dans le cas d'un champ électrique de faible intensité. Cependant, son obtention sous forme analytique est difficile en raison de la non-linéarité du système d'équations, ce qui nous a conduit à la recherche d'une solution auto-similaire. L'étude de la stabilité linéaire d'un écoulement stationnaire a plusieurs objectifs : 1) identifier les mécanismes physiques de base du comportement du système ; 2) identifier le mécanisme responsable de l'apparition d'une instabilité ; 3) obtenir des données pour la vérification d'expériences numériques.Pour discrétiser le problème stationnaire et le problème de stabilité linéaire le long d'une coordonnée spatiale dirigée perpendiculairement à la surface de la frontière, la méthode de Galerkin a été appliquée.Parmi les résultats les plus significatifs, il convient de mentionner le détachement de la couche de diffusion sous un angle θ = θ 0 (ce résultat a été prédit par Levich pour d'autres problèmes) et la découverte de l'instabilité d'une solution stationnaire en couches.La modélisation numérique a montré un bon accord avec les résultats analytiques.Sur la face frontale de la particule (la région des ions entrants), une zone de dessalement est formée, où K ≈ 0. Comme le prévoit la théorie, une mince couche de diffusion électroniquement-neutre se forme à la limite de la zone de dessalement. L'analyse asymptotique prédit la séparation de la couche limite avec un angle θ = θ 0 ≈ 70 -75 • . Il a également été montré que, dans la région du flux de cations entrant, une zone de dessalement se formait, que la charge d'espace se déformait et finissait par atteindre un maximum à une certaine distance de la surface de la particule. Une telle charge est appelée une charge sortante en électrochimie. La distribution de la fonction de courant Ψ montre la formation du très connu vortex électroconvectif de Dukhin-Mishchuk à θ ≈ θ 0 . L'intensité du vortex augmente avec l'augmentation de E ∞ . Les résultats sont en bon accord qualitatif avec les expériences de Mishchuk et Takhistov. La cause de la formation d'un vortex toroïdal est l'inhomogénéité de la répartition de la charge à la surface et, par conséquent, l'inhomogénéité de la vitesse de glissement à la surface de la particule. La modélisation numérique a confirmé le phénomène intéressant et jusqu'alors inconnu pour les petites particules de la perte de stabilité de la solution à zéro nombre de Reynolds. Lorsque le champ électrique externe dépasse une certaine valeur critique, E ∞ > E * ∞ , la solution numérique du problème non stationnaire n'atteint pas une solution stationnaire lorsque t → ∞. De petites perturbations superposées aux données initiales simulant l'augmentation du bruit externe entraînent une solution non stationnaire : l'instabilité électrocinétique se manifeste alors. Dans les membranes planes, cette instabilité a été découverte par Rubinstein et Zalzmann. Pour une membrane échangeuse de cations, cette instabilité apparaît du côté de la surface par laquelle les ions pénètrent dans la particule. Le côté opposé de la surface est stable. Les résultats de nos simulations numériques et des analyses semi-analytiques ont été comparés aux données expérimentales. Les données expérimentales des travaux de Barani et al., Mishchuk et Takhistova, Mishchuk et Dukhin, Mishchuk et Barinova ont été utilisés. Dans ces travaux, le rayon de la particule échangeuse de cations varie entre 0.5 et 600 µm ; la concentration de la solution de NaCl, pour laquelle des comparaisons ont été faites, est de 10 -4 mol/l. L'intensité du champ électrique externe est comprise entre 1 V /m et 100 kV /m. De plus, les particules échangeuses d'ions utilisées dans les expériences étaient fabriquées à partir de divers types de matériaux. La dépendance de la vitesse des particules avec le champ électrique possède un grand intérêt pratique, de sorte que la comparaison est effectuée pour la vitesse électrophorétique et la mobilité de ces particules. Un bon accord a été trouvé pour les forces de champs électriques de faibles et fortes intensités. En particulier, pour un champ électrique de forte intensité, la formule de Dukhin, U ∞ = 2κE 2 ∞ , a montré un écart avec les données expérimentales. Dans le travail présent, on a montré que U ∞ ∼ E 4/et est en très bon accord avec les données expérimentales. "Dedicated to my mother Anna Frants" Acknowledgement My long PhD journey was the most interesting and difficult period in my life. I owe this opportunity to two wonderful people -my advisors -Professor Sakir Amiroudin and Evgeny Demekhin. It is not possible to express how lucky I am to have two advisors, two brilliant, but very different people. Professor Demekhin has shown me how passionate one can be about science. His energy and constant movement forward are always amazing to me. He thinks big and prefers not fo focus on detail. Professor Amiroudine has shown me another side. His appreciation of a very deep understanding and paying great attention to detail reveal how it is important to know every small detail in the whole picture. I am deeply grateful to both of them for guiding and supporting me all this time, as well as for their patience with me. This work would not be possible without the financial support of the French Government and the Embassy of France in Russia through the "Vernadski" program. I am grateful to the Campus France team and especially Mr. Farid Saadoun, who organized my stay in France and did everything to make my life there comfortable. Colleagues have also greatly influenced me and my work. I want to thank all collective of I2M Lab for the warm welcome and good attitude. Special thanks to Madam Valerie Thouard for helping me with administrative procedures and to Deewakar and Vibhor for helping me with my first steps in a new place. I would like to express sincere thanks to my Russian colleagues Dr. Georgy Ganchenko and Dr. Vladimir Shelistov for their interest and concern for this work.I would further like to expand my sincere thanks to a few, but very true friends: Katya and her sister Linda and their relatives for treating me as one of their own family. I also want to thank my parents, especially my mother who has always supported me in all phases of my life's journey. She knows how to give me strength in the most difficult circumstances and always believes in me. Thank to my elder brother Vladimir who has encouraged me to go through this PhD journey.Many thanks to my husband Aleksey, who had been waiting for me for the last three years and had courageously endured time without me during my staying in France.

Finally, I would

  like to thank Dr. E. Sobolev from Financial University, Dr. E. Stroganova and Dr. I. Rayushkina from KubSU who helped me a lot in different aspects of my PhD program.

  1.1). This is because cation exchange membranes are only permeable for cations. The anion exchange membranes perform the opposite way compared to cation exchange membranes (see Fig.1.1).

Figure 1 . 1 - 2 An overview of electrokinetic phenomena 2 . 1

 11221 Figure 1.1 -Scheme of cation and anion exchange membrane.

Figure 2 . 1 -

 21 Figure 2.1 -Scheme of EDL for (a) flat surface and (b) spherical surface.

Figure 2 . 2 -

 22 Figure 2.2 -Three types of particles in an electrolyte solution under the influence of an external electric field: (a) dielectric particle, (b) metal particle and (c) ionselective particle.

a 2

 2 /D: characteristic time; D/a: characteristic velocity; Φ 0 = RT /F : thermal potential; c ∞ : equilibrium concentration of ions away from the particle; µD/a 2 : characteristic pressure;

Fig. 3 .

 3 Fig. 3.1 illustrates the example of staggered grid. A staggered mesh is used to represent the mesh functions corresponding to continuous functions: electric potential, two velocity components, and two concentrations. The potential-nodes and concentration-nodes are in the geometrical centers of the cells. The velocity-nodes are at the centers of faces of the cell, shifted from the center-nodes at the half-cell distance towards the r or θ directions.

Figure 3 . 1 -

 31 Figure 3.1 -Staggered grid for numerical scheme. The electrostatic potential is defined in nodes i -1, i and i + 1. The velocity component is defined in center nodes i -1/2 and i + 1/2.

  and multiply Eq. (3.15) by -

  By adding Eq. (3.18) and (3.19), us introduce the stream function in spherical coordinates, U = -1 r sin θ ∂Ψ ∂r ; V = 1 r 2 sin θ ∂Ψ ∂θ and substitute the stream function into ω, By multiplying Eq. (3.20) by r 3 sin θ gives, r 2 ∂ 2 Ψ ∂r 2 -sin θ ∂ ∂θ 1 sin θ ∂Ψ ∂θ = ω r 3 sin θ Substitute ω into Stokes equation, and after some simplification it leads to a fourth-order equation in Ψ,

  [START_REF] Yariv | Electro-osmotic flows over highly polarizable dielectric surfaces[END_REF]) is a O(τ ) perturbation to the right-hand side of(3.24), which leads to a O(τ 2 ) variation in w * and in F * . Then, the right-hand side of (3.28) is a O(τ 2 ) perturbation to the right-hand side of (3.25), which leads to a O(τ 3 ) variations in w * * and in F * * . At last, the right-hand side of (3.29) is a O(τ 3 ) perturbation to the right-hand side of (3.26), causing a O(τ 4 ) variation in w n+1 .

Figure 4 . 1 -

 41 Figure 4.1 -Structure of electrical double layer in weak electric field.

(

  or overlap) region, where O(ν) z +∞. This condition is that the outer limit of the inner solution has to match the inner limit of the outer solution,

  44) -(4.48), it is necessary to specify four boundary conditions (two for the function Φ and one for c + and c -), but the boundary conditions (4.47) -(4.49) are only three. The fourth boundary condition will be determined from the matching condition of the outer expansion with the inner.

  is the function of p, which can be calculated numerically. It is possible to give an approximation of the function δ(p) using a simple analytical expression, δ ≈ 1 -0.11 ln(p) = 1 -0.11 ζ. (4.80)

Figure 4 . 2 -

 42 Figure 4.2 -The geometry of the projections of normal and tangential vectors to the particle surface r = 1.

  .93) After substituting (4.84), (4.91) -(4.93) into (4.90) and after some transformations, one can obtain the expression for viscous and Maxwell stresses,

Fig. 4 .

 4 Fig.4.3(a) shows the space charge distribution ρ(r, θ) = c + -c -at E ∞ = 0.05. This value of E ∞ is considered to be small, therefore ρ should correspond to the analytical

Figure 4 . 3 -

 43 Figure 4.3 -(a) Space charge distribution, ρ(r, θ) = c + -c -, for a weak electric field, E ∞ = 0.05 and p = 3. (b) The cross-section of ρ(r) for the angle θ = 0 (solid line) and θ = 180 • (dashed line); 1 is the solution obtained from first-order expansion.

Figure 4 . 4 -

 44 Figure 4.4 -(a) Electrical conductivity, K(r, θ) = c + + c -, for a weak electric field, E ∞ = 0.05 and p = 3. (b) The cross-section of K(r) for the angle θ = 0 (solid line) and θ = 180 • (dashed line); curve 1 is the solution obtained by first-order expansion. The behavior inside the EDL is shown in the insets.

  4.6 and 4.5 are shown the distribution of the space charge density ρ and the concentration of salt K near the surface of the ion-selective microparticle. The charge density distribution (or EDL thickness) becomes asymmetric.

Figure 4 . 5 -

 45 Figure 4.5 -(a) Space charge distribution, ρ(r, θ), for moderate electric field strength, E ∞ = 3, p = 3. (b) The section ρ(r) for θ = 0 • (solid line) and θ = 180 • (dashed line).

Figure 4 . 6 -

 46 Figure 4.6 -(a) Electrical conductivity, K(r, θ) = c + + c -, for moderate electric field strength E ∞ = 3, p = 3. (b) The cross-section K(r) for θ = 0 (solid line) and θ = 180 • (dashed line). The behavior inside the EDL is shown in the insets.

  4.3). In the case of higher values of electric field the tangential component of the electric field E s = ∂Φ ∂θ at the boundary of the EDL is subject to the equation (4.71) and it is proportional to sin θ. The quasi-equilibrium slip velocity is directed clockwise for angle 0 • < θ < 180 • , reaches a maximum at θ = 90 • and disappears at θ = 0 • and θ = 180 • (see Fig. 4.7), where the tangential component of the electric field is zero (see equation (4.79)).

Figure 4 . 7 -

 47 Figure 4.7 -Direction of electroosmotic slip velocity along the particle.

Fig. 4 .

 4 Fig. 4.8 represents the dependence of the electrophoretic velocity on the external electric field for various values of the parameter p. The curves for different values of p approach each other with increasing electric field strength and ultimately collapse into one (at the value E ∞ ≈ 5). Consequently, electrokinetic phenomena are now completely

Figure 4 . 8 -

 48 Figure 4.8 -Electrophoretic mobility vs the electric field strength E ∞ for different values of p.

Figure 4 . 9 -

 49 Figure 4.9 -The dependence of the particle velocity on the cations concentration p on its surface for electric field strength E ∞ = 5.

Fig. 4 .

 4 Fig. 4.9 complements the picture of Fig. 4.8. The particle velocity U ∞ dependence on parameter p at the fixed electric field strength E ∞ = 5 becomes almost independent on p

Figure 5 . 2 -

 52 Figure 5.2 -I, II, and III correspond to thin boundary layers embedded in each other: EDL, space charge region, and diffusion layer, respectively. IV corresponds to the electrically neutral bulk of the electrolyte, where K = 2. The points 1 and 1 separate the areas of incoming and outgoing ion fluxes. Inset A: Typical structure of the incoming ion flux region. Inset B: Typical structure of the outgoing ion flux region.

Figure 5 . 3 -

 53 Figure 5.3 -Schematic diagrams of bifurcation. (a) The solution of the equation (5.11)has an imperfect pitchfork bifurcation, which is structurally unstable[START_REF] Arnold | Geometrical Methods in the Theory of Ordinary Differential Equations[END_REF] (b)In the solution of the equation (5.10) for small but finite values of ν, this type of bifurcation is destroyed. The solid line corresponds to the solution that has physical meaning.

  [START_REF] Van Der Put | Electrical conductivity of dilute and concentrated aqueous dispersions of monodisperse polystyrene particles. influence of surface conductance and double-layer polarization[END_REF] into(5.20), one can obtain the normal velocity component on the edge of SCR,V m = -U m ∂y m

. 23 )

 23 At y = y m at the boundary of the desalination zone, the salt concentration is assumed to be zero, and the ion flux normal to the particle surface and determined by diffusion is, y = y m : K = 0, ∂K ∂y = j.(5.24)At the outer boundary of the diffusion layer, the salt concentration tends to an equilibrium value, y -y m δ : K = 2.(5.25)In a thin diffusion layer, the tangential velocity can be set equal to a constant with respect to y, equal to the electroosmotic slip velocity, U = U m (X), but it varies along the surface of the particle. To find the normal component of the velocity V , one needs to use the equation of mass balance (2.31), 26), one can get, y ym(X,t) ∂ ∂y V sin X dy -y ym(X,t) ∂ ∂X U sin X dy = 0, (5.27)

  .30) even for the stationary case, ∂/∂t = 0, would greatly complicate the model if the problem did not have an analytic self-similar solution. One can look for a solution to the equation (5.29) in the following self-similar form, K = K(Y ) where Y = η -X, t) is an arbitrary self-similar function. After substituting (5.31) into (5.29) and separating b-terms and δ-terms, one can obtain a one-dimensional partial differential equation, sin X = W m sin X. (5.33) The third algebraic equation can be obtained from the solution of the ODE (5.31), K = const 1 erf (Y ) + const 2 by substituting the boundary condition, (5.30) in (5.31),

µm 10 4 × 10 -

 410

  k+1 P k (cos s)P k (cos θ) sin sds(5.43)

Fig. 5 .

 5 Fig. 5.4 shows the distribution of K at E ∞ = 5 (moderate electric field). In a thin EDL, K(r) decays exponentially, but its change slows down outside the EDL.

Figure 5 . 4 -

 54 Figure 5.4 -(a) The electric conductivity K(r, θ) for the moderate electric field strength E ∞ = 5. The dashed line stands for the outer edge of the diffusion layer. (b) Cross-section K(r) at fixed angles θ = 0 • and θ = 180 • . The behavior inside the EDL (the depleted and enriched regions correspondingly) is shown in the insets.

Figure 5

 5 Figure 5.5 -(a) The electric conductivity K(r, θ) for E ∞ = 10. The dashed line stands for the outer edge of the diffusion layer. The values have been clipped from the actual maximum down to K = 5 in order to achieve contrast with the diffusion layer. (b) Cross-section of K(r) at fixed angles θ = 0 • and θ = 180 • . The behavior inside the EDL (the depleted and enriched regions correspondingly) is shown in the insets.

Figure 5 . 6 -

 56 Figure 5.6 -(a) The electric conductivity K(r, θ) for E ∞ = 20. The dashed line stands for the outer edge of the diffusion layer. The values have been clipped from the actual maximum down to K = 5 in order to achieve contrast with the diffusion layer. (b) Cross-section of K(r) at fixed angles θ = 0 • and θ = 180 • .The behavior inside the EDL (the depleted and enriched regions correspondingly) is shown in the insets.

  (2.34)). The numerical estimation of this angle is θ 0 ≈ 76 • . At moderate and large values of E ∞ , the non-equilibrium effects become dominant and are responsible for the formation of space charge. The space charge consists of a thin double layer (EDL) and the space charge region itself with a maximum located far from the particle surface (see Fig. 5.7(b)).

  a) The charge density ρ(r, θ) for the moderate electric field strength E ∞ = 5. (b) Cross-section of ρ(r) at fixed angles θ = 0 • and θ = 180 • near the particle surface.

Figure 5 . 8 -Figure 5 . 9 -

 5859 Figure 5.8 -(a) The charge density ρ(r, θ) for E ∞ = 10. (b) Cross-section of ρ(r) at fixed angles θ = 0 • and θ = 180 • near the particle surface.

Figure 5 . 10 -

 510 Figure 5.10 -Theoretical stream function distribution for (a) E ∞ = 5, (b) E ∞ = 10, (c) E ∞ = 20. (d) The picture of electrokinetic flow from [69] around a 1-mm cation selective particle at Ẽ∞ = 100 V/cm, which corresponds to dimensionless E ∞ = 200 and Debye number ν ≈ 8.68 • 10 -5 .

Figure 5 .

 5 Figure 5.11 -(a) The charge density ρ = c + -c -distribution and (b) electric conductivity K = c + +c -distribution at E ∞ = 30. I and II -regions of the standing and travelling waves, respectively. Arrows show the direction of propagation.

Fig. 5 .

 5 Fig. 5.11 shows the charge density ρ and the salt concentration K far from the instability threshold for a strongly non-linear regime, E ∞ = 30. Periodic oscillations of the charge density and salt concentration arise in the region I (see Fig. 5.11). These oscillations are perpendicular to the surface of the microparticle and are quite strong. They completely destroy the classical structure of the SCR, the depletion region of the electrolyte, and the diffusion layer, but these oscillations do not affect the EDL. They can be interpreted as highly non-linear standing waves. At the boundaries of the region I, as well as in the region II, these standing waves cause non-linear traveling waves. Fig. 5.11(a) shows that

θ

  s is close to 120 • , this angle did not change in other numerical calculations with different ~ 120F igure 5.12 -(a) Blow-up of region II of Fig. 5.11(a). (b) Spikes for a flat membrane [45, 72]. The normal to the membrane coordinate is strongly compressed in comparison to the direction along the membrane. values of E ∞ . In the region I in Fig. 5.11, all unknowns change periodically with time, forming a standing wave near the pole of the particle. The charge density ρ and the salt concentration K at θ = 180 • as a function of the radius r for different times are shown in Fig. 5.13.

Figure 5 .

 5 Figure 5.13 -(a) Cross-sections of ρ(r) at θ = 180 • at different time instants. (b) Crosssections of K(r) at θ = 180 • at different time instants. The wave is bouncing back from the outer edge of the diffusion layer, which thickness defined by δ.

N

  aCl the dimensional period changes from 0.05 sec. up to 5 sec. (the corresponding typical oscillation frequency ranges from f = 20 Hz to f = 0.2 Hz). The velocity field around the particle completes the picture. An image of the stream function for an unstable regime is shown in Fig. 5.14. A standing wave at θ = 180 • creates Rubinstein-Saltzman microvortices in the region I, which move towards a decrease in the angle θ. As for the flat membrane [44, 72], these microvortices are located in the diffusion layer. When the angle θ = θ 0 is reached, the Rubinstein-Zaltzman microvortices ultimately merge with the Dukhin-Mischuk (region II) vortices and completely disappear. Microvortices in region I cannot pass through the Dukhin-Mischuk vortices (II), and therefore the region of the outgoing ion flux 0 < θ < θ 0 is completely motionless.

Figure 5 . 14 - 1 E ¥ = 5 E ¥ = 10 E ¥ = 30 E ¥ = 20 Dt

 51415103020 Figure 5.14 -Stream function distribution for the unstable regime at E ∞ = 30. Here I corresponds to the Rubinstein-Zaltzman micro-vortices and II stands for the Dukhin-Mishchuk vortices. The arrows show the direction of the Rubinstein-Zaltzman micro-vortices propagation.

Figure 5 . 15 -

 515 Figure5.15 -The electric current normalized to its maximum value along the particle surface 0 < θ < 180 • for different strengths of the external electric field E ∞ . The inset shows the variation of the ion flux j/j max in time for the unstable case, which shows a travelling wave. In this inset, all subsequent curves are for t 0 + n∆t, n = 0, 1, 2, 3 (red arrow shows the direction of wave propagation).

Figure 5 . 16 -

 516 Figure 5.16 -Time records of the electrophoretic velocity U ∞ at the different strengths of the electric field E ∞ .

Figure 5 . 17 -

 517 Figure 5.17 -Diffusion boundary layer separation angle θ 0 vs χ at κ = 0.26. The solid line stands for the semi-analytical solution and the triangles stand for numerical results.

Figure 5 . 18 -Figure 5 . 19 -Figure 5 . 20 -

 518519520 Figure 5.18 -The normalized drop of the electrical potential F inside the region 0 < y < y m as a function of χ at θ = 180 • , F = F m , at κ = 0, 26. The solid line stands for the semi-analytical solution and the triangles stand for the numerical results. The inset shows F m along the particle surface at χ = 1.

Figure 5 . 21 - 1

 5211 Figure 5.21 -The electrophoretic velocity according to different theoretical approaches as a function of χ. Curves 1 and 2 are the predictions of our semi-analytical theory: 1 correspond to κ = 0.1 and 2 correspond to κ = 1. Curve 3 stands for the DNS with κ = 0.26 and ν = 0.002. The dashed line corresponds to Dukhin's relation (5.55).

Figure 6 . 1 -Figure 6 . 2 -

 6162 Figure 6.1 -Dependence of the particle velocity on the cations concentration p on its surface at ν = 0.0086, which correspond to the particle radius a = 5µm (a) for E ∞ = 0.05 and (b) for E ∞ = 0.8.

6 . 3 .Figure 6 . 3 -

 6363 Figure 6.3 -Slip velocity normalized to E ∞ along the particle surface (triangles correspond to the analytical solution).

Figure 6 . 4 -

 64 Figure 6.4 -The solid line corresponds to the results of direct numerical simulations, the triangular markers correspond to the experimental work [39]. The dimensionless parameter χ is (a) χ = 0.035, (b) χ = 0.088, and (c) χ = 0.151.Here, the angle is counted backwards to match the experiments.

Fig. 6 .Figure 6 . 5 -

 665 Fig. 6.5 shows the comparison of the velocity u ∞ with the experimental data. For better clarity, all experimental points are presented in dimensionless logarithmic coordinates, u ∞ vs E ∞ .

F = Q 2 Figure 6 . 6 -

 266 Figure 6.6 -Comparison of the theoretical and experimental electrophoretic velocity in the universal variables, U ∞ vs χ.The markers for the experimental points are given in Table6.1; the solid line stands for the results of our DNS at ν = 0.002, and the dashed line stands for the semi-analytic results. Inset A: behavior of u ∞ for small E ∞ , the solid line stands the DNS, and the dashed line stands for the analytic solution. Inset B: behavior of u ∞ for large E ∞ .

the

  Nernst-Planck-Poisson-Stokes system. The presented solution of this problem shows a qualitative and quantitative picture of the ion concentration and charge density fields, electric field, and hydrodynamic field in the vicinity of the microparticle. The loss of stability of space charge near the particle surface at sufficiently high electric field strengths and the subsequent transition to the chaotic motion has been revealed. It should be emphasized that the instability was obtained at zero Reynolds numbers. The analytical and numerical values of the electrophoretic velocity, electroosmotic velocity, and electrokinetic potential have been successfully compared with the experimental data available in the literature.

1. 1 4 / 3 ∞c 5 ( 1 ) = c 5 ( 1

 1435151 Scheme of cation and anion exchange membrane. . . . . . . . . . . . . . . 2.1 Scheme of EDL for (a) flat surface and (b) spherical surface. . . . . . . . . 2.2 Three types of particles in an electrolyte solution under the influence of an external electric field: (a) dielectric particle, (b) metal particle and (c) ion-selective particle. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Schematic of the flow near the granule under the external electric field, E ∞ and the fluid flow velocity U ∞ at infinity which are directed along thex-coordinate. The electrophoretic velocity of the particle has an opposite direction. In the spherical polar system, x = r cos θ and y = r sin θ, U = U θ is the tangential velocity at the particle surface and V = U r is the normal velocity at the particle surface. At the co-moving reference frame, the far-field velocity condition is U → -U ∞ sin θ and V → U ∞ cos θ. . . . . 3.1 Staggered grid for numerical scheme. The electrostatic potential is defined in nodes i -1, i and i + 1. The velocity component is defined in center nodes i -1/2 and i + 1/2. . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1 Structure of electrical double layer in weak electric field. . . . . . . . . . . 4.2 The geometry of the projections of normal and tangential vectors to the particle surface r = 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3 (a) Space charge distribution, ρ(r, θ) = c + -c -, for a weak electric field, E ∞ = 0.05 and p = 3. (b) The cross-section of ρ(r) for the angle θ = 0 (solid line) and θ = 180 • (dashed line); 1 is the solution obtained from first-order expansion. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4 (a) Electrical conductivity, K(r, θ) = c + + c -, for a weak electric field, E ∞ = 0.05 and p = 3. (b) The cross-section of K(r) for the angle θ = 0 (solid line) and θ = 180 • (dashed line); curve 1 is the solution obtained by first-order expansion. The behavior inside the EDL is shown in the insets.4.5 (a) Space charge distribution, ρ(r, θ), for moderate electric field strength, E ∞ = 3, p = 3. (b) The section ρ(r) for θ = 0 • (solid line) and θ = 180 • (dashed line). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.6 (a) Electrical conductivity, K(r, θ) = c + + c -, for moderate electric field strength E ∞ = 3, p = 3. (b) The cross-section K(r) for θ = 0 (solid line) and θ = 180 • (dashed line). The behavior inside the EDL is shown in the insets. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.7 Direction of electroosmotic slip velocity along the particle. . . . . . . . . . 4.8 Electrophoretic mobility vs the electric field strength E ∞ for different values of p. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.9 The dependence of the particle velocity on the cations concentration p on its surface for electric field strength E ∞ = 5. . . . . . . . . . . . . . . . . . 5.1 The structure of a non-equilibrium EDL near a cation-selective particle in a strong electric field. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2 I, II, and III correspond to thin boundary layers embedded in each other: EDL, space charge region, and diffusion layer, respectively. IV corresponds to the electrically neutral bulk of the electrolyte, where K = 2. The points 1 and 1 separate the areas of incoming and outgoing ion fluxes. Inset A: Typical structure of the incoming ion flux region. Inset B: Typical structure of the outgoing ion flux region. . . . . . . . . . . . . . . . . . . . . . . . . 5.3 Schematic diagrams of bifurcation. (a) The solution of the equation (5.11) has an imperfect pitchfork bifurcation, which is structurally unstable [68] (b) In the solution of the equation (5.10) for small but finite values of ν, this type of bifurcation is destroyed. The solid line corresponds to the solution that has physical meaning. . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.4 (a) The electric conductivity K(r, θ) for the moderate electric field strength E ∞ = 5. The dashed line stands for the outer edge of the diffusion layer. (b) Cross-section K(r) at fixed angles θ = 0 • and θ = 180 • . The behavior inside the EDL (the depleted and enriched regions correspondingly) is shown in the insets. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.5 (a) The electric conductivity K(r, θ) for E ∞ = 10. The dashed line stands for the outer edge of the diffusion layer. The values have been clipped from the actual maximum down to K = 5 in order to achieve contrast with the diffusion layer. (b) Cross-section of K(r) at fixed angles θ = 0 • and θ = 180 • . The behavior inside the EDL (the depleted and enriched regions correspondingly) is shown in the insets. . . . . . . . . . . . . . . . . . . . . 5.6 (a) The electric conductivity K(r, θ) for E ∞ = 20. The dashed line stands for the outer edge of the diffusion layer. The values have been clipped from the actual maximum down to K = 5 in order to achieve contrast with the diffusion layer. (b) Cross-section of K(r) at fixed angles θ = 0 • and θ = 180 • .The behavior inside the EDL (the depleted and enriched regions correspondingly) is shown in the insets. . . . . . . . . . . . . . . . . . . . . 5.7 (a) The charge density ρ(r, θ) for the moderate electric field strength E ∞ = 5. (b) Cross-section of ρ(r) at fixed angles θ = 0 • and θ = 180 • near the particle surface. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.8 (a) The charge density ρ(r, θ) for E ∞ = 10. (b) Cross-section of ρ(r) at fixed angles θ = 0 • and θ = 180 • near the particle surface. . . . . . . . . . 5.9 (a) The charge density ρ(r, θ) for E ∞ = 20. (b) Cross-section of ρ(r) at fixed angles θ = 0 • and θ = 180 • near the particle surface. . . . . . . . . . 5.10 Theoretical stream function distribution for (a) E ∞ = 5, (b) E ∞ = 10, (c) E ∞ = 20. (d) The picture of electrokinetic flow from [69] around a 1-mm cation selective particle at Ẽ∞ = 100 V/cm, which corresponds to dimensionless E ∞ = 200 and Debye number ν ≈ 8.68 • 10 -5 . . . . . . . . . 5.11 (a) The charge density ρ = c + -c -distribution and (b) electric conductivity K = c + + c -distribution at E ∞ = 30. I and II -regions of the standing and travelling waves, respectively. Arrows show the direction of propagation. 5.12 (a) Blow-up of region II of Fig. 5.11(a). (b) Spikes for a flat membrane [45, 72]. The normal to the membrane coordinate is strongly compressed in comparison to the direction along the membrane. . . . . . . . . . . . . . 5.13 (a) Cross-sections of ρ(r) at θ = 180 • at different time instants. (b) Crosssections of K(r) at θ = 180 • at different time instants. The wave is bouncing back from the outer edge of the diffusion layer, which thickness defined by δ. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.14 Stream function distribution for the unstable regime at E ∞ = 30. Here I corresponds to the Rubinstein-Zaltzman micro-vortices and II stands for the Dukhin-Mishchuk vortices. The arrows show the direction of the Rubinstein-Zaltzman micro-vortices propagation. . . . . . . . . . . . . . . 5.15 The electric current normalized to its maximum value along the particle surface 0 < θ < 180 • for different strengths of the external electric field E ∞ . The inset shows the variation of the ion flux j/j max in time for the unstable case, which shows a travelling wave. In this inset, all subsequent curves are for t 0 + n∆t, n = 0, 1, 2, 3 (red arrow shows the direction of wave propagation). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.16 Time records of the electrophoretic velocity U ∞ at the different strengths of the electric field E ∞ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.17 Diffusion boundary layer separation angle θ 0 vs χ at κ = 0.26. The solid line stands for the semi-analytical solution and the triangles stand for numerical results. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.18 The normalized drop of the electrical potential F inside the region 0 < y < y m as a function of χ at θ = 180 • , F = F m , at κ = 0, 26. The solid line stands for the semi-analytical solution and the triangles stand for the numerical results. The inset shows F m along the particle surface at χ = 1. 5.19 The thickness Y m of the SCR at different points of the particle surface. The solid line stands for the semi-analytical solution and the triangles joined by the dashed line stand for the numerical results. . . . . . . . . . . . . . . 5.20 The ion flux J on the surface versus the angle θ for different values of χ. The dashed line stands for the DNS results and the solid line stands for the analytical solution. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.21 The electrophoretic velocity according to different theoretical approaches as a function of χ. Curves 1 and 2 are the predictions of our semi-analytical theory: 1 correspond to κ = 0.1 and 2 correspond to κ = 1. Curve 3 stands for the DNS with κ = 0.26 and ν = 0.002. The dashed line corresponds to Dukhin's relation (5.55). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.1 Dependence of the particle velocity on the cations concentration p on its surface at ν = 0.0086, which correspond to the particle radius a = 5µm (a) for E ∞ = 0.05 and (b) for E ∞ = 0.8. . . . . . . . . . . . . . . . . . . . . . 6.2 Comparison of the theoretical particle velocity (a) and its mobility (b) with experiments. The solid line corresponds to the numerical solution, the dashed line to the analytical solution (4.98). Experimental data from [40] are labeled with round markers; triangular markers stand for data from [41]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.3 Slip velocity normalized to E ∞ along the particle surface (triangles correspond to the analytical solution). . . . . . . . . . . . . . . . . . . . . . . . 6.4 The solid line corresponds to the results of direct numerical simulations, the triangular markers correspond to the experimental work [39]. The dimensionless parameter χ is (a) χ = 0.035, (b) χ = 0.088, and (c) χ = 0.151. Here, the angle is counted backwards to match the experiments. . . 6.5 The experimental dimensionless velocity u ∞ as a function of E ∞ . Regions I, II and III stands for the regions of the electrophoresis of the first kind, transition region and the region of the electrophoresis of the second kind respectively. Lines 1, 2 and 3 stands for u∞ ∼ const 1 E ∞ , u ∞ ∼ const 2 E and u ∞ ∼ const 2 E 2 ∞ ,respectively. . . . . . . . . . . . . . . . . . . . . . . . 6.6 Comparison of the theoretical and experimental electrophoretic velocity in the universal variables, U ∞ vs χ. The markers for the experimental points are given in Table 6.1; the solid line stands for the results of our DNS at ν = 0.002, and the dashed line stands for the semi-analytic results. Inset A: behavior of u ∞ for small E ∞ , the solid line stands the DNS, and the dashed line stands for the analytic solution. Inset B: behavior of u ∞ for large E ∞ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . i = 2 : a 5 (2) ψ(0, j) + b 5 (2) ψ(1, j) + c 5 (2) ψ(2, j) + d 5 (2) ψ(3, j)+ +e 5 (2) ψ(4, j) = R(2, j) From the boundary condition at r = 1 the stream function ψ(0, j) = 0. We can write the following relation in order to express unknown ψ(-1, j), We need to correct the coefficient c 5 (1) write the formulas for the coefficients a 5 , b 5 , c 5 , d 5 , e 5 . The first derivative for the function is defined in nodes.

2 U ∞ r 2 2 U ∞ rn 2 (

 2222 = U ∞ sin θ, V = -U ∞ cos θ ∂Ψ ∂r r=Rmax = -U r sin θ = -U ∞ r sin 2 θ Ψ| r=Rmax = ∂Ψ ∂r dr = -U ∞ sin 2 θ rdr = -1 sin 2 θ + CWe assume the integration constant equal to zero. From the boundary conditions atr = R max one can get, ψ(Im, j) = -1 Im) stn 2 (j) ψ (Im + 1, j) = -U ∞ rn(Im) stn 2 (j) For i = Im -2 : a 5 (i) ψ(Im -4, j)+b 5 (i) ψ(Im -3, j)+c 5 (i) ψ(Im -2, j)+d 5 (i) ψ(Im -1, j)+ e 5 (i) ψ(Im, j) = R(Im -2, j)We need to correct the right part of equation in the following way, R(Im -2, j) → R(Im -2, j)e 5 (i) ψ(Im, j).For i = Im -1 : a 5 (i) ψ(Im -3, j) + b 5 (i) ψ(Im -2, j) + c 5 (i) ψ(Im -1, j) + d 5 (i) ψ(Im, j) + e 5 (i) ψ(Im + 1, j) = R(Im -1, j)We can write the following relation in order to express ψ(Im + 1, j), , j) -ψ(Im -1, j) rm1(Im)+ ψ(Im + 1, j) -ψ(Im, j) rm1(Im + 1) ψ(Im + 1, j) = rm1(Im + 1) rm1(Im) ψ(Im -1, j) + 1 -rm1(Im + 1) rm1(Im) ψ(Im, j)+ +2 rm1(Im + 1) ∂Ψ ∂r r=Rmax = rm1(Im + 1) rm1(Im) ψ(Im -1, j)+ + 1 -rm1(Im + 1) rm1(Im) ψ(Im, j) + 2 rm1(Im + 1) ψ (Im + 1, j)We need to correct the coefficient c5(Im -1) and the right part of the equation,i = Im -1 : -1, j) to R(Im -1, j) -d5(i) ψ(Im, j)--e5(i) 1 -rm1(Im + 1) rm1(Im) ψ(Im, j) + 2 rm1(Im + 1) ψ (Im + 1, j) -For the final five-diagonal matrix a(i) and e(i) remain unchanged. b 5 (i) → b 5 (i) + λ j • b 3 (i) c 5 (i) → c 5 (i) + λ j • c 3 (i) + λ 2 j d 5 (i) → d 5 (i) + λ j • d 3 (i)

Table 2 .

 2 1 gives dimensional physical quantities that are closest to the different experimental values.

Table 6 .

 6 respectively.

	Marker	Radius, µm Papers	E ∞
	+	600	Baran et al. [78], Barany et al. [71]	2.5 -8 kV/m
	•	375	Baran et al. [78], Barany et al. [71, 41],	2.5 -30 kV/m
			Mishchuk et al. [79]
		250	Baran et al. [78], Barany et al. [71, 41],	0.5 -100 kV/m
			Mishchuk et al. [79]
	×	112.5	Baran et al. [78]	10 -90 kV/m
		100	Baran et al. [78], Barany et al. [71, 41],	0.5 -100 kV/m
			Mishchuk et al. [79]
		62.5	Baran et al. [78]	15 -90 kV/m
		50	Baran et al. [78], Barany et al. [71, 41],	0.5 -100 kV/m
			Mishchuk et al. [40, 79]
		25	Baran et al. [78], Barany et al. [71, 41],	2.5 -90 kV/m
			Mishchuk et al. [40, 79]
		5	Baran et al. [78], Barany et al. [71, 41],	1 -10 kV/m
			Mishchuk et al. [40, 79]
	♦	2.5	Baran et al. [78], Barany et al. [71, 41],	1 -10 kV/m
			Mishchuk et al. [40]
		0.5	Baran et al. [78], Barany et al. [71, 41],	0.2 -10 kV/m
			Mishchuk et al. [40]

1 -Experimental data from different papers (indicated). The region I of electrophoresis of the first kind is clearly seen since the velocity here is linearly proportional to E ∞ . After the transition region II, the region III of strongly non-linear electrophoresis, predicted by Dukhin, appears. According to Dukhin's prediction (see Eq. (5.54)) u

Appendix A. Calculation of matrix operator for Poisson equation

Let us show how to determine these coefficients. Function Φ is defined in center nodes in r and θ directions. The first derivative with respect to r is defined as,

multiplying by r 2 we get,

and finally, the second derivative takes the form,

Let us group the coefficients for each Φ,

For the Poisson equation, the matrix of the operator with regards to the variable r has three-diagonal form. The system of equations is the following for fixed index j:

... ... ... ... ... ... ... ...

where Φ 0 and Φ Im+1 have to be defined from the boundary conditions.

Function Φ is defined in the center nodes, but for the boundary conditions, it must be defined in nodes. The cosine and sine functions at nodes and center nodes from angle θ are denoted as follows, stn(j) = sin(tn(j))

We need to correct the coefficients b 1 and b Im according to the boundary conditions.

where R is the right hand side of Poisson equation. From the boundary condition,

That means that we need to reassign the coefficient b(1) in the following way b(1)

In the last layer of nodes at i = Im we have,

After substitution of Φ(Im + 1, j) from Eq. (A.10) in the Eq. (A.14) one can get,

Then coefficient at Φ(Im, j) needs to be changed as well as the

Let us discretize the second term of Poisson equation.

• stn(j) Φ(i, j + 1) -Φ(i, j) tn1(j) stn(j -1) Φ(i, j) -Φ(i, j -1) tn1(j -1) (A.18)

Φ(i, j) : -stn(j -1) stm(j) tm1(j) tn1(j -1) -stn(j) stm(j) tm1(j) tn1(j) = e j (A.20)

Consider the matrix operator in the equation (3.14). The coefficients a i and c i remain unchanged, and b i must be corrected as follows: b i = b i + λ in the outer loop over j and the inner loop over i.

In this case, the coefficients d(j), e(j), f (j) are not involved in the process of solving equation, they needed only for defining the operator T θ and its eigenvectors and eigenvalues. Eventually, we end up with the tridiagonal matrix, which is solved by the tridiagonal matrix algorithm for each column of the right-hand matrix.

Appendix B. Calculation of matrix operator for Stokes equation

Consider the numerical scheme for Eq. (3.23). The four-order term r 4 ∂ 4 ∂r 4 in this equation can be approximated by a five-point stencil,

The second-order term r 4 ∂ 2 ∂r 2

1

Ψ can be approximated by a threepoint stencil,

Function Ψ defined in nodes for radial direction and the boundary conditions must be defined in nodes.