Mathematical Modeling of the Influence of Electrical Heterogeneity on the Processes of Salt Ion Transfer in Membrane Systems with Axial Symmetry Taking into Account Electroconvection

1. Introduction

Operation of electromembrane systems (EMS) in overlimit current modes has a significant potential for a significant reduction in investment costs for desalination plants. The study by Wessling M., Stockmeier F. et al. [1] is devoted to the analysis of electroconvective vortices using the surface patterning method or geometric and electrical heterogeneity. In [1], a comparison of electroconvection development on two membranes modified using patterns was performed. An analysis was carried out of the increase in the area of electroconvective vortices, their structural stability in a steady state, as well as the possibility of controlling electroconvective vortices by introducing certain structural heterogeneities of the membrane surface to reduce energy losses in the current-voltage characteristic (CVC) and more economical operation in EMS in overlimit modes, significantly increasing their efficiency.

As the analysis of works [2,3,4,5] shows, the phenomenon of overlimit mass transfer is currently being intensively studied, however, many issues are still poorly understood. For example, it is necessary to study what effect modifications of the membrane surface have on the development of electroconvection. Electrical inhomogeneity arises due to the inhomogeneous distribution of charges on the membrane surface, which can be caused by both the geometric features of the system [6,7] and the properties of the membrane material. In work [8], it is shown that such inhomogeneities lead to local zones with increased electric field strength, which, in turn, affects the rate of ion migration, and, accordingly, the emergence and development of electroconvective vortices. This, in turn, leads to an increase in the rate of ion transfer through the membrane. In works [9,10,11,12], the importance of taking into account electrical inhomogeneity in modeling ion transport is also emphasized.

In recent years, there has been significant progress in understanding the role of electrical heterogeneity and electroconvection in EMS. A new model was proposed in [13] that takes into account nonlinear effects in electroconvection, which allowed for a more accurate description of mass transfer processes under high current density conditions.

Experimental studies conducted in [14] confirmed that taking into account electrical heterogeneity allows for a more accurate description of mass transfer processes in membrane systems. Theoretical models developed on the basis of these data take into account both ion migration and convective flows, which makes them more realistic.

Thus, the study of electrical heterogeneity and its effect on ion transport with electroconvection taken into account is an important area of research into increasing the efficiency of separation and purification of solutions and is of great practical importance. Further theoretical research in this area should be aimed at developing more accurate models. The advantage of constructing 3D models with axial symmetry is the property of equal accessibility at sublimiting currents, i.e. the fact that the thickness of the diffusion layer in such systems is constant up until the moment of the onset of electroconvection. At the same time, we note the insufficient development of such mathematical models and numerical methods for solving boundary value problems for overlimit current modes.

This article is a continuation of the work [15], where the basic 3D model of electroconvection with axial symmetry was formulated and investigated. In this article, a modification of the basic model in the form of an annular membrane disk is considered, associated with the addition of non-conductivity areas on the cation-exchange membrane (CEM) and their influence on the transfer process is investigated.

2. Methods

2.1. Mathematical Model of Transport for Annular Membrane Disk

Let us consider the transfer of salt ions during rotation of a cation-exchange membrane disk, with non-conducting and conducting regions, inside a vertically standing cylindrical cell around the central axis, taking into account electroconvection, and formulate the corresponding mathematical model.

2.1.1. Solution Area

When formulating a mathematical model, as well as the subsequent numerical solution due to axial symmetry (Figure 1a), it is sufficient to describe half of the cross-section of the cylindrical region, where it is necessary to determine the equations and boundary conditions (Figure 1b). To do this, it is necessary to move from a rectangular coordinate system $(x,y,z)$ to a cylindrical one $(r,\phi ,z)$. However, the velocities in the angular direction dffer from zero, so the model must include all three components of velocity (radial, azimuthal and axial), flows, electric field strength, although they depend only on two arguments $(r,z)$, i.e. we have a 2D model. Moreover, with axial symmetry, the derivatives with respect to the component $\phi$ are zero, so all equations in the cylindrical coordinate system can be simplified in this case.

2.1.2. System of Equations

To model the transfer processes, a coupled system of Navier – Stokes equations with volumetric electric force $\overrightarrow{f}$ and Nernst – Planck – Poisson equations is used, which in a cylindrical coordinate system $(r,\phi ,z)$ takes the form [15]:

$$\rho {\displaystyle \frac{\partial u}{\partial t}}+\rho \left(u{\displaystyle \frac{\partial u}{\partial r}}-{\displaystyle \frac{v^2}{r}}+w{\displaystyle \frac{\partial u}{\partial z}}\right)+{\displaystyle \frac{\partial p}{\partial r}}=\eta \left[{\displaystyle \frac{1}{r}}{\displaystyle \frac{\partial }{\partial r}}\left(r{\displaystyle \frac{\partial u}{\partial r}}\right)-{\displaystyle \frac{u}{r^2}}+{\displaystyle \frac{{\partial }^{2}u}{\partial z^2}}\right]+f_r$$

$$\rho {\displaystyle \frac{\partial v}{\partial t}}+\rho \left(u{\displaystyle \frac{\partial v}{\partial r}}-{\displaystyle \frac{uv}{r}}+w{\displaystyle \frac{\partial v}{\partial z}}\right)=\eta \left[{\displaystyle \frac{1}{r}}{\displaystyle \frac{\partial }{\partial r}}\left(r{\displaystyle \frac{\partial v}{\partial r}}\right)-{\displaystyle \frac{v}{r^2}}+{\displaystyle \frac{{\partial }^{2}v}{\partial z^2}}\right]+f_{\phi }$$

$$\rho {\displaystyle \frac{\partial w}{\partial t}}+\rho \left(u{\displaystyle \frac{\partial w}{\partial r}}+w{\displaystyle \frac{\partial w}{\partial z}}\right)+{\displaystyle \frac{\partial p}{\partial z}}=\eta \left[{\displaystyle \frac{1}{r}}{\displaystyle \frac{\partial }{\partial r}}\left(r{\displaystyle \frac{\partial w}{\partial r}}\right)+{\displaystyle \frac{{\partial }^{2}w}{\partial z^2}}\right]+f_z$$

