1. Introduction
Electrospray atomization is based on the use of an electric field to produce a liquid jet, which ultimately breaks into charged droplets (Cloupeau & Prunet-Foch Reference Cloupeau and Prunet-Foch1989; Fernández de La Mora Reference Fernández de la Mora2007). There exist several electrospraying modes, among which the so-called cone-jet has the unique ability to generate droplets with narrow size distributions and controllable average diameter ranging from hundreds of micrometres down to a few nanometres (Gañán-Calvo et al. Reference Gañán-Calvo, López-Herrera, Herrada, Ramos and Montanero2018; Rosell-Llompart, Grifoll & Loscertales Reference Rosell-Llompart, Grifoll and Loscertales2018). Typically, a liquid is fed through a capillary needle towards its tip while applying a potential difference with respect to a facing electrode. The liquid meniscus at the tip adopts a conical shape referred to as a Taylor cone (Taylor Reference Taylor1964), which emits a thin jet from its apex. The fluid in the jet accelerates forced by the electric field, undergoes Rayleigh instability and breaks into charged droplets (Rayleigh Reference Rayleigh1882). The consideration of different balances in the cone-jet and experimental data has led to scaling laws for the current emitted in the form of charged droplets and for the radius of the jet. For example, Gañán-Calvo et al. (Reference Gañán-Calvo, Barrero and Pantano-Rubiño1993), postulating that the surface charge must be in a state of electrostatic quasiequilibrium everywhere in the cone-jet, derive the following scaling law for the current (Gañán-Calvo Reference Gañán-Calvo1999, Reference Gañán-Calvo2004):
\begin{equation} I\cong 2.47 (\gamma K Q)^{1/2}, \end{equation}
while Fernández de La Mora & Loscertales (Reference Fernández de La Mora and Loscertales1994) propose that the surface charge must depart from equilibrium near the tip of the Taylor cone, leading to
\begin{equation} I\cong f(\varepsilon ) (\gamma K Q)^{1/2}. \end{equation}
Here
$K$
,
$\gamma$
,
$\varepsilon$
and
$Q$
are the electrical conductivity, surface tension, dielectric constant and the volumetric flow rate of the liquid, respectively. The scaling law (1.1) fits well the experimental data for a large number of liquids (Gañán-Calvo et al. Reference Gañán-Calvo, López-Herrera, Herrada, Ramos and Montanero2018). It is remarkable that Gañán-Calvo et al. (Reference Gañán-Calvo, Barrero and Pantano-Rubiño1993) obtained the accurate factor of 2.47 as early as 1993. Experimental (Fernández de La Mora & Loscertales Reference Fernández de La Mora and Loscertales1994; Chen & Pui Reference Chen and Pui1997) and numerical (Gamero-Castaño & Magnani Reference Gamero-Castaño and Magnani2019b
) studies have found a dependence of the current on the dielectric constant, as suggested by (1.2). We note that the dielectric constant survives in the equations describing the electrohydrodynamics of cone-jets, specifically in the jump condition for the electric field across the surface, and therefore it must have an effect on the solution. However, this effect should be small if the surface charge is not far from equilibrium. The diameter of the jet and droplets follows the scaling law (Gañán-Calvo et al. Reference Gañán-Calvo, Lasheras, Dávila and Barrero1994; Gañán-Calvo Reference Gañán-Calvo1999, Reference Gañán-Calvo2004)
\begin{equation} d \approx \left ( \frac {\rho \varepsilon _{0}Q^{3}}{\gamma K} \right )^{1/6}, \end{equation}
where
$\varepsilon _0$
is the permittivity of vacuum and
$\rho$
is the density of the liquid. This result has been validated by experiments and numerical calculations (Gañán-Calvo et al. Reference Gañán-Calvo, Davila and Barrero1997, Reference Gañán-Calvo, López-Herrera, Herrada, Ramos and Montanero2018; Gamero-Castaño & Hruby Reference Gamero-Castaño and Hruby2002; Gamero-Castaño & Magnani Reference Gamero-Castaño and Magnani2019b
).
Taylor and Melcher’s leaky dielectric (LD) model (Melcher & Taylor Reference Melcher and Taylor1969; Melcher Reference Melcher1981; Saville Reference Saville1997) has become the default first-principles formulation for the study of cone-jets. The LD model assumes that the volumetric charge density is zero. Furthermore, it considers a zero-thickness surface where all free charge is concentrated in the form of a surface charge variable with an associated conservation equation. These assumptions simplify the model: charge neutrality in the bulk eliminates the need to track the concentration of ionic species, and their conservation equations are replaced by a simple Ohm’s law with a constant electrical conductivity; the electric potential in the bulk only needs to fulfil the Laplace equation, instead of the Poisson equation; and the surface charge approximation eliminates the need to resolve the Debye layer that naturally forms under the surface. Despite these simplifying assumptions, numerical solutions based on the LD model have successfully predicted the behaviour of cone-jets. Hartman et al. (Reference Hartman, Brunner, Camelot, Marijnissen and Scarlett1999) reproduced experimental data, including the scaling law for the electrospray current. Higuera (Reference Higuera2003) used Taylor’s potential (Taylor Reference Taylor1964) as a far field boundary condition to obtain universal solutions independent of the geometry of the electrodes; the use of Taylor’s potential as a far-field boundary condition appears for the first time in an analytical solution by Gañán-Calvo (Reference Gañán-Calvo1997). Herrada et al. (Reference Herrada, López-Herrera, Gañán-Calvo, Vega, Montanero and Popinet2012) simulated cone-jets including the geometry of the electrodes, and were able to reproduce measurements of the geometry of cone-jets and of the emitted current. Gamero-Castaño & Magnani (Reference Gamero-Castaño and Magnani2019a , Reference Gamero-Castaño and Magnanib ) employed the far field Taylor potential to find accurate solutions for cone-jets in a wide range of dimensionless flow rates, electric Reynolds numbers and dielectric constants, and to study the minimum flow rate. Recently, Magnani & Gamero-Castaño (Reference Magnani and Gamero-Castaño2024) have extended their LD formulation to account for ohmic and viscous dissipation and the associated self-heating, which is key for liquids with high electrical conductivities (Gamero-Castaño Reference Gamero-Castaño2010, Reference Gamero-Castaño2019). In addition to the inherently steady cone-jet, the LD model has also been employed to simulate time-dependent tip streaming problems (Collins et al. Reference Collins, Jones, Harris and Basaran2008, Reference Collins, Sambath, Harris and Basaran2013; Gawande, Mayya & Thaokar Reference Gawande, Mayya and Thaokar2020; Wagoner et al. Reference Wagoner, Vlahovska, Harris and Basaran2021).
The volume of fluid (VOF) and phase-field techniques have been used to account for non-zero volumetric charge density in cone-jets. In the simplest case the concentration fields of ionic species are ignored while using a volumetric charge density that fulfils a conservation of charge equation with convection and conduction terms (Herrada et al. Reference Herrada, López-Herrera, Gañán-Calvo, Vega, Montanero and Popinet2012), and a constant electrical conductivity (for a detailed discussion see Gañán-Calvo et al. Reference Gañán-Calvo, López-Herrera, Herrada, Ramos and Montanero2018). The VOF solution by Gañán-Calvo et al. (Reference Gañán-Calvo, López-Herrera, Rebollo-Munoz and Montanero2016) shows that the assumptions of the LD model fail in ultrafast ejection and tip streaming from Taylor cones, i.e. in time-dependent problems. Herrada et al. (Reference Herrada, López-Herrera, Gañán-Calvo, Vega, Montanero and Popinet2012) have shown that the LD model and the VOF technique agree well in steady-state calculations of cone-jets. Recently, several authors have used the VOF technique to time-march their algorithms to obtain the steady-state cone-jet solution (Dastourani, Jahannama & Eslami-Majd Reference Dastourani, Jahannama and Eslami-Majd2018; Huh & Wirz Reference Huh and Wirz2022; Mai et al. Reference Mai, Vu, Dinh, Vu, Tran, Dau and Ngo2023; Suo et al. Reference Suo, Zhang, Huang, Wang, Jia, Yang, Zhang, Li, Tu and Song2023). The VOF and phase field techniques have been used to study other electrohydrodynamic problems (Tomar et al. Reference Tomar, Gerlach, Biswas, Alleborn, Sharma, Durst, Welch and Delgado2007; López-Herrera et al. Reference López-Herrera, Popinet and Herrada2011; Ferrera et al. Reference Ferrera, López-Herrera, Herrada, Montanero and Acero2013; Gañán-Calvo et al. Reference Gañán-Calvo, López-Herrera, Rebollo-Munoz and Montanero2016; Misra & Gamero-Castaño Reference Misra and Gamero-Castaño2022, Reference Misra and Gamero-Castaño2023).
There are few studies that retain the full electrokinetic details to study steady cone-jet or unsteady tip-streaming. The extent to which the bulk can be assumed neutral is governed by the ratio between the Debye length and the characteristic length scale. When this ratio is not sufficiently small the electrical conductivity in the bulk cannot be regarded as constant but must be expressed in terms of the concentrations and mobilities of ionic species, and the LD model should fail. The VOF technique combined with the standard Poisson–Nernst–Plank (PNP) electrokinetic formulation has been employed to study the breakup of electrified jets (López-Herrera et al. Reference López-Herrera, Gañán-Calvo, Popinet and Herrada2015) and the evolution of electrified droplets (Pillai et al. Reference Pillai, Berry, Harvie and Davidson2015, Reference Pillai, Berry, Harvie and Davidson2016). These studies have shown that electrokinetic details become relevant when the Debye length is comparable to the characteristic length scale.
The electrokinetic and LD approach were unified by Schnitzer & Yariv (Reference Schnitzer and Yariv2015) and Mori & Young (Reference Mori and Young2018). Mori & Young (Reference Mori and Young2018) consider the weak electrolyte limit; Schnitzer & Yariv (Reference Schnitzer and Yariv2015) consider the strong electrolyte limit and, following Baygents & Saville (Reference Baygents and Saville1990), assume the existence of an imbalance in surface concentrations of anions and cations, thereby solving an additional surface charge conservation equation despite fully resolving the Debye layer. Mori & Young (Reference Mori and Young2018) claim that this conservation equation is an arbitrary requirement without justification. Despite these differences, both analyses employ rigorous scaling arguments to extend the earlier work by Baygents & Saville (Reference Baygents and Saville1990) and Saville (Reference Saville1997), and demonstrate that the LD model can be derived from the more general electrokinetic formulation under certain limiting conditions (see § 2). However, López-Herrera et al. (Reference López-Herrera, Herrada and Gañán-Calvo2023) have recently pointed out that cone-jets are well-described by the LD model despite not fulfilling one of these conditions. To the best of our knowledge, the work by López-Herrera et al. (Reference López-Herrera, Herrada and Gañán-Calvo2023) is the first investigation of cone-jets that includes a detailed electrokinetic formulation for both weak and strong electrolytes. This study compares LD and PNP solutions of cone-jets finding that, in the case of strong electrolytes, both electrokinetic effects and the polarity of the cone-jet strongly influence the total current. However, for these strong electrolytes, experiments agree much better with the LD solution than with the PNP solution. The authors hypothesize that the poor performance of the PNP model is due to the lack of a surface adsorption process in the model, as postulated by Schnitzer & Yariv (Reference Schnitzer and Yariv2015).
Highly concentrated solutions such as ionic liquids exhibit overscreening and overcrowding of ions near the interface (Kilic, Bazant & Ajdari Reference Kilic, Bazant and Ajdari2007a ; Kornyshev Reference Kornyshev2007; Bazant, Storey & Kornyshev Reference Bazant, Storey and Kornyshev2011; Storey & Bazant Reference Storey and Bazant2012). These effects lead to complex behaviours like charge density oscillations (Hayes, Warr & Atkin Reference Hayes, Warr and Atkin2010; Perkin Reference Perkin2012; Li, Endres & Atkin Reference Li, Endres and Atkin2013; Smith, Lee & Perkin Reference Smith, Lee and Perkin2016; Noh & Jung Reference Noh and Jung2019). These phenomena occur because the charges at the interface are overcompensated by counter-ions, which in turn are overcompensated by additional ions, creating a non-monotonic charge distribution close to the interface (Bazant et al. Reference Bazant, Storey and Kornyshev2011). Several modifications to the PNP model, which is strictly applicable to dilute solutions, have been proposed (Kornyshev Reference Kornyshev2007; Kilic et al. Reference Kilic, Bazant and Ajdari2007a , Reference Kilic, Bazant and Ajdarib ; Bazant et al. Reference Bazant, Storey and Kornyshev2011; Storey & Bazant Reference Storey and Bazant2012). These modifications account for the finite size of the ion and short-range Coulomb forces. These models have been used to describe the behaviour of ionic liquids (Lee et al. Reference Lee, Kondrat, Vella and Goriely2013, Reference Lee, Im and Kang2015; Stout & Khair Reference Stout and Khair2014; Wang et al. Reference Wang, Bao, Pan and Sun2017), but have not been applied to cone-jets of ionic liquids, which have only been studied with the LD formulation (Coffman et al. Reference Coffman, Martínez-Sánchez, Higuera and Lozano2016, Reference Coffman, Martínez-Sánchez and Lozano2019; Gallud & Lozano Reference Gallud and Lozano2022; Magnani & Gamero-Castaño Reference Magnani and Gamero-Castaño2023, Reference Magnani and Gamero-Castaño2024).
Building on the work by López-Herrera et al. (Reference López-Herrera, Herrada and Gañán-Calvo2023), this article revisits the electrokinetic modelling of cone-jets of fully dissociated salts. This is partially due to our own research, which focuses on cone-jets of highly conducting liquids that necessarily involve strong electrolytes. An important objective is to analyse the conditions under which the electrokinetic and LD formulations coincide and depart. Section 2 analyses the ranges of ion diffusivities and Debye lengths realized in experiments. Section 3 introduces the PNP and modified electrokinetic (MEK) models and the numerical implementation. The models account for self-heating and non-isothermal behaviour. Section 4 analyses the solutions of the LD, PNP and MEK models for dilute electrolytes and ionic liquids. Finally, a summary of the findings is presented in § 5.
2. Applicability of the LD model to cone-jets
The implementation of charge transport is the main difference between LD and electrokinetic formulations. While the LD model does not consider a volumetric charge density and employs Ohm’s law in the bulk
\begin{equation} \tilde {\boldsymbol{\nabla }}\boldsymbol{\cdot }(K \boldsymbol{\tilde E})=0, \end{equation}
the PNP formulation tracks the concentrations of anions and cations, defines a volumetric charge density
$\tilde {\rho }_{e}$
with their balance and requires the electric potential in the bulk to fulfil the Poisson equation,
\begin{align} \frac {\partial \tilde {n}^{\pm }}{\partial \tilde {t}}+\boldsymbol{\tilde {\nabla }} \boldsymbol{\cdot } (\tilde {n}^{\pm } \boldsymbol{\tilde {u}})&=\boldsymbol{\tilde {\nabla }}\boldsymbol{\cdot }\left (D^{\pm }\boldsymbol{\tilde {\nabla }}\tilde {n}^{\pm }\pm \frac {D^{\pm }}{k_{B}\tilde {T}}e\tilde {n}^{\pm }\boldsymbol{\tilde {\nabla }}\tilde {\phi }\right ), \notag\\ \varepsilon \varepsilon _{0}\tilde {\nabla }^{2}\tilde {\phi }&=-\tilde {\rho }_{e}, \quad \tilde {\rho }_{e}=e(\tilde {n}^{+}-\tilde {n}^{-}), \end{align}
where
$e$
,
$k_B$
,
$\tilde {T}$
are the elementary charge, Boltzmann constant and the temperature;
$\boldsymbol{\tilde E}$
,
$\tilde {\phi }$
and
$\varepsilon$
stand for the electric field, electric potential and the dielectric constant of the liquid; and
$D^{\pm }$
and
$\tilde {n}^{\pm }$
are the diffusivity and the concentration of the anion (−) and cation (+). Throughout the article dimensional variables are capped with a tilde, while dimensionless variables are uncapped. Baygents & Saville (Reference Baygents and Saville1990) and Schnitzer & Yariv (Reference Schnitzer and Yariv2015) have shown that the LD model is a limit of the more general electrokinetic formulation under the following conditions:
-
(i) the characteristic length (
$r_c$
) of the liquid domain is much larger than the Debye length (
$\lambda _D$
):
$\Lambda _{D}=({r_{c}}/{\lambda _{D}})\gg 1$
; -
(ii) the local electric potential at the interface of the liquid domain (
$\phi _s$
) is much larger than the thermal potential (
$\phi _{Th}=({k_{B}\tilde {T}_{o}})/{e}$
):
$\Gamma =({\phi _{s}})/({\phi _{Th}})\gg 1$
; -
(iii) the electric field outside the fluid domain (
$\phi _s/r_c)$
is much smaller than the electric field in the Debye layer (
$\phi _{Th}/\lambda _{D}$
):
$\Lambda _{D}/\Gamma \gg 1$
.
We next analyse whether these conditions, referred to as the Baygent–Saville limit (Schnitzer & Yariv Reference Schnitzer and Yariv2015; Mori & Young Reference Mori and Young2018; López-Herrera et al. Reference López-Herrera, Herrada and Gañán-Calvo2023) are fulfilled in typical cone-jets. For simplicity we assume a dilute electrolyte of a fully dissociated
$1:1$
salt. The Debye length is given by
\begin{equation} \lambda _{D}=\left (\frac {\varepsilon \varepsilon _{0}k_{B}\tilde {T}}{\tilde {n} e^2}\right )^{1/2}=\left (\frac {\varepsilon \varepsilon _{0}(D^{+}+D^{-})}{K}\right )^{1/2}, \end{equation}
where
$\tilde {n}$
is the salt concentration. The electrical conductivity is proportional to the ion concentrations and diffusivities,
\begin{equation} K=\frac {e^2}{k_{B}\tilde {T}}(\tilde {n}^{+}D^{+}+\tilde {n}^{-}D^{-}). \end{equation}
While the diffusion coefficients for various salts dissolved in water and for some ionic liquids are known (Lide Reference Lide2004; Tokuda et al. Reference Tokuda, Tsuzuki, Susan, Hayamizu and Watanabe2006), extrapolating these values to other solvents is not generally possible. However, for strong electrolytes frequently employed in electrosprays, we can estimate the sum of the diffusion coefficients with the reported bulk electrical conductivity and concentration,
\begin{equation} D^{+}+D^{-}\approx \frac {k_{B}\tilde {T}}{e^2}\left (\frac {K}{\tilde {n}}\right ). \end{equation}
Figure 1. Conductivity versus molar concentration for solutions of lithium chloride (LiCl) and 1-ethyl-3-methylimidazolium bis((trifluoromethyl)sulphonyl)imide (EMI-Im) salts in water (W), ethylene glycol (EG), propylene carbonate (PC) and tributyl phosphate (TBP). Here
$N_A$
refers to Avogadro’s number.
Figure 1 shows the bulk electrical conductivity versus the bulk salt concentration for several electrolytes, and table 1 lists the average diffusion coefficients obtained from linear fittings of these data. The diffusivity estimates in table 1 are rough approximations, based on the assumption of full salt dissociation, an assumption that may not hold in low dielectric solvents. The diffusion coefficients of the ionic liquids EMI-Im, 1-ethyl-3-methylimidazolium tetrafluoroborate (EMI-BF4) and ethylammonium nitrate (EAN) in table 1 are reported by Noda, Hayamizu & Watanabe (Reference Noda, Hayamizu and Watanabe2001), Tokuda et al. (Reference Tokuda, Tsuzuki, Susan, Hayamizu and Watanabe2006) and Filippov et al. (Reference Filippov, Alexandrov, Gimatdinov and Shah2021), respectively.
Table 1. Estimates of diffusion coefficients for solutions of lithium chloride and EMI-Im salts in water (W), PC, EG and TBP along with pure ionic liquids EMI-Im, EMI-BF4 and EAN.
The appropriate length and electric potential needed to evaluate
$\Lambda _{D}$
and
$\Gamma$
are the radius of the jet and the radius times the normal component of the outward electric field on the surface, both of which vary along the axis. We will next use the scaling laws of cone-jets to express these quantities (Fernández de La Mora & Loscertales Reference Fernández de La Mora and Loscertales1994; Gañán-Calvo et al. Reference Gañán-Calvo, López-Herrera, Herrada, Ramos and Montanero2018). These scaling laws are strictly applicable to isothermal cone-jets, and therefore the physical properties defining the different dimensionless groups are constant. However, our cone-jet model incorporates self-heating and temperature variations, and considers temperature-dependent electrical conductivity, viscosity (
$\mu$
) and ion diffusivities (all other physical properties are considered constant). Thus, we will use the values of these physical properties at the inlet temperature
$\tilde {T}_o$
, that is
$\mu _o$
,
$K_o$
and
$D^{\pm }_o$
, when defining the different dimensionless groups. The characteristic scales for the length and electric potential are
\begin{align}&\qquad\qquad\qquad\qquad\qquad r_{c} = d_0\Pi _{Q}^{1/2}, \end{align}
\begin{align}& \phi _{s}= \frac {\sigma }{\varepsilon _{0}}r_{c} = \frac {I r_{c}^2}{2 \varepsilon _{0} Q} = \left (\frac {\rho ^2 \gamma K_o Q^3}{2^6 \varepsilon _{0}^4} \right )^{1/6} =\left (\frac { \gamma ^4}{2^6 \varepsilon _{0} \rho K_o^2 } \right )^{1/6} \Pi _Q^{1/2}, \end{align}
where
$\sigma$
stands for the surface charge density, and
$d_0=(\gamma \varepsilon _{0}^{2}/(\rho K{_o}^2))^{1/3}$
is a characteristic length that only depends on physical properties. Note that
$r_c$
and
$\phi _s$
are expressed in terms of the dimensionless flow rate
$\Pi _{Q}$
,
\begin{equation} \Pi _{Q}=\frac {\rho K_o Q}{\varepsilon _{0}\gamma }. \end{equation}
We do this because in the traditional description of cone-jets, i.e. excluding electrokinetic effects, the state of an isothermal cone-jet is specified by three dimensionless numbers. The typical choices for this triplet are
$\Pi _{Q}$
,
$\varepsilon$
and the so-called electric Reynolds number
${Re}$
,
\begin{equation} {Re}=\left (\frac {\varepsilon _{0}\rho \gamma ^{2}}{ \mu {_o}^{3}K_o}\right )^{1/3}. \end{equation}
Finally,
$\Lambda _D$
and
$\Gamma$
can be expressed in terms of these dimensionless numbers as
\begin{align}& \Lambda _{D}= \left [\frac {\varepsilon _{0}\gamma ^2}{ \rho ^2K_{o}\left(D_{o}^{+}+D_{o}^{-}\right)^3}\right ]^{1/6} \left (\frac {\Pi _{Q}}{\varepsilon }\right )^{1/2}, \end{align}
\begin{align}&\qquad\quad \Gamma = \frac {e\, \gamma ^{2/3}}{2 k_{B}\tilde {T}_{o}\,\varepsilon _{0}^{1/6} K_{o}^{1/3}\rho ^{1/6}} \Pi _Q^{1/2}. \end{align}
These expressions are excellent approximations for
$\Lambda _{D}$
and
$\Gamma$
at the point where the conduction current through the bulk of the liquid equals the surface current, i.e. in the initial region of the jet. The local values of
$\Lambda _{D}$
and
$\Gamma$
decrease downstream as the jet thins down.
Figure 2. Cone-jet states reproduced in experiments (TBP, PC and EG-based solutions with various electrical conductivities, and ionic liquids EMI-Im and EAN) in
$\Lambda _{D}-\Gamma$
space. The label next to each data series indicates the Reynolds number of the solution.
Figure 2 shows states of experimental cone-jets in
$\Lambda _{D} - \Gamma$
space, for several liquids including TBP (
$\varepsilon =8.91$
), EG (
$\varepsilon =37.6$
), PC (
$\varepsilon =64.9$
) (Gamero-Castaño & Magnani Reference Gamero-Castaño and Magnani2019a
; Gamero-Castaño Reference Gamero-Castaño2019) and the pure ionic liquids EMI-Im (
$\varepsilon =13.05$
) and EAN (
$\varepsilon =27$
) (Caballero-Pérez & Gamero-Castaño Reference Caballero-Pérez and Gamero-Castaño2025). States along a data series correspond to cone-jets of a solution with fixed conductivity and varying flow rate, with a Reynolds number indicated by a label. Although all these cone-jets fulfil the first two Baygent–Saville conditions,
$\Lambda _{D} \gg 1$
and
$\Gamma \gg 1$
, none of them satisfies the third condition
$\Lambda _{D}/\Gamma \gg 1$
. Therefore, and as already pointed out by López-Herrera et al. (Reference López-Herrera, Herrada and Gañán-Calvo2023), the LD model should not be suitable for these cone-jets. However, numerical solutions of the LD model accurately reproduce observable features such as the current and the jet diameter (Hartman et al. Reference Hartman, Brunner, Camelot, Marijnissen and Scarlett1999; Higuera Reference Higuera2003; Herrada et al. Reference Herrada, López-Herrera, Gañán-Calvo, Vega, Montanero and Popinet2012; Gamero-Castaño Reference Gamero-Castaño2019; Gamero-Castaño & Magnani Reference Gamero-Castaño and Magnani2019a
). Section 4 will show that the criterion
$\varepsilon \Lambda _{D}^{2}/\Gamma \gt 1$
instead of
$\Lambda _{D}/\Gamma \gg 1$
is the correct condition for the validity of the LD model in cone-jets.
2.1. Highly concentrated electrolytes
The PNP equations (2.2) are only valid for dilute solutions and do not account for short-range Coulomb interactions and the finite size of ions. Highly concentrated electrolytes such as ionic liquids exhibit overcrowding and overscreening effects leading to oscillations in the volumetric charge density (Bazant et al. Reference Bazant, Storey and Kornyshev2011; Gebbie et al. Reference Gebbie, Dobbs, Valtiner and Israelachvili2015; Smith et al. Reference Smith, Lee and Perkin2016). Due to these effects, the Debye layer in concentrated electrolytes is substantially larger than predicted by (2.3). For example, (2.3) yields
$\lambda _D \approx 0.1$
nm for EMI-Im, which is much smaller than the average ion diameter, approximately 0.75 nm (Sun et al. Reference Sun, Zhuo, Chen, Du, Zhang and Wang2022), and therefore unphysical. On the contrary, the Debye length of EMI-Im evaluated in experiments is approximately 7 nm (Smith et al. Reference Smith, Lee and Perkin2016).
The PNP equations and (2.3) are only applicable when the separation between the ions is much greater than the Bjerrum length
$l_{B}$
(Lee et al. Reference Lee, Perez-Martinez, Smith and Perkin2017),
\begin{equation} \delta _{l_{B}}=l_{B}^{3}\tilde {n}\lt \lt 1, \end{equation}
where
$l_{B}=e^{2}/(4\pi \varepsilon \varepsilon _{0} k_{B}\tilde {T})$
is the distance at which the ion–ion interaction energy equals the thermal energy. Ionic liquids such as EMI-Im and EAN, with
$\delta _{l_{B}}\sim {\mathcal{O}}(10^2)$
, do not fulfil (2.12).
Kilic et al. (Reference Kilic, Bazant and Ajdari2007a , Reference Kilic, Bazant and Ajdarib ) and Bazant et al. (Reference Bazant, Storey and Kornyshev2011) have proposed modifications to the PNP equations to account for crowding and screening effects. These equations capture the spatial oscillations of the volumetric charge density observed near charged interfaces in ionic liquids (Lee, Im & Kang Reference Lee, Im and Kang2013; Wang et al. Reference Wang, Bao, Pan and Sun2017). The MEK read
\begin{align} \frac {\partial \tilde {n}^{\pm }}{\partial \tilde {t}}+\boldsymbol{\tilde {\nabla }} \boldsymbol{\cdot } (\tilde {n}^{\pm } \boldsymbol{\tilde {u}})&=\boldsymbol{\tilde {\nabla }}\boldsymbol{\cdot }\left [D^{\pm }\boldsymbol{\tilde {\nabla }}\tilde {n}^{\pm }\pm \frac {D^{\pm }}{k_{B}T}e\tilde {n}^{\pm }\boldsymbol{\tilde {\nabla }}\tilde {\phi }+\frac {a^{3}D^{\pm }n^{\pm }\boldsymbol{\tilde {\nabla }}(\tilde {n}^{+}+\tilde {n}^{-})}{1-a^{3}(\tilde {n}^{+}+\tilde {n}^{-})}\right ], \notag\\&\quad \varepsilon \varepsilon _{0}\big(1-l_{c}^{2}\tilde {\nabla }^{2}\big)\tilde {\nabla }^{2}\tilde {\phi }=-\tilde {\rho }_{e}, \quad \quad \tilde {\rho }_{e}=e(\tilde {n}^{+}-\tilde {n}^{-}). \end{align}
The last term in the ion conservation equation accounts for the finite ion radius
$a$
, and its role on the overcrowding near the charged interface. The Poisson equation for the electric potential is replaced by a fourth-order differential equation where ion–ion correlation effects are accounted for by the electrostatic correlation length
$l_{c}$
. Here
$l_{c}$
is bounded by the ion radius and the Bjerrum length,
$a \leqslant l_{c}\leqslant l_{B}$
(Bazant et al. Reference Bazant, Storey and Kornyshev2011), and we will use
$l_{c}=l_{B}$
in our calculations, unless stated otherwise.
Figure 3. Schematic of the computational domain. The cone-jet profile
$R(z)$
is fixed and computed with the LD model (Magnani & Gamero-Castaño Reference Magnani and Gamero-Castaño2024).
3. Governing equations and numerical implementation
Figure 3 illustrates the axisymmetric computational domain. The radius of the cone-jet,
$R(z)$
, is taken from the numerical solution of the LD model for a given
$ \{\varepsilon , \Pi _Q, {Re} \}$
state (Magnani & Gamero-Castaño Reference Magnani and Gamero-Castaño2024), and is kept fixed. The surface separates the liquid phase
$\Omega ^{i}$
from the surrounding vacuum
$\Omega ^{o}$
. The electrostatic Taylor potential
$\phi _{T}$
is used as the far-field boundary condition on the surface
$\Gamma _{1}^{o}\cup \Gamma _{2}^{i}$
, which allows modelling the physics without accounting for the geometry of the electrodes (Gañán-Calvo Reference Gañán-Calvo1997; Higuera Reference Higuera2003; Gamero-Castaño & Magnani Reference Gamero-Castaño and Magnani2019b
). The velocity vector
$\boldsymbol{u}$
, the pressure
$p$
, the temperature
$T$
, the electric potentials in the liquid
$\phi ^{i}$
and gas
$\phi ^{o}$
phases, and the concentrations of anions and cations
$n^{\pm }$
are the dependent variables of the MEK model. The governing equations are made dimensionless with characteristic scales for length
$r_{c}$
, velocity
$u_{c}$
, pressure
$p_{c}$
, temperature
$T_{c}$
, current
$I_{c}$
and concentration
$n_{c}$
, with derived scales for the electric field
$E_{c}$
, the electric potential
$\phi _{c}$
and power
$P_{c}$
,
\begin{align}& r_{c}=\left (\frac {\rho \varepsilon _{0}Q^{3}}{\gamma K_{o}}\right )^{1/6}\!, \!\quad u_{c}=\frac {Q}{\pi r_{c}^{2}}, \quad p_{c}=\frac {\gamma }{r_{c}}, \!\quad T_{c}=\frac {1}{\pi c}\left (\frac {\gamma K_{o}}{\varepsilon _{0}\rho }\right )^{2/3}\!, \!\quad I_{c}=(\gamma K_{o} Q)^{1/2}, \notag\\ &\qquad\qquad n_{c}=\frac {k_{B}\tilde {T}_{o}}{e^{2}}\frac {K_{o}}{D^{+}_{o}+D^{-}_{o}}, \quad E_{c}=\frac {I_{c}}{\pi r_{c}^{2}K_{o}}, \quad \phi _{c}=r_{c}E_{c}, \quad P_{c}=\phi _{c}I_{c}, \end{align}
where
$c$
stands for the specific heat capacity of the liquid. It should be noted that the electric potential scale
$\phi _c$
differs from
$\phi _s$
introduced in § 2. Specifically,
$\phi _c$
is based on the tangential electric field, while
$\phi _s$
is based on the outward normal electric field, as defined in (2.7). In the non-dimensionalization of the model equations, the scale
$\phi _c$
helps to make the dimensionless variables of order one. The dimensionless numbers parametrizing the model are listed in table 2. The model solves the continuity equation, the equations of conservation of momentum and energy, the modified Poisson’s equation for the electric potential in the liquid, the Laplace’s equation for the potential in the surrounding vacuum and the modified ion conservation equations,
\begin{align}&\qquad\qquad\qquad\qquad\qquad\qquad\quad \boldsymbol{\nabla } \boldsymbol{\cdot} \boldsymbol{u}=0, \end{align}
\begin{align}&\,\,\, \boldsymbol{u} \cdot \boldsymbol{\nabla }\boldsymbol{u} =-\frac {\pi ^{2}\boldsymbol{\nabla }p}{\sqrt {\Pi _{Q}}} + \frac {\pi \mu \boldsymbol{\nabla } \cdot (\boldsymbol{\nabla u} + \boldsymbol{\nabla u}^{T})}{\mu _{o}{Re}\sqrt {\Pi _{Q}}} - \frac { \pi \varepsilon \Lambda _{D}^{2}}{2\Gamma \sqrt {\Pi _{Q}}}(n^{+}-n^{-})\boldsymbol{\nabla }\phi ^{i},\\[-10pt] \nonumber \end{align}
\begin{align}&\qquad\quad\,\,\, \boldsymbol{u} \cdot \nabla T =\frac {\nabla ^{2} T}{Pe_{Th}}+\frac {\mu (\boldsymbol{\nabla u} + \boldsymbol{\nabla u}^{T}) : \nabla \boldsymbol{u}}{\mu _{o}{Re}\sqrt {\Pi _{Q}}} + \frac {K\boldsymbol{\nabla }\phi ^{i}\cdot \boldsymbol{\nabla }\phi ^{i}}{K_{o}\sqrt {\Pi _{Q}}}, \end{align}
\begin{align}&\qquad (1-\delta _{c}^{2}\nabla ^{2})\nabla ^{2}\phi ^{i}=-\frac {\pi \Lambda _{D}^{2}}{2\Gamma }\sqrt {\Pi _{Q}}(n^{+}-n^{-}), \quad \quad \quad \nabla ^{2}\phi ^{o}=0, \end{align}
\begin{align}& \boldsymbol{u} \cdot \boldsymbol{\nabla }n^{\pm }=\frac {1}{Pe_{D}^{\pm }}\boldsymbol{\nabla } \cdot \frac {D^\pm }{D^{\pm }_o}\left (\boldsymbol{\nabla }n^{\pm }\pm \frac {2\Gamma }{\pi \sqrt {\Pi _{Q}}} n^{\pm }\boldsymbol{\nabla }\phi ^{i}+\frac {\nu _{a}n^{\pm }\boldsymbol{\nabla }(n^{+}+n^{-})}{2-\nu _{a}(n^{+}+n^{-})}\right ). \end{align}
The last two terms in (3.4) are the viscous and ohmic power dissipation densities,
\begin{align} \tilde {P}_{\mu }^{\prime \prime \prime }&=\mu (\boldsymbol{\nabla \tilde {u}} + \boldsymbol{\nabla \tilde {u}}^{T}) : \nabla \boldsymbol{\tilde {u}}, \notag \\ \tilde {P}_{\Omega }^{\prime \prime \prime }&=K\boldsymbol{\tilde {E}}^{i}\cdot \boldsymbol{\tilde {E}}^{i}. \end{align}
Table 2. Dimensionless numbers governing non-isothermal cone-jets. The dimensionless numbers are represented in terms of the state variables
$\Pi _{Q}$
and
$\varepsilon$
, which are typically used to characterize the state of cone-jets. Dimensionless groups are defined with the viscosity, electrical conductivity and ion diffusivities at the inlet temperature,
$\mu _{o}$
,
$K_{o}$
,
$D^{\pm }_{o}$
,
$\tilde {T}_{o}$
. Here
$\kappa$
stands for the thermal conductivity of the liquid.
In order to study the effect of varying anion/cation diffusivities at constant bulk electrical conductivity, we define the relative cation diffusivity,
$D^{+}_{r}$
,
\begin{equation} D^{+}_{r} = \frac {D^{+}}{D^{+} + D^{-}}, \end{equation}
where
$D^{+}_{r}=0.5$
corresponds to equal cation and anion diffusivities. The temperature-dependent ion diffusivity and viscosity are exponential functions of the temperature given by
\begin{align} \frac {D^{\pm }(\tilde {T})}{D^{\pm }_{o}}&=Y_{D}\tilde {T}e^{\frac {B_{D}}{\tilde {T}-T_{D}}}, \notag \\ \frac {\mu (\tilde {T})}{\mu _{o}} &=Y_{\mu }e^{\frac {B_{\mu }}{\tilde {T}-T_{\mu }}}, \end{align}
where the constants
$Y_{D}$
,
$B_{D}$
,
$T_{D}$
,
$Y_{\mu }$
,
$B_{\mu }$
,
$T_{\mu }$
are specific to the liquid and are given in table 3. We compute them using published data of the dependence of the viscosity and the electrical conductivity on temperature, (2.4) and fixing
$D_{r}^{+} = 0.5$
. The reason for using
$D^{\pm }\propto \tilde {T}\exp (B_{D}/(\tilde {T}-T_{D}))$
is the temperature dependence in (2.4) and the fact that electrical conductivity is well approximated by an exponential function of temperature (Magnani & Gamero-Castaño Reference Magnani and Gamero-Castaño2024; Tokuda et al. Reference Tokuda, Tsuzuki, Susan, Hayamizu and Watanabe2006). The coefficients for EG and PC are obtained from Magnani & Gamero-Castaño (Reference Magnani and Gamero-Castaño2024), while Appendix A describes the calculations of the coefficients for the ionic liquids.
Table 3. Coefficients for computing temperature-dependent viscosities and ion diffusivities, (3.9). The reference temperatures are 298 K for EG and PC, and 294 K for EMI-Im and EAN.
At the surface we enforce the surface kinematic condition, the balance of tangential stresses, the continuity of the electric potential and the vanishing of the ion correlation term, the customary jump of the electric displacement field across the interface between two dielectrics and the zero heat and ion fluxes,
\begin{align}&\qquad\qquad\qquad\qquad\qquad \boldsymbol{u} \cdot \boldsymbol{n}=0, \end{align}
\begin{align}&\qquad\qquad\qquad\quad \boldsymbol{t}\cdot (\boldsymbol{\nabla u}+\boldsymbol{\nabla u}^{T})\boldsymbol{n}=0, \end{align}
\begin{align}&\qquad\qquad\quad \phi ^{i}=\phi ^{o}\!, \quad \boldsymbol{n}\cdot \delta _{c}^{2}\boldsymbol{\nabla }(\nabla ^{2}\phi ^{i})=0, \end{align}
\begin{align}&\qquad\qquad\qquad\qquad\,\,\, E_n^{o}-\varepsilon E_{n}^{i}=0, \end{align}
\begin{align}&\qquad\qquad\qquad\qquad\quad\, \boldsymbol{n}\cdot \boldsymbol{\nabla }T=0, \end{align}
\begin{align}& \left (\boldsymbol{\nabla }n^{\pm }\pm \frac {2\Gamma }{\pi \sqrt {\Pi _{Q}}} n^{\pm }\boldsymbol{\nabla }\phi ^{i}+\frac {\nu _{a}n^{\pm }\boldsymbol{\nabla }(n^{+}+n^{-})}{2-\nu _{a}(n^{+}+n^{-})}\right ) \boldsymbol{\cdot} \boldsymbol{n}=0. \end{align}
Note that (3.13) implies zero surface charge, in agreement with Mori & Young (Reference Mori and Young2018) and in contrast with the formulations by Saville (Reference Saville1997) and Schnitzer & Yariv (Reference Schnitzer and Yariv2015). Furthermore, we do not require fulfilment of the balance of normal stresses across the surface because the position of the surface is fixed, adopted from the LD solution. The main goal of this article is to determine the regions in
$\Lambda _D - \Gamma$
space in which the electrokinetic and LD formulations coincide or diverge, and this is best accomplished by computing and comparing solutions at varying
$\Lambda _D$
and
$ \Gamma$
while keeping all other parameters, including the shape of the cone-jet, constant. Future work will focus on solving self-consistent electrokinetic formulations, i.e. models that also fulfil the balance of normal stresses on the surface. The Taylor potential is given by
\begin{equation} \phi _{T}=- 4.228\Pi _{Q}^{1/4}\big(z^{2}+r^{2}\big)^{1/4}P_{1/2}\left (\frac {z}{\sqrt {z^{2}+r^{2}}}\right ), \end{equation}
where
$P_{1/2}$
is the Legendre function of the first kind of degree
$1/2$
. At the inlet surface
$\Gamma _{1}^{i}$
we prescribe a sink flow with net flow rate
$Q$
, and values for the ion densities, the electric potential and temperature,
\begin{align} \boldsymbol{u} ( \Gamma _{1}^{i})=-\frac {1}{2\left(1-\cos\theta _{T}\right)\left(z^{2}+r^{2}\right)^{3/2}}\boldsymbol{r}, \, n^{\pm } (\Gamma _{1}^{i}) =1, \, \phi ^{i} (\Gamma _{1}^{i})=0,\, T(\Gamma _{1}^{i})=T_{o}, \end{align}
where
$\boldsymbol{r}$
is the position vector. Note that
$Q$
does not appear explicitly in the dimensionless boundary condition due to its presence in the definition of the characteristic velocity. At the outlet boundary
$\Gamma _{2}^{i}$
we impose zero velocity, ion and heat fluxes and the Taylor potential
\begin{align} \boldsymbol{n}\cdot \boldsymbol{\nabla }\boldsymbol{u}(\Gamma _{2}^{i})=0, \quad \boldsymbol{n}\cdot \boldsymbol{\nabla }n^{\pm }(\Gamma _{2}^{i})=0, \quad \boldsymbol{n}\cdot \boldsymbol{\nabla }T (\Gamma _{2}^{i})=0, \quad \phi ^{i} (\Gamma _{2}^{i})=\phi _{T}. \end{align}
The system of equations is solved with COMSOL Multiphysics. For the mass and momentum conservation equations, we use the inbuilt steady-state laminar flow module. The conservation of energy is solved using the steady-state heat transfer module, with additional source terms accounting for viscous and ohmic dissipation. We use the partial differential equation interface to define our own modules for solving the Poisson/Laplace equations for the electric potential and the conservation equations for ionic species. The electric potential is discretized using cubic-order Lagrange elements due to the modified fourth-order Poisson equation (3.5) in the formulation. The higher-order nature of the equation benefits from the increased accuracy using cubic elements, particularly near sharp potential gradients within the Debye layer. The velocity, temperature and ionic concentrations are discretized using quadratic-order Lagrange elements. Quadratic elements are well-suited for resolving the gradients of ionic concentrations, while avoiding excessive computational cost. The pressure is discretized with linear Lagrange elements. The resulting system of equations is solved using the MUMPS (multifrontal massively parallel sparse) solver, which is well-suited for large, sparse matrices. A relative tolerance of
$10^{- 4}$
was imposed on the residuals of all dependent variables: the electric potential, ionic concentrations, velocity, temperature and pressure. To accurately resolve the Debye layer near the liquid–gas interface, an adaptive mesh refinement strategy is employed, as depicted in figure 4. The refinement is based on a gradient criterion of the form
$|\nabla n^{\pm }|h\lt \epsilon$
, where
$h$
is the local mesh size, and
$\epsilon$
is a user-defined threshold. This criterion ensures that regions exhibiting strong gradients in the ion concentration, particularly near the interface, are locally refined. Convergence was also verified by monitoring the total current along the axis. The total current, which includes contributions from conduction, convection and diffusion, was observed to remain constant, indicating that the solution was well resolved and grid-independent.
Figure 4. Example of adaptive grid used to resolve the Debye layer.
We use the model presented above to simulate cone-jets of both ionic liquids and dilute electrolytes. For the latter, we make the electrostatic correlation length
$\delta _{c}$
and the ion volume fraction
$\nu _{a}$
equal to zero, thus recovering the standard PNP formulation (2.2). Although the article presents several LD solutions, the focus is on the electrokinetic modelling of cone-jets. Thus, we skip a detailed description of the LD model which can be found elsewhere (Gamero-Castaño & Magnani Reference Gamero-Castaño and Magnani2019b
; Magnani & Gamero-Castaño Reference Magnani and Gamero-Castaño2024), and simply point out that although the convergence of the solving scheme depends on the initial position of the surface, the errors in all equations of the LD model tend to zero. This supports the expectation that the numerical LD solution, which includes the position of the surface, corresponds to the unique physical solution.
4. Results
4.1. Liquids with low and medium electrical conductivity,
$\delta _{c} \ll 1$
,
$\nu _{a} \ll 1$
We investigate the PNP solution,
$\delta _{c}=\nu _{a}=0$
in (3.5) and (3.6), and compare it with the LD solution. We study two liquids with different dielectric constants, EG and PC, at varying Reynolds number, dimensionless flow rate and dimensionless Debye length. Unless otherwise specified, we assume equal cation and anion diffusivities,
$D^{+}_{r}=0.5$
. Table 4 lists the physical properties and dimensionless parameters.
Table 4. Physical properties and dimensionless parameters as a function of
${Re}$
,
$\Pi _{Q}$
and
$\Lambda _{D}$
for EG and PC solutions (
$\tilde {T}_{o}=298$
K).
Figure 5. Comparison between the PNP and LD solutions: (a) EG at
${Re}=0.095$
,
$\Pi _{Q}=64$
and
$\Lambda _{D}=35.0$
; and (b) PC at
${Re}=0.38$
,
$\Pi _{Q}=100$
and
$\Lambda _{D}=24.5$
. The value of the ion diffusivity for both cases is
$D^{+}_{o}+D^{-}_{o}=2\times 10^{-9}\,\rm m^{2}\,s^-{^1}$
. For plotting purposes the origin of the axial coordinate is placed at the maximum of
$R''(z)$
(Gamero-Castaño & Magnani Reference Gamero-Castaño and Magnani2019b
).
Figure 6. Cone-jet, velocity streamlines (cyan), equipotential lines (rainbow) and variation of ionic concentrations and relative electrical conductivity at three different axial positions (
$I_{b}=0.95I$
,
$I_{b}=I_{s}$
and
$I_{s}=0.95I$
) for: (a) EG,
${Re}=0.095$
,
$\Pi _{Q}=64$
,
$\Lambda _{D}=35.0$
; and (b) PC,
${Re}=0.38$
,
$\Pi _{Q}=100$
,
$\Lambda _{D}=24.5$
. The value of ion diffusivity for both cases is
$D^{+}_{o}+D^{-}_{o}=2\times 10^{-9}\,\rm m^{2}\,s^-{^1}$
. The electrical conductivity is normalized with its value at the axis,
$K_{b}$
.
Figure 5(a) compares the axial variation of several features of the numerical solution for EG at
${Re}=0.095$
,
$\Pi _{Q}=64$
and
$\Lambda _{D}=35.0$
(
$\Gamma =2575$
,
$Pe_{Th} = 39.9$
and
$Pe_{D}^{\pm } = 3675$
;
$\Lambda _{D}/\Gamma =0.0136$
), while figure 5(b) shows the comparison for PC,
${Re}=0.38$
,
$\Pi _{Q}=100$
and
$\Lambda _{D}=24.5$
(
$\Gamma =1864$
,
$Pe_{Th} = 30$
and
$Pe_{D}^{\pm } = 2521$
;
$\Lambda _{D}/\Gamma =0.013$
). The ratios
$\Lambda _{D}/\Gamma$
are much smaller than unity. Since this violates condition (iii) of the Baygent–Saville limit, the PNP and LD solutions are expected to be significantly different. The conduction and convection currents,
$I_{b}$
and
$I_{s}$
, are given by
\begin{align}& \tilde {I}_{b}(\textit{PNP})\!=2\pi\!\! \int _{0}^{\tilde {R}(z)}\!\!\frac {e^2}{k_{B}T}(D^{+}\tilde {n}^{+}+D^{-}\tilde {n}^{-})\tilde {E}_{z}^{i}\tilde {r}\,{\rm d}\tilde {r}, \,\, \tilde {I}_{b}(LD)=2\pi\!\! \int _{0}^{\tilde {R}(z)}\!K\tilde {E}_{z}^{i}\tilde {r}\,{\rm d}\tilde {r}, \end{align}
\begin{align}&\qquad\quad \tilde {I}_{s}(PNP)=2\pi \int _{0}^{\tilde {R}(z)}e(\tilde {n}^{+}-\tilde {n}^{-})\tilde {u}_{z}\tilde {r}\,{\rm d}\tilde {r}, \quad \tilde {I}_{s}(LD)=2\pi \tilde {R}(z)\tilde {\sigma }\tilde {u}_{s}, \end{align}
where
$\tilde {\sigma }$
and
$\tilde {u}_s$
are the surface charge density and the surface velocity in the LD solution. The diffusion current
$I_{d}$
only exists in the PNP formulation and is given by
\begin{equation} \tilde {I}_{d} (PNP)=2\pi \int _{0}^{\tilde {R}(z)} e\left (D^{-}\frac {\partial \tilde {n}^-}{\partial \tilde {z}}-D^{+}\frac {\partial \tilde {n}^+}{\partial \tilde {z}}\right )\tilde {r}\,{\rm d}\tilde {r}. \end{equation}
Conservation of charge requires the total current
$I=I_{b}(z)+I_{s}(z)+I_{d}(z)$
to be constant along the axis. The diffusion current is negligible for the conditions studied in figure 5, and the agreement between the conduction and the convection currents for both models is excellent. It is important to note that the conduction current is dominant in the Taylor cone, the convection current is dominant far downstream in the jet, and that the importance of each term abruptly reverses in a relatively short section at the base of the jet referred to as the transition region; the axial location where both currents are equal is known as the current crossover point. The electrokinetic model does not include the net surface charge postulated by Saville (Reference Saville1997) and Schnitzer & Yariv (Reference Schnitzer and Yariv2015). In our solutions the surface charge density in the LD model matches the integral of the PNP volumetric charge density
\begin{equation} \tilde {\sigma } (PNP)=\frac {1}{\tilde {R}(z) \sqrt {\tilde {R}'^2+1}} \int _{0}^{\tilde {R}(z)} e(\tilde {n}^{+}-\tilde {n}^{-})\tilde {r}\,{\rm d}\tilde {r}. \end{equation}
This is illustrated in the second panel, where
$\sigma (PNP)$
is indistinguishable from the surface charge density in the LD solution. This second panel also shows that both solutions have the same electric field in the outer side of the interface. The agreement in the normal electric field also highlights that the surface charge density in the LD model is equivalent to
$\sigma (PNP)$
and thus, for dilute solutions and within the assumptions of this work, there is no need to introduce an artificial surface charge in the electrokinetic formulation. While López-Herrera et al. (Reference López-Herrera, Herrada and Gañán-Calvo2023) hypothesized that surface charge accumulation or adsorption may play a role in certain regimes, none is required to explain the behaviour observed here. However, this does not rule out the possibility that interfacial charge adsorption may be important under conditions not captured by the present model. The third panel in figure 5 compares the ohmic
$P_{\Omega }^{\prime }$
and viscous
$P_{\mu }^{\prime }$
linear power densities
\begin{align} P_{\Omega }^{\prime }=2\pi \int _{0}^{R(z)}P_{\Omega }^{\prime \prime \prime }r\,{\rm d}r,\notag\\ P_{\mu }^{\prime }=\frac {2\pi }{{Re}}\int _{0}^{R(z)}P_{\mu }^{\prime \prime \prime }r\,{\rm d}r, \end{align}
and the total change in temperature
$\Delta T$
relative to the upstream temperature. The linear power densities and the total change in temperature computed by the LD and PNP models are again in excellent agreement. The ohmic power density has a broader distribution than the viscous power density (the latter is confined to a shorter section of the cone-jet where the radius of the surface rapidly changes). The maximum temperature increase in the jet is
$\Delta \tilde {T}= 1.0$
K for EG, and
$\Delta \tilde {T}= 2.5$
K for PC. Self-heating for these solutions is not significant due to their relatively large
${Re}$
numbers (Magnani & Gamero-Castaño Reference Magnani and Gamero-Castaño2024). In summary, figure 5 shows that the LD model is a valid electrokinetic limit for these cone-jets, despite the violation of the Baygent–Saville condition
$ \Lambda _{D}/\Gamma \gg 1$
.
Figure 6 shows the streamlines and equipotential lines for the PNP solutions in figure 5, together with the radial variation of the ionic concentrations and the electrical conductivity at three axial locations: upstream,
$I_b=0.95I$
; at the current crossover,
$I_{b}=I_{s}$
; and downstream,
$I_{s}=0.95I$
. The choice of radial variable,
$1-r/R(z)$
, highlights variations near the surface. In the upstream region the majority of the cross-section remains electrically neutral, and departure from neutrality only occurs within a thin Debye layer with a thickness of
$0.02R(z)$
for EG and
$0.03R(z)$
for PC. The thickness of the layer increases to
$0.11R(z)$
at the current crossover, and it becomes
$0.25R(z)$
for EG and
$0.3R(z)$
for PC at the downstream location. Relative to its value on the axis, the electrical conductivity shows a maximum variation of roughly
$40\,\%$
and
$10\,\%$
for EG and PC, respectively, but this occurs in the thin Debye layer and therefore its effect on the conduction current is negligible. Despite these variations, the PNP model exhibits excellent agreement with the LD model. The concentration and electrical conductivity profiles reported by López-Herrera et al. (Reference López-Herrera, Herrada and Gañán-Calvo2023) exhibit significantly larger changes in the Debye layer. For instance, under certain conditions, the authors observe relative electrical conductivities exceeding 100.
Figure 7. Flow recirculation pattern for PC,
$\Pi _Q = 16$
and
${Re} = 1.67$
. The sink flow boundary condition at the inlet is compatible with flow recirculation.
The streamlines shown in figure 6 displaying a sink flow may seem to indicate that the velocity boundary condition at the inlet (3.17) does not allow for the recirculation patterns often observed in cone-jets (Herrada et al. Reference Herrada, López-Herrera, Gañán-Calvo, Vega, Montanero and Popinet2012). However, the sink flow boundary conditions can lead to recirculation zones in the cone when the conditions require it. Like all other aspects of the LD solution, the formation of a recirculation zone is a function of the three dimensionless numbers parameterizing it, i.e.
$\Pi _Q, \, {Re} \text{ and } \varepsilon$
(the dependence of recirculation on
$\varepsilon$
is obviously negligible). For example, figure 7 shows a recirculation cell for PC,
$\Pi _Q = 16$
and
${Re} = 1.67$
, a case solved in the present article. On the other hand, the dimensionless flow rates and Reynolds numbers for the cone-jets in figure 6 are too high and too low, respectively, to induce recirculation. Note also that even if there were recirculation in the solutions shown in figure 6, it would not likely be observable because the plotted region does not extend far enough upstream. The particular form of the inlet velocity profile for a given flow rate does not affect significantly the cone-jet characteristics determined locally at the base of the jet and the transition region, as long as the inlet is placed sufficiently upstream. There is a large, quasihydrostatic region in the Taylor cone that can accommodate any particular inlet velocity condition.
Figure 8. Anion density, currents, electric field on the surface and surface charge, and concentration and conductivity profiles for EG,
${Re}=0.095$
,
$\Pi _{Q}=64$
and
$\Lambda _{D}=11.1, 14.5 \text{ and } 35.0$
(
$\varepsilon \Lambda _{D}^{2}/\Gamma =1.799, 3.07 \text{ and} 17.88$
). The corresponding ion diffusivities
$D^{+}_{o}+D^{-}_{o}$
are
$2\times 10^{-8}\,\rm m^{2}\,s^-{^1}$
,
$1.2\times 10^{-8}\,\rm m^{2}\,s^-{^1}$
and
$2\times 10^{-9}\,\rm m^{2}\,s^-{^1}$
.
Figure 8 compares the solution of the PNP model for three cone-jets with equal dielectric constant, dimensionless flow rate and Reynolds number, and different
$\Lambda _{D}$
ratios: EG,
${Re}=0.095$
,
$\Pi _{Q}=64$
and
$\Lambda _{D}=$
11.1, 14.5 and 35.0. In all cases
$\Gamma =2575$
, while the
$\varepsilon \Lambda _{D}^{2}/\Gamma$
ratios are 1.80, 3.07 and 17.9. Note that the three cone-jets share the same LD solution, which coincides with the PNP solution for
$\Lambda _{D}=$
35.0 (see figure 5
a). These three cone-jets can be recreated in experiments by dissolving different salts in EG, and adjusting the concentration inversely to the diffusion coefficients so that the bulk electrical conductivity remains constant. Figure 8(a i) shows contour plots of the anion concentration. For
$\Lambda _{D}=11.1$
the jet is depleted of anions along
$z\gtrsim 38$
. This possibility was discussed by Fernández de La Mora (Reference Fernández de la Mora2007), who suggested that complete charge separation could trigger the onset of the minimum flow rate. Figure 8(a ii) compares the conduction, convection and diffusion currents. Figure 5(a) showed agreement between the PNP and LD solutions for
$\Lambda _{D}=35.0$
, which is not the case for the two new and lower values of
$\Lambda _{D}$
. As
$\Lambda _{D}$
decreases the total current decreases, and the diffusion current becomes more significant in the upstream region. The latter trend was also noted by López-Herrera et al. (Reference López-Herrera, Herrada and Gañán-Calvo2023). The crossover between the conduction and convection currents occurs farther upstream for
$\Lambda _{D}=11.1$
, and immediately precedes complete charge separation. Figure 8(a iii) shows that, as
$\Lambda _{D}$
decreases, the maximum normal electric field decreases, while the maximum of the tangential electric field increases. Figure 8(b) shows ion concentration profiles along with the electrical conductivity. Interestingly, for
$\Lambda _{D} = 11.1$
and despite the small difference with respect to the
$\Lambda _{D} = 14.5$
case, the jet is depleted of anions downstream of the current crossover point. In both cases, the electrical conductivity exhibits a much stronger increase near the interface than for
$\Lambda _{D}=35.0$
.
Figure 9. Total current as a function of the dimensionless flow rate and the electric Reynolds number for PC cone-jets. Comparison between PNP solution (
$D^{+}_{o}+D^{-}_{o}=5.62\times 10^{-10}\,\rm m^{2}\,s^-{^1}$
for all cases), LD solution and experimental data (Gamero-Castaño Reference Gamero-Castaño2019).
Figure 9 compares the total currents computed with the PNP and LD models with experimental values (Gamero-Castaño Reference Gamero-Castaño2019), for cone-jets of PC as a function of the dimensionless flow rate and for two different Reynolds numbers,
${Re}=1.67$
and
${Re}=0.38$
. For the PNP model we use the values of diffusivity in table 1. When the current is normalized with the scale
$I_{o}=\sqrt {\varepsilon _{0}\gamma ^{2}/\rho }$
, the data follow the scaling law
$\tilde {I}/I_{o}\simeq 2.5\sqrt {\Pi _{Q}}$
(Fernández de La Mora & Loscertales Reference Fernández de La Mora and Loscertales1994; Gañán-Calvo et al. Reference Gañán-Calvo, López-Herrera, Herrada, Ramos and Montanero2018). The PNP solution coincides with the LD solution, and both are in good agreement with experimental results.
Figure 10. Ratio between the total currents obtained with the PNP and LD models as a function of
$\Lambda _{D}/\Gamma$
(a), and
$\varepsilon \Lambda _{D}^{2}/\pi \Gamma$
(b).
Figure 10 quantifies the conditions needed for the agreement between the LD and PNP solutions. We use the ratio of the total currents,
$I_{PNP}/I_{LD}$
, as figure of merit, which is computed for several values of
$\Lambda _{D}$
,
$\Pi _Q$
,
${Re}$
and
$\varepsilon$
. Figure 10(a) shows that the LD and PNP solutions coincide for values of
$\Lambda _{D}/\Gamma$
much smaller than unity, in disagreement with the third condition of the Baygent–Saville limit. On the other hand, figure 10(b) indicates that the LD model coincides with the PNP formulation for
\begin{equation} \frac {\varepsilon }{\pi }\frac {\Lambda _{D}^{2}}{\Gamma } \gtrsim 1. \end{equation}
Furthermore, regardless of the value, the
$I_{PNP}/I_{LD}$
ratios collapse relatively well on a single-valued function of
$({\varepsilon }/{\pi })({(\Lambda _{D}^{2})}/{\Gamma })$
. Therefore,
$({\varepsilon }/{\pi })({(\Lambda _{D}^{2})}/{\Gamma })\gtrsim 1$
, instead of
$({\Lambda _{D}}/{\Gamma }) \gg 1$
, is one of the required conditions for the application of the LD model in cone-jets. A more illuminating form of criterion (4.6a
) is obtained by writing
$\Lambda _{D}$
and
$\Gamma$
in terms of the dimensionless flow rate and dielectric constant,
\begin{equation} \frac {2\varepsilon _0 ^{1/2} k_{B}\tilde {T}}{\pi e \rho ^{1/2}(D^{+}+D^{-})}\Pi _{Q}^{1/2} \gtrsim 1. \end{equation}
Most dilute electrolytes fulfil this condition, and therefore their cone-jets are accurately described by the simpler LD model: the dimensionless flow rates of cone-jets are always of order unity or larger; and for typical electrosprays operated at room temperature, the fraction preceding
$\Pi _{Q}^{1/2}$
is larger than one except for electrolytes with very large ionic diffusivities. For example, the values of this factor for the TBP, EG and PC solutions in table 1 are
$27.2$
,
$15.5$
and
$2.5$
, respectively, while it is relatively small,
$0.56$
, only in the case of W + LiCl, which has the largest combined diffusivity. It is worth noting that
$\Lambda _{D}$
and
$\Gamma$
appear in the model (3.3) and (3.5) as
$\Lambda _{D}^2/\Gamma$
. Only (3.6), where the electric current term is proportional to
$\Gamma$
, exhibits a different dependency, but this term is divided by a large diffusion Peclet number. This explains the main
$\Lambda _{D}^2/\Gamma$
dependency of the
$I_{PNP}/I_{LD}$
ratio observed in figure 10(b). Note also that
$\Lambda _{D}^2/\Gamma$
is the ratio between the normal derivatives of the electric field in the Debye layer and in the surrounding vacuum. Thus, in cone-jets of dilute electrolytes, it is the ratio of the derivatives of the electric field rather than the electric fields what determines whether the LD model is an accurate alternative to the electrokinetic formulation.
Figure 11. Ratio between the surface charge and its equilibrium value for the cone-jets in figure 10 (obtained from the LD solution).
Gañán–Calvo and collaborators have postulated since the early 1990s (Gañán-Calvo et al. Reference Gañán-Calvo, Barrero and Pantano-Rubiño1993; Gañán-Calvo Reference Gañán-Calvo1999, Reference Gañán-Calvo2004) that the surface charge density in cone-jets is always in electrostatic quasiequilibrium. They argue that the electrical relaxation time
$t_e = \varepsilon \varepsilon _o/K$
is much smaller than the residence time of the flow
$t_f$
, and as a result any potential volumetric charge density quasi-instantaneously vanishes, net charge only existing at the surface. Gañán-Calvo et al. (Reference Gañán-Calvo, López-Herrera, Herrada, Ramos and Montanero2018) also conclude that surface charge in quasiequilibrium, or equivalently
$t_e/t_f\ll 1$
, is a sufficient condition for the application of the LD model in cone-jets, due to the explicit assumption of zero volumetric charge density. The cone-jets modelled in the present study are appropriately reproduced by the LD model despite not fulfilling
$t_e/t_f\ll 1$
, suggesting that this sufficient condition may be overly conservative. In the transition region of cone-jets the ratio
$t_e/t_f$
is well approximated by
$\varepsilon /(\pi \Pi _Q^{1/2})$
(Gamero-Castaño & Magnani Reference Gamero-Castaño and Magnani2019b
), resulting on
$t_e/t_f \approx 1.50$
for the cone-jet in figure 5(a) and
$t_e/t_f \approx 2.10$
for the cone-jet in figure 5(b). Although the surface charge is far from electrostatic equilibrium under these conditions,
$t_e/t_f = {O}(1)$
, the LD and PNP solutions in figure 5 coincide. Figure 11 further demonstrates that the surface charge can be far from electrostatic equilibrium in the transition region. This figure plots the local ratio between the surface charge and its equilibrium value,
$\sigma _{eq} = \varepsilon _0E_n^o$
, for the six cone-jets analysed in figure 10. All flow rates are well above the minimum flow rates recorded in experiments (Gamero-Castaño Reference Gamero-Castaño2019). Although the surface charge is in quasiequilibrium sufficiently upstream and downstream from the origin, its value is significantly smaller than
$\sigma _{eq}$
in the central transition region, the divergence increasing with
$\varepsilon /(\pi \Pi _Q^{1/2})$
(Gamero-Castaño & Magnani Reference Gamero-Castaño and Magnani2019b
). Despite the significant departure of the surface charge from equilibrium, figures 5 and 10 show that the LD and PNP solutions for these cone-jets coincide as long as the necessary condition
$\varepsilon \Lambda _{D}^{2}/(\pi \Gamma ) \gtrsim 1$
is fulfilled. Finally, it is important to note that, regardless of whether
$t_e/t_f\ll 1$
is fulfilled, the volumetric charge density in a cone-jet is zero as long as electrical conduction in the bulk is governed by Ohm’s law with a constant value of the electrical conductivity, which is a fundamental assumption in the LD model. This is due to the fact that fluid particle lines in the bulk of the transition region do not traverse upstream regions containing a net volumetric charge density (Melcher & Taylor Reference Melcher and Taylor1969).
4.1.1. Ion diffusivity effects
The previous section considered equal cation and anion diffusivities,
$D^{+}_{r}=0.5$
. It is possible to design dilute electrolytes by varying the concentration and nature of the salt dissolved in a given liquid, so that they have identical physical properties (
$\rho$
,
$\gamma$
,
$K$
,
$\mu$
and
$\varepsilon$
) and different cation and anion diffusivities. Although the traditional scaling laws and the numerical solution of the LD model for these electrolytes are identical, the different diffusivities may influence the cone-jets (López-Herrera et al. Reference López-Herrera, Herrada and Gañán-Calvo2023). Figure 12 plots the total current computed with the PNP model as a function of
$D^{+}_{r}$
, for cone-jets of EG and PC at several dimensionless flow rates, Reynolds numbers and
$({\varepsilon \Lambda _{D}^{2}})/({\pi \Gamma })$
ratios. Different anion/cation diffusivities influence the total current, but this effect is strong only when
$({\varepsilon \Lambda _{D}^{2})}/({\pi \Gamma })\lt 1$
. For example, for EG operating at
${Re}=0.095$
,
$\Pi _{Q}=64$
and
${(\varepsilon \Lambda _{D}^{2})}/{(\pi \Gamma)} = 5.7$
, the dimensionless current slightly increases from
$2.20$
to
$2.39$
as
$D^{+}_{r}$
increases from 0.1 to 0.9; however, for the same values of the Reynolds number and the dimensionless flow rate but
${(\varepsilon \Lambda _{D}^{2})}/{(\pi \Gamma)} = 0.9$
, the current varies from 1.61 to 2.77. Cone-jets of PC exhibit a similar trend.
Figure 12. Total current as a function of the relative cation diffusivity
$D^{+}_{r}$
and for several
${Re}$
,
$\Pi _{Q}$
and
${(\varepsilon \Lambda _{D}^{2})}/{(\pi \Gamma)}$
values, for (a) EG and (b) PC.
The ion density and electrical conductivity profiles in figure 13 for EG,
${Re}=0.95$
and
$\Pi _{Q}=64$
, help explaining the influence of disparate diffusivities on the total current. The profiles correspond to the current crossover point,
$I_{b}=I_{s}$
, and cross-sections upstream and downstream of this point,
$I_b=0.95I$
and
$I_{s}=0.95I$
. The profiles for
${(\varepsilon \Lambda _{D}^{2})}/{(\pi \Gamma)} = 5.7$
exhibit a large central region where the liquid is electrically neutral and the local electrical conductivity is equal to the bulk conductivity. The Debye layers are thin, and the ionic concentrations inside the layers are affected by the relative diffusivity: for
$D^{+}_{r}=0.1$
, that is when the anion is significantly more mobile than the cation, the depletion of the anion in the Debye layer is larger than for
$D^{+}_{r}=0.9$
; meanwhile, the cation concentration in the Debye layer is slightly smaller for
$D^{+}_{r}=0.1$
than for
$D^{+}_{r}=0.9$
. This translates into a smaller electrical conductivity in the Debye layer and a smaller total current for
$D^{+}_{r}=0.1$
than for
$D^{+}_{r}=0.9$
. But the net effect of the varying conductivity on the total current is small, because the Debye layer occupies a small fraction of the cross-section. On the other hand, the profiles for
${(\varepsilon \Lambda _{D}^{2})}/{(\pi \Gamma)} = 0.9$
are very different: the Debye layers are more extensive, occupying the whole cross-section at the current crossover point and farther downstream. The relative diffusivity has the same effect as before, lowering significantly the electrical conductivity when the anion is much more mobile than the cation, and increasing significantly the conductivity when the cation is much more mobile than the anion. The effect of
$D^{+}_{r}$
on the electrical conductivity has now a much larger effect on the total current, see figure 12, because the Debye layer extends throughout the cross-section of the jet.
Figure 13. Concentration profiles and electrical conductivity for EG,
${Re}=0.095$
and
$\Pi _{Q}=64$
. We consider two extreme values of the cation relative diffusivity
$D^{+}$
, and three different axial positions,
$I_{b}=0.95I$
,
$I_{b}=I_{s}$
and
$I_{s}=0.95I$
: (a)
${(\varepsilon \Lambda _{D}^{2})}/{(\pi \Gamma)}=5.7$
; (b)
${(\varepsilon \Lambda _{D}^{2})}/{(\pi \Gamma)}=0.9$
. The electrical conductivity is normalized with its value at the axis,
$K_{b}$
.
Although the possibility for a strong influence of
$D^{+}_{r}$
on the operation of cone-jets exists, it seems that this effect is not important in most cases considered in the current study. The main reason for this is that most dilute electrolytes operate in the region
${(\varepsilon \Lambda _{D}^{2})}/{(\pi \Gamma)} \gtrsim 1$
. In addition, relative diffusivities do not seem to exhibit extreme values: an extreme case such as a HCl-based electrolyte, with the H
$^{+}$
ion having a mobility roughly
$4.65$
times higher than Cl
$^{-}$
, has a relative diffusivity
$D^{+}_{r}\approx 0.8$
; in the case of the ionic liquid EMI-Im, the diffusion coefficient of EMI
$^{+}$
is only
$1.2$
times larger than that of Im
$^{-}$
(Tokuda et al. Reference Tokuda, Tsuzuki, Susan, Hayamizu and Watanabe2006); and the relative diffusivity in a LiCl-water electrolyte is
$D^{+}_{r}\approx 0.33$
. On the other hand, disparate anion/cation diffusivities may have a substantial effect on applications and processes that intrinsically depend on the ion concentrations in the Debye layer, such as electrospray mass spectrometry and ion field emission, regardless of the value of
${(\varepsilon \Lambda _{D}^{2})}/{(\pi \Gamma)}$
. The effect would be more accentuated in water-based solutions near the minimum flow rate (water is a common medium in electrospray mass spectrometry) since, as discussed in the previous section,
${(\varepsilon \Lambda _{D}^{2})}/{(\pi \Gamma)} \lesssim 1$
under such conditions.
4.2. Ionic and highly conducting liquids,
$\delta _{c}\sim {O}(10^{-1})$
,
$\nu _{a} \sim {O}(10^{-1})$
In this section we apply the PNP and MEK models to the ionic liquids EMI-Im and EAN, and compare these solutions with the LD model. Table 5 lists the physical properties and the dimensionless numbers appearing in the governing equations. The radius of the cone-jet
$R(z)$
for a given
$ \{\varepsilon , \Pi _Q, {Re} \}$
state is obtained from the LD solution (Magnani, Caballero-Pérez & Gamero-Castaño Reference Magnani, Caballero-Pérez and Gamero-Castaño2025). Given that the diffusivities of the cation and the anion are similar for each liquid (Noda et al. Reference Noda, Hayamizu and Watanabe2001; Tokuda et al. Reference Tokuda, Tsuzuki, Susan, Hayamizu and Watanabe2006; Filippov et al. Reference Filippov, Alexandrov, Gimatdinov and Shah2021), we simply use
$D^{+}_{r}=0.5$
.
Table 5. Physical properties and dimensionless numbers for ionic liquid EMI-Im and EAN (
$\tilde {T}_{o}=294$
K).
Figure 14 compares the solutions of the three models for EMI-Im at a flow rate
$\Pi _{Q}=200$
. Figure 14(a) shows the radius of the cone-jet and the variation of the different contributions to the total current. The MEK model includes the additional ion correlation current (Kilic et al. Reference Kilic, Bazant and Ajdari2007b
),
\begin{equation} \tilde {I}_{c}=2\pi \int _{0}^{\tilde {R}(z)}\frac {ea^{3}(D^{-}\tilde {n}^{-}-D^{+}\tilde {n}^{+})}{1-a^{3}(\tilde {n}^{+}+\tilde {n}^{-})}\frac {\partial }{\partial \tilde {z}}(\tilde {n}^{+}+\tilde {n}^{-})\tilde {r}\,{\rm d}\tilde {r}, \end{equation}
which is negligible in all cases considered in this study. The diffusion current is substantial only in the MEK solution, negative and approximately
$4.5\,\%$
of the total current far upstream in the cone, and negligible along the jet. Despite this, the MEK solution yields a slightly larger total current (
$167$
nA) than the PNP solution (
$163$
nA), and the latter is larger than the current of the LD solution (
$156.5$
nA). These values are in good agreement with the 177 nA recorded in experiments (Caballero-Pérez & Gamero-Castaño Reference Caballero-Pérez and Gamero-Castaño2025). The higher electric current yielded by the electrokinetic models is due to more intense ohmic dissipation, see figure 14(b), which increases self-heating, the temperature and the electrical conductivity, figure 14(c). Despite differences in the total current, all three solutions yield very similar electric fields on the surface, see figure 14(d). While the differences in the currents, temperature and ohmic dissipation profiles between the MEK solution and the other two solutions could have been expected, it is worth noting that the PNP and LD solutions do not coincide despite the fulfilment of
$\varepsilon \Lambda _{D}^{2}/(\pi \Gamma ) \gg 1$
(the value of this ratio is 250). This is due to the substantial self-heating taking place in the cone-jets, and the strong dependence of self-heating on the electrical conductivity, which includes positive feedback (Magnani & Gamero-Castaño Reference Magnani and Gamero-Castaño2024). The more accurate description and the higher values of the electrical conductivity associated with the electrokinetic formulations lead to higher self-heating and total currents.
Figure 14. Comparison between the MEK, PNP and LD solutions for ionic liquid EMI-Im operating at
$\Pi _{Q}=200$
: (a) cone-jet radius and conduction, convection, diffusion and ion correlation currents; (b) ohmic and viscous linear power densities; (c) electrical conductivity and change in temperature; (d) electric field on the surface.
Figure 15. Ion concentration and electrical conductivity profiles obtained with the MEK and PNP models for EMI-Im and
$\Pi _{Q}=200$
, at three different locations:
$I_{b}=0.95I$
,
$I_{b}=I_{s}$
and
$I_{s}=0.95I$
. The electrical conductivity is normalized with its value at the axis,
$K_{b}$
.
Figure 15 shows the ion concentration profiles at the cross-sections where
$I_b=0.95I$
,
$I_{b}=I_{s}$
and
$I_{s}=0.95I$
. The Debye layer resolved by the PNP model has a thickness of approximately
$0.2{-}0.4$
nm, clearly unphysical in view of the diameter of the EMI-Im molecule, approximately 0.75 nm. The thickness of the layer where the variations of the ion concentrations are significant is much larger in the MEK solution, approximately
$3.5{-}5$
nm. Moreover, the variation of the volumetric charge density,
$e(n^+ - n^-)$
, is non-monotonic, with a region of negative charge density closer to the axis followed by a much more intense positive charge density near the surface. Furthermore, near the surface, the concentration of the cation is much larger in the MEK solution than in the PNP solution. Such charge density oscillations are observed in the surface force measurements of ionic liquids confined between charged surfaces (Hayes et al. Reference Hayes, Warr and Atkin2010; Perkin Reference Perkin2012; Li et al. Reference Li, Endres and Atkin2013; Smith et al. Reference Smith, Lee and Perkin2016). Despite the drastically different concentration profiles, the total current yielded by the PNP and MEK models nearly coincide. Note also that the electrical conductivity is nearly constant along the cross-sections in both electrokinetic solutions, as a result of the constancy of the total ion concentration,
$n^+ + n^-$
, and the equal ion diffusivities,
$D^{+}_{r}=0.5$
. The near constancy of
$n^+ + n^-$
along the radial coordinate is due to the very small radial components of the electric field and the velocity inside the slender jet. Adding (3.6) for the cation and the anion and integrating once yields
\begin{align} &\frac {D^+}{Pe_{D}^{+} D^{+}_o}\left (\boldsymbol{\nabla }({n^+}+{n^-})+\frac {2\Gamma }{\pi \sqrt {\Pi _{Q}}} ({n^+}-{n^-})\boldsymbol{\nabla }\phi ^{i}+\frac {\nu _{a}({n^+}+{n^-})\boldsymbol{\nabla }(n^{+}+n^{-})}{2-\nu _{a}(n^{+}+n^{-})}\right )\notag\\&\quad -\boldsymbol{u} ({n^+}+{n^-}) = 0, \end{align}
and, since
$u_r \cong 0$
and
$E^i_r \cong 0$
(note that
$\Gamma$
is almost one order of magnitude smaller for EMI-Im than for the EG and PC solutions analysed in figure 6), the radial component of this equation simplifies to
\begin{equation} \left [1+\frac {\nu _{a}({n^+}+{n^-})}{2-\nu _{a}(n^{+}+n^{-})} \right ]\frac {\partial }{\partial r} ({n^+}+{n^-})= 0\, \Longrightarrow \, {n^+}+{n^-} = f(z). \end{equation}
Figure 16. Influence of the ratio between the electrostatic correlation and Bjerrum lengths,
$l_{c}/l_{B}$
, on the MEK solution for EMI-Im at
$\Pi _{Q}=200$
: (a) the total emitted current and temperature increase (evaluated at
$I_s=0.95I$
); (b) concentration profiles at three different axial positions,
$I_{b}=0.95I$
,
$I_{b}=I_{s}$
and
$I_{s}=0.95I$
.
As discussed in § 2, the analysis so far has assumed that the electrostatic correlation length and the Bjerrum length are equal,
$l_c = l_B$
. Figure 16 explores the impact of varying this ratio using EMI-Im at a flow rate of
$\Pi _Q=200$
as case study. Figure 16(a) shows the total emitted current and the temperature increase, the latter evaluated in the downstream region where
$I_{s}=0.95I$
. When
$l_{c}/l_{B}=0$
, the MEK model recovers the classical PNP behaviour with no electrostatic correlations, whereas
$l_{c}/l_{B}=1$
represents an approximated upper bound for the electrostatic correlation length, and is the case considered in figures 14 and 15. As
$l_{c}/l_{B}$
increases, the temperature intensifies, enhancing ion mobility and resulting in a higher total current. This leads to a maximum increase in
$\Delta \tilde {T}$
and
$\tilde {I}$
of
$13.8\,\%$
and
$2.9\,\%$
, respectively. Figure 16(b) compares the concentration profiles of cations and anions at three axial positions:
$I_b=0.95I$
,
$I_{b}=I_{s}$
and
$I_{s}=0.95I$
. As
$l_{c}/l_{B}$
increases, the effective Debye length expands, resulting in more pronounced volume charge density oscillations near the interface. The effective thickness of the charge oscillation layer ranges from approximately 1.2–1.75 nm for
$l_{c}/l_{B}=0.25$
, to 3.5–5 nm for
$l_{c}/l_{B}=1$
. Furthermore, higher values of
$l_{c}/l_{B}$
amplify the intensity of both the negative and positive charge densities, indicating stronger electrostatic correlation.
Figure 17 shows the total change in temperature along the cone-jet, evaluated at the current crossover and at the point where
$I_{s}=0.95I$
, as a function of the dimensionless flow rate. The upstream temperature is 294 K in all cases. The position of the current crossover is important because it is representative of the centre of the transition region, where variations in physical properties have the strongest effects. The total temperature increases in typical cone-jets for ionic liquids exceeds several hundred degrees. The electrokinetic models produce temperature differences significantly larger than the LD model, and the temperatures in the MEK solution are always higher than in the PNP solution. The latter is due to the non-monotonic and wider Debye layer of the MEK model, its effect on the radial component of the electric field and the associated increase in the ohmic dissipation; different cation/anion diffusivities,
$D^{+}_{r} \neq 0.5$
, would have further increased the difference between the MEK and PNP solutions.
Figure 17. Total change in temperature at the current crossover and at
$I_{s}=0.95I_{T}$
for ionic liquids EMI-Im and EAN. Comparison between the MEK, PNP and LD solutions.
It is worth noting that the temperature increases for these ionic liquids are much more substantial than in the cone-jets of EG and PC shown in figure 5. This is primarily due to the much larger electrical conductivity of ionic liquids, and its effect on the radius of the jet, the axial length of the transition region and directly on ohmic dissipation, which is proportional to
$K$
. The following expression is an approximated scaling law for the temperature increase (Magnani & Gamero-Castaño Reference Magnani and Gamero-Castaño2024):
\begin{equation} {\Delta \tilde {T}}_{est} \approx \frac {7.74}{{Re}^{0.11}{\Pi _Q}^{0.11}} \frac {1}{\pi c} \left (\frac {\gamma K}{\varepsilon _0\rho }\right )^{2/3} \approx \frac {7.74}{\pi c} \left (\frac {\gamma K}{\varepsilon _0\rho }\right )^{2/3}, \end{equation}
which clearly indicates the importance of the electrical conductivity on
$\Delta \tilde {T}$
(the electrical conductivity is the only physical property in (4.10) that can vary by many orders of magnitude in electrosprayed solutions). The main assumptions in this law are negligible evaporation and heat flux on the surface of the cone-jet, conditions which are fulfilled in cone-jets of ionic liquids electrosprayed in vacuum. For example, for a temperature increase of 10 °C in a water cone-jet, the electrical conductivity needs to be approximately 0.27 S m−1, which is a substantial value. For this conductivity, the characteristic radius of the jet would be
$r_c = 4.2$
nm for a dimensionless flow rate of unity, and 42 nm for
$\Pi _Q = 100$
. Strong thermal effects in electrosprays of water under atmospheric conditions may not be easily observable not only because of the high electrical conductivity required, the smallness of the jets and droplets and the speed with which the droplets would evaporate, but also because of heat transfer with the surrounding gas and the high latent heat of vaporization and vapour pressure of water, which enhance evaporative cooling.
Finally, figure 18 shows the total current emitted by cone-jets of EMI-Im and EAN as a function of the dimensionless flow rate. We compare the numerical solutions of the three models with experimental data. It is worth noting that the current yielded by the MEK model is always slightly larger than that of the PNP model, and both are larger than the LD solution. The electrokinetic solution for EMI-Im is in good agreement with the experimental results, slightly underpredicting the total current. The electrokinetic solution for EAN is lower than the experimental current. We suspect that this is due to inaccuracies in the
$D^{+}(\tilde {T})$
function employed in the calculations, which is derived from measurements of the electrical conductivity in the temperature range between
$256$
K and
$391$
K (Bouzón-Capelo et al. Reference Bouzón-Capelo, Méndez-Morales, Carrete, López Lago, Vila, Cabeza, Rodríguez, Turmine and Varela2012). The electrokinetic model predicts temperatures in excess of 400 K for all flow rates simulated, and as high as 827 K at the lowest flow rates, i.e. well above the experimental temperature range for which we have conductivity values. This large extrapolation combined with the exponential
$K(\tilde {T})$
function, the positive feedback between self-heating and the electrical conductivity, and the strong dependence of the solution on the electrical conductivity, are the likely cause for the disagreement. Additional factors that may contribute to the discrepancy include the use of a fixed interface provided by the LD solution, rather than an surface self-consistently obtained with the electrokinetic models; the assumption of a fully dissociated ionic liquid, whereas in reality partial dissociation or ion pairing may occur; and evaporation of the ionic liquid, which is not included in the model and may be significant for EAN.
Figure 18. Total current versus dimensionless flow rate for ionic liquids EMI-Im and EAN. Comparison between experimental data and the MEK, PNP and LD solutions. Here MCP and MGC refer to experimental results by Caballero-Pérez & Gamero-Castaño (Reference Caballero-Pérez and Gamero-Castaño2025) and Gamero-Castaño & Cisquella-Serra (Reference Gamero-Castaño and Cisquella-Serra2020), respectively.
In summary, because of the strong temperature variations and the positive feedback between self-heating and the electrical conductivity, the improved description of both the electrical conductivity and the ohmic dissipation made possible by an electrokinetic model increases the accuracy of the solution compared with the LD formulation. Furthermore, the MEK model is superior to the PNP formulation because, as figure 15 illustrates, the Debye layer resolved by the latter is unphysical, and therefore it does not improve the fidelity of the physics compared with the LD formulation. The accurate description of the ion concentration fields by the MEK model is also important to investigate associated phenomena. For example, ion-field emission (Iribarne & Thomson Reference Iribarne and Thomson1976), which we do not consider presently but can become important in charged nanodroplets (Loscertales & Fernández De La Mora Reference Loscertales and Fernández De La Mora1995; Misra & Gamero-Castaño Reference Misra and Gamero-Castaño2023), cone-jets (Gamero-Castaño Reference Gamero-Castaño2002) and Taylor cones (Gallud & Lozano Reference Gallud and Lozano2022; Magnani & Gamero-Castaño Reference Magnani and Gamero-Castaño2023), is a kinetic phenomenon in which both the ion concentration and the temperature near the surface play crucial roles.
5. Conclusions
We have modelled electrosprays operating in the cone-jet mode using both the PNP equations and a MEK formulation that accounts for overscreening and overcrowding of ions (Kilic et al. Reference Kilic, Bazant and Ajdari2007b ; Bazant et al. Reference Bazant, Storey and Kornyshev2011). While the former is strictly applicable to dilute electrolytes (liquids with low and medium electrical conductivities), the latter is designed for concentrated electrolytes such as ionic liquids and other liquids with high electrical conductivities. The results from the PNP and MEK models are compared with the traditional LD model and experimental data.
The LD and PNP solutions coincide for isothermal cone-jets as long as the required condition
$\varepsilon \Lambda _{D}^{2}/(\pi \Gamma )\gtrsim 1$
is fulfilled. This is less restrictive than the Baygent–Saville condition
$ \Lambda _{D}/\Gamma \gg 1$
(Baygents & Saville Reference Baygents and Saville1990; Saville Reference Saville1997; Schnitzer & Yariv Reference Schnitzer and Yariv2015). In practice, most cone-jets of liquids with low and moderate electrical conductivities fulfil
$\varepsilon \Lambda _{D}^{2}/(\pi \Gamma )\gtrsim 1$
. In agreement with Mori & Young (Reference Mori and Young2018), the electrokinetic formulation does not require the inclusion of the surface charge postulated by Saville (Reference Saville1997) and Schnitzer & Yariv (Reference Schnitzer and Yariv2015), to reproduce either experimental cone-jets or the LD solution. Disparate ion diffusivities can significantly affect the current of the cone-jet, but only when
$\varepsilon \Lambda _{D}^{2}/(\pi \Gamma ) \lt 1$
.
The MEK model provides the best results for cone-jets of ionic liquids, predicting higher total currents and temperature increases than the PNP and LD models. This is due to a better description of the Debye layer, which results in larger ohmic dissipation and self-heating. The MEK model resolves a non-monotonic charge density within the Debye layer, similar to experimental observations in ionic liquids confined between charged surfaces (Smith et al. Reference Smith, Lee and Perkin2016). In contrast, the PNP model predicts an unphysically thin Debye layer. Self-heating in cone-jets of ionic liquids is significant, increasing at decreasing flow rate. Temperature variations of several hundred degrees are common in EMI-Im and EAN cone-jets. This causes large variations in physical properties such as viscosity, electrical conductivity and ion diffusivity, which need to be accounted for when modelling these cone-jets, or analysing experimental data.
As a continuation to this work, we plan to investigate the applicability of the
$\varepsilon \Lambda _{D}^{2}/(\pi \Gamma )\gtrsim 1$
condition to the modelling of time-dependent electrohydrodynamic processes like tip streaming. A second interesting problem is the modelling of phenomena such as ion-field emission that require the accurate description of the Debye layer made possible by the MEK formulation.
Acknowledgments
The authors would like thank Professor Gañán-Calvo and Professor López-Herrera for multiple discussions on the application of electrokinetics to cone-jets. The authors would also like to thank M. Caballero-Pérez for providing the experimental measurements of ionic liquids used in this study.
Funding
This work was funded by the Air Force Office of Scientific Research, award no. FA9550-21-1-0200.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Regression laws for the electrical conductivity, viscosity and ion diffusivity
The coefficients
$Y_D$
,
$B_D$
,
$T_D$
,
$Y_\mu$
,
$B_\mu$
and
$T_\mu$
listed in table 3 for the regression laws of the ion diffusivity and viscosity are obtained from experimental measurements of the electrical conductivity and viscosity. For the two ionic liquids we have collected experimental data from the ionic liquid database managed by the National Institute of Standards and Technology (Kazakov et al. Reference Kazakov, Magee, Chirico, Paulechka, Diky, Muzny, Kroenlein and Frenkel2022). This database includes data from different publications. At the moment of writing this article, the database was last updated on 4 June 2024. For EMI-Im, the database contains 176 measurements for the electrical conductivity and 687 for the viscosity. For EAN, the database includes 110 measurements for the electrical conductivity and 95 for the viscosity. The data were harvested from the database using the Python library pyILT2.
Figure 19. Experimental data and regression laws for the electrical conductivity and viscosity of EMI-Im.
Figure 20. Experimental data and regression laws for the electrical conductivity and viscosity of EAN.
For each property, the experimental data includes a measured value
$y_i$
, temperature
$T_i$
and an uncertainty
$\pm \delta _i$
. The fittings are obtained by minimizing
\begin{equation} \epsilon =\sum _{i=1}^N\left [\frac {1}{1+\delta _i/y_i}\left (\frac {f(\boldsymbol{p},T_i)}{y_i}-1\right )\right ]^2, \end{equation}
where
$N$
is the total number of points and
$f(\boldsymbol{p},T_i)$
is the fitting function selected. We use
\begin{equation} f(\boldsymbol{p},T)=p_a\exp \left (\frac {p_b}{T-p_c}\right ), \end{equation}
where
$p_a$
,
$p_b$
and
$p_c$
correspond to
$Y_K$
,
$B_K$
,
$T_K$
for the electrical conductivity and
$Y_\mu$
,
$B_\mu$
,
$T_\mu$
for the viscosity. Figures 19 and 20 show the electrical conductivity and viscosity fittings for EMI-Im and EAN. The fittings for the electrical conductivity are used directly in the LD model (Magnani & Gamero-Castaño Reference Magnani and Gamero-Castaño2024). We also use the fitting for the electrical conductivity to evaluate the coefficients inside the exponential for the average ion diffusivity,
$B_D=B_K$
and
$T_D=T_K$
, while the
$Y_D$
coefficient is obtained from the experimental values shown in table 1. This process ensures that the electrical conductivity employed in the LD model can be made to coincide with the bulk conductivity in the electrokinetic formulations, enabling consistent comparison of the solutions. The coefficients for the PC and EG mixtures are taken from Magnani & Gamero-Castaño (Reference Magnani and Gamero-Castaño2024), which employ a similar procedure to compute
$Y_K$
,
$B_K$
,
$T_K$
,
$Y_\mu$
,
$B_\mu$
and
$T_\mu$
.
Open access
