1. Introduction
Turbulence is ubiquitous in both natural and engineering systems, and its control remains a fundamental challenge in turbulence. This study examines Rayleigh–Bénard convection (RBC), a classical model for investigating thermal convection (Ahlers, Grossmann & Lohse Reference Ahlers, Grossmann and Lohse2009; Lohse & Shishkina Reference Lohse and Shishkina2024). Rayleigh–Bénard convection occurs in a fluid layer confined between two horizontal plates, where the lower plate is heated and the upper plate is cooled. The resulting temperature difference generates convection driven by density variations, providing a widely used framework for exploring the fundamental physics of natural convection (Lohse & Xia Reference Lohse and Xia2010; Chillà & Schumacher Reference Chillà and Schumacher2012; Sun & Zhou Reference Sun and Zhou2014; Shishkina Reference Shishkina2021).
Early studies on thermally driven RBC primarily examined the collective influence of the Rayleigh number, Prandtl number and geometric configuration on heat transfer efficiency and flow structure, particularly in terms of scaling relations (Ahlers et al. Reference Ahlers, Grossmann and Lohse2009). A key feature of RBC is large-scale circulation (LSC) (see figure 1 a), a dominant convective loop driven by the temperature difference between the top and bottom plates. The LSC’s flow direction and intensity are influenced by system parameters and external perturbations (Brown & Ahlers Reference Brown and Ahlers2006; Zhou, Sun & Xia Reference Zhou, Sun and Xia2007). In recent years, studies have increasingly focused on the transition of LSC from steady-state heat transfer and flow characteristics to unsteady behaviours, particularly phenomena such as LSC cessation and reversal (Benzi Reference Benzi2005; Assaf, Angheluta & Goldenfeld Reference Assaf, Angheluta and Goldenfeld2011; Vasil’ev & Frick, Reference Vasil’ev and Frick2011; Wagner & Shishkina Reference Wagner and Shishkina2013; Ni, Huang & Xia Reference Ni, Huang and Xia2015; Xia et al. Reference Xia, Wan, Liu, Wang and Sun2016). These dynamic behaviours not only reveal a complex nonlinear dynamics within the system but also relate to natural phenomena such as circulation reversals in the Earth’s atmosphere and oceans (Ahlers & Xu Reference Ahlers and Xu2001; Xi & Xia Reference Xi and Xia2007; Sugiyama et al. Reference Sugiyama, Ni, Stevens, Chan, Zhou, Xi, Sun, Grossmann, Xia and Lohse2010). Flow reversal refers to a phenomenon where the dominant flow direction undergoes a sudden or gradual reversal, typically accompanied by energy redistribution and flow reorganisation. Flow reversal can generally be classified into two categories (Bahr et al. Reference Bahr, Gerbatsch, Bahr and Weiss1995; Cioni, Ciliberto & Sommeria Reference Cioni, Ciliberto and Sommeria1997; Sun, Xi & Xia Reference Sun, Xi and Xia2005). The first is reversal driven by azimuthal rotation; in cylindrical geometries, due to azimuthal symmetry, the LSC slowly rotates around the vertical axis. When the rotation angle reaches 180°, the circulation direction reverses. This process typically occurs under high Rayleigh number conditions and is influenced by system noise and external perturbations. The second is reversal induced by cessation; in this case, the momentum of the original flow gradually decays until it vanishes completely, after which the flow direction re-establishes due to random perturbations or system reorganisation.
Schematic of LSC in a square-cavity cell: (a) LSC in thermally driven RBC; (b) LSC in solutally driven RBC within an electrochemical system under a negligible EC scenario; (c) EC vortex structure. As EC intensity increases, small-scale EC vortices develop near the membrane–solution surface, gradually disrupting the LSC structure. (d) Ion concentration distribution near the bottom membrane–solution surface. Here,
$c_{+}$
and
$c_{-}$
represent the concentrations of cations and anions, respectively. The inset illustrates the ESC layer, which exhibits a peak structure with space charge density (
$c_{+}-c_{-}$
).
The effects of boundary conditions, cell geometry and fluid properties on flow reversal have been extensively studied. Huang et al. (Reference Huang, Wang, Xi and Xia2015) found that applying a fixed heat flux boundary condition to the bottom plate of a convective cell significantly reduces the frequency of flow reversal, highlighting the critical role of boundary conditions in regulating reversal behaviour. Zhang et al. (Reference Zhang, Xia, Zhou and Chen2020) demonstrated that local constant temperature control on the sidewall can effectively suppress or enhance flow reversal in two-dimensional RBC. In studies on non-Oberbeck–Boussinesq effects, Xia et al. (Reference Xia, Wan, Liu, Wang and Sun2016) discovered that the growth of corner vortices can induce flow reversal, revealing a new dynamic process distinct from traditional mechanisms. Wang et al. (Reference Wang, Lai, Song and Tong2018) found that tilting the convection cell suppresses flow reversal characterised by a single main vortex, whereas increasing the aspect ratio of the tilt promotes reversal associated with a double main vortex. Chen, Wang & Xi (Reference Chen, Wang and Xi2020) conducted comparative experiments on the flow reversal behaviour of LSC in right cylindrical and cornerless units, showing that flow reversal can still occur in cornerless units, with the reversal frequency closely linked to Rayleigh number variations. Beyond geometry and boundary conditions, other strategies have been explored to control LSC. Emran & Shishkina (Reference Emran and Shishkina2020) investigated the impact of wall roughness on the LSC dynamics. Sterl, Li & Zhong (Reference Sterl, Li and Zhong2016) found that external periodic rotational forcing modulates flow reversal. Hu, Zhang & Xia (Reference Hu, Zhang and Xia2024) demonstrated that adjusting the relative position of partial isothermal plates can either accelerate or completely stabilise flow reversal. Additionally, alternative approaches, such as applying magnetic fields, have been explored for controlling the flow dynamics (Zurner et al. Reference Zurner, Schindler, Vogt, Eckert and Schumacher2020).
However, these studies have primarily focused on classical thermally driven RBC. Recent research (de Valença et al. Reference de Valença, Kurniawan, Wagterveld, Wood and Lammertink2017) has revealed that RBC can also occur in small-scale electrochemical systems, such as those with ion-selective membranes separated by a 2 mm gap. Under an external vertical electric field, selective salt transport near a perfectly ion-selective membrane creates a non-neutral concentration distribution (see black regions in figures 1 b and 1 c), known as the extended space charge (ESC) layer (Liu, Zhou & Shi Reference Liu, Zhou and Shi2020; Rubinstein & Zaltzman Reference Rubinstein and Zaltzman2000; Zaltzman & Rubinstein Reference Zaltzman and Rubinstein2007). In the ESC layer, the cation concentration is high while the anion concentration is nearly zero, resulting in a localised hump-like structure in the space charge density (see figure 1 d). The electric field volumetric force within the ESC layer can overcome viscous forces, leading to the formation of vortex structures near the membrane–solution surface, known as electroconvection (EC) (Li, Archer & Koch Reference Li, Archer and Koch2019; Mani & Wang Reference Mani and Wang2020; Chen et al. Reference Chen, Zhu, Wang and Liu2024). The ESC layer generates a fluid density gradient related to salt concentration. When buoyancy forces are sufficiently large, the interaction between electric and buoyancy forces drives the coupled transport of solutally driven RBC and EC. This process is crucial for optimising ion transport in electrochemical systems, as it frequently occurs in practical applications such as electroplating (Debruyn Reference Debruyn1995; Marshall et al. Reference Marshall, Mocskos, González, Dengra, Molina and Iemmi2006), electrodeposition (Huth et al. Reference Huth, Swinney, McCormick and Kuhn1995; Mocskos et al. Reference Mocskos, González, Molina and Marshall2011) and electrodialysis (de Valença et al. Reference de Valença, Wagterveld, Lammertink and Tsai2015). As the external voltage increases, the ESC layer induces a strong electric field, generating EC near the bottom membrane surface (Rubinstein et al. Reference Rubinstein, Manukyan, Staicu, Rubinstein, Zaltzman, Lammertink, Mugele and Wessling2008; Chang, Demekhin & Shelistov Reference Chang, Demekhin and Shelistov2012; Kwak et al. Reference Kwak, Pham, Lim and Han2013; Demekhin, Nikitin & Shelistov Reference Demekhin, Nikitin and Shelistov2014; Kang & Kwak Reference Kang and Kwak2020). When the electrochemical device is placed vertically, the directions of the electric field ( E ) and gravitational field ( g ) align, which can simultaneously induce RBC and EC. This indicates the possibility of RBC–EC coupling. Karatay et al. (Reference Karatay, Andersen, Wessling and Mani2016) were the first to numerically predict the coupling effects between RBC and EC, demonstrating that this interaction significantly influences the flow structure. Later, de Valença et al. (Reference de Valença, Kurniawan, Wagterveld, Wood and Lammertink2017) experimentally confirmed this coupling effect, showing that RBC stretches EC, forming a salt-finger structure.
Based on the aforementioned discussion, evidently, in electrochemical systems, EC and RBC are inevitably coupled. This interaction not only affects mass transport but also influences flow behaviour. In RBC (Huang et al. Reference Huang, Wang, Xi and Xia2015; Sterl et al. Reference Sterl, Li and Zhong2016; Wang et al. Reference Wang, Lai, Song and Tong2018; Chen et al. Reference Chen, Wang and Xi2020; Emran & Shishkina Reference Emran and Shishkina2020; Zhang et al. Reference Zhang, Xia, Zhou and Chen2020; Zurner et al. Reference Zurner, Schindler, Vogt, Eckert and Schumacher2020; Hu et al. Reference Hu, Zhang and Xia2024), scalar field transport, specifically the temperature and concentration fields, plays a crucial role, raising the question of whether these processes contribute to the formation of LSC structures and flow reversal in solutally driven RBC. Thus, investigating the coupled transport of EC and RBC not only deepens our understanding of novel mechanisms for controlling LSC and flow reversal in electrochemical systems but also provides a theoretical foundation for optimising the design of electrochemical devices (Mühlenhoff et al. Reference Mühlenhoff, Eckert, Heinze and Uhlemann2007).
Herein, we performed direct numerical simulations (DNS) of solutal convection in a two-dimensional electrochemical cell bounded on both horizontal sides by a cation-selective membrane. For the first time, we reproduced the flow reversal of LSC observed in thermally driven RBC within an electrochemical system and explored the flow pattern (see figure 1 b,c)). By adjusting the voltage and electrolyte properties, the EC can be enhanced, resulting in the disappearance of the LSC reversal phenomenon, which is similar to the stabilisation of the LSC (figure 1 c). Unlike conventional approaches (Huang et al. Reference Huang, Wang, Xi and Xia2015; Xia et al. Reference Xia, Wan, Liu, Wang and Sun2016; Zhang et al. Reference Zhang, Xia, Zhou and Chen2020; Zurner et al. Reference Zurner, Schindler, Vogt, Eckert and Schumacher2020) that achieve LSC stabilisation, our study reveals a mechanism based on the coupling between EC and RBC. These EC vortices interact dynamically with the LSC and eliminate reversals via time-dependent forces within the ESC, without requiring any geometric or material asymmetry. This mechanism can be viewed as a form of flow–flow interaction, where one flow controls the behaviour of another.
2. Theory
This study examined the geometric configuration shown in figure 1(b), comprising a binary electrolyte layer confined in a square cavity, bounded by two perfectly cation-selective membranes and two electrically neutral sidewalls. Both the gravitational and electric fields were aligned downward. A square cavity was chosen owing to its relevance in electrochemical device design and its frequent use in studies on flow reversal in thermally driven RBC (Mühlenhoff et al. Reference Mühlenhoff, Eckert, Heinze and Uhlemann2007). The computational domain measured 1 × 1.
2.1. Governing equations
The solutal convection in electrochemical systems involves the spatio-temporal evolution of electric potential, charge density and fluid motion. Thus, it can be described by the Poisson–Nernst–Planck and Navier–Stokes equations. For simplicity, this study considers the effects of a monovalent (
$\tilde{z}_{+}=-\tilde{z}_{-}=1$
) salt solution. Assuming the fluid flow of salt solution is incompressible and applying the Oberbeck–Boussinesq approximation, the flow is governed by the incompressible Navier–Stokes and continuity equations. Although Joule dissipation plays a crucial role in thermally driven RBC, significantly affecting flow stability and the formation of LSC, in this study, we neglect the effects of Joule dissipation. This simplification is commonly adopted in electrochemical studies (Karatay et al. Reference Karatay, Andersen, Wessling and Mani2016). Thus the coupled governing equations are as follows:
\begin{align}&\qquad \frac{\partial \tilde{c}_{\pm }}{\partial \tilde{t}}=-\tilde{\boldsymbol{\nabla }}\boldsymbol{\cdot }\left(-\tilde{D}_{\pm }\tilde{\boldsymbol{\nabla }}\tilde{c}_{\pm }-\frac{\tilde{z}_{\pm }\tilde{D}_{\pm }\tilde{F}}{\tilde{R}\tilde{T}}\tilde{c}_{\pm }\tilde{\boldsymbol{\nabla }}\tilde{V}+\tilde{\boldsymbol{U}}\tilde{c}_{\pm }\right)\! , \end{align}
\begin{align}&\qquad\qquad\qquad\qquad\qquad \tilde{\varepsilon }_{0}\tilde{\varepsilon }_{\mathrm{r}}\tilde{{\nabla} }^{2}\tilde{V}=-\tilde{\rho }_{e} , \end{align}
\begin{align}& \tilde{\rho }_{0}\frac{\partial \tilde{\boldsymbol{U}}}{\partial \tilde{t}}+\tilde{\rho }_{0}\tilde{\boldsymbol{U}}\boldsymbol{\cdot }\tilde{\boldsymbol{\nabla }}\tilde{\boldsymbol{U}}=-\tilde{\boldsymbol{\nabla }}\tilde{P}+\tilde{\eta }\tilde{{\nabla} }^{2}\tilde{\boldsymbol{U}}+\tilde{\rho }_{e}\tilde{\boldsymbol{E}}+\tilde{\rho }_{0}\boldsymbol{\cdot }\tilde{\beta }\boldsymbol{\cdot }\Delta \tilde{c}\boldsymbol{\cdot }\tilde{\boldsymbol{g}} , \end{align}
\begin{align}&\qquad\qquad\qquad\qquad\qquad\quad \tilde{\boldsymbol{\nabla }}\boldsymbol{\cdot }\tilde{\boldsymbol{U}}=0 . \end{align}
In the Nernst–Planck equation,
$\tilde{t}$
represents the time,
$\tilde{F}$
is the Faraday constant,
$\tilde{R}$
is the gas constant,
$\tilde{T}$
is the temperature,
$\tilde{V}$
is the electric potential,
$\tilde{\boldsymbol{U}}$
is the fluid velocity and
$\tilde{c}_{\pm }$
are the concentrations of cations and anions, respectively. Also,
$\tilde{D}_{\pm }$
is the diffusion coefficients of the cations and anions. In (2.1),
$-\tilde{D}_{\pm }\tilde{\boldsymbol{\nabla }}\tilde{c}_{\pm }$
represents the ion diffusion flux,
$-\tilde{c}_{\pm }\tilde{z}_{\pm }\tilde{D}_{\pm }\tilde{F}\tilde{\boldsymbol{\nabla }}\tilde{V}/(\tilde{R}\tilde{T})$
represents the electrophoretic flux and
$\tilde{\boldsymbol{U}}\tilde{c}_{\pm }$
represents the convective flux. Equation (2.1) describes the coupling between diffusion, electric fields and fluid flow. In the Poisson equation,
$\tilde{\rho }_{e}=\tilde{F}(\tilde{c}_{+}-\tilde{c}_{-})$
represents the space charge density,
$\tilde{\varepsilon }_{0}$
is the vacuum permittivity and
$\tilde{\varepsilon }_{\mathrm{r}}$
is the relative permittivity of the salt solution. Local variations in ion concentration directly affect the electric potential distribution.
In the Navier–Stokes equation,
$\tilde{\eta }$
is the dynamic viscosity,
$\tilde{P}$
is the pressure,
$\tilde{\boldsymbol{E }}=-\tilde{\boldsymbol{\nabla }}\tilde{V}$
is the electric field and
$\tilde{\rho }_{e}\tilde{\boldsymbol{E}}$
is the electric body force. The buoyancy force
$\tilde{\rho }_{0}\tilde{\beta }\Delta \tilde{c}\tilde{\boldsymbol{g}}$
arises from the density variation of the salt solution near the membrane surface, where
$\tilde{\boldsymbol{g}}$
is the gravitational acceleration, directed along the negative normal direction. Under the Oberbeck–Boussinesq approximation, it is assumed that the fluid density variation is assumed to be a linear function (
$\tilde{\rho }-\tilde{\rho }_{0}=\tilde{\rho }_{0}\tilde{\beta }\Delta \tilde{c}$
) of the average salt concentration (Shi & Liu Reference Shi and Liu2018). where
$\tilde{\rho }$
is the fluid density,
$\tilde{\rho }_{0}$
is the reference density at the initial concentration,
$\Delta \tilde{c}$
is the change in salt concentration relative to the initial concentration and
$\tilde{\beta }$
is the volumetric expansion coefficient. In this study, we consider a salt solution, where the volumetric expansion coefficient is
$\tilde{\beta }\sim O(10^{-3})[1\mathrm{mM}^{-1}]$
at room temperature (de Valença et al. Reference de Valença, Kurniawan, Wagterveld, Wood and Lammertink2017).
As mentioned earlier, the governing equations for ion transport, potential distribution and fluid motion are described by the Nernst–Planck equation, Poisson equation and Navier–Stokes equation, respectively. To facilitate analysis, we adopt their dimensionless forms. The coordinates, velocity, pressure, time, concentration and electric potential are non-dimensionalised using the distance between the two ion-selective membranes
$\tilde{H}$
, the diffusion velocity
$\tilde{D}_{0}/\tilde{H}$
, the diffusion time
$\tilde{H}^{2}/\tilde{D}_{0}$
, the initial concentration
$\tilde{c}_{0}$
, the osmotic pressure
$\tilde{\eta }\tilde{D}_{0}/\tilde{H}^{2}$
and thermal voltage
$\tilde{V}_{T}=\tilde{k}_{B}\tilde{T}/\tilde{z}_{+}\tilde{e}$
, where
$\tilde{D}_{0}$
is the average diffusion coefficient of the ions,
$\tilde{k}_{B}$
is the Boltzmann constant and
$\tilde{e}$
is the elementary charge. In this study, the characteristic time scale is normalised using the diffusion time rather than the free-fall time, as ion diffusion fundamentally governs the flow (Mani & Wang Reference Mani and Wang2020). Based on these dimensionless references, the following dimensionless Poisson–Nernst–Planck and Navier–Stokes equations are obtained:
\begin{align}&\qquad\qquad\qquad\quad \frac{\partial c_{\pm }}{\partial t}=\boldsymbol{\nabla }\boldsymbol{\cdot }\left(D_{\pm }\boldsymbol{\nabla }c_{\pm }+z_{\pm }D_{\pm }c_{\pm }\boldsymbol{\nabla }V-\boldsymbol{U}\!c_{\pm }\right)\! ,\end{align}
\begin{align}&\qquad\qquad\qquad\qquad\qquad\qquad\quad \lambda _{D}^{2}{\nabla} ^{2}V=-\rho _{e} , \end{align}
\begin{align}& \frac{1}{Sc}\left(\frac{\partial \boldsymbol{U}}{\partial t}+\boldsymbol{U}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{U}\right)+\boldsymbol{\nabla }P-{\nabla} ^{2}\boldsymbol{U}=-\frac{\kappa }{{\lambda _{D}}^{2}}\left(c_{+}-c_{-}\right)\boldsymbol{\nabla }V-Ra\boldsymbol{\cdot }\frac{\left(c_{+}+c_{-}\right)}{2}\boldsymbol{\cdot }\boldsymbol{e}_{\boldsymbol{y}} , \end{align}
\begin{align}&\qquad\qquad\qquad\qquad\qquad\qquad\qquad \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{U}=0 . \end{align}
Here,
$t$
represents the dimensionless time, and
$c_{\pm }$
represent the dimensionless concentrations of cations and anions, respectively. Also,
$D_{+}$
and
$D_{-}$
represent the diffusion coefficients of cations and anions, respectively. In this paper,
$D_{+}=D_{-}=1$
is selected. Further,
$z_{\pm }$
represent the valences of the cation and anion, respectively,
$V$
is the dimensionless electric potential,
$\boldsymbol{U}$
is the dimensionless velocity vector field,
$\rho _{e}=c_{+}-c_{-}$
is the dimensionless space charge density,
$P$
is the dimensionless pressure and
$\boldsymbol{e}_{\boldsymbol{y}}$
is a unit vector directed downward. The parameter
$\lambda _{D}=\sqrt{\tilde{\varepsilon }_{0}\tilde{\varepsilon }_{\mathrm{r}}\tilde{k}_{B}\tilde{T}/(\tilde{z}_{+}\tilde{e})^{2}\tilde{c}_{0}}/\tilde{H}$
is the dimensionless Debye length, which describes the ratio of the dimensional Debye length
$\Big(\sqrt{\tilde{\varepsilon }_{0}\tilde{\varepsilon }_{\mathrm{r}}\tilde{k}_{B}\tilde{T}/(\tilde{z}_{+}\tilde{e})^{2}\tilde{c}_{0}}\Big)$
to the distance between two selective membranes (
$\tilde{H}$
). The EC exhibits weak sensitivity; thus, we select
$\lambda _{D}=10^{-3}$
(Karatay et al. Reference Karatay, Andersen, Wessling and Mani2016). The Schmidt number
$Sc=\tilde{\eta }/(\tilde{\rho }_{0}\tilde{D}_{0})$
is a dimensionless parameter that characterises the relationship between momentum transport and mass diffusion, defined as the ratio of the kinematic viscosity to the solute diffusion coefficient. Here, the Schmidt number is assumed to be a typical value (
$Sc=1000$
) (Druzgalski, Andersen & Mani Reference Druzgalski, Andersen and Mani2013; Karatay et al. Reference Karatay, Andersen, Wessling and Mani2016). Also,
$-\kappa (c_{+}-c_{-})\boldsymbol{\nabla }V/{\lambda _{D}}^{2}$
represents the electric field force, which controls the EC intensity, while
$Ra\boldsymbol{\cdot }({(c_{+}+c_{-})}/{2})\boldsymbol{\cdot }\boldsymbol{e}_{\boldsymbol{y}}$
represents the buoyant force, which controls the RBC intensity.
Building on the strategies outlined in Karatay et al. (Reference Karatay, Andersen, Wessling and Mani2016) regarding dimensionless parameters, this study focuses on exploring the combined effects of the hydrodynamic coupling coefficient, voltage difference and Rayleigh number on flow reversal. The first key dimensionless parameter is the hydrodynamic coupling coefficient
\begin{equation} \kappa =\big(\tilde{\varepsilon }_{0}\tilde{\varepsilon }_{\mathrm{r}}/\tilde{\eta }\tilde{D}_{0}\big)\big(\tilde{k}_{B}\tilde{T}/\tilde{z}_{+}\tilde{e}\big)^{2} .\end{equation}
The hydrodynamic coupling coefficient describes the impact of the physical properties of the solution on the flow. Increasing the hydrodynamic coupling coefficient enhances the EC strength.
The second key dimensionless parameter is the dimensionless voltage difference
\begin{equation} \Delta V=\Delta \tilde{V}/\tilde{V}_{T} .\end{equation}
The dimensionless voltage difference describes the impact of the voltage difference on EC. This is a crucial dimensionless parameter for initiating the EC and creating the ESC layer. Increasing the dimensionless voltage difference can enhance the strength of the EC flow.
The third key dimensionless parameter is the Rayleigh number
\begin{equation} Ra=\tilde{g}\tilde{\rho }_{0}\tilde{H}^{3}\tilde{\beta }\tilde{c}_{0}/(\tilde{\eta }\tilde{D}_{0}) .\end{equation}
This represents the Rayleigh number related to salt concentration. The Rayleigh number is sufficient to trigger the RBC.
2.2. Boundary conditions
As the top (
$y=1$
) is a perfectly cation-selective membrane with cation-selective permeability, only cations can pass through the surface. Therefore, the anion flux on the cation-selective membrane surface was considered zero, and a fixed cation concentration value was applied. The cation-selective membrane surface was assumed to follow no-slip conditions. A voltage was applied to the upper surface. The corresponding boundary equations are as follows (Mani & Wang Reference Mani and Wang2020):
\begin{equation} \left(D_{-}\boldsymbol{\nabla }c_{-}+z_{-}D_{-}\boldsymbol{\nabla }V\right)\boldsymbol{\cdot }\boldsymbol{n}=0 , \boldsymbol{U}=\boldsymbol{0} , c_{+}=2 , V=\Delta V ,\end{equation}
where
$\boldsymbol{n}$
is the unit outward normal vector at the boundary surface.
Similarly, for the bottom (
$y=0$
) perfectly cation-selective membrane, the boundary conditions are similar to those described above (Mani & Wang Reference Mani and Wang2020)
\begin{equation} \left(D_{-}\boldsymbol{\nabla }c_{-}+z_{-}D_{-}\boldsymbol{\nabla }V\right)\boldsymbol{\cdot }\boldsymbol{n}=0 , \boldsymbol{U}=\boldsymbol{0} , c_{+}=2 , V=0 .\end{equation}
For the fixed walls on the left and right sides (
$x=0$
and
$x=1$
), we assume that they are electrically neutral, with no-slip conditions and no flux of cations and anions. The corresponding boundary equations are as follows:
\begin{equation} \left(D_{-}\boldsymbol{\nabla }c_{-}+z_{-}D_{-}\boldsymbol{\nabla }V\right)\boldsymbol{\cdot }\boldsymbol{n}=0 , \left(D_{+}\boldsymbol{\nabla }c_{+}+z_{+}D_{+}\boldsymbol{\nabla }V\right)\boldsymbol{\cdot }\boldsymbol{n}=0 , \boldsymbol{U}=\boldsymbol{0} .\end{equation}
2.3. Critical selection of flow pattern
When EC dominates, increasing the hydrodynamic coupling coefficient and applied voltage enhances the electric body force, forming a double EC vortex. Conversely, when RBC dominates, increasing Ra amplifies the buoyancy force, reinforcing the LSC. The competition between RBC and EC determines the dominant flow regime; when EC surpasses RBC, EC vortex structure dominates the flow. Analysing the critical balance between LSC vorticity and EC vorticity is expected to yield a scaling law for the flow pattern. At the critical vortex height
$d_{EC,c}$
, the criterion can be expressed as
\begin{equation} \frac{u_\textit{LSC}}{1-d_{EC,c}}\sim \frac{u_\textit{EC}}{d_{EC,c}} ,\end{equation}
where
$u_\textit{LSC}$
is the velocity of LSC,
$u_\textit{EC}$
is the velocity of EC and
$d_{EC,c}$
is the critical height of the EC vortex. Here,
$(1-d_{EC,c})/d_{EC,c}$
can be treated as a constant-fitting parameter as
$d_{EC,c}$
is a fixed constant value (Kwak, Pham & Han Reference Kwak, Pham and Han2017). Therefore, scaling (2.15) can be expressed as
\begin{equation} u_\textit{LSC}\sim u_\textit{EC} .\end{equation}
Scaling (2.16) highlights the interplay between EC and RBC in governing the flow structure. When the EC velocity exceeds the LSC velocity, a double EC vortex is formed; conversely, when the LSC velocity surpasses the EC velocity, the flow becomes RBC dominated, and LSC persists.
To derive the LSC velocity, we apply the Grossmann–Lohse theory (Ahlers et al. Reference Ahlers, Grossmann and Lohse2009), a well-validated framework for analysing transport phenomena in RBC systems (Grossmann & Lohse Reference Grossmann and Lohse2000, Reference Grossmann and Lohse2001, Reference Grossmann and Lohse2002; Yang et al. Reference Yang, Verzicco and Lohse2015, Reference Yang, van der Poel, Ostilla-Mónico, Sun, Verzicco, Grossmann and Lohse2016; Shishkina, Grossmann & Lohse Reference Shishkina, Grossmann and Lohse2016; Wang et al. Reference Wang, Lohse and Shishkina2021a , Reference Wang, Liu, Verzicco, Shishkina and Lohseb ). Based on the exact global balance equations for kinetic dissipation, the theory provides the following formulation:
\begin{equation} \tilde{\varepsilon }_{u}=\tilde{\nu }\big\langle \left(\partial _{i}\tilde{u}_{j}\right)^{2}\big\rangle _{V}\sim \frac{\tilde{\nu }^{3}}{\tilde{H}^{4}}\boldsymbol{\cdot }Ra^{3/2} ,\end{equation}
where
$\tilde{\nu }=\tilde{\eta }/\tilde{\rho }_{0}$
denotes kinematic viscosity. According to the Grossmann–Lohse theory, the dissipation rate can be decomposed into the contributions of the boundary layer (BL) and bulk region (bulk) as follows:
\begin{equation} \tilde{\varepsilon }_{u}=\tilde{\varepsilon }_{u}^\textit{BL}+\tilde{\varepsilon }_{u}^{bulk} .\end{equation}
The dissipation rates are determined by contributions from the two BLs (figure 2
a). The BL thickness, denoted as
$\tilde{\lambda }_{u}$
, scales as
$\tilde{\lambda }_{u}\sim \tilde{H}/\sqrt{{Re}}$
(Ahlers et al. Reference Ahlers, Grossmann and Lohse2009). The Reynolds number, Re, is defined as
${Re}=\tilde{U}_{rms}\tilde{H}/\tilde{\nu }$
, where
$\tilde{U}_{rms}$
represents the root-mean-squared velocity, obtained through a time–space average of the velocity field. Through DNS, we found that the BL thickness follows the scaling
$\tilde{\lambda }_{u}/\tilde{H}\sim 1/\sqrt{{Re}}\sim 1/Ra^{-3/10}$
(figure 2
b), confirming the applicability of the Grossmann–Lohse theory to electrochemical systems.
(a) Schematic representation of the splitting of kinetic dissipation rates. The thickness of the kinetic BL, denoted as
$\tilde{\lambda }_{u}$
, is defined by the intersection of the tangents to the velocity profiles of the BL and bulk layers (Zhu, Fu & De Paoli Reference Zhu, Fu and De Paoli2024). (b) Scaling behaviour of the kinetic thickness.
The dissipation rate within the BL can be estimated as
\begin{equation} \tilde{\varepsilon }_{u}^\textit{BL}=\tilde{\nu }\big\langle \left(\partial _{i}\tilde{u}_{j}\right)^{2}\big\rangle _{V}\sim \tilde{\nu }\boldsymbol{\cdot }\left(\frac{\tilde{U}}{\tilde{\lambda }_{u}}\right)^{2}\boldsymbol{\cdot }\frac{\tilde{\lambda }_{u}}{\tilde{H}}=\frac{\tilde{\nu }^{3}}{\tilde{H}^{4}}\boldsymbol{\cdot }{Re}^{5/2} .\end{equation}
Considering the balance of dissipation rates from (2.17) and (2.19) and neglecting the contributions of dissipation in the bulk region, we have
\begin{equation} {Re}\sim Ra^{3/5} .\end{equation}
We further assumed that the root-mean-squared velocity of the RBC approximately corresponded to the LSC velocity. Thus, we have
\begin{equation} u_\textit{LSC}\sim Ra^{3/5} .\end{equation}
Next, we derived the relationship for the rolling velocity of the EC vortices near the membrane–solution surface. For a given Ra, as the applied voltage increases, the EC strength increases. By considering the balance among the viscous term, the local pressure gradient and the electric body force near the membrane–solution surface, we obtain
\begin{equation} {\nabla} ^{2}\boldsymbol{U}\sim -\boldsymbol{\nabla }P-\kappa {\nabla} ^{2}V\boldsymbol{\cdot }\boldsymbol{\nabla }V .\end{equation}
Considering the form of the Maxwell stress near the membrane–solution interface, this yields (Rubinstein & Zaltzman Reference Rubinstein and Zaltzman2000)
\begin{equation} \frac{\partial ^{2}u}{\partial y^{2}}/\kappa \sim \frac{1}{2}\frac{\partial }{\partial x}\left(\frac{\partial V}{\partial y}\right)^{2}-\left(\frac{\partial ^{2}V}{\partial y^{2}}\right)\frac{\partial V}{\partial x} .\end{equation}
In (2.23), the first term on the right-hand side,
$[\partial (\partial V/\partial y)^{2}]/(2\partial x)$
, describes the contribution of the pressure-driven force to the EC velocity, whereas the second term,
$-(\partial ^{2}V/\partial y^{2})(\partial V/\partial x)$
, represents the contribution of the tangential electric body force within the ESC layer to the EC velocity. Based on the Poisson–Nernst–Planck equation, we can obtain the voltage within the ESC layer (Rubinstein & Zaltzman Reference Rubinstein and Zaltzman2000)
\begin{equation} V\sim \frac{2\sqrt{2}}{3}\lambda _{D}^{-1}j^{1/2}y_\textit{ESC}^{\frac{3}{2}}-\frac{2\sqrt{2}}{3}\lambda _{D}^{-1}j^{1/2}\left(y_\textit{ESC}-y\right)^{\frac{3}{2}} .\end{equation}
In (2.24),
$y_\textit{ESC}$
represents the thickness of the ESC layer, which can be identified by the location corresponding to the peak of the ESC. Substituting (2.24) into (2.23) and performing two integrations, we obtain
\begin{equation} u_\textit{EC}/\kappa =u\left(y=y_\textit{ESC}\right)/\kappa \sim -\frac{1}{9}\lambda _{D}^{-2}y_\textit{ESC}^{3}\frac{\partial j}{\partial x}-\frac{4}{9}\lambda _{D}^{-2}y_\textit{ESC}^{3}\frac{\partial j}{\partial x}-\frac{4}{3}\lambda _{D}^{-2}jy_\textit{ESC}^{2}\frac{\partial y_\textit{ESC}}{\partial x}.\end{equation}
Assuming that the voltage (
$V_\textit{ESC}$
) within the ESC layer is approximately equal to the total voltage difference (
$\Delta V$
) and satisfies the relationship
$V(y=y_\textit{ESC})-V(y=0)\sim \Delta V$
, the expression for the length of the ESC layer can be derived from (2.24) as follows:
\begin{equation} y_\textit{ESC}\sim \frac{\sqrt[3]{9}}{2}\lambda _{D}^{2/3}\Delta V^{2/3}j^{-1/3} .\end{equation}
By substituting (2.26) in (2.25), the velocity of the EC vortex can be obtained as follows:
\begin{equation} u_\textit{EC}/\kappa \sim -\frac{1}{8}V^{2}j^{-1}\frac{\partial j}{\partial x}-V\frac{\partial V}{\partial x} .\end{equation}
Using the first-order approximation
$1/\partial x\sim 1/\partial y\sim 1/d_\textit{EC}$
and the fact that
$d_\textit{EC}\sim \Delta V\sqrt{\kappa }$
(Kwak et al. Reference Kwak, Pham, Lim and Han2013), we can obtain
\begin{equation} u_\textit{EC}/\Delta V\sim \sqrt{\kappa } .\end{equation}
For a given hydrodynamic coupling coefficient scenario, we observed that the velocity has a linear relationship with the applied voltage, which is consistent with the previous results (Rubinstein et al. Reference Rubinstein, Manukyan, Staicu, Rubinstein, Zaltzman, Lammertink, Mugele and Wessling2008). Considering the velocity balance relations from (2.16), we can derive the scaling relation for the critical Rayleigh number
\begin{equation} Ra_{c}\sim \Delta V^{5/3}\kappa ^{5/6} .\end{equation}
Scaling (2.29) reveals the role of the relation between EC and RBC in regulating flow reversal. When the scaling factor
$\Delta V^{5/3}\kappa ^{5/6}$
exceeds the critical Rayleigh number, the flow pattern is dominated by the EC. When the scaling factor
$\Delta V^{5/3}\kappa ^{5/6}$
is lower than the critical Rayleigh number, the flow pattern is dominated by LSC. Subsequently, DNS was used to explore the flow behaviour within the parameter spaces of
$Ra , \Delta V$
and
$\kappa$
, which facilitated a comprehensive analysis of the flow pattern.
3. Discussion
For the simulation parameters, variations in the applied voltage do not affect the concentration difference and the hydrodynamic coupling coefficient. As a result, within our dimensionless parameter space, the Rayleigh number, the hydrodynamic coupling coefficient and the applied voltage are treated as independent parameters (Karatay et al. Reference Karatay, Andersen, Wessling and Mani2016). Based the existing approach (Karatay et al. Reference Karatay, Andersen, Wessling and Mani2016), we consider
$0\leq \kappa \leq 1 , 20\leq \Delta V\leq 100$
and
$10^{5}\leq Ra\leq 2.5\times 10^{6}$
. In the absence of an applied electric field, the velocity, pressure and electric potential fields remain undeveloped and are initialised as
$c_{+}(x,y,t=0)=c_{-}(x,y,t=0)=1 , \boldsymbol{U}(x,y,t=0)=\boldsymbol{0} , V(x,y,t=0)=0$
and
$P(x,y,t=0)=0$
. Direct numerical simulation of the Navier–Stokes and Poisson–Nernst–Planck equations was performed using the finite element method.
Previous studies (Karatay et al. Reference Karatay, Andersen, Wessling and Mani2016) indicated that a two-dimensional model can effectively capture the flow characteristics in coupled RBC and EC systems. Therefore, this study employs a two-dimensional computational domain to reduce the computational costs associated with three-dimensional models. Accurate resolution of the ESC layer requires careful grid refinement. We used a non-uniform structured quadrilateral grid with local refinement near the membrane, dividing the computational domain into a 300(x) × 300(y) grid. To precisely capture the ESC layers, we set the minimum mesh size near the ion-selective membrane surface to O(
$\lambda _{D}$
), while the maximum mesh size in the ESC layer was set to O(
$\lambda _{D}^{2/3}$
). The longitudinal grid cell ratio was set to 103 with a symmetric geometric distribution, while an asymmetric grid set-up was used in the horizontal direction. This fine mesh configuration effectively captures the flow induced by the ESC. Grid configurations of 600(x) × 700(y), 500(x) × 600(y) and 300(x) × 300(y) were analysed sequentially. The numerical results confirmed that the 300(x) × 300(y) grid configuration effectively captures the EC, and RBC. For further details on the simulation set-up, please refer to our previous studies (Liu et al. Reference Liu, Zhou and Shi2020, Reference Liu, Zhou and Shi2022; Xu et al. Reference Xu, Gu, Liu, Huo, Zhou, Rubinstein, Bazant, Zaltzman, Rubinstein and Deng2020).
For pure EC (
$Ra=0$
), the voltage–current curve serves as a key benchmark, exhibiting three distinct stages. Initially, the current increases linearly with voltage; this is followed by a saturation phase. Beyond a critical threshold, EC emerges, where EC vortices enhance cation transport toward the membrane (Mani & Wang Reference Mani and Wang2020). A parametric voltage scan confirms that the present DNS results are consistent with existing DNS (Druzgalski et al. Reference Druzgalski, Andersen and Mani2013; Mani & Wang Reference Mani and Wang2020) regarding the overall trends (figure 3
a).
Benchmark results. (a) Voltage–current (
$\Delta V$
-
$j$
) characteristics for pure EC, with parameters
$\kappa =0.5$
and
$Ra=0$
. (b) Relationship between the time- and space-averaged flux (
$j$
) and the Rayleigh number (Ra), with parameters
$\Delta V=20$
and
$\kappa =0$
. The simulation data closely follow the scaling law, confirming the accuracy of our DNS in capturing transport scaling behaviour. The ion flux is quantified by computing the difference between cation and anion fluxes, followed by averaging over both temporal and spatial domains.
For pure RBC (
$\Delta V=20$
and
$\kappa =0$
), we found that the ion flux satisfies
$Ra^{1/3}$
, which is similar to the heat flux scaling in thermally driven RBC (Malkus Reference Malkus1954; Priestley Reference Priestley1954; Yang et al. Reference Yang, van der Poel, Ostilla-Mónico, Sun, Verzicco, Grossmann and Lohse2015) (figure 3b
). This law characterises the relationship between the transport efficiency and Rayleigh number, a fundamental criterion in RBC. Our simulation results across different Ra values align with the scaling results, demonstrating the accuracy of our DNS in capturing transport behaviour. These two validations enable further exploration of the flow behaviour in solutally driven RBC.
3.1. Flow characteristics
When EC is weak or negligible, we first demonstrate the LSC and flow reversal in an electrochemical system, which is similar to the thermally driven RBC (figure 4 a) (Chandra & Verma Reference Chandra and Verma2013; Chen et al. Reference Chen, Huang, Xia and Xi2019). This similarity arises because the solutally driven RBC presented in this paper shares fundamental characteristics with thermally driven RBC – both involve the transport of scalar fields (concentration and temperature fields), with buoyancy as the driving force (buoyancy caused by a temperature gradient and buoyancy caused by a salt concentration gradient). Both thermally and solutally driven RBC exhibit analogous flow structures, including LSC and corner vortices (Chandra & Verma Reference Chandra and Verma2013; Chen et al. Reference Chen, Huang, Xia and Xi2019). The LSC, also known as the ‘wind’ of convection, emerges from the background turbulent flow when the Rayleigh number is sufficiently high. In confined convection, the size of the LSC is primarily determined by the size of the convection cell. In electrochemical systems, the LSC size is significantly smaller than in thermally driven RBC, owing to the much smaller size of typical electrochemical cells (de Valença et al. Reference de Valença, Kurniawan, Wagterveld, Wood and Lammertink2017). To understand the mechanism of flow reversal, we adopt the Brown and Ahlers model (Brown & Ahlers Reference Brown and Ahlers2007), where the system has two quasi-stable solutions, and quasi-random turbulent forcing induces transitions between them. The clockwise and counterclockwise flow states correspond to the two quasi-stable solutions. Flow reversal depends on disturbances caused by quasi-random turbulent forcing in the system. Small irregular disturbances lead to random effects, triggering slight deviations in the LSC and the periodic reconnection of the corner vortices, ultimately resulting in flow reversal (figures 4 a and 4 c).
Flow characteristics. (a) The LSC and flow reversal in an electrochemical system when EC is negligible. The ESC forms near the bottom perfect cation-selective membrane surface, creating a concentration gradient, generating a density difference, thus forming RBC. The periodic reconstruction of corner vortices leads to LSC reversal. Here,
$Ra=5\times 10^{5} , \Delta V=20$
and
$\kappa =0.05$
. (b) The EC vortex pattern. Here,
$Ra=5\times 10^{5} , \Delta V=80$
and
$\kappa =0.5$
. As soon as the double vortex formed by the EC dominates the flow, the global angular momentum is zero when EC vortices are symmetric. Temporal and spatial fluctuations of salt concentration at the cutline (y = 0.5), corresponding to those shown in (c) LSC and (d) EC scenarios, respectively. The transient characteristics of (e) angular momentum and (f) ion flux. For more flow videos, see the supporting materials are available at https://doi.org/10.1017/jfm.2025.10678.
When the EC exceeds the RBC, an EC vortex occurs. A quasi-periodically oscillating vertical transport is observed along the y-axis between the two counter-rotating EC vortices, extending from the bottom to the top of the cell. This motion reflects a dynamic exchange of ion concentration between cation-selective membranes through the region between the vortices. It can be closely associated with the long-wave oscillations observed, as shown in figure 4(d) or figure 7(c). The EC vortex pairs in the bulk zone drive small-scale vortices in the ESC zone, extending oscillations across the entire cell along the region near the
$x=0.5$
centreline, and causing high-frequency oscillations.
Notably, the similar stabilisation results induced by the EC vortex, along with factors such as the non-Oberbeck–Boussinesq effect (Xia et al. Reference Xia, Wan, Liu, Wang and Sun2016), boundary conditions (Zhang et al. Reference Zhang, Xia, Zhou and Chen2020), magnetic field (Zurner et al. Reference Zurner, Schindler, Vogt, Eckert and Schumacher2020) and boundary conditions (Huang et al. Reference Huang, Wang, Xi and Xia2015), can stabilise LSC, although their underlying physical mechanisms differ significantly. The non-Oberbeck–Boussinesq effect primarily introduces asymmetry by altering the fluid’s thermal expansion properties, internal heating modifies the flow structure through the distribution of heat sources, boundary condition variations directly impose constraints on the flow and magnetic fields exert control by generating Lorentz forces. In contrast, in this study, the electric field force induces the formation of EC vortices near the bottom boundary. When the influence of EC surpasses that of RBC, LSC disappears, which is similar to a phenomenon of LSC stabilisation in thermally driven RBC. This strategy of using one type of flow to control another is fundamentally different from conventional approaches (Huang et al. Reference Huang, Wang, Xi and Xia2015; Xia et al. Reference Xia, Wan, Liu, Wang and Sun2016; Zhang et al. Reference Zhang, Xia, Zhou and Chen2020; Zurner et al. Reference Zurner, Schindler, Vogt, Eckert and Schumacher2020).
To quantify the flow reversal, we calculate the global angular momentum using (3.1). This formula has been successfully used for quantitative characterisation of reversal (Sugiyama et al. Reference Sugiyama, Ni, Stevens, Chan, Zhou, Xi, Sun, Grossmann, Xia and Lohse2010). To quantify the ion transport efficiency, we calculated the ion flux (Druzgalski et al. Reference Druzgalski, Andersen and Mani2013). The instantaneous global angular momentum and instantaneous flux are defined as follows:
\begin{align}&\qquad\qquad\qquad\qquad\quad L\left(t\right)=\left\langle \left(x-0.5\right)v-\left(y-0.5\right)u\right\rangle , \end{align}
\begin{align}& j\left(t\right)=-\left\langle D_{+}\boldsymbol{\nabla }c_{+}+z_{+}D_{+}c_{+}\boldsymbol{\nabla }V-vc_{+}-\left(D_{-}\boldsymbol{\nabla }c_{-}+z_{-}D_{-}c_{-}\boldsymbol{\nabla }V-vc_{-}\right)\right\rangle\! , \end{align}
where
$\langle \cdot \rangle$
denotes the surface average. Figures 4(e) and 4(f) show the changes in the angular momentum and ion flux during the flow reversal process. During the initial stage of activation, there is a reversal of the LSC as the sign of
$L(t)$
undergoes periodic changes, indicating a flow reversal. However, after t = 0.008, the EC gradually affects the bulk region, influencing the development of the local LSC, thereby eliminating the LSC. The blue dotted line in figure 4(e) shows period oscillations when L(t = 0.07) = 0, indicating that the flow field temporarily reaches a nearly symmetric state. When the flow reaches a steady fluctuation state, the sign of
$L(t)$
no longer changes, indicating that LSC is eliminated. In addition, EC significantly enhances the ion transport efficiency in the RBC by approximately five times (figure 4
f).
3.2. Mechanism of pattern selection
Next, we investigated the mechanism underlying the pattern transitions. In electrochemical systems, ion-selective transport at the membrane–solution surface leads to a non-electroneutral ESC layer near the membrane surface (figure 5 a). Within the ESC layer, the inhomogeneous space charge distribution induces an intense local electric field (Mani & Wang Reference Mani and Wang2020), which not only drives ion motion along the membrane surface, generating EC vortex structures, but also directly acts on the LSC, introducing significant electric volume forces (figures 5 b and 5 c).
Driving force analysis of EC. (a) Space charge density distribution. (b) Electric field volume force distribution. (c) Electric field body force profile, which is calculated by taking the line average and time average along the x direction. The strong electric field volume force within the ESC layer can suppress the development of corner vortices and LSC structures, thereby stabilising LSC. Here,
$Ra=5\times 10^{5} , \Delta V=80$
and
$\kappa =0.5$
.
This strong electric field establishes quasi-stable EC vortex structures (Mani & Wang Reference Mani and Wang2020), which eliminate the LSC. In other words, the strong electric forcing within the ESC not only disrupts the original LSC structure but also reorganises the flow pattern, transitioning the system from a buoyancy-driven regime (LSC structure) to an electrically driven regime (EC structure) dominated by electric body forces. Unlike thermally driven RBC, where reorientation, reversal or cessation of flow may occur (Huang et al. Reference Huang, Wang, Xi and Xia2015; Xia et al. Reference Xia, Wan, Liu, Wang and Sun2016; Zhang et al. Reference Zhang, Xia, Zhou and Chen2020; Zurner et al. Reference Zurner, Schindler, Vogt, Eckert and Schumacher2020), this EC control strategy uniquely alters the convective behaviour without inducing such effects.
Local perturbations within the ESC – such as space charge distribution, electric body forces and vorticity – create asymmetry (Mani & Wang Reference Mani and Wang2020), leading to a non-zero global angular momentum. This periodic asymmetry results from strong electric forcing and uneven ion transport, and leads to the continuous ejection of small-scale EC vortices from the bottom surface. These interactions drive period oscillations in global angular momentum. To confirm this mechanism, we analysed the spatial distribution of vorticity and statistically compared local flow properties. Figure 6 presents the vorticity distribution within the ESC (figure 6 a,b) alongside the electric body force distribution (figure 6 c). The findings indicate that the strong electric field in the ESC creates localised inhomogeneities, which in turn generate vortex pair asymmetry. The interactions between these local small-scale EC vortices within the ESC layer extend the oscillation to the entire cell.
Asymmetry caused by EC vortices within ESC. (a) Spatial distribution of vorticity. (b) Distribution of vorticity and (c) electric field force along the section line (y = 0.06). Asymmetric distribution of vortices in the ESC layer due to strong electric field forces and local disturbances. Here,
$Ra=5\times 10^{5} , \Delta V=80$
and
$\kappa =0.5$
.
Effect of voltage on flow reversal. Effect of voltage on the spatio-temporal fluctuations of (a) global angular momentum, (b) ion flux and (c) salt concentration. As the voltage gradually increases, the EC strength increases, which leads to LSC stabilisation.
Effect of the hydrodynamic coupling coefficient on flow characteristics. Effect of hydrodynamic coupling coefficient on the fluctuations of (a) global angular momentum, (b) ion flux and (c) salt concentration.
3.3. Statistical characteristics in parameter space
To investigate EC regulation of flow reversal, we analysed the parameter space defined by the Rayleigh number, the hydrodynamic coupling coefficient and the applied voltage. Figure 7 illustrates the impact of voltage on flow reversal. As previously discussed, increasing the voltage strengthens EC, leading to vortex pair; the flow reversal has disappeared. When the applied voltage is
$\Delta V=50$
, the flow initially undergoes periodic reversals due to the reconnection of corner vortices. However, as EC fully develops and overcomes buoyancy forces, the LSC is eliminated. Figure 8 demonstrates the effect of the hydrodynamic coupling coefficient on flow reversal. At
$\kappa =0.5$
, LSC is completely eliminated. The enhancement of the electrostatic force also led to the EC vortex. As the hydrodynamic coupling coefficient increased, the global angular momentum, which alternated between positive and negative changes, was confined to below zero, locking the flow and making the EC vortex pair dominant. Therefore, when the LSC is eliminated, the salt concentration fluctuation centre exhibits a low-concentration strip structure. Additionally, increasing the hydrodynamic coupling coefficient improves the ion transport efficiency, as EC vortices enhance ion movement along the membrane surface. Figure 9 shows the influence of the Rayleigh number. As the Rayleigh number increases, RBC intensifies, and the frequency of flow reversal accelerates. This behaviour aligns with classical thermally driven RBC, where higher Rayleigh numbers promote more frequent reversals.
Effect of Rayleigh number on flow characteristics. Effect of voltage on the spatio-temporal fluctuations of (a) global angular momentum, (b) ion flux and (c) salt concentration. As the Rayleigh number increases, flow reversals become more frequent. These three points correspond to (99.6, 7.5
$\times 10^{5}$
), (99.6, 1
$\times 10^{6}$
) and (99.6, 2.5
$\times 10^{6}$
) in figure 10(b).
Critical selection of different flows. (a) For fixed
$Ra=5\times 10^{5}$
, the effect of the hydrodynamic coupling coefficient and voltage on flow pattern. The critical curve is
$(Ra/\alpha )^{3/5}=\Delta V\kappa ^{1/2}$
. The hollow points correspond to EC vortex, while solid points correspond to the presence of LSC. (b) Critical scaling law. Critical selection of the flow pattern for the phase diagram replotted in dimensionless numbers:
$\Delta V^{5/3}\kappa ^{5/6}$
vs.
$Ra$
. The hydrodynamic coupling coefficient and applied voltage describe the effect of the EC, while the Rayleigh number describes the effect of the RBC. Competition between RBC and EC dominates the flow pattern. For strong RBC, the corresponding Ra is larger, requiring a larger
$\Delta V^{5/3}\kappa ^{5/6}$
to form an EC vortex. For convenience, two area colours are used to distinguish flow regimes. Grey areas represent the EC vortex, whereas white areas represent the presence of LSC.
3.4. Critical scaling law
Finally, a comprehensive analysis of the flow pattern is conducted by systematically exploring the broad parameter space defined by
$Ra , \Delta V$
and
$\kappa$
. Figure 10(a) presents a phase diagram illustrating the combined effects of these parameters on the flow pattern under a fixed
$Ra=5\times 10^{5}$
. For a fixed
$\kappa$
, an increase in the applied voltage is required to strengthen the EC. For a fixed voltage, the EC can be enhanced by increasing the hydrodynamic coupling coefficient, thereby forming an EC vortex structure. These findings were further extended to establish a generalised framework for predicting the critical parameter values necessary for the flow pattern across different Rayleigh numbers. The critical scaling law is given by
\begin{equation} Ra_{c}=\alpha \Delta V^{5/3}\kappa ^{5/6} .\end{equation}
In our parameter space (
$0\leq \kappa \leq 1 , 20\leq \Delta V\leq 100$
and
$10^{5}\leq Ra\leq 2.5\times 10^{6}$
), the fitting coefficient was
$\alpha \approx 1751$
. When scaling
$\Delta V^{5/3}\kappa ^{5/6}$
on the x-axis and the Rayleigh number on the y-axis, the normalised distribution of the scattered data was obtained (figure 10
b). At
$Ra=5\times 10^{5}$
, we derived the critical EC intensity required for the double-roll EC pattern, which is jointly determined by the hydrodynamic coupling coefficient and applied voltage. The critical threshold was determined using
$\Delta V^{5/3}\kappa ^{5/6}\approx 286$
. When the EC intensity (
$\Delta V^{5/3}\kappa ^{5/6}$
) surpasses a critical value, the EC dominated the flow. In this regime, the flow was characterised by an EC vortex pair.
We further extended the results to scenarios with different Rayleigh numbers, as indicated by the black pentagram markers in figure 10(b). For a fixed voltage (
$\Delta V=80$
), we investigated the interplay between RBC and EC by systematically varying the hydrodynamic coupling coefficient (from 0 to 1). As the Rayleigh number increases, the buoyancy-driven flow becomes more dominant, and thus a stronger electric field force is required to form the EC vortex structure. This implies that a higher EC intensity is required to effectively compete with LSC. When the EC intensity surpasses a critical threshold, electrostatic forces disrupt the LSC structure, preventing its reversal and transitioning of the system into an EC-dominated regime characterised by a double-roll EC vortex pair.
4. Conclusion
This study reveals the flow reversal phenomenon of LSC induced by the periodic reconstruction of corner vortices in electrochemical systems and proposes a strategy for controlling the flow pattern through EC. The EC flow alters the LSC dynamics by inducing a strong electric field force within the ESC layer. The results demonstrate that enhancing EC significantly improves the ion transport efficiency, replacing LSC with EC vortex-dominated flow patterns. By establishing a critical relation based on Grossmann–Lohse theory, we clarify the joint regulatory effects of the Rayleigh number, voltage and hydrodynamic coupling coefficient on RBC and EC. Our findings offer insights into the control of flow pattern in electrochemical systems, with implications for electrochemical applications such as desalination, energy harvesting and sensing technologies.
Inspired by thermally driven RBC (Ahlers et al. Reference Ahlers, Grossmann and Lohse2009; Lohse & Xia Reference Lohse and Xia2010; Lohse & Shishkina Reference Lohse and Shishkina2024), exploring the Reynolds number scaling (Ahlers et al. Reference Ahlers, Grossmann and Lohse2009), transport scaling (Ahlers et al. Reference Ahlers, Grossmann and Lohse2009) and cell aspect ratio (Ahlers et al. Reference Ahlers, Bodenschatz, Hartmann, He, Lohse and Reiter2022) in solution-driven RBC would be interesting. We will extend the Grossmann–Lohse theory to analyse the coupled interactions between EC and RBC in electrochemical systems. In addition, the Joule heating and a varying Schmidt number are not considered in the present study, which represent limitations of the current work. A more comprehensive investigation on these factors deserves to be conducted in future work.
Supplementary movies
Supplementary movies are available at https://doi.org/10.1017/jfm.2025.10678.
Acknowledgements
We thank Professor D.-S. Deng of Fudan University and the three reviewers for their helpful comments on the draft.
Funding
This work was financially supported by the China Postdoctoral Science Foundation (No. 2023M741970) and the Postdoctoral Fellowship Program of the CPSF (No. GZB20230361).
Declaration of interests
The authors report no conflict of interest.