$${\overrightarrow{j}}_i=({\displaystyle \frac{F}{RT}}{z}_i{D}_i{C}_i{E}_r-D_i{\displaystyle \frac{\partial C_i}{\partial r}}+C_{i}u){\overrightarrow{e}}_r+C_{i}v{\overrightarrow{e}}_{\phi }+({\displaystyle \frac{F}{RT}}{z}_i{D}_i{C}_i{E}_z-D_i{\displaystyle \frac{\partial C_i}{\partial z}}+C_{i}w){\overrightarrow{e}}_z$$

(1)

$${\displaystyle \frac{\partial C_i}{\partial t}}=-{\displaystyle \frac{1}{r}}{\displaystyle \frac{\partial }{\partial r}}(r{j}_{i,r})-{\displaystyle \frac{\partial j_{i,z}}{\partial z}}, i=1,2$$

(2)

$${\displaystyle \frac{{\partial }^2\mathsf{\Phi }}{\partial r^2}}+{\displaystyle \frac{{\partial }^2\mathsf{\Phi }}{\partial z^2}}=-{\displaystyle \frac{1}{r}}{\displaystyle \frac{\partial \mathsf{\Phi }}{\partial r}}-{\displaystyle \frac{F}{\epsilon }}(z_1{C}_1+z_2{C}_2)$$

(3)

$$\overrightarrow{I}=F(z_1{\overrightarrow{j}}_1+z_2{\overrightarrow{j}}_2)$$

(4)

Where

$$f_r=\epsilon \left({\displaystyle \frac{{\partial }^2\mathsf{\Phi }}{\partial r^2}}+{\displaystyle \frac{{\partial }^2\mathsf{\Phi }}{\partial z^2}}\right){\displaystyle \frac{\partial \mathsf{\Phi }}{\partial r}}+{\displaystyle \frac{\epsilon }{r}}{\left({\displaystyle \frac{\partial \mathsf{\Phi }}{\partial r}}\right)}^2$$

$$f_{\phi }=0$$

$$f_z=\epsilon \left({\displaystyle \frac{{\partial }^2\mathsf{\Phi }}{\partial r^2}}+{\displaystyle \frac{{\partial }^2\mathsf{\Phi }}{\partial z^2}}\right){\displaystyle \frac{\partial \mathsf{\Phi }}{\partial z}}+{\displaystyle \frac{\epsilon }{r}}{\displaystyle \frac{\partial \mathsf{\Phi }}{\partial r}}{\displaystyle \frac{\partial \mathsf{\Phi }}{\partial z}}$$

$\overrightarrow{u}=(u,v,w)$, where $u,v,w$ – radial, azimuthal and axial components of the solution flow velocity $\overrightarrow{u}$; index 1 refers to cations, index 2 to anions; ${\overrightarrow{j}}_1,{\overrightarrow{j}}_2$ – flows; $C_1,C_2$ - concentration in solution; – charge numbers; $\overrightarrow{I}$ – current density; $D_1, D_2$ – diffusion coefficients; $\mathsf{\Phi }$ – electric field potential; $\overrightarrow{E}=-\nabla \mathsf{\Phi }$ – electric field strength; $\epsilon$ – dielectric permeability of the electrolyte; $F$ – Faraday constant; $R$ – universal gas constant; $T$ – absolute temperature; $t$ – time; $\rho$– density; $\eta$ – dynamic viscosity and $p$ – pressure.

2.1.3. Boundary Conditions

In the studied section of the cylindrical region (Figure 1), boundary 2 models the part of space infinitely distant from CEM, in which the condition of electroneutrality is satisfied and the concentration of the solution is constant ($C_0$). Boundary 2 is also considered an anode and an open boundary (inlet) for the solution, for the velocity here the condition of the absence of normal voltage is set, and the pressure is considered equal to zero. This boundary is considered an equipotential surface, and $\mathsf{\Phi }(t,r,0)=d_{\phi }$. Boundary 1 is the axis of symmetry. Boundary 5 is considered the outlet for the solution, therefore the condition of removal of ions only by convective flow is set ${\overrightarrow{j}}_i=-\overrightarrow{u}\cdot C_i, i=1,2$, and for the velocity the condition of absence of normal stress. In addition, for the potential the boundary condition is used $-\overrightarrow{n}\cdot {(r{\displaystyle \frac{\partial \mathsf{\Phi }}{\partial r}},{\displaystyle \frac{\partial \mathsf{\Phi }}{\partial z}})}^T=0$, here is $\overrightarrow{n}$ the normal vector. At boundary 4 the CEM surface is considered equipotential: $\mathsf{\Phi }(t,r,H)=d_0=0$. Boundaries 3, 4 correspond to the rotating CEM, for the azimuthal velocity the condition is used here: $v(t,r,H)=\omega r$, where $\omega$ is the angular velocity of rotation of the disk. For anions the condition of impermeability (absence of flow) is set here: $-\overrightarrow{n}\cdot {\overrightarrow{j}}_2=0$. The conductivity region is located at boundary 4, so the cation exchange membrane forms a ring (Figure 1b). At boundary 3, corresponding to the non-conducting region, there is a condition of impermeability (no flow): $-\overrightarrow{n}\cdot {\overrightarrow{j}}_i=0$, $i=1,2$. To derive the condition for the potential at the non-conducting boundary 3, the equality of the current density $\overrightarrow{I}$ to zero at this boundary is used. The normal component of the current density at this boundary is equal to zero due to the conditions on the flows. We determine the radial component of the potential gradient from the condition that the radial component of the current density is equal to zero from the formula for the current (4):

$$\begin{array}{l}\overrightarrow{I}=F(-{\displaystyle \frac{F}{RT}}(z_1^2{D}_1{C}_1+z_2^2{D}_2{C}_2){\displaystyle \frac{\partial \mathsf{\Phi }}{\partial r}}-(z_1{D}_1{\displaystyle \frac{\partial C_1}{\partial r}}+z_2{D}_2{\displaystyle \frac{\partial C_2}{\partial r}})+(z_1{C}_1+z_2{C}_2)u){\overrightarrow{e}}_r+\\ +((z_1{C}_1+z_2{C}_2)v){\overrightarrow{e}}_{\phi }+\\ +F(-{\displaystyle \frac{F}{RT}}(z_1^2{D}_1{C}_1+z_2^2{D}_2{C}_2){\displaystyle \frac{\partial \mathsf{\Phi }}{\partial z}}-(z_1{D}_1{\displaystyle \frac{\partial C_1}{\partial z}}+z_2{D}_2{\displaystyle \frac{\partial C_2}{\partial z}})+(z_1{C}_1+z_2{C}_2)w){\overrightarrow{e}}_z\end{array}$$

