1. Introduction
Electrokinetic phenomena (Stone, Stroock & Ajdari Reference Stone, Stroock and Ajdari2004; Squires & Quake Reference Squires and Quake2005) are ubiquitous in nature, such as porous media (Revil & Linde Reference Revil and Linde2006; Revil et al. Reference Revil, Linde, Cerepi, Jougnot, Matthäi and Finsterle2007), charged microchannels (Xuan, Sinton & Li Reference Xuan, Sinton and Li2004; Pascall & Squires Reference Pascall and Squires2011), electrochemical systems (Xiao et al. Reference Xiao, Wang, Zhao, Zhang, Song and Vikesland2024) and polarisable interfaces (Bazant & Squires Reference Bazant and Squires2010; Revil et al. Reference Revil, Ghorbani, Su, Cai and Hu2024). In this paper we focus on electrokinetic effects near polarisable (metal) surfaces. Early theoretical models (Dukhin Reference Dukhin1991) of electro-osmotic flows around polarisable surfaces were developed by Dukhin. The mechanisms were later formalised under the term induced-charge electro-osmosis (ICEO) by Ramos et al. (Reference Ramos, Morgan, Green and Castellanos1998, Reference Ramos, Morgan, Green and Castellanos1999). When an external electric field is applied to a metal surface, the surface becomes polarised, inducing a redistribution of ions in the surrounding electrolyte. This leads to the formation of an asymmetric electric double layer (EDL) adjacent to the metal. The accumulated induced charges interact with the tangential component of the electric field within the EDL, generating electro-osmotic slip and driving vortical flows in the bulk. This phenomenon is referred to as ICEO (Ramos et al. Reference Ramos, Morgan, Green and Castellanos1998, Reference Ramos, Morgan, Green and Castellanos1999). A canonical asymptotic model has been developed (Bazant & Squires Reference Bazant and Squires2004, Reference Bazant and Squires2010; Squires & Bazant Reference Squires and Bazant2004; Pascall & Squires Reference Pascall and Squires2010), especially under assumptions of low potentials and a thin EDL, which is
$u_{s}\sim ({\varepsilon _{w}E_{0}^{2}a}/{\mu })$
, where
$u_{s}$
is the induced velocity along the external electrical field
$E_{0}$
,
$\varepsilon _{w}=80\varepsilon _{0}$
is the dielectric constant of the aqueous solution,
$\varepsilon _{0}$
is the vacuum dielectric constant,
$a$
is the effective radius of the surface and
$\mu$
is the solution viscosity.
In weak-field conditions, this mechanism is well captured by the canonical asymptotic model, whereas under strongly polarised surfaces, extended space-charge (ESC) layers and nonlinear charging currents substantially modify the induced flows. The term ‘strongly polarised surfaces’ is used to describe the strongly ion-depleted state under high applied voltages, where capacitive charging of the EDL leads to nonlinear ion transport and the emergence of ESC structures. The classical model has been shown to significantly overestimate the ICEO velocity (Levitan et al. Reference Levitan, Devasenathipathy, Studer, Ben, Thorsen, Squires and Bazant2005; Harnett et al. Reference Harnett, Templeton, Dunphy-Guzman, Senousy and Kanouff2008). Strong nonlinearities arise under strong electric fields, which may lead to flow instabilities (Bazant, Thornton & Ajdari Reference Bazant, Thornton and Ajdari2004). Consequently, the discrepancies between predictions and experimental results can be as large as several orders of magnitude. Due to surface conduction (Messinger & Squires Reference Messinger and Squires2010; Liu et al. Reference Liu, Li, Song, Pan, Li and Xuan2017), a significant amount of counter-ions was required to neutralise field-induced surface charges, resulting in the formation of a high-concentration EDL on the upper and lower surfaces of the cylinder. To theoretically match experimental data of ICEO velocities, some models introduced an adjustable parameter for the ionic radius. However, the actual ion radius is only a few angstroms, which prompted further corrections for finite ion radius effects in the EDL to explain deviations from the canonical asymptotic model (Bazant & Squires Reference Bazant and Squires2010; Bazant, Storey & Kornyshev Reference Bazant, Storey and Kornyshev2011). In this framework, the effective ionic radius modifies the internal structure of the EDL, and such structural changes can in turn alter the slip velocity at the interface, providing one possible mechanism for the observed deviations. In addition, the Debye length may significantly affect velocity prediction (Gregersen et al. Reference Gregersen, Andersen, Soni, Meinhart and Bruus2009; Schnitzer & Yariv Reference Schnitzer and Yariv2012).
The study to date by direct numerical simulation (DNS) (Davidson, Andersen & Mani Reference Davidson, Andersen and Mani2014) of ICEO under intense electric fields demonstrated that the localised flow instability did not depend on adjustable parameters such as ion size (Bazant et al. Reference Bazant, Storey and Kornyshev2011) or inequality of ionic diffusion coefficients (Khair & Balu Reference Khair and Balu2020), but rather on the ESC layer (Rubinstein & Zaltzman Reference Rubinstein and Zaltzman2000; Gu et al. Reference Gu, Xu, Huo, Rubinstein, Bazant and Deng2019). When the external voltage was much higher than the thermal voltage, numerous small-scale vortices formed on the metallic cylinder, causing localised flow disruption and dissipating flow energy. These chaotic vortices were sometimes also called electroconvection (Rubinstein & Zaltzman Reference Rubinstein and Zaltzman2000; Gu et al. Reference Gu, Xu, Huo, Rubinstein, Bazant and Deng2019; Kang & Kwak Reference Kang and Kwak2020; Guan et al. Reference Guan, He, Wang, Song, Zhang and Wu2021; Traoré & Wu Reference Traoré and Wu2013; Zaltzman & Rubinstein Reference Zaltzman and Rubinstein2007). In recent years, Kim et al. have made advances through DNS in elucidating the formation of ESC layers near electrodes and the resultant chaotic vortices (Kim, Davidson & Mani Reference Kim, Davidson and Mani2019). Their work has substantially deepened our understanding of the interplay between ESC structures and the flow dynamics, providing a crucial foundation for both experimental observations and theoretical interpretations of chaotic vortices (Davidson et al. Reference Davidson, Andersen and Mani2014; Kim et al. Reference Kim, Davidson and Mani2019). However, the effects of the hydrodynamic coupling coefficient and the dimensionless Debye layer length on the flow have not been explored (Davidson et al. Reference Davidson, Andersen and Mani2014; Kim et al. Reference Kim, Davidson and Mani2019). Although previous studies have demonstrated that the ESC layer can reduce the ICEO velocity, the underlying mechanism governing ESC formation remains insufficiently understood, particularly how EDL charging triggers the development of the ESC (Davidson et al. Reference Davidson, Andersen and Mani2014; Lazo et al. Reference Lazo, Peng, Xiang, Shiyanovskii and Lavrentovich2014; Peng et al. Reference Peng, Lazo, Shiyanovskii and Lavrentovich2014; Kim et al. Reference Kim, Davidson and Mani2019). Moreover, the canonical asymptotic models do not account for the influence of the ESC, and therefore cannot capture the scaling behaviour of the velocity under high-voltage conditions, limiting their applicability in analysing ICEO phenomena (Bazant & Squires Reference Bazant and Squires2004; Pascall & Squires Reference Pascall and Squires2010; Davidson et al. Reference Davidson, Andersen and Mani2014).
In this paper, our study focuses on an ICEO system devoid of externally imposed fluid disturbances and selective membranes, where flow arises purely from the nonlinear coupling between electric fields and ionic migration near the metal surface. We observe a velocity scaling transition driven by the strong EDL charging under strong electric fields, which differs qualitatively from the shear flow-dominated systems (Kang & Kwak Reference Kang and Kwak2020; Liu et al. Reference Liu, Zhou and Shi2020b , Reference Liu, Zhou and Shi2022). To uncover the underlying physics, we combine theoretical analysis with DNS to develop a framework centred on the local breakdown of electroneutrality, which leads to a velocity scaling transition. This framework not only advances the mechanistic understanding of ESC formation but also systematically explains the physical origin of the velocity scaling shift in unstable flow near strongly polarised surfaces. The strong EDL charging, along with localised flow instabilities, is responsible for the discrepancies between the theoretical and experimental data. Our work thus extends and complements the existing work (Bazant & Squires Reference Bazant and Squires2004; Pascall & Squires Reference Pascall and Squires2010; Davidson et al. Reference Davidson, Andersen and Mani2014; Kim et al. Reference Kim, Davidson and Mani2019), offering new theoretical insights that enrich the understanding of strong EDL charging.
2. Theory
2.1. Dimensionless governing equations
The ICEO system involves multiscale spatio-temporal evolution and multiphysics coupling of fluid flow, ionic transport and electric fields, spanning from nanometres to micrometres. This dynamics can be described by the coupled Poisson–Nernst–Planck and incompressible Navier–Stokes equations (Davidson et al. Reference Davidson, Andersen and Mani2014; Kim et al. Reference Kim, Davidson and Mani2019). For clarity and analytical convenience, the governing equations are presented in their non-dimensional form. The characteristic reference parameters for non-dimensionalisation are selected as follows: the cylinder radius
$a$
, diffusion velocity
$D_{0}/a$
, diffusion time
$a^{2}/D_{0}$
, bulk concentration
$c_{0}$
, osmotic pressure
$\mu D_{0}/a^{2}$
, diffusion flux
$D_{0}c_{0}/a$
and thermal voltage
$k_{B}T/e$
. Here,
$D_{0}$
is the diffusion coefficient,
$k_{B}$
is the Boltzmann constant,
$T$
is the temperature and
$e$
is the elementary charge. The dimensionless forms of governing equations for transports of mass, momentum, species and electrical fields are
\begin{align}\frac{\partial \tilde{c}_{1}}{\partial \tilde{t}}+\tilde{\boldsymbol{\nabla }}\boldsymbol{\cdot }{\boldsymbol{\tilde{J}}}_{1}=0, \\[-28pt] \nonumber \end{align}
\begin{align}{\boldsymbol{\tilde{J}}}_{1}=-\tilde{D}_{1}\tilde{\boldsymbol{\nabla }}\tilde{c}_{1}-\tilde{D}_{1}\tilde{c}_{1}\tilde{\boldsymbol{\nabla }}\tilde{\varphi }+\tilde{\boldsymbol{U}}\boldsymbol{\cdot }\tilde{c}_{1},\\[-28pt] \nonumber \end{align}
\begin{align}\frac{\partial \tilde{c}_{2}}{\partial \tilde{t}}+\tilde{\boldsymbol{\nabla }}\boldsymbol{\cdot }{\boldsymbol{\tilde{J}}}_{2}=0,\\[-28pt] \nonumber \end{align}
\begin{align}{\boldsymbol{\tilde{J}}}_{2}=-\tilde{D}_{2}\tilde{\boldsymbol{\nabla }}\tilde{c}_{2}+\tilde{D}_{2}\tilde{c}_{2}\tilde{\boldsymbol{\nabla }}\tilde{\varphi }+\tilde{\boldsymbol{U}}\boldsymbol{\cdot }\tilde{c}_{2},\\[-28pt] \nonumber \end{align}
\begin{align}-\tilde{\lambda }_{D}^{2}\tilde{\boldsymbol{\nabla} }^{2}\tilde{\varphi }=\tilde{c}_{1}-\tilde{c}_{2},\\[-28pt] \nonumber \end{align}
\begin{align}\widetilde{Re}\frac{\partial \tilde{\boldsymbol{U}}}{\partial \tilde{t}}+\widetilde{Re}\tilde{\boldsymbol{U}}\boldsymbol{\cdot }\boldsymbol{\nabla }\tilde{\boldsymbol{U}}-\tilde{\boldsymbol{\nabla} }^{2}\tilde{\boldsymbol{U}}=-\tilde{\boldsymbol{\nabla }}\tilde{P}-\tilde{\kappa }\frac{\left(\tilde{c}_{1}-\tilde{c}_{2}\right)}{\tilde{\lambda }_{D}^{2}}\tilde{\boldsymbol{\nabla }}\tilde{\varphi },\\[-28pt] \nonumber \end{align}
\begin{align}\tilde{\boldsymbol{\nabla }}\boldsymbol{\cdot }\tilde{\boldsymbol{U}}=0,\\[-6pt] \nonumber\end{align}
where
$\tilde{c}_{1}$
and
$\tilde{c}_{2}$
represent the dimensionless concentrations of cations and anions, respectively,
${\boldsymbol{\tilde{J}}}_{1}$
and
${\boldsymbol{\tilde{J}}}_{2}$
are the dimensionless fluxes of cations and anions, respectively,
$\tilde{D}_{1}$
and
$\tilde{D}_{2}$
are the dimensionless diffusion coefficients of the cations and anions, respectively;
$\tilde{D}_{1}=1.13$
and
$\tilde{D}_{2}=2.03$
(NaCl) are considered here. Also,
$\tilde{\varphi }$
,
$\tilde{\boldsymbol{U}}$
and
$\tilde{t}$
denote the dimensionless concentration, velocity and time, respectively. The Reynolds number, denoted as
$\widetilde{Re}$
, is given by
$\rho D_{0}/\mu$
, where
$\rho$
is the density and
$\mu$
is the viscosity. We consider an aqueous solution with a reference ion diffusion coefficient of
$D_{0}=2\times10^{-9}\mathrm{m}^{2}\,\mathrm{s}^{-1}$
. Consequently, the Reynolds number is a fixed value of 0.001 (Davidson et al. Reference Davidson, Andersen and Mani2014).
Three key dimensionless parameters can be distinguished from the governing equations. The first one is the hydrodynamic coupling coefficient
\begin{equation}\tilde{\kappa }=\frac{V_{T}^{2}\varepsilon _{w}}{\mu D_{0}},\end{equation}
where
$V_T$
is the thermal voltage, which describes the influence of the physical properties of the salt solution on the ICEO. By altering the properties of the salt solution, the flow state of ICEO can be adjusted.
The second key dimensionless parameter is the dimensionless Debye length
\begin{equation}\tilde{\lambda }_{D}=\frac{\sqrt{\varepsilon _{w}k_{B}T/[(ze)^{2}c_{0}]}}{a},\end{equation}
where
$z$
is the chemical valence, which describes the ratio of the Debye length (
$\sqrt{\varepsilon _{w}k_{B}T/[(ze)^{2}c_{0}]}$
) to the radius of the metallic cylinder (
$a$
). Adjusting the radius of the cylinder and the bulk concentration can also influence the local flow field and charge layer characteristics.
The last is the commonly used dimensionless electric potential difference
\begin{equation}\Delta \tilde{\varphi}=\frac{\Delta \varphi }{k_{B}T/e},\end{equation}
which represents the ratio of the electric potential difference between the upper and lower electrodes
$(\Delta \varphi )$
at the two ends. The dimensionless voltage above is the most critical parameter for triggering the electrokinetic response.
This study focuses on a canonical ICEO configuration (Davidson et al. Reference Davidson, Andersen and Mani2014), in which a two-dimensional metallic cylinder is immersed in an aqueous binary electrolyte subjected to a DC electric field. A square domain is considered, with a fixed metal cylinder (
$\tilde{a}=1$
) placed at its centre. The domain size (
$\tilde{L}=50$
) is large enough to ensure that the results are independent of the boundary. At the cylinder surface (
$\tilde{a}=1$
), boundary conditions include no flux, no slip, no penetration and a fixed electric potential (Davidson et al. Reference Davidson, Andersen and Mani2014). In the far field, given concentrations (
$\tilde{c}_{1}=\tilde{c}_{2}=1$
) and a uniform electric field (
$\Delta \tilde{\varphi }/\tilde{L}$
) are imposed (Davidson et al. Reference Davidson, Andersen and Mani2014).
2.2. Theoretical analysis of ICEO velocity
We begin by deriving the velocity of ICEO. Two distinct mechanisms are considered: (i) at low voltages, the ICEO velocity is driven by the EDL, where the current induced by EDL charging can be neglected (Bazant & Squires Reference Bazant and Squires2004); and (ii) at high voltages, the ICEO velocity is governed by the ESC, where the charging current of the EDL plays a central role (Davidson et al. Reference Davidson, Andersen and Mani2014). Ignoring the angular variations, the velocity in the standard asymptotic model can be described by (2.11) (Squires & Bazant Reference Squires and Bazant2004). Based on the dimensionless reference values used in the main text, the dimensionless velocity
$\tilde{u}_{s}$
in the canonical asymptotic model is given by
\begin{equation}\tilde{u}_{s}\sim (\Delta \tilde{\varphi }/\tilde{L})^{2}\;\text{for the EDL}\hbox{-}\text{triggered ICEO}.\end{equation}
Equation (2.11) describes the EDL-triggered velocity. As the voltage increases, ions accumulate on the metallic cylinder surface. This leads to enhanced capacitive charging of the EDL, generating an ion current. The ion current is induced on the upper surface of the polarised metal primarily due to the cation flux, while, another additional ion current appears on the lower part of the polarised metal primarily due to the anion flux. This analysis will be confirmed later based on our DNS.
Here, we derive the scaling of velocity at the ESC. The transition from EDL to ESC primarily depends on the magnitude of the ion flux generated during EDL charging. By direct simulations, the critical flux is approximately O(1) for our cases. As analysed in Davidson et al. (Reference Davidson, Andersen and Mani2014), we focus more on the statistical behaviour of velocities at the inflow points and outflow points (0° and 180°). In our theoretical analysis, the angular variations are similarly ignored, and a local Cartesian coordinate system is used to analyse the statistical behaviour of these inflow and outflow points, which is established with the centre of the circle at (
$\tilde{x}_{0}, \tilde{y}_{0}$
). For convenience, we consider
$\tilde{x}=\tilde{x}_{0}-0$
and
$\tilde{y}=\tilde{y}_{0}-1$
.
We draw on existing theories (Rubinstein & Zaltzman Reference Rubinstein and Zaltzman2001; Zaltzman & Rubinstein Reference Zaltzman and Rubinstein2007; Xu et al. Reference Xu, Gu, Liu, Huo, Zhou and S.2020) and apply them to the ICEO analysis, and finally derive the velocity scaling law in ESC-triggered ICEO.
When the flow is in a steady state and ICEO is at low Reynolds numbers, the inertial and transient terms can be neglected. Thus, the Navier–Stokes equation is simplified to
\begin{equation}-\tilde{{\nabla} }^{2}\tilde{\boldsymbol{U}}=-\tilde{\boldsymbol{\nabla }}\!\tilde{P}-\tilde{\kappa }\frac{(\tilde{c}_{1}-\tilde{c}_{2})}{\tilde{\lambda }_{D}^{2}}\tilde{\boldsymbol{\nabla }}\tilde{\varphi }.\end{equation}
It is well known that the vortex roll is triggered by the slip velocity, expected to yield a quantitative understanding of the flow. The Stokes equation in the ESC layer can be simplified to the following form (Rubinstein & Zaltzman Reference Rubinstein and Zaltzman2001):
\begin{equation}\frac{\partial ^{2}\tilde{u}}{\partial \tilde{y}^{2}}=\frac{1}{2}\tilde{\kappa }\frac{\partial }{\partial \tilde{x}}\left(\frac{\partial \tilde{\varphi }}{\partial \tilde{y}}\right)^{2}-\tilde{\kappa }\left(\frac{\partial ^{2}\tilde{\varphi }}{\partial \tilde{y}^{2}}\right)\frac{\partial \tilde{\varphi }}{\partial \tilde{x}}.\end{equation}
Equation (2.13) shows that the flow is dominated by tangential variations, allowing the neglect of the second-order derivative of velocity along the tangential direction. The flow is in a quasi-steady state with a low Reynolds number, dominated by viscous, Maxwell pressure and electrostatic forces. The first term on the right-hand side,
$(1/2)\tilde{\kappa }({\partial }/{\partial \tilde{x}})({\partial \tilde{\varphi }}/{\partial \tilde{y}})^{2}$
, describes the contribution of the pressure driving force to the slip velocity (Rubinstein & Zaltzman Reference Rubinstein and Zaltzman2001). The second term,
$-\tilde{\kappa}({\partial ^{2}\tilde{\varphi }}/{\partial \tilde{y}^{2}})(\partial \tilde{\varphi }/\partial \tilde{x})$
, describes the contribution of the electrostatic driving force within the ESC layer. The third term,
${\partial ^{2}\tilde{u}}/{\partial \tilde{y}^{2}}$
, describes the viscous force.
The velocity primarily depends on the electric potential distribution within the ESC layer. Integrating (2.13) twice over the range
$0\leq \tilde{y}\leq \tilde{y}_{\textit{esc}}$
yields the solution for the slip velocity associated with ESC-triggered ICEO, where
$\tilde{y}_{\textit{esc}}$
is the ESC thickness. In the steady-state Nernst–Planck equation, neglecting the transient terms, we obtain
\begin{equation}\tilde{\boldsymbol{\nabla }}\boldsymbol{\cdot }{\boldsymbol{\tilde{J}}}_{1}=0,\tilde{\boldsymbol{\nabla }}\boldsymbol{\cdot }{\boldsymbol{\tilde{J}}}_{2}=0.\end{equation}
Our DNS shows that the accumulation of cations on the surface of the metallic cylinder, strong capacitive charging occurs in the upper EDL, resulting in a local cation current (
$\tilde{j}_{1}$
) near the metal surface, while the local anion current is nearly zero. In addition, we neglected the convection term in the Nernst–Planck equation. This assumption is justified because the ESC region is extremely thin, of the order of the Debye length, such that the influence of fluid motion on the potential distribution within the ESC is negligible. Therefore, (2.14) can be simplified to a one-dimensional transport form:
$\tilde{j}_{1}=({{\rm d}\tilde{c}_{1}}/{{\rm d}\tilde{y}})+\tilde{c}_{1}({{\rm d}\tilde{\varphi }}/{{\rm d}\tilde{y}})$
and
$0=({{\rm d}\tilde{c}_{2}}/{{\rm d}\tilde{y}})-\tilde{c}_{2}({{\rm d}\tilde{\varphi }}/{{\rm d}\tilde{y}})$
. Adding these two equations yields
\begin{equation}\tilde{j}_{1}=\frac{{\rm d}\left(\tilde{c}_{1}+\tilde{c}_{2}\right)}{{\rm d}\tilde{y}}+\frac{{\rm d}\tilde{\varphi }}{{\rm d}\tilde{y}}(\tilde{c}_{1}-\tilde{c}_{2}).\end{equation}
Consider the one-dimensional form of the Poisson equation, we have
\begin{equation}-\frac{{\rm d}^{2}\tilde{\varphi }}{{\rm d}\tilde{y}^{2}}=\frac{(\tilde{c}_{1}-\tilde{c}_{2})}{\tilde{\lambda }_{D}^{2}}.\end{equation}
Substituting (2.16) into (2.15), we can obtain
\begin{equation}\tilde{j}_{1}=\frac{{\rm d}\left(\tilde{c}_{1}+\tilde{c}_{2}\right)}{{\rm d}\tilde{y}}-\tilde{\lambda }_{D}^{2}\frac{{\rm d}\tilde{\varphi }}{{\rm d}\tilde{y}}\frac{{\rm d}^{2}\tilde{\varphi }}{{\rm d}y^{2}}.\end{equation}
The integration limits in (2.17) are from
$\tilde{y}_{\textit{esc}}$
to 1, and an expression for the salt concentration can be obtained as
\begin{equation}\tilde{c}_{1}+\tilde{c}_{2}=\tilde{j}_{1}\left(\tilde{y}-\tilde{y}_{\textit{esc}}\right)+\frac{1}{2}\tilde{\lambda }_{D}^{2}\left(\frac{{\rm d}\tilde{\varphi }}{{\rm d}\tilde{y}}\right)^{2},\end{equation}
where
$\tilde{y}_{\textit{esc}}$
is an integration constant, corresponding to the thickness of the ESC layer. Furthermore, based on the equations
$\tilde{j}_{1}=({{\rm d}\tilde{c}_{1}}/{{\rm d}\tilde{y}})+\tilde{c}_{1}({{\rm d}\tilde{\varphi }}/{{\rm d}\tilde{y}})$
and
$0=({{\rm d}\tilde{c}_{2}}/{{\rm d}\tilde{y}})-\tilde{c}_{2}({{\rm d}\tilde{\varphi }}/{{\rm d}\tilde{y}})$
, the transport for the space charge can also be derived
\begin{equation}\tilde{j}_{1}=\frac{{\rm d}\left(\tilde{c}_{1}-\tilde{c}_{2}\right)}{{\rm d}\tilde{y}}+\frac{{\rm d}\tilde{\varphi }}{{\rm d}\tilde{y}}(\tilde{c}_{1}+\tilde{c}_{2}).\end{equation}
Substituting (2.16) and (2.18) into (2.19), we obtain
\begin{equation}\tilde{j}_{1}=\frac{{\rm d}\left(-\tilde{\lambda }_{D}^{2}\frac{{\rm d}^{2}\tilde{\varphi }}{{\rm d}\tilde{y}^{2}}\right)}{{\rm d}\tilde{y}}+\frac{{\rm d}\tilde{\varphi }}{{\rm d}\tilde{y}}\left(\tilde{j}_{1}\left(\tilde{y}-\tilde{y}_{\textit{esc}}\right)+\frac{1}{2}\tilde{\lambda }_{D}^{2}\left(\frac{{\rm d}\tilde{\varphi }}{{\rm d}\tilde{y}}\right)^{2}\right).\end{equation}
Multiplying both sides of (2.20) by the dimensionless Debye length, yields
\begin{equation}\tilde{\lambda }_{D}\tilde{j}_{1}=\tilde{\lambda }_{D}\frac{{\rm d}\left(-\tilde{\lambda }_{D}^{2}\frac{{\rm d}^{2}\tilde{\varphi }}{{\rm d}\tilde{y}^{2}}\right)}{{\rm d}\tilde{y}}+\tilde{\lambda }_{D}\frac{{\rm d}\tilde{\varphi }}{{\rm d}\tilde{y}}\left(\tilde{j}_{1}\left(\tilde{y}-\tilde{y}_{\textit{esc}}\right)+\frac{1}{2}\tilde{\lambda }_{D}^{2}\left(\frac{{\rm d}\tilde{\varphi }}{{\rm d}\tilde{y}}\right)^{2}\right).\end{equation}
Due to the small value of
$\tilde{\lambda }_{D}$
and the cation flux
$\tilde{j}_{1}$
being
$O(10)$
, thus
$\tilde{\lambda }_{D}\tilde{j}_{1}\approx 0$
. For
$\tilde{\lambda }_{D}({{\rm d}(-\tilde{\lambda }_{D}^{2}({{\rm d}^{2}\tilde{\varphi }}/{{\rm d}\tilde{y}^{2}}))}/{{\rm d}\tilde{y}})=-\tilde{\lambda }_{D}^{3}({{\rm d}({{\rm d}^{2}\tilde{\varphi }}/{{\rm d}\tilde{y}^{2}})}{{\rm d}\tilde{y}})$
, taking the first-order approximation, we obtain
$-\tilde{\lambda }_{D}^{3}({{\rm d}({{\rm d}^{2}\tilde{\varphi }}/{{\rm d}\tilde{y}^{2}})}{{\rm d}\tilde{y}})=-\tilde{\lambda }_{D}^{3}({\tilde{\varphi }}/{(\tilde{y}_{\textit{esc}})^{3}})$
, where
$\tilde{y}_{\textit{esc}}\approx \tilde{\lambda }_{D}^{2/3}\tilde{\varphi }^{2/3}$
. The analytical expression for
$\tilde{y}_{\textit{esc}}$
will be derived later (2.16). Therefore,
$\tilde{\lambda }_{D}({{\rm d}(-\tilde{\lambda }_{D}^{2}({{\rm d}^{2}\tilde{\varphi }}/{{\rm d}\tilde{y}^{2}})})/{{\rm d}\tilde{y}})\approx 0$
. Equation (2.21) can be simplified to
\begin{equation}\tilde{\lambda }_{D}\frac{{\rm d}\tilde{\varphi }}{{\rm d}\tilde{y}}\left(\tilde{j}_{1}\left(\tilde{y}-\tilde{y}_{\textit{esc}}\right)+\frac{1}{2}\tilde{\lambda }_{D}^{2}\left(\frac{{\rm d}\tilde{\varphi }}{{\rm d}\tilde{y}}\right)^{2}\right)=0.\end{equation}
Therefore, (2.22) has a physically meaningful solution given by
\begin{equation}\tilde{\lambda }_{D}\frac{{\rm d}\tilde{\varphi }}{{\rm d}\tilde{x}}=-\sqrt{2\tilde{j}_{1}\left(\tilde{y}_{\textit{esc}}-y\right)}.\end{equation}
When the cation flux vanishes, no ESC layer forms. In contrast, the presence of a finite ion flux induces a strong electric field within the ESC, which is essential for its formation. Integrating (2.23) twice in the ESC layer (
$0\leq \tilde{y}\leq \tilde{y}_{\textit{esc}}$
) yields
\begin{equation}\tilde{\varphi }=\left(\frac{2}{3}\right)\big(\tilde{\lambda }_{D}^{-1}\sqrt{2\tilde{j}_{1}}\big)\big[\tilde{y}_{\textit{esc}}^{\frac{3}{2}}-\left(\tilde{y}_{\textit{esc}}-\tilde{y}\right)^{\frac{3}{2}}\big]+{C},\end{equation}
where
$\tilde{\varphi }$
describes the electric potential distribution within the ESC layer, and
${C}$
is the integration constant. It is important to note that
${C}$
does not affect the derivation of the velocity and is therefore not solved here. The polarisation potential on the metal column surface is
$\tilde{\varphi }_{\textit{pmi}}$
. Then, the electric potential difference within the ESC layer is given by
\begin{equation}\tilde{\varphi }_{\textit{esc}}-\tilde{\varphi }_{\textit{pmi}}=\frac{2}{3}\big(\tilde{\lambda }_{D}^{-1}\sqrt{2\tilde{j}_{1}}\big)(\tilde{y}_{\textit{esc}})^{\frac{3}{2}}\!,\end{equation}
where
$\tilde{\varphi }$
describes the electric potential distribution within the ESC layer, and
${C}$
is the integration constant. It is important to note that
${C}$
does not affect the derivation of the velocity and is therefore not solved here. Denoting
$\tilde{\varphi }(\tilde{y}=0)=\tilde{\varphi }_{\textit{pmi}}$
and
$\tilde{\varphi }(\tilde{y}=\tilde{y}_{\textit{esc}})=\tilde{\varphi }_{\textit{esc}}$
, the electric potential difference within the ESC layer is given by
\begin{equation}\tilde{y}_{\textit{esc}}=\frac{\sqrt[3]{9}}{2}\tilde{\lambda }_{D}^{\frac{2}{3}}\big(\tilde{j}_{1}\big)^{-\frac{1}{3}}({\tilde{\varphi }_{\textit{esc}}}-{\tilde{\varphi }_{\textit{pmi}}})^{2/3}.\end{equation}
To obtain the expression of the velocity profile, based on the electric potential distribution within the ESC and the expression for the thickness of the ESC layer, we substitute (2.24) into (2.23) and note that both the cation flux
$\tilde{j}_{1}$
and the ESC layer thickness
$\tilde{y}_{\textit{esc}}$
are functions of
$\tilde{y}$
and thus obtain
\begin{align}\frac{\partial ^{2}\tilde{u}}{\partial \tilde{y}^{2}}&=\frac{1}{2}\tilde{\kappa }\frac{\partial }{\partial \tilde{x}}\left(\frac{\partial \big[\frac{2}{3}\big(\tilde{\lambda }_{D}^{-1}\sqrt{2\tilde{j}_{1}}\big)\big[-\big(\tilde{y}_{\textit{esc}}-\tilde{y}\big)^{\frac{3}{2}}\big]\big]}{\partial \tilde{y}}\right)^{2} \nonumber \\ &\quad -\tilde{\kappa }\left(\frac{\partial }{\partial \tilde{y}}\frac{\partial \big[\frac{2}{3}\big(\tilde{\lambda }_{D}^{-1}\sqrt{2\tilde{j}_{1}}\big)\big[-\big(\tilde{y}_{\textit{esc}}-\tilde{y}\big)^{\frac{3}{2}}\big]\big]}{\partial \tilde{y}}\right)\frac{\partial \big[\frac{2}{3}\big(\tilde{\lambda }_{D}^{-1}\sqrt{2\tilde{j}_{1}}\big)\big[-\big(\tilde{y}_{\textit{esc}}-\tilde{y}\big)^{\frac{3}{2}}\big]\big]}{\partial \tilde{x}}. \end{align}
The cation flux
$\tilde{j}_{1}$
and the ESC thickness
$\tilde{y}_{\textit{esc}}$
are functions of
$\tilde{y}$
, and (2.27) simplifies to
\begin{equation}\tilde{\kappa }^{-1}\frac{\partial ^{2}\tilde{u}}{\partial \tilde{y}^{2}}=\frac{2}{3}\tilde{\lambda }_{D}^{-2}\left(\tilde{y}_{\textit{esc}}-\tilde{y}\right)\frac{\partial \tilde{j}_{1}}{\partial \tilde{x}}+\frac{1}{3}\frac{1}{\sqrt{\tilde{y}_{\textit{esc}}-\tilde{y}}}\tilde{\lambda }_{D}^{-2}\tilde{y}_{\textit{esc}}^{\frac{3}{2}}\frac{\partial \tilde{j}_{1}}{\partial \tilde{x}}+\tilde{\lambda }_{D}^{-2}\tilde{j}_{1}\frac{\tilde{y}_{\textit{esc}}^{\frac{1}{2}}}{\sqrt{\tilde{y}_{\textit{esc}}-\tilde{y}}}\frac{\partial \tilde{y}_{\textit{esc}}}{\partial \tilde{x}}.\end{equation}
By integrating (2.28) once over
$0\leq \tilde{y}\leq \tilde{y}_{\textit{esc}}$
, and considering that vortex streamlines cannot penetrate the metal surface but slide tangentially along it (Davidson et al. Reference Davidson, Andersen and Mani2014), we conclude that
${\partial \tilde{u}}/{\partial \tilde{y}}=0$
and
${\partial \tilde{u}}/{\partial \tilde{x}}\neq 0$
at ESC. Thus, we have
\begin{align}\tilde{\kappa }^{-1}\frac{\partial \tilde{u}}{\partial \tilde{y}}&=-\frac{1}{3}\tilde{\lambda }_{D}^{-2}\left(\tilde{y}_{\textit{esc}}-\tilde{y}\right)^{2}\frac{\partial \tilde{j}_{1}}{\partial \tilde{x}} \nonumber \\ &\quad-\frac{2}{3}\left(\tilde{y}_{\textit{esc}}-\tilde{y}\right)^{\frac{1}{2}}\tilde{\lambda }_{D}^{-2}\tilde{y}_{\textit{esc}}^{\frac{3}{2}}\frac{\partial \tilde{j}_{1}}{\partial \tilde{x}}-\tilde{\lambda }_{D}^{-2}\tilde{j}_{1}\tilde{y}_{\textit{esc}}^{\frac{1}{2}}\left(\tilde{y}_{\textit{esc}}-\tilde{y}\right)^{1/2}\frac{\partial \tilde{y}_{\textit{esc}}}{\partial \tilde{x}}.\end{align}
Integrating (2.29) and considering the no-slip boundary condition, we obtain
\begin{equation}\tilde{\kappa }^{-1}\tilde{u}(\tilde{y}=\tilde{y}_{\textit{esc}})=-\frac{5}{9}\left(\tilde{y}_{\textit{esc}}\right)^{3}\tilde{\lambda }_{D}^{-2}\frac{\partial \tilde{j}_{1}}{\partial \tilde{x}}-\frac{4}{3}\tilde{\lambda }_{D}^{-2}\tilde{j}_{1}\left(\tilde{y}_{\textit{esc}}\right)^{2}\frac{\partial \tilde{y}_{\textit{esc}}}{\partial \tilde{x}}.\end{equation}
Considering that the local ESC layer satisfies the following relations (Rubinstein & Zaltzman Reference Rubinstein and Zaltzman2001):
\begin{equation}\frac{2\sqrt{2\tilde{j}_{1}}}{3\tilde{\lambda }_{D}}\tilde{y}_{\textit{esc}}^{3/2}=\tilde{\varphi }+\ln\left(\tilde{\lambda }_{D}\sqrt{\frac{2\tilde{j}_{1}}{\tilde{y}_{\textit{esc}}}}\right),\end{equation}
and multiplying both sides of (2.31) by the dimensionless Debye length and neglecting the lower-order terms, we obtain
\begin{align} \frac{2\sqrt{2\tilde{j}_{1}}}{3}\tilde{y}_{\textit{esc}}^{3/2}\approx -\tilde{\lambda }_{D}ln\sqrt{\tilde{y}_{\textit{esc}}} \end{align}.
Taking the derivative of (2.32), we obtain
\begin{equation}\frac{\sqrt{2}}{3}\frac{\partial \tilde{j}_{1}/\partial \tilde{x}}{\sqrt{\tilde{j}_{1}}}\boldsymbol{\cdot }\tilde{y}_{\textit{esc}}^{3/2}+\sqrt{2}\sqrt{\tilde{j}_{1}}\boldsymbol{\cdot }\frac{\partial \tilde{y}_{\textit{esc}}}{\partial \tilde{x}}\boldsymbol{\cdot }\sqrt{\tilde{y}_{\textit{esc}}}\approx -\frac{1}{2}\tilde{\lambda }_{D}\frac{\partial \tilde{y}_{\textit{esc}}/\partial \tilde{x}}{\tilde{y}_{\textit{esc}}}.\end{equation}
Dividing both sides of (2.33) by
$\tilde{y}_{\textit{esc}}^{3/2}/\sqrt{\tilde{j}_{1}}$
yields
\begin{equation}\frac{\partial \tilde{j}_{1}/\partial \tilde{x}}{\tilde{j}_{1}}+3\frac{\partial \tilde{y}_{\textit{esc}}/\partial \tilde{x}}{\tilde{y}_{\textit{esc}}}\approx -\frac{3}{2\sqrt{2}}\tilde{\lambda }_{D}\frac{{\tilde{j}_{1}}^{1/2}}{\tilde{y}_{\textit{esc}}^{5/2}}\frac{\partial \tilde{y}_{\textit{esc}}}{\partial \tilde{x}}.\end{equation}
Ignoring the lower-order terms in (2.34), we get
\begin{equation}\frac{\partial \tilde{j}_{1}/\partial \tilde{x}}{\tilde{j}_{1}}\approx -3\frac{\partial \tilde{y}_{\textit{esc}}/\partial \tilde{x}}{\tilde{y}_{\textit{esc}}}.\end{equation}
Taking the derivative of (2.32), we obtain
\begin{equation}\tilde{\kappa }^{-1}\tilde{u}(\tilde{y}=\tilde{y}_{\textit{esc}})=-\frac{1}{9}\tilde{\lambda }_{D}^{-2}\left(\tilde{y}_{\textit{esc}}\right)^{3}\frac{\partial \tilde{j}_{1}}{\partial \tilde{x}}.\end{equation}
Substituting (2.26) into (2.36), yields
\begin{equation}\tilde{u}(\tilde{y}=\tilde{y}_{\textit{esc}})=-8^{-1}\tilde{\kappa }({\tilde{\varphi }_{\textit{esc}}}-{\tilde{\varphi }_{\textit{pmi}}})^{2}\frac{\partial \tilde{j}_{1}}{\tilde{j}_{1}\partial \tilde{x}}.\end{equation}
Using (2.35), (2.37) simplifies to
\begin{equation}\tilde{u}(\tilde{y}=\tilde{y}_{\textit{esc}})\approx \frac{3}{8}\tilde{\kappa }({\tilde{\varphi }_{\textit{esc}}}-{\tilde{\varphi }_{\textit{pmi}}})^{2}\frac{\partial \tilde{y}_{\textit{esc}}/\partial \tilde{x}}{\tilde{y}_{\textit{esc}}}.\end{equation}
Noting the vortex velocity at the edge of the ESC, we use a first-order approximation (Kwak et al. Reference Kwak, Pham, Lim and Han2013),
$1/\partial \tilde{x}\sim 1/\partial \tilde{y}\sim (1/\tilde{y}_{\textit{esc}})\sim 1/(\tilde{\lambda }_{D}^{(2/3)}({\tilde{\varphi }_{\textit{esc}}}-{\tilde{\varphi }_{\textit{pmi}}})^{2/3})$
. This is because the size of the vortex on the surface of the metal column is comparable to the thickness of the ESC. This approximation comes from Davidson et al. (Reference Davidson, Andersen and Mani2014). Therefore, we can get
\begin{equation}\tilde{u}_{s}=\tilde{u}\left(\tilde{y}=\tilde{y}_{\textit{esc}}\right)\sim \tilde{\kappa }\tilde{\lambda }_{D}^{-\frac{2}{3}}\left(\tilde{\varphi }_{\textit{esc}}-\tilde{\varphi }_{\textit{pmi}}\right)^{\frac{4}{3}}.\end{equation}
We assume that
$\tilde{\varphi }_{\textit{esc}}$
satisfies a linear relation with the anode electric potential (Kwak et al. Reference Kwak, Pham, Lim and Han2013),
$\Delta \tilde{\varphi }\sim \tilde{\varphi }_{\textit{esc}}$
, and adopt the approximation
$\tilde{\varphi }_{\textit{esc}}-\tilde{\varphi }_{\textit{pmi}}\sim \Delta \tilde{\varphi }/\tilde{L}\boldsymbol{\cdot }\tilde{a}=\Delta \tilde{\varphi }/\tilde{L}$
(Davidson et al. Reference Davidson, Andersen and Mani2014), where
$\tilde{a}=1$
refers to the radius of the metal cylinder. The tangential velocity
$\tilde{u}_{s}$
of the vortex and the normal velocity
$\tilde{v}$
are approximately linearly related (Li et al. Reference Li, Townsend, Archer and Koch2022). This allows one to use the normal velocity
$\tilde{v}$
to analyse the tangential velocity
$\tilde{u}_{s}$
. We can obtain the scaling of the vortex velocity
\begin{equation}\tilde{v}\sim \tilde{\kappa }\tilde{\lambda }_{D}^{-\frac{2}{3}}\big(\Delta \tilde{\varphi }/\tilde{L}\big)^{\frac{4}{3}}\text{for the ESC}\hbox{-}\text{triggered ICEO}.\end{equation}
Scaling (2.40) represents the key result of this study, as it reveals the scaling law in ESC-triggered ICEO. Equations (2.11) and (2.40) describes the scaling behaviour of velocity. To distinguish between the critical points of transition, we observe that the vortex formation near the ESC edges signifies the onset of the velocity scaling transition, which can be predicted by linear instability theory. Rubinstein & Zaltzman (Reference Rubinstein and Zaltzman2000, Reference Rubinstein and Zaltzman2003) and Demekhin, Shelistov & Polyanskikh (Reference Demekhin, Shelistov and Polyanskikh2011) have established an implicit marginal stability relation that jointly regulates voltage, Debye length, hydrodynamic coupling coefficient and wavenumber
$k : f(k)=4({\sinh k\cosh k-k}/{\sinh k\cosh k+k-2k^{2}\coth k})=({\tilde{\kappa }(\Delta \tilde{\varphi }/\tilde{L})^{2}}/{8})$
. We further consider that the critical wavenumber
$k_{c}$
depends on the dimensionless Debye length
$\tilde{\lambda }_{D}$
, satisfying
$k_{c}\sim | ln (\tilde{\lambda }_{D})| /3$
(Rubinstein & Zaltzman Reference Rubinstein and Zaltzman2003). Hence, the corrected critical stability is given by
\begin{equation}f(k_{c})=\frac{\tilde{\kappa }\big(\Delta \tilde{\varphi }/\tilde{L}\big)^{2}}{8}.\end{equation}
When
${\tilde{\kappa }(\Delta \tilde{\varphi }/\tilde{L})^{2}}/{8}\gt f(k_{c})$
, the flow is unstable. When
${\tilde{\kappa }(\Delta \tilde{\varphi }/\tilde{L})^{2}}/{8}\lt f(k_{c})$
, the flow is stable. The scaling region with a 4/3 exponent corresponds to the strong-field regime (flow unstable), while the scaling region with an exponent of 2 corresponds to the weak-field regime.
3. Results and discussion
The parameter space considered is (Davidson et al. Reference Davidson, Andersen and Mani2014; Karatay et al. Reference Karatay, Andersen, Wessling and Mani2016; Kim et al. Reference Kim, Davidson and Mani2019; Mani & Wang Reference Mani and Wang2020):
$0\leq \tilde{\kappa }\leq 1$
,
$10^{-4}\leq \tilde{\lambda }_{D}\leq 10^{-2}$
and
$0\leq \Delta \tilde{\varphi }\leq 100$
. Dimensionless equations are solved numerically in a two-dimensional Cartesian coordinate system via the finite element method with extremely fine grids. For details of the DNS (finite element method), and grid, please refer to our previous work (Xu et al. Reference Xu, Gu, Liu, Huo, Zhou and S.2020; Liu, Zhou & Shi Reference Liu, Zhou and Shi2020a
).
Validation. (a) Ion concentrations. The inset shows the distribution results of cation and anion concentrations from Davidson et al. (Reference Davidson, Andersen and Mani2014), demonstrating that the simulation in this study aligns with the existing DNS. (b) Grid independence. The structure of the space-charge layer (obtained by subtracting the cation concentration from the anion concentration) is used to verify grid independence. Here, the parameters are
$\tilde{\lambda }_{D}=0.0014$
,
$\tilde{\kappa }=0.5$
and
$E_{0}a/V_{T}=30$
.
Strong EDL capacitive charging triggers flow instability. Experimental (Peng et al. Reference Peng, Lazo, Shiyanovskii and Lavrentovich2014) and simulation results. The DNS shows an instantaneous snapshot of the dimensionless salt concentration and streamlines.
$\Delta \tilde{\varphi}/\tilde{L}=16$
corresponding to
$\Delta \tilde{\varphi}=800$
. The simulation parameters are
$\Delta \tilde{\varphi}/\tilde{L}=16$
,
$\tilde{\lambda}_{D}=5\boldsymbol{\cdot }10^{-4}$
and
$\tilde{\kappa}=0.5$
.
The concentration curve serves as a key validation. Note that solving ICEO requires significantly refined meshes under strong electric fields, which is a computational challenge for ICEO. Specifically, a local super-fine mesh refinement is necessary within the EDL and the ESC to accurately capture the ESC induced by strong charge accumulation in the EDL. This aspect was often overlooked in previous ICEO simulations. For weak electric fields, a radial arrangement of O(1) mesh points within the EDL is typically sufficient. However, under strong electric fields, a large number of cations accumulate in the EDL, inducing the formation of the ESC. Consequently, O(10) mesh points are required within the EDL under strong-field conditions to capture the dynamics of ICEO effectively. The verification of grid independence ensures the consistency of simulation results across different mesh resolutions, thereby confirming the robustness and reliability of the model in numerical computations. Figure 1(a) presents a comparison with existing literature, demonstrating that the current DNS results align well with those reported in previous studies. Additionally, figure 1(b) investigates the impact of the minimum mesh size within the EDL on the space-charge layer, revealing that capturing EDL charging behaviour requires a fine mesh resolution. A minimum radial grid size of 0.02
$\tilde{\lambda }_{D}$
within the EDL effectively captures ICEO spanning a wide range of voltages.
3.1. Space-charge layer transition
Figure 2 shows a snapshot of salt concentration and streamlines, where EDL capacitive charging leads to the ESC formation. The low-concentration area along the bulk zone is similar to capacitive deionisation (Rubin et al. Reference Rubin, Suss, Biesheuvel and Bercovici2016; Oyarzun et al. Reference Oyarzun, Hemmatifar, Palko, Stadermann and Santiago2018; Ramachandran et al. Reference Ramachandran, Hawks, Stadermann and Santiago2018). Additionally, it is noteworthy that the net flux in this system is maintained by the EDL capacitance and the salt jet. Ions are stored within the EDL, and vortices drive cations out of the EDL, forming enrichment jets at 0° and 180°. When the electric field volumetric forces within the ESC overcome the viscous forces, local small-scale vortex structures are formed, disrupting the symmetrical flow pattern. The background flow (comprising two pairs of vortices) pushes the small-scale vortices to detach at 0° and 180°, resulting in a four-way vortex shedding phenomenon around the cylinder and a localised salt jet.
The formation mechanism of ESC. (a) The schematic of EDL strong capacitive charging illustrates the subtle processes involved in EDL charging development at the metal–cylinder interface. It highlights the contribution of the electric field to the arrangement of ions. (b) ion concentration, and (c) electric field strength at the cut line (dashed line in figure 2). The EDL structure transitions to an ESC structure, in which a small number of cations are present, while the anion concentration is nearly zero. (d) Ion current induced by EDL capacitive charging. This time current is computed spatially by the radial cation flux on the upper surface of the metal cylinder minus the anion flux. (e) Variation of the cation and anion fluxes with the external electric field. The upper surface EDL exhibits cation selectivity. This ion flux is computed temporally and spatially by the radial cation flux on the upper surface of the metal cylinder minus the anion flux. In (b)–(c), the DNS parameters are
$\Delta \tilde{\varphi }/\tilde{L}=16$
,
$\tilde{\lambda }_{D}=5 \times 10^{-4}$
and
$\tilde{\kappa }=0.5$
.
The EDL capacitive charging triggers scaling transition of velocity. (a) Fully developed concentration field and streamline distribution. Here,
$A_{1}$
and
$A_{2}$
are the velocities at two different points in the inflow, and the positions of
$A_{1}$
and
$A_{2}$
are
$\tilde{y}=\sqrt{3}$
and
$\tilde{y}=\tilde{\delta }_{\mathrm{esc}}$
+ 1, respectively. Also,
$B_{1}$
and
$B_{2}$
are the velocities at two different points in the outflow, and the positions of
$B_{1}$
and
$B_{2}$
are
$\tilde{x}=-\sqrt{3}$
and
$\tilde{x}=-\tilde{\delta }_{\mathrm{esc}}-1$
, respectively. Here, we choose
$\tilde{\delta }_{\mathrm{esc}}=0.05$
. The positions of
$A_{1}$
and
$B_{1}$
are consistent with those in Davidson et al. (Reference Davidson, Andersen and Mani2014). (b) Velocity scaling. To highlight the influence of vortices, the statistics points (
$A_{2}$
and
$B_{2}$
) are considered. The DNS results in this paper are basically consistent with those in Davidson et al. (Reference Davidson, Andersen and Mani2014). The simulation parameters are
$\Delta \tilde{\varphi }/\tilde{L}=30$
,
$\tilde{\lambda }_{D}=0.0014$
and
$\tilde{\kappa }=0.5$
.
Figure 3(a) schematically illustrates the formation of the ESC layer induced by strong capacitive charging of the EDL near an ideally polarisable interface. When a strong electric field is applied to a metal cylinder, surface conduction generates opposite charges on its upper and lower surfaces, which attract counter-ions and establish EDLs of opposite polarity. The applied field further enhances ion accumulation within the EDLs, leading to strong capacitive charging current. This charging drives cations to penetrate beyond the Debye layer into the diffusion region, thereby giving rise to the ESC and the emergence of local vortices.
To highlight the structural changes of the charge layer, figure 3(b–c) shows the increase in cation concentration within the EDL due to capacitive charging, followed by the transition from EDL to ESC. This transition is driven by the ionic current induced by capacitive charging, see figure 3(d), which produces a strong negative electric field on the upper surface. The resulting field pushes cations into the diffuse layer while excluding anions, creating a localised region of high cation concentration. Consequently, an ESC layer emerges between the EDL and the diffuse layer, characterised by a finite charge density. As shown in figure 1(e), the EDL thus acquires an ion-selective transport property: on the upper surface, the cation current increases with the applied field while the anion flux remains negligible, whereas by symmetry the lower surface exhibits anion-selective transport.
3.2. Velocity scaling transition
We focus on the statistical characteristics of the radial velocity at the inlet and the outlet locations, as shown in figure 4(a). Here, we observe a transition in velocity scaling. Under weak electric fields, the flow velocity follows the canonical quadratic scaling predicted by the asymptotic model. However, as the field strength increases, vortices and ESC structures emerge due to strong capacitive charging, causing the velocity to deviate from the canonical prediction (figure 2 b). This transition reflects a fundamental change in the charge layer. In the EDL-dominated regime, the quadratic scaling arises because the slip velocity is determined solely by the induced potential, with no concentration depletion and negligible charging current. By contrast, in the ESC regime, ion depletion and a finite charging current govern the dynamics, giving rise to the 4/3 exponent. Therefore, in our derivation, we emphasise the contribution of ion flux, see (2.15). The absence of saturation can be attributed to the persistent strong local electric field maintained by the ESC, which continues to drive the vortex. At still higher applied fields, additional processes such as electrochemical reactions are expected to occur, but these lie beyond the present scope and will be addressed in future work.
To eliminate the influence of electric field frequency,
$S(\omega )$
was considered (Pascall & Squires Reference Pascall and Squires2010). By normalising the velocity by
$S(\omega )$
, the velocities across different frequencies can be unified. Therefore, the velocity scaling can be corrected
\begin{equation}\tilde{v}=\alpha \tilde{\kappa }\tilde{\lambda }_{D}^{-2/3}\big(\Delta \tilde{\varphi }/\tilde{L}\big)^{4/3}(1+{\omega ^{2}}\tau _{c}^{2})^{-1}\;\text{for the ESC}\hbox{-}\text{triggered ICEO}.\end{equation}
In the parameter space in this paper, the fitting coefficient is
$\alpha \approx -0.0007$
. In scaling (61),
$\tilde{\kappa }$
describes the contribution of the electrolyte properties,
$\tilde{\lambda }_{D}$
describes the contribution of the space-charge layer, and
$\Delta \tilde{\varphi }/\tilde{L}$
is the contribution of the electric field. Here,
$S(\omega )=1/(1+\omega ^{2}\tau _{c}^{2})$
is the contribution of the electric field frequency. Using scaling, all DNS data and experimental points collapse approximately on the theoretical straight line, as shown in figure 5.
Velocity scaling in a broad parameter space. (a) The effect of dimensionless Debye length, hydrodynamic coupling coefficient and electric field on the velocity. Under the action of strong electric field, the velocity satisfies the 4/3 power law. (b) Universal scaling laws of velocity. The experimental setup described in Peng et al. (Reference Peng, Lazo, Shiyanovskii and Lavrentovich2014) involves an Au sphere with a diameter of
$50\,{\unicode{x03BC}} \mathrm{m}$
, placed between two metal plates spaced
$10\,\mathrm{mm}$
apart. The solution concentration is
$0.1\mathrm{mM}$
, corresponding to a dimensionless Debye length of
$\tilde{\lambda }_{D}\approx 0.005$
. To normalise the experimental data, we adopt the dimensionless reference values used in this study: a reference length of
$10\,{\unicode{x03BC}} \mathrm{m}$
, a reference voltage of
$25.6\mathrm{mV}$
and a diffusion velocity of
$2 \times 10^{-4}\,\mathrm{ms}^{-1}$
. In the scaling factor
$S(\omega )=1/(1+\omega ^{2}\tau _{c}^{2})$
, where
$\omega$
is the electric field frequency which is set to
$1\mathrm{kHZ}$
in the experiment. The characteristic resistance–capacitance charging time
$\tau _{c}$
can be calculated using
$\tau _{c}=\varepsilon _{w}(\pi | x_{e}| )\sigma \lambda _{D}$
(Pascall & Squires Reference Pascall and Squires2010), with
$x_{e}$
being the position where velocity is statistically measured–
$20\,{\unicode{x03BC}} \mathrm{m}$
in the experiment (Peng et al. Reference Peng, Lazo, Shiyanovskii and Lavrentovich2014). Here,
$\sigma$
is the conductivity of the aqueous solution, and
$\lambda _{D}$
is the Debye length. Note that we do not consider the non-monotonic effects of frequency on velocity, as studied in Davidson et al. (Reference Davidson, Andersen and Mani2014) and Peng et al. (Reference Peng, Lazo, Shiyanovskii and Lavrentovich2014).
The scaling law presented in (3.1) deviates fundamentally from the predictions of standard asymptotic models (Bazant & Squires Reference Bazant and Squires2004). In the regime of EDL, the flow near a polarisable interface is governed by a balance between the electrical body force within the EDL and viscous stresses, resulting in a quadratic dependence of velocity on the applied electric field. However, as the external electric field is increased, strong charging of the EDL near the polarisable surface generates intense local ionic currents, which in turn induce the formation of an ESC layer. In this regime, the potential gradients are no longer confined within the thin EDL, but act over a much broader ESC region. Consequently, the dominant force balance shifts to one between the electrical body force in the ESC and viscous stresses. This reorganisation of electrohydrodynamic force balance leads to a velocity scaling transition. Besides, unlike previous studies involving externally driven flows and ion-selective membranes, the present work investigates an ICEO system that is entirely free of externally imposed flow disturbances and selective membranes. The flow here is spontaneously driven by nonlinear EDL charging alone, underscoring a fundamental difference in the underlying physics compared with shear flow systems (Kwak et al. Reference Kwak, Pham, Lim and Han2013; Kang & Kwak Reference Kang and Kwak2020; Liu et al. Reference Liu, Zhou and Shi2020a ).
3.3. Flow pattern
In figure 6(a–c), we plot the flow patterns in a broad parameter space. It is found that the dimensionless Debye length and the hydrodynamic coupling coefficient have a significant impact on the flow pattern, rather than the electric field, which is common sense. As the dimensionless Debye length decreases, small-scale vortices appear near the cylinder, leading to an irregular concentration boundary layer structure. Under the influence of the background flow, the flow becomes increasingly unstable and is accompanied by enhanced velocity. Physically, decreasing the dimensionless Debye length reduces the thickness of the ESC layer and enhances the local electrostatic driving force, leading to the burst of small-scale vortices. The background flow drives these small-scale vortices to roll along the surface of the metal cylinder and shed at 0° and 180°. At 0° and 180°, the vortices are compressed by the background flow (comprising two pairs of vortices), resulting in the salt jets. Meanwhile, vortices located away from 0° and 180° continue to maintain an approximately elliptical structure.
Flow pattern. Effect of electric field on (a) flow pattern and (d) velocity distribution at the cutline. Here, the parameters are set to
$\tilde{\lambda }_{D}=0.005$
and
$\tilde{\kappa }=0.5$
. Effect of the hydrodynamic coupling coefficient on (b) flow pattern and (e) velocity distribution. Here, the parameters are set to
$\Delta \tilde{\varphi }/\tilde{L}=20$
and
$\tilde{\lambda }_{D}=10^{-3}$
. Effect of the dimensionless Debye layer thickness on (c) flow pattern and (f) velocity distribution. In the diffuse layer, the radial velocity and distance satisfy the decay power law of −3. Here,
$\Delta \tilde{\varphi }/\tilde{L}=30$
and
$\tilde{\kappa }=0.5$
.
Under different hydrodynamic coupling coefficients, we observe a third pathway that enhances flow instability. Increasing the hydrodynamic coupling coefficient strengthens the coupling between fluid and ion transport, leading to the vortices bursting. This bursting pattern causes the flow field structure to differ from the effects of adjusting the dimensionless Debye length and voltage. As the hydrodynamic coupling coefficient increases, vortices transition from semi-elliptical structures to elongated band-like structures, forming long salt jets at 0° and 180°.
In figure 6(d–f), the radial velocity distribution spanning a broad parameter space of dimensionless parameters is plotted, and extracted from the cutline (white dashed lines) in figure 6(a–c). It was found that the radial velocity distribution in the diffuse layer conforms to the predictions of the standard asymptotic model, with the radial velocity decaying according to a power law of −3 with distance. However, in the ESC layer, due to the influence of local small-scale vortices, the velocity field single-peak structure transitions to a double-peak structure, and even to a triple-peak structure.
3.4. Flow state transition
Finally, we reveal how different dimensionless Debye lengths, electric fields and hydrodynamic coupling coefficients jointly regulate the flow state and identify the flow state shown in the Supplementary Materials. Since unstable flow can drive spatio-temporal fluctuations in the salt concentration field, the flow regime is identified by examining the space–time evolution of salt concentration, as illustrated in figure 7. As reported in the literature (Davidson et al. Reference Davidson, Andersen and Mani2014), the concentration field and the velocity field exhibit a strong positive correlation in microfluidic systems. Therefore, analysing the concentration field provides a reliable proxy for determining the flow state. Specifically, we define the stable four-vortex configuration, accompanied by a steady salt-jet structure, as the stable flow state. However, when the four-vortex structure is accompanied by localised small-scale vortices, the system enters an unstable flow state, characterised by plume-like, unsteady salt-jet structures; see figure 7. Using DNS, we can get
$k_{c}(\tilde{\lambda }_{D})=| ln (\tilde{\lambda }_{D})| /3-1.36$
. The critical stability equation of the flow transition can be expressed as
The spatio-temporal evolution of salt concentration (
$\tilde{c}_{1}$
+
$\tilde{c}_{2}$
)/2. At the cutline (
$\tilde{x}=1.1$
,
$-0.5\leq \tilde{y}\leq 0.5$
), the spatio-temporal evolution of salt concentration is shown in (a) case 1 (
$\tilde{\lambda}_{D}=0.0005$
,
$\Delta \tilde{\varphi}/\tilde{L}=16$
and
$\tilde{\kappa}=0.5$
) and (b) case 2 (
$\tilde{\lambda}_{D}=0.0014$
,
$\Delta \tilde{\varphi}/\tilde{L}=22$
and
$\tilde{\kappa}=0.5$
). (c) The velocity fluctuation (
$\tilde{u}^{\prime}$
) curves. For case 1, the flow is unstable, which is consistent with our theoretical analysis. Local small-scale vortices erupt, forming a fascinating feather texture structure and irregular fluctuations in velocity. However, for case 2, the flow is stable, forming a long band-like salt-jet structure. Case 1 and case 2 correspond to point 1 (unstable) and point 2 (stable) in figure 8, respectively.
The EDL capacitive charging triggers flow state transition. When the dimensionless parameters exceed the critical stability curve, the flow is in an unstable state; otherwise, it remains stable. The flow state is identified through the spatio-temporal evolution of salt concentration. Physically, a decrease in
$\tilde{\lambda }_{D}$
corresponds to the volumetric force of the ESC layer becoming stronger, triggering small-scale vortex bursts.
\begin{equation}4\frac{\sinh k_{c}(\tilde{\lambda }_{D})\boldsymbol{\cdot }\cosh k_{c}(\tilde{\lambda }_{D})-k_{c}(\tilde{\lambda }_{D})}{\sinh k_{c}(\tilde{\lambda }_{D})\boldsymbol{\cdot }\cosh k_{c}(\tilde{\lambda }_{D})+k_{c}(\tilde{\lambda }_{D})-2k_{c}^{2}(\tilde{\lambda }_{D})\boldsymbol{\cdot }\coth k_{c}(\tilde{\lambda }_{D})}=\tilde{\kappa }\big(\Delta \tilde{\varphi }/\tilde{L}\big)^{2}/8.\end{equation}
In figure 8, the critical stability curve predicts the influence of the dimensionless parameter space on the flow. As
$\tilde{\lambda }_{D}$
increases, the values of
$\tilde{\kappa }$
and
$\Delta \tilde{\varphi }/\tilde{L}$
required for flow instability also increase. For
$\tilde{\lambda }_{D}=0.0014$
and
$\tilde{\kappa }=0.5$
, the critical electric field is approximately
$(\Delta \tilde{\varphi}/\tilde{L})\approx 24$
, which is in good agreement with the results reported in Davidson et al. (Reference Davidson, Andersen and Mani2014). Additionally, we conducted further calculations to discuss scenarios where
$\tilde{\lambda }_{D}$
takes on larger (
$\tilde{\lambda }_{D}=0.01$
) and smaller (
$\tilde{\lambda }_{D}=0.0002$
) values. When
$\tilde{\lambda }_{D}=0.0002$
, the critical
$\Delta \tilde{\varphi }/\tilde{L}$
required for flow instability is approximately greater than 12, which is consistent with our DNS.
4. Conclusions
This study employs direct simulations and theoretical analysis to investigate capacitive charging, ESC, scaling and instability near an ideally polarisable metal cylinder, elucidating the mechanisms by which the velocity deviates from predictions of the canonical asymptotic model for a strongly polarised scenario. The key conclusions are as follows.
Mechanism of the strong capacitive charging triggers instability transition and ESC formation. With increasing external voltage, the amount of charge stored in the EDL rises, leading to significant capacitive charging and the generation of ionic currents, which induce a strong electric field. Under the influence of this strong electric field, cations extend over a finite region, forming the ESC layer with a small amount of cations.
The velocity of ICEO can be described by a power-law combination of the dimensionless Debye length, the hydrodynamic coupling constant and the external electric potential difference. This scaling is supported by both DNS and experimental data, establishing the transitional behaviour from weak to strong fields. It identifies the impact of the ESC layer on the velocity scaling regime and extends the canonical asymptotic model.
A phase diagram of the flow state is presented, constructed from the dimensionless Debye length, the hydrodynamic coupling constant and the external electric field. The flow stability and instability are described by the critical stability relation. Our DNSs confirm this analysis, offering a theoretical perspective for stability control.
In conclusion, this study reveals the transition behaviour of ICEO through DNS and theoretical analysis, clarifying the critical role of EDL capacitive charging in the ICEO and departing from both the classical asymptotic theories (Bazant & Squires Reference Bazant and Squires2004; Squires & Bazant Reference Squires and Bazant2004) and the shear sheltering (Kang & Kwak Reference Kang and Kwak2020; Liu et al. Reference Liu, Zhou and Shi2020b , Reference Liu, Zhou and Shi2022). Our approach can be extended to consider the analysis of ICEO in diffusiophoresis (Gupta, Rallabandi & Stone Reference Gupta, Rallabandi and Stone2019) and polyelectrolyte solutions (Chen, Perazzo & Stone Reference Chen, Perazzo and Stone2020).
Funding
We thank the four reviewers for their valuable comments on the draft. The authors express their sincere gratitude for the financial support provided by the NSF of China (No. 12272207, 12432013), China Postdoctoral Science Foundation (No. 2023M741970) and the Postdoctoral Fellowship Program of CPSF (No. GZB20230361).
Declaration of interests
The authors declare no conflict of interest.