Therefore, from $I_r=0$, we get:

$$F(-{\displaystyle \frac{F}{RT}}(z_1^2{D}_1{C}_1+z_2^2{D}_2{C}_2){\displaystyle \frac{\partial \mathsf{\Phi }}{\partial r}}-(z_1{D}_1{\displaystyle \frac{\partial C_1}{\partial r}}+z_2{D}_2{\displaystyle \frac{\partial C_2}{\partial r}})+(z_1{C}_1+z_2{C}_2)u)=0$$

Since the radial component of the velocity is $u=0$, we obtain:

$$-{\displaystyle \frac{F}{RT}}(z_1^2{D}_1{C}_1+z_2^2{D}_2{C}_2){\displaystyle \frac{\partial \mathsf{\Phi }}{\partial r}}-(z_1{D}_1{\displaystyle \frac{\partial C_1}{\partial r}}+z_2{D}_2{\displaystyle \frac{\partial C_2}{\partial r}})=0$$

From here: ${\displaystyle \frac{\partial \mathsf{\Phi }}{\partial r}}=-{\displaystyle \frac{RT(z_1{D}_1\frac{\partial C_1}{\partial r}+z_2{D}_2\frac{\partial C_2}{\partial r})}{F(z_1^2{D}_1{C}_1+z_2^2{D}_2{C}_2)}}$.

Thus, at boundary 3, the following condition must be set for the potential:

$${\displaystyle \frac{\partial \mathsf{\Phi }}{\partial r}}=-{\displaystyle \frac{RT(z_1{D}_1\frac{\partial C_1}{\partial r}+z_2{D}_2\frac{\partial C_2}{\partial r})}{F(z_1^2{D}_1{C}_1+z_2^2{D}_2{C}_2)}}$$

It is assumed that the cell is initially completely filled with an ideally mixed solution of sodium chloride with a concentration $C_0$ and the same ideally mixed solution is fed into it through the boundary 2. Therefore, a constant concentration is taken as the initial condition $C_1(0,r,z)=C_2(0,r,z)=C_0$.

2.2. Algorithm for Hybrid Numerical-Analytical Solution

2.2.1. Problems of Numerical Solution of Boundary Value Problem

In order to clarify the problems of numerical solution of the boundary value problem, we will move to a dimensionless form using characteristic quantities.

The problem has two characteristic quantities $H$- the cell height and $r_0$ is the outer radius of the annular membrane. Since in practice $H=a\cdot r_0$, where proportionality coefficient $a=1÷5$, therefore they have the same order, and is taken $r_0$ as a characteristic quantity, since it is convenient to express the characteristic speed through it. As the characteristic velocity $V_0$, we will use the azimuthal component of the velocity V₀ = ω₀r₀.

Let us assume that:

$$\begin{array}{l}{r}^{(u)}={\displaystyle \frac{r}{r_0}}; z^{(u)}={\displaystyle \frac{z}{r_0}}; t^{(u)}={\displaystyle \frac{t{V}_0}{r_0}}; {\overrightarrow{u}}^{(u)}={\displaystyle \frac{\overrightarrow{u}}{V_0}}; {С}_i^{(u)}={\displaystyle \frac{C_i}{C_0}};{\overrightarrow{j}}_i^{(u)}={\displaystyle \frac{{\overrightarrow{j}}_i}{j_0}};{\displaystyle \frac{D_0{C}_0}{r_0}}{\overrightarrow{j}}_i^{(u)}={\overrightarrow{j}}_i;\\ D_i^{(u)}={\displaystyle \frac{D_i}{D_0}};{\overrightarrow{E}}^{(u)}={\displaystyle \frac{r_{0}F}{RT}}\overrightarrow{E};Pe={\displaystyle \frac{r_0{V}_0}{D_0}};{\epsilon }^{(u)}={\displaystyle \frac{\epsilon R{T}_0}{r_0^2{C}_0{F}^2}};{\mathsf{\Phi }}^{(u)}={\displaystyle \frac{F}{R{T}_0}}\mathsf{\Phi };p_0=\rho V_0^2;\\ {\eta }^{(u)}={\displaystyle \frac{\eta }{D_0{C}_{0}F}};K_{el}={\displaystyle \frac{RT{C}_0}{\rho V_0^2}};\mathrm{Re}={\displaystyle \frac{r_0{V}_0}{\nu }}.\end{array}$$

Next, we will make the transition in the equations from dimensional quantities to dimensionless ones. Then we will get:

$$\begin{array}{l}{\overrightarrow{j}}_i^{(u)}=(z_i{D}_i^{(u)}{C}_i^{(u)}{E}_r^{(u)}-D_i^{(u)}\frac{\partial C_i^{(u)}}{\partial r^{(u)}}+Pe{C}_i^{(u)}{u}^{(u)}\left)\right){\overrightarrow{e}}_r+\\ +Pe{C}_i^{(u)}{v}^{(u)}{\overrightarrow{e}}_{\phi }+(z_i{D}_i^{(u)}{C}_i^{(u)}{E}_z^{(u)}-D_i^{(u)}\frac{\partial C_i^{(u)}}{\partial z^{(u)}}+Pe{C}_i^{(u)}{w}^{(u)}){\overrightarrow{e}}_z+\\ Pe\frac{\partial C_i^{(u)}}{\partial t^{(u)}}=-\frac{1}{r^{(u)}}\frac{\partial }{\partial r^{(u)}}\left(r^{(u)}{j}_{i,r}^{(u)}\right)-\frac{\partial j_{i,z}^{(u)}}{\partial z^{(u)}},i=1,2\\ {\epsilon }^{(u)}(\frac{1}{r^{(u)}}\frac{\partial }{\partial r^{(u)}}(r^{(u)}{E}_r^{(u)})+\frac{\partial }{\partial z^{(u)}}{E}_z^{(u)})=z_1{C}_1^{(u)}+z_2{C}_3^{(u)}\\ {\overrightarrow{I}}^{(u)}=z_1{\overrightarrow{j}}_1^{(u)}{z}_2{\overrightarrow{j}}_2^{(u)}\\ \frac{\partial u^{(u)}}{\partial t^{(u)}}+(u^{(u)}\frac{\partial u^{(u)}}{\partial r^{(u)}}\frac{{(v^{(u)})}^2}{r^{(u)}}+w^{(u)}\frac{\partial u^{(u)}}{\partial z^{(u)}})+\frac{\partial p^{(u)}}{\partial r^{(u)}}=\\ =\frac{1}{\mathrm{Re}}\left[\frac{1}{r^{(u)}}\frac{\partial }{\partial r^{(u)}}\right(r^{(u)}\frac{\partial u^{(u)}}{\partial r^{(u)}})-\frac{u^{(u)}}{{(r^{(u)})}^2}+\frac{{\partial }^2{u}^{(u)}}{\partial z^{(u)2}}]+K_{\mathit{ek}}{f}_r^{(u)}\\ \frac{\partial v^{(u)}}{\partial t^{(u)}}+(u^{(u)}\frac{\partial v^{(u)}}{\partial r^{(u)}}-\frac{u^{(u)}{v}^{(u)}}{\partial r^{(u)}}+w^{(u)}\frac{\partial v^{(u)}}{\partial z^{(u)}})=\frac{1}{\mathrm{Re}}\left[\frac{1}{r^{(u)}}\frac{\partial }{\partial r^{(u)}}\right(r^{(u)}\frac{\partial v^{(u)}}{\partial r^{(u)}})-\frac{v^{(u)}}{{(r^{(u)})}^2}+\frac{{\partial }^2{v}^{(u)}}{\partial z^{(u)2}}]\\ \frac{\partial w^{(u)}}{\partial t^{(u)}}+(u^{(u)}\frac{\partial w^{(u)}}{\partial r^{(u)}}+w^{(u)}\frac{\partial w^{(u)}}{\partial z^{(u)}})+\frac{\partial p^{(u)}}{\partial z^{(u)}}=\\ =\frac{1}{\mathrm{Re}}\left[\frac{1}{r^{(u)}}\frac{\partial }{\partial r^{(u)}}\right(r^{(u)}\frac{\partial w^{(u)}}{\partial r^{(u)}})+\frac{{\partial }^2{w}^{(u)}}{\partial z^{(u)2}}]+K_{\mathit{ek}}{f}_r^{(u)}\\ \frac{1}{r^{(u)}}\frac{\partial }{\partial r^{(u)}}(r^{(u)}{u}^{(u)})+\frac{\partial w^{(u)}}{\partial z^{(u)}}=0\end{array}$$

where

$$K_{ek}=K_{el}{\epsilon }^{(u)} \mathrm{and} {\overrightarrow{f}}^{(u)}=E^{(u)}({\displaystyle \frac{1}{r^{(u)}}}{\displaystyle \frac{\partial }{\partial r^{(u)}}}(r^{(u)}{E^{(u)}}_r)+{\displaystyle \frac{\partial {E_z}^{(u)}}{\partial z^{(u)}}})$$

Let us evaluate and clarify the meaning of dimensionless parameters. The system of Nernst – Planck – Poisson and Navier – Stokes equations in a cylindrical coordinate system obtained in the process of non-dimensionalization includes 4 dimensionless parameters:

1. The Peclet number shows the ratio of the coefficients of kinematic viscosity and diffusion. With a characteristic channel height and an outer radius of the membrane disk of 1 mm, we obtain an estimate of the Peclet number (Table 1):

$$Pe={\displaystyle \frac{r_0{V}_0}{D_0}}={\displaystyle \frac{{\omega }_0{r}_0^2}{D_0}}={\displaystyle \frac{{\omega }_0[rad/s]\cdot {10}^{-6}[m^2]}{1.33\cdot {10}^{-9}[m/s^2]}}=0.75\cdot {10}^3\cdot {\omega }_0$$

The Peclet number can be considered a large parameter, that is, in this problem there will be a developed diffusion layer, and the convective transfer of salt ions prevails over the diffusion transfer.

2. Dimensionless number ${\epsilon }^{(u)}$:

$${\epsilon }^{(u)}={\displaystyle \frac{\epsilon RT}{r_0^2{C}_0{F}^2}}={\displaystyle \frac{8.854\cdot {10}^{-12}\cdot 8.314\cdot 293}{{10}^{-6}\cdot {96485}^2}}\cdot \left[{\displaystyle \frac{1}{C_0}}\right]\approx 2.3\cdot {10}^{-12}\cdot \left[{\displaystyle \frac{1}{C_0}}\right]$$

From this formula it follows that when changing $C_0$ from 0.01 mol/m³ to 100 mol/m³ the parameter ${\epsilon }^{(u)}$ changes from $2,3\cdot {10}^{-10}$ to $2,3\cdot {10}^{-14}$. Therefore, ${\epsilon }^{(u)}$ it can be considered a small parameter, which has a meaning in the form of the ratio of the square of the thickness of the region of equilibrium space charge to the square of the intermembrane distance [16]: ${\epsilon }^{(u)}={\displaystyle \frac{RT\epsilon }{r_0^2{C}_0{F}^2}}=2{\left[{\displaystyle \frac{l_d}{r_0}}\right]}^2$, where $l_d=\sqrt{{\displaystyle \frac{RT\epsilon }{2C_0{F}^2}}}$ – Debye length. Thus, the space charge region (SCR) located as part of the double electric boundary layer near the CEM is of the order of $\sqrt{{\epsilon }^{(u)}}$. It is easy to show that in this region the gradients of concentrations and electric field strengths are of the order of ${\displaystyle \frac{1}{\sqrt{{\epsilon }^{(u)}}}}$, as a result of which the boundary value problem is related to stiff problems and therefore serious problems arise in the numerical solution.

3. Reynolds number Re, which is the ratio of the inertial force $F_{in}={\rho }_0{H}^2{V}_0^2$ to the viscous friction force $F_{tr}=\nu {\rho }_0{V}_{0}H$.

Table 2 presents the results of estimating the Reynolds number for the desalination chamber. If we take the outer radius of the membrane disk $r_0$ of about 1 mm, and consider a liquid with the same viscosity as water $\nu =1.006\cdot {10}^{-6}$ m²/s, then $\mathrm{Re}={\displaystyle \frac{V_0{r}_0}{\nu }}={\displaystyle \frac{{\omega }_0{r}_0^2}{\nu }}={\displaystyle \frac{{\omega }_0\cdot {10}^{-6}[m^2]}{1.006\cdot {10}^{-6}\left[\raisebox{1ex}{$m^2$}\!\left/ \!\raisebox{-1ex}{$s$}\right.\right]}}=0.994{\omega }_0$.

From this table it is clear that the forces of inertia and viscous friction are approximately the same, therefore it is necessary to use the Navier-Stokes equation without any simplifications, in addition, from the adhesion conditions it follows that a hydrodynamic layer (friction layer) arises, but from the value of the Reynolds number it follows that it is not developed.

4. The dimensionless parameter $K_{el}$ is the ratio of electrical force to inertial force. Let's calculate it for a membrane radius of $r_0={10}^{-3}$ m (Table 3):

$$K_{el}={\displaystyle \frac{RT{C}_0}{\rho r_0^2{\omega }_0^2}}={\displaystyle \frac{8.314\cdot 293}{1.002\cdot {10}^3\cdot {10}^{-6}}}\cdot {\displaystyle \frac{C_0}{{\omega }_0^2}}={\displaystyle \frac{2.436\cdot {10}^3}{1.002\cdot {10}^{-3}}}\cdot {\displaystyle \frac{C_0}{{\omega }_0^2}}=2.43\cdot {10}^6\cdot {\displaystyle \frac{C_0}{{\omega }_0^2}}$$

From Table 3 and the formula it is evident that with increasing speed ${\omega }_0$ the value $K_{el}$ quickly decreases, and with increasing concentration it increases; in the Navier – Stokes equation this value occurs in the complex $K_{ek}=K_{el}\cdot {\epsilon }^{(u)}={\displaystyle \frac{RT{C}_0}{\rho r_0^2{\omega }_0^2}}{\epsilon }^{(u)}$, hence it is evident that electroconvection practically does not depend on the initial concentration, but strongly depends on the angular velocity of rotation, namely: with increasing angular velocity electroconvection quickly dies out.

The presence of a small parameter ${\epsilon }^{(u)}$ at the derivative in the Poisson equation shows that the problem belongs to the type of singularly perturbed problems for a quasilinear system of partial differential equations.

2.2.2. The Structure of the Diffusion Layer Near the CEM and the General Idea of the Numerical-analytical Solution

As mentioned above, the Peclet number shows that a diffusion layer arises near the CEM, which in turn consists of an electroneutrality region and a space charge region directly adjacent to the membrane [17]. The structure of the space charge region is determined by the current density, which in turn is determined by the potential jump. At currents below a certain critical value, called the limiting value [17], the diffusion layer consists of an electroneutrality region and a quasi-equilibrium space charge region (quasi-equilibrium SCR), which is several orders of magnitude smaller than the electroneutrality region. At overlimiting currents, in addition to the quasi-equilibrium space charge region, an extended region of the spatial layer arises, which has a size significantly smaller, but already comparable with the electroneutrality region [9]. In the extended SCR, electromigration prevails over diffusion, and in the electroneutrality region, the migration and diffusion flows are equal. In a quasi-equilibrium SCR, the migration flow is equal to the diffusion flow in the first approximation, but they are opposite in direction, which leads to the current becoming zero in the first approximation. In the quasi-equilibrium region of the space charge, the concentrations and electric field strength increase exponentially and therefore their gradients have very large values compared to the values in the region of electroneutrality and the extended region of the space charge. Thus, the boundary value problems of membrane electrochemistry that take into account the quasi-equilibrium region of the space charge are hard problems and are difficult to solve numerically, require an exponential change in the grid step in the quasi-equilibrium region of the space charge and, accordingly, very small time steps, which leads to the accumulation of errors [18]. At the same time, a numerical study of the properties of a quasi-equilibrium SCR made it possible to establish its main properties:

1. quasi-equilibrium SCR is formed almost instantly, in about 10^-5 seconds;

2. the thickness of quasi-equilibrium SCR does not depend on the radius (r) of the membrane disk, except for the vicinity of $r=0$;

3) quasi-equilibrium SCR is also quasi-stationary, i.e. it is practically independent of time;

4) the axial and radial velocities in quasi-equilibrium SCR are close to zero, and the azimuthal can be considered practically constant and equal to $\omega r$, i.e. the electrolyte solution near the membrane in quasi-equilibrium SCR rotates as a whole.

In this connection, a natural idea arises to use an analytical solution in the quasi-equilibrium region in combination with a numerical solution in the remaining region. To implement this idea, three problems must be solved:

1. find an analytical solution in the quasi-equilibrium SCR;

2. determine the boundary conditions on the left boundary of the quasi-equilibrium SCR for the numerical solution;

3. splice the analytical and numerical solutions, i.e. find the constants included in the analytical solution and determine the boundary between the analytical and numerical solutions.

To solve the last two problems, note that at sub-limit currents, the cation concentration decreases in the electroneutrality region and increases in the quasi-equilibrium SCR, and at overlimit currents, the concentration also decreases in almost the entire extended SCR, and the region of increasing cations (RICC) includes a small part of the extended SCR and the entire quasi-equilibrium SCR. Thus, the region of increasing cations (RICC) at sublimiting currents coincides with the quasi-equilibrium SCR, and at overlimiting currents it almost coincides with the quasi-equilibrium SCR, since it also includes a small part of the extended SCR. In addition, RICC has all the above-listed properties of the quasi-equilibrium SCR. Therefore, in what follows, we will use RICC instead of the extended SCR. The boundary between RICC and the rest of the solution region is the points at which the cation concentration has a local minimum.

2.2.3. Analytical Solution of Boundary Value Problem in RICC

The analytical solution in RICC depends on the magnitude of the potential jump, namely, sublimit or overlimit. In the sublimit regime, RICC coincides with the quasi-equilibrium SCR, and in the overlimit regime, these regions are very close. As was noted above, the RICC thickness does not depend on the radial coordinate r, the radial ($u$) and axial ($w$) velocities are close to zero, and the azimuthal ($v$) velocity is almost constant in RICC. Therefore, only the third component remains in equation (1):

$$j_{i,z}={\displaystyle \frac{F}{RT}}{z}_i{D}_i{C}_i{E}_z-D_i{\displaystyle \frac{\partial C_i}{\partial z}},i=1,2$$

In addition, in RICC all unknown functions are practically independent of time t, therefore the left-hand side of equation (2) is equal to zero, and in the right-hand side, due to independence from r, the derivative of the flow with respect to z is equal to zero, in equations (3) and (4) only components dependent on z will also remain and, thus, to find unknown functions in RICC, we obtain a system of equations:

$$j^{(u)}i,z=z_i{D}_i^{(u)}{C}_i^{(u)}{E}_z^{(u)}-D_i^{(u)}{\displaystyle \frac{\partial C_i^{(u)}}{\partial z^{(u)}}},i=1,2$$

(5)

$${\displaystyle \frac{\partial j_{{i,}_z}^{(u)}}{\partial z^{(u)}}}=0,i=1,2$$

(6)

$${\epsilon }^{(u)}{\displaystyle \frac{\partial E_z^{(u)}}{\partial z^{(u)}}}=z_1{C}_1^{(u)}+z_2{C}_2^{(u)}$$

(7)

$$I_z^{(u)}=z_1{j}_{1,z}^{(u)}+z_2{j}_{2,z}^{(u)}$$

(8)

Thus, the system of equations (5-8), for example, for a 1:1 solution of electrolyte with ideally selective CEM takes the form (the index “u” is omitted for simplicity):

$${\displaystyle \frac{d{C}_1}{dz}}=C_1{E}_z-I_z$$

(9)

$${\displaystyle \frac{d{C}_2}{dz}}=-C_2{E}_z$$

(10)

$$\epsilon {\displaystyle \frac{d{E}_z}{dz}}=C_1-C_2$$

(11)

where $E(z,\epsilon )=-{\displaystyle \frac{d\mathsf{\Phi }}{dz}}$ – electric field strength, and $\mathsf{\Phi }(z,\epsilon )$– is its potential; $C_1(z,\epsilon )$, $C_2(z,\epsilon )$ – respectively, the desired concentrations of cations and anions, $I_z$ – current, $\epsilon >0$ – small parameter.

As boundary conditions in dimensionless form for $z={\displaystyle \frac{{\overline{z}}_m}{H}}$, the conditions of merging with the numerical solution are set, and for $z=1$:

$$C_1(1,\epsilon )=C_{1,m}$$

(12)

$$\left(-C_2{\displaystyle \frac{\partial \mathsf{\Phi }}{\partial z}}+{\displaystyle \frac{\partial C_2}{\partial z}}\right)(1,\epsilon )=0$$

(13)

$$\mathsf{\Phi }(1,\epsilon )=0$$

(14)

The system of equations (9-11) represents a singularly perturbed boundary value problem.

1. Solution in the pre-limit case

Let us make the substitution $\xi ={\displaystyle \frac{z-1}{\sqrt{\epsilon }}}$, $E_z={\displaystyle \frac{\tilde{E}(\xi ,\epsilon )}{\sqrt{\epsilon }}}$, $C_1(z,\epsilon )={\tilde{C}}_1(\xi ,\epsilon )$, $C_2(z,\epsilon )={\tilde{C}}_2(\xi ,\epsilon )$, then for small $\epsilon$, in the initial approximation, we obtain the classical system of Boltzmann-Debye equations, which describes the quasi-equilibrium SCR:

$${\displaystyle \frac{d{\tilde{C}}_1}{d\xi }}={\tilde{C}}_1\tilde{E}$$

(15)

$${\displaystyle \frac{d{\tilde{C}}_2}{d\xi }}=-{\tilde{C}}_2\tilde{E}$$

(16)

$${\displaystyle \frac{d\tilde{E}}{d\xi }}={\tilde{C}}_1-{\tilde{C}}_2$$

(17)

with the corresponding boundary conditions (12-14). Therefore, in the pre-limit regime, RICC coincides with the quasi-equilibrium SCR.

The boundary value problem has an approximate analytical solution:

$$E(z,\epsilon )={\displaystyle \frac{1}{\sqrt{\epsilon }}}{\displaystyle \frac{4\sqrt{\beta }{e}^{\sqrt{-\alpha }\frac{z-1}{\sqrt{\epsilon }}}}{1-\beta e^{\sqrt{-4\alpha }\frac{z-1}{\sqrt{\epsilon }}}}}(\sqrt{-\alpha })$$

(18)

$$C_1={\displaystyle \frac{1}{2}}\epsilon {\displaystyle \frac{dE}{dz}}+{\displaystyle \frac{1}{4}}\epsilon E^2-{\displaystyle \frac{1}{2}}\alpha$$

(19)

$$C_2=-{\displaystyle \frac{1}{2}}\epsilon {\displaystyle \frac{dE}{dz}}+{\displaystyle \frac{1}{4}}\epsilon E^2-{\displaystyle \frac{1}{2}}\alpha$$

(20)

where $\alpha =-(C_1({\overline{z}}_m,\epsilon )+C_2({\overline{z}}_m,\epsilon ))\approx -(C_1(1,\epsilon )+C_2(1,\epsilon ))<0$, $\beta$ – some positive number that is determined from the condition $\tilde{E}(0,\epsilon )=\sqrt{2(C_{1,m}+\alpha )}$ equivalent to the condition $C_1(1,\epsilon )=C_{1,m}$.

$E({\overline{z}}_m-0,\epsilon )$ – the value of the numerical solution at the point ${\overline{z}}_m$ must be finite and, accordingly, $E({\overline{z}}_m+0,\epsilon )$ (the value of the analytical solution in RICC) must also be finite. This condition is satisfied if ${\overline{z}}_m=1-k\sqrt{\epsilon }\left|\ln \epsilon \right|$, then we obtain $E({\overline{z}}_m+0,\epsilon )=4\sqrt{\beta }\sqrt{-\alpha }$, for $\epsilon \to 0$, if we take $k={\displaystyle \frac{1}{2\sqrt{-\alpha }}}$. The equality $E({\overline{z}}_m-0,\epsilon )=E({\overline{z}}_m+0,\epsilon )$ is satisfied using a higher approximation.

The points ${\overline{z}}_m$ obtained from the numerical solution by the independent finite element method ${\overline{z}}_m\approx 1.976$ mm and the analytical solution in RICC ${\overline{z}}_m\approx 1.999$ mm coincide with an accuracy of 1.1%.

Concentrations are merged by virtue of the choice $\alpha$ and formulas (19-20).

2. Solution in the overlimit case

The appearance of an extended SCR in the overlimit mode, where ion concentrations are low and the electric field strength tends to infinity as${\displaystyle \frac{1}{\sqrt{\epsilon }}}$, requires a different procedure.

At overlimit currents, there are no anions in the region $({\overline{z}}_m,1]$, i.e. $C_2(z,\epsilon )\equiv 0$, $z\in ({\overline{z}}_m,1]$.

Therefore, the system of equations (9)-(11) takes the form:

$${\displaystyle \frac{d{C}_1}{dz}}=C_1{E}_z-I_z$$

$$\epsilon {\displaystyle \frac{d{E}_z}{dz}}=C_1$$

This system has an approximate analytical solution:

$$E(z,\epsilon )={\displaystyle \frac{1}{\sqrt{\epsilon }}}{\displaystyle \frac{4\sqrt{\beta }{e}^{b\frac{z-1}{\sqrt{\epsilon }}}}{1-\beta e^{2b\frac{z-1}{\sqrt{\epsilon }}}}}b$$

(21)

$$C_1(z,\epsilon )={\displaystyle \frac{4b^2\sqrt{\beta }{e}^{b\frac{z-1}{\sqrt{\epsilon }}}(1+\beta e^{2b\frac{z-1}{\sqrt{\epsilon }}})}{{(1-\beta e^{2b\frac{z-1}{\sqrt{\epsilon }}})}^2}}$$

(22)

where $b>0$ is the integration constant.

3. Splicing solutions and determining constants

Let's put it this way ${\overline{z}}_k=1-k_1\sqrt{\epsilon }+\mathrm{...}\to 1,\epsilon \to 0+$.

For the splicing procedure $E(z,\epsilon )$:

$$E({\overline{z}}_k+0,\epsilon )=E({\overline{z}}_k-0,\epsilon )$$

we use $E({\overline{z}}_k+0,\epsilon )$ from the analytical solution in the interval $\left({\overline{z}}_k,1\right]$, and for the calculation of $E({\overline{z}}_k-0,\epsilon )$ we use the function $E(z,\epsilon )={\displaystyle \frac{\sqrt{2(I_{z}z-I_{пр})}}{\sqrt{\epsilon }}}$, which is a continuation into interval $({\overline{z}}_m,{\overline{z}}_k)$ the solution in the extended SCR. Therefore:

$${\displaystyle \frac{1}{\sqrt{\epsilon }}}{\displaystyle \frac{4\sqrt{\beta }{e}^{-k_{1}b}}{1-\beta e^{-2k_{1}b}}}b={\displaystyle \frac{\sqrt{2(I_z-I_{пр})}}{\sqrt{\epsilon }}}$$

Thus $k_1=-{\displaystyle \frac{1}{b}}\ln \left({\displaystyle \frac{-4b\sqrt{\beta }+\sqrt{8\beta (2b^2+I_z-I_{пр})}}{2\beta \sqrt{2(I_z-I_{пр})}}}\right)$.

Let ${\overline{z}}_m=1-k_2\sqrt{\epsilon }\left|\ln \epsilon \right|+\mathrm{...}\to 1,\epsilon \to 0+$, then:

$$C_1({\overline{z}}_m-0,\epsilon )=C_1({\overline{z}}_m+0,\epsilon )$$

where for the calculation $C_1({\overline{z}}_m-0,\epsilon )$ we use the function $C_1(z,\epsilon )={\displaystyle \frac{I_z}{\sqrt{2(I_{z}z-I_{пр})}}}\sqrt{\epsilon }$ solution in RICC (part of the extended SCR), and in $C_1({\overline{z}}_m+0,\epsilon )$ we use the solution in the interval $\left({\overline{z}}_k,1\right]$, extended to the interval $({\overline{z}}_m,{\overline{z}}_k)$.

Where do we get: $k_2={\displaystyle \frac{1}{2b}}$ and $\beta ={\left({\displaystyle \frac{I_z}{4b^2\sqrt{2(I_z-I_{пр})}}}\right)}^2$.

To find b we use the condition $C_1(1,\epsilon )=C_{1,m}$, from which $b={\displaystyle \frac{1}{2}}\sqrt{{\displaystyle \frac{C_{1,m}{(1-\beta )}^2}{\sqrt{\beta }(1+\beta )}}}$.

2.2.4. Algorithm for Hybrid Numerical-Analytical Solution

1. We numerically solve the boundary value problem of the model without RICC and find $C_1(1,\epsilon )$, $C_2(1,\epsilon )$.

2. We find the potential jump numerically. Next, we find the potential jump for the basic model using the relation:

$$\begin{array}{l}{\mathsf{\Phi }}_0={\displaystyle \underset{0}{\overset{1}{\int }}Е(z,\epsilon )dz}={\displaystyle \underset{0}{\overset{{\overline{z}}_m}{\int }}Е(z,\epsilon )dz}+{\displaystyle \underset{{\overline{z}}_m}{\overset{{\overline{z}}_k}{\int }}Е(z,\epsilon )dz}+{\displaystyle \underset{{\overline{z}}_k}{\overset{1}{\int }}Е(z,\epsilon )dz}=\\ ={\displaystyle \underset{0}{\overset{{\overline{z}}_m}{\int }}Е(z,\epsilon )dz}+{\displaystyle \underset{{\overline{z}}_m}{\overset{{\overline{z}}_k}{\int }}Е(z,\epsilon )dz}-{\displaystyle \underset{{\overline{z}}_k}{\overset{1}{\int }}\frac{d{C}_2}{C_2}=-}\ln {\displaystyle \frac{C_2(1,\epsilon )}{C_2({\overline{z}}_k,\epsilon )}}+{\displaystyle \underset{0}{\overset{{\overline{z}}_m}{\int }}Е(z,\epsilon )dz}+{\displaystyle \underset{{\overline{z}}_m}{\overset{{\overline{z}}_k}{\int }}Е(z,\epsilon )dz}\end{array}$$

Taking into account that ${\overline{z}}_k\approx {\overline{z}}_m\approx 1$, we obtain:

$${\mathsf{\Phi }}_0\approx -\ln {\displaystyle \frac{C_2(1,\epsilon )}{C_2({\overline{z}}_k,\epsilon )}}+{\displaystyle \underset{0}{\overset{1}{\int }}Е(z,\epsilon )dz}$$

$${\mathsf{\Phi }}_0={\mathsf{\Phi }}_{RICC}+{\mathsf{\Phi }}_{NS}$$

Here, the first term ${\mathsf{\Phi }}_{RICC}$– is the potential jump in the region of increasing cations (near CEM), and the second term ${\mathsf{\Phi }}_{NS}$ – is the potential jump calculated numerically.

3. We find the analytical solution in RICC using formulas (18-20) for the sub-limit current mode and using formulas (21-22) for the overlimit current mode.

4. Using steps 1) and 3), we obtain the solution to the boundary value problem of the basic model.

3. Results

3.1. Computational Experiments for Verification of Hybrid Numerical Method

To verify the results of the proposed hybrid numerical-analytical method, their comparison with the results of an independent numerical solution by the finite element method is used. A large number of computational experiments were carried out with different initial concentrations, angular velocities, etc. An independent numerical solution was performed using the finite element method with a non-uniform grid. When performing computational experiments, the number of elements in the grid varied from 6305 to 147415 elements. It turned out that a grid consisting of 9581 elements is quite sufficient, since with a larger number of elements, the numerical results do not change, and with a significantly smaller number, they begin to change. A time step of 0.01 seconds was used.

3.1.1. Comparison of Results for the Annular Disk Model Obtained by the Hybrid Method and the Independent Finite Element Method

Let us consider the results of the numerical study at the initial concentration $C_0=0.01 mol/m^3$. Figure 2a, b show the graphs of the cation concentration and the solution flow lines for the annular disk model obtained using the hybrid method, and Figure 2c, d show those obtained using the independent finite element method. Figure 3 a-d) similarly show the results of the numerical study, but already at the initial concentration $C_0=0.1 mol/m^3$. Figure 2a, c and Figure 3a, c show that the maximum concentration occurs near the boundary of the non-conductivity and conductivity regions, and the concentration also increases on the symmetry axis. Figure 2b, d and Figure 3b, d show that vortices are formed in the conductivity region, and the number of vortices increases with increasing concentration.

An animation of the flow (streamlines) near the membrane disk at the initial concentration $C_0=0.01 mol/m^3$ when solved by the hybrid method is given in [Video1], and when solved by the finite element method in [Video2].

An animation of the flow (streamlines) near the membrane disk at the initial concentration $C_0=0.1 mol/m^3$ when solved by the hybrid method is given in [Video3], and when solved by the finite element method in [Video4].

Thus, at low concentrations of the order of 0.1, 0.01 mol/m³ and less, the finite element method is preferable, and the hybrid method works better with increasing initial concentration of 0.1, 1, 10 mol/m³, etc. In the general area of applicability, the methods lead to results that coincide with an accuracy of about 1%.

3.2. Basic Laws of Occurrence and Development of Electroconvection and Mass Transfer for the Model with an Annular Membrane

To identify the main patterns, a number of numerical experiments were conducted, the initial concentration varied from 0.01 mol/m³ to 100 mol/m³, the angular velocity from ${\displaystyle \frac{\pi }{2}}$ to $4\pi$ rad/s, and the potential jump from 0.1 to 2 V. The most typical results are given below.

Figure 4 shows the solution flow lines at an angular velocity of ${\displaystyle \frac{\pi }{2}}$ rad/s, an initial concentration of 10 mol/m³, and a potential jump of 0.1 V. It is shown that at small potential jumps, vortices do not appear, and the solution flow lines correspond to logarithmic spirals. Also, at a small potential jump, a diffusion layer is formed, the thickness of which is approximately constant, with the exception of a small neighborhood near the junction of the non-conductivity and conductivity regions.

Figure 5 shows the results of the study with a potential jump of 0.3 V (Figure 5a-c) and a potential jump of 1.5 V (Figure 5d-f), initial concentration $C_0=10 mol/m^3$ and variation of the angular velocity: ${\displaystyle \frac{\pi }{2}}$ rad/s (Figure 5a, d), $\pi$ (Figure 5b, e) and $2\pi$ rad/s (Figure 5c, f). It is shown that the greater the angular velocity, the greater the potential jump required for the occurrence of electroconvection.

In real devices, concentrations are of the order of 10 or more mol/m³, so the results for 1 and 10 mol/m³ will be presented below.

Let us consider the case when the initial angular velocity is $\omega ={\displaystyle \frac{\pi }{2}}$rad/s, and the potential jump is $d_{\phi }=1.5$V. The distribution of cation concentration is shown in Figure 6a, where it is evident that the maximum concentration is achieved at the junction of the non-conductivity and conductivity regions. Based on Figure 6b-d), it can be concluded that relaxation associated with the matching of boundary, initial conditions, etc. occurs within 2 seconds. Vortices arise at the junction of the non-conductivity and conductivity regions and flow downstream, gradually decreasing in size (Figure 6d-e). This is due to the fact that the azimuthal component of the velocity $v=\omega r$ increases with increasing r.

An animation of the flow (streamlines) near the membrane disk at the initial concentration $C_0=1 mol/m^3$ is shown in [Video5], and at $C_0=10 mol/m^3$ in [Video6].

4. Discussion

The article considers the influence of heterogeneity areas on the transport of salt ions in systems with axial symmetry taking into account electroconvection. The process of formation of electroconvective vortices in the potentiostatic mode is shown with the addition of conductivity and non-conductivity regions to the cation-exchange rotating membrane disk, a modification in the form of an annular disk is considered. An algorithm for a hybrid numerical-analytical solution is formulated, which allows performing calculations at a high initial concentration, potential jump and disk rotation speed. Using the hybrid solution method, the main patterns of occurrence and development of electroconvection and mass transfer for the model with an annular membrane are derived, and a comparison of the results for the annular disk model obtained by the hybrid method and the independent finite element method is carried out. The vortex formation process can be enhanced by adding non-conductivity regions to the CEM, while if the annular membrane model is used, the vortices are formed at the junction of the non-conductivity and conductivity regions. At a small potential jump and/or at a high angular velocity, a diffusion layer is formed, the thickness of which is approximately constant, with the exception of a small area near the junction of the regions of non-conductivity and conductivity.