1. Introduction
The electrospray technique is based on the generation of charged droplets and ions from a liquid under the influence of an electric field. These particles are important for several applications, ranging from mass spectrometry to materials science and medicine. In mass spectrometry, electrospray ionisation has become a critical tool for the analysis of biomolecules, enabling the study of complex proteins and nucleic acids (Fenn et al. Reference Fenn, Mann, Meng, Wong and Whitehouse1989; De Vijlder et al. Reference De Vijlder, Valkenborg, Lemière, Romijn, Laukens and Cuyckens2018; Prabhu et al. Reference Prabhu, Williams, Wilm and Urban2023). Electrospray is also useful in the fabrication of nanomaterials (Loscertales et al. Reference Loscertales, Barrero, Guerrero, Cortijo, Marquez and Gañán-Calvo2002; Wang et al. Reference Wang, Zheng, He, Wei, Liu, Lin, Zheng and Sun2013; Xu, Tao & Lozano Reference Xu, Tao and Lozano2018; Li et al. Reference Li, Lin, Han, Mu and Yu2021), drug delivery systems (Sill & Von Recum Reference Sill and von Recum2008; Steipel et al. Reference Steipel, Gallovic, Batty, Bachelder and Ainslie2019; Alfatama, Shahzad & Choukaife Reference Alfatama, Shahzad and Choukaife2024) and electric propulsion for space applications (Gamero-Castaño & Hruby Reference Gamero-Castaño and Hruby2001; Krejci et al. Reference Krejci, Mier-Hicks, Thomas, Haag and Lozano2017; Jenkins, Krejci & Lozano Reference Jenkins, Krejci and Lozano2018).
A crucial aspect of electrospray technology is the role of the liquid’s electrical conductivity
$K$
. Unlike other relevant physical properties like the surface tension
$\gamma$
, the dielectric constant
$\varepsilon$
, the density
$\rho$
and the viscosity
$\mu$
, the electrical conductivity can be varied by many orders of magnitude, and its value affects all features of an electrospray such as the electric current and the size of the droplets. Liquids with high conductivities (Gamero-Castaño & Fernandez de la Mora Reference Gamero-Castaño and Fernandez de la Mora2000; Larriba et al. Reference Larriba2007) produce nanometric-size droplets with high charge-to-mass ratios. At the same time, they are also prone to significant temperature increases due to energy dissipation and self-heating (Gamero-Castaño Reference Gamero-Castaño2010, Reference Gamero-Castaño2019; Magnani & Gamero-Castaño Reference Magnani and Gamero-Castaño2024). This alters the liquid properties, thereby affecting the electrospray process. Understanding the physics underlying electrospraying is essential for optimising its applications. Researchers have dedicated significant effort to both experimental and theoretical investigations in order to understand the mechanisms of electrospray atomisation. Experimentally, advances in high-speed imaging (Parvin et al. Reference Parvin, Galicia, Gauntt, Carney, Nguyen, Park, Heffernan and Vertes2005; Alexander, Paine & Stark Reference Alexander, Paine and Stark2006; Hsu et al. Reference Hsu, Prabhu, Chang, Hsu, Buchowiecki and Urban2023) and measuring techniques (Gamero-Castaño & Hruby Reference Gamero-Castaño and Hruby2002; Gamero-Castaño & Cisquella-Serra Reference Gamero-Castaño and Cisquella-Serra2021; Thuppul et al. Reference Thuppul, Collins, Wright, Uchizono and Wirz2021) have provided insights into the dynamics of droplet formation and charge distribution, while theoretical models, including computational fluid dynamics, have helped predict electrospray behaviour under various conditions. Traditionally, these models have been based on the leaky-dielectric assumption under isothermal conditions (Melcher & Taylor Reference Melcher and Taylor1969; Basaran et al. Reference Basaran, Patzek, Benner and Scriven1995; Gañán-Calvo Reference Gañán-Calvo1997, Reference Gañán-Calvo1999; Saville Reference Saville1997; Higuera Reference Higuera2003; Fernandez de la Mora Reference Fernández de la Mora2007; Herrada et al. Reference Herrada, López-Herrera, Gañán-Calvo, Vega, Montanero and Popinet2012; Gamero-Castaño & Magnani Reference Gamero-Castaño and Magnani2019; López-Herrera et al. Reference López-Herrera, Herrada, Gamero-Castaño and Gañán-Calvo2023). Magnani & Gamero-Castaño (Reference Magnani and Gamero-Castaño2024) extended this approach to non-isothermal cone jets by including dissipation, self-heating and temperature-dependent physical properties in the model. Thermal effects are magnified by high conductivities and low flow rates, which are required for electrospray applications such as electrospray propulsion and energetic nanodroplet impact.
Ionic liquids are often employed in electrospray research due to their distinctive properties (Prado & Weber Reference Prado and Weber2016; Kaur, Kumar & Singla Reference Kaur, Kumar and Singla2022). In particular, their high electrical conductivities and extremely low vapour pressures make them ideal for vacuum applications requiring the generation of charged nanodroplets. For example, ionic liquids are the propellants of choice in electrospray propulsion (Romero-Sanz, De Carcer & Fernandez de la Mora Reference Romero-Sanz, De Carcer and Fernandez de la Mora2005; Legge & Lozano Reference Legge and Lozano2011; Prince et al. Reference Prince, Fritz and Chiu2012; Grustan-Gutiérrez & Gamero-Castaño Reference Grustan-Gutiérrez and Gamero-Castaño2017; Krejci et al. Reference Krejci, Mier-Hicks, Thomas, Haag and Lozano2017; Cisquella-Serra et al. Reference Cisquella-Serra, Galobardes-Esteban and Gamero-Castaño2022). In this article we improve an existing cone-jet model (Gamero-Castaño & Magnani Reference Gamero-Castaño and Magnani2019; Magnani & Gamero-Castaño Reference Magnani and Gamero-Castaño2024) to study four ionic liquids: 1-ethyl-3-methylimidazolium bis((trifluoromethyl)sulfonyl)imide (EMI-Im, CAS 174899-82-2), 1-butyl-3-methylimidazolium tricyanomethane (BMI-TCM, CAS 878027-73-7), 1-ethyl-3-methylimidazolium trifluoroacetate (EMI-TFA, CAS 174899-65-1) and ethylammonium nitrate (EAN, CAS 22113-86-6). The model employs the leaky-dielectric formulation with the addition of dissipation, self-heating and temperature-dependent physical properties. Magnani & Gamero-Castaño (Reference Magnani and Gamero-Castaño2024) first included thermal effects to study solutions of propylene carbonate and ethylene glycol with low and moderate electrical conductivities, for which self-heating ranges from negligible to noticeable (the maximum temperature increases computed were 15.7 °C for propylene carbonate and 1.80 °C for ethylene glycol). This choice of liquids made it possible to vary gradually the electrical conductivity while keeping constant all other physical properties, so that the transition from negligible to noticeable self-heating could be studied both numerically and experimentally. Magnani & Gamero-Castaño (Reference Magnani and Gamero-Castaño2024) also used the model to study one flow rate of EMI-Im, to verify the much higher effects of self-heating in liquids with sufficiently high electrical conductivity,
$K\gtrsim 0.5$
S m−1. The present study extends Magnani & Gamero-Castaño (Reference Magnani and Gamero-Castaño2024) by:
-
(i) extensively analysing liquids in which self-heating largely drives the experimental behaviour and the numerical solution;
-
(ii) improving the numerical model by implementing a new optimisation method and two strategies to eliminate unphysical oscillations in the position of the free surface, which speed up the convergence and reduce excursions of the iterative algorithm used to solve the model equations;
-
(iii) evaluating in Appendix B the importance of including other temperature-dependent physical properties in addition to the electrical conductivity and viscosity (Magnani & Gamero-Castaño Reference Magnani and Gamero-Castaño2024), which due to their exponential dependence on temperature are the obvious candidates to couple self-heating and the electrohydrodynamics of the flow;
-
(iv) extending the model to include ion field emission from the surface of the steady cone jet, see Appendix C.
An important goal of the present article is to highlight the elevated temperature increases taking place in cone jets of ionic liquids, and the strong departure from the commonly assumed isothermal behaviour. Section 2 describes in detail the model, the discretisation of the equations and the solution scheme. Section 3 describes the numerical solutions, including the comparison between model results, experimental data and the isothermal solution; and the analysis of key features such as temperature increases and dissipation densities that are only accessible through numerical solutions. Finally, § 4 summarises the findings.
2. Cone-jet model
We use the traditional leaky-dielectric model (Melcher & Taylor Reference Melcher and Taylor1969; Saville Reference Saville1997), with the inclusion of self-heating resulting from ohmic and viscous dissipation. Self-heating is coupled to the electrohydrodynamics through the viscosity and the electrical conductivity, which are strong functions of temperature (Stokes & Mills Reference Stokes and Mills1966; FogelSon & Likhachev Reference FogelSon and Likhachev2001); all other physical properties are kept constant. The model is steady state and axisymmetric, and includes a liquid phase (the cone jet) separated from a surrounding vacuum by a free surface
$R(x)$
. Inside the liquid the velocity
$\boldsymbol{v}$
, pressure
$p$
, temperature
$T$
and electric potential
$\varPhi _i$
fields fulfil the equations of conservation of mass, momentum and energy, and a simplified form of charge conservation based on Ohm’s law, while the electric potential
$\varPhi _o$
in the vacuum fulfils the Laplace equation
\begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{v}&=0, \end{align}
\begin{align} \sqrt {\varPi _Q}\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{v}+\pi ^2\boldsymbol{\nabla }p&=\frac {\pi \left [\mu {\nabla} ^2\boldsymbol{v}+\left (\boldsymbol{\nabla }\boldsymbol{v}+{\boldsymbol{\nabla }\boldsymbol{v}}^T\right )\boldsymbol{\cdot }\boldsymbol{\nabla }\mu \right ]}{\mu _0 Re}, \end{align}
\begin{align} \sqrt {\varPi _Q}\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{\nabla }T-\frac {\sqrt {\varPi _Q}}{\textit{Pe}}{\nabla} ^2T&=\frac {\mu \left (\boldsymbol{\nabla }\boldsymbol{v}+{\boldsymbol{\nabla }\boldsymbol{v}}^T\right )\boldsymbol{:}\boldsymbol{\nabla }\boldsymbol{v}}{\mu _0 Re}+\frac {K{\boldsymbol{\nabla }\varPhi _i}\boldsymbol{\cdot }{\boldsymbol{\nabla }\varPhi _i}}{K_0}, \end{align}
\begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }\left (\frac {K \boldsymbol{E}_i}{K_0}\right )&=0 \rightarrow K{\nabla} ^2\varPhi _i+\boldsymbol{\nabla }\varPhi _i\boldsymbol{\cdot }\boldsymbol{\nabla }K=0, \end{align}
\begin{align} {\nabla} ^2\varPhi _o&=0, \end{align}
where
$\mu _0$
and
$K_0$
are the viscosity and the electrical conductivity at 294.15 K. E stands for the electric field. The two terms in the right-hand side of the temperature equation (2.1c
) are the viscous dissipation density
$P^{\prime\prime\prime}_\mu$
and the ohmic dissipation density
$P^{\prime\prime\prime}_\varOmega$
. Equations (2.1a
)–(2.3f
) are in dimensionless form, hereafter all dimensional variables are capped with a tilde. Tables 1 and 2 show the characteristic scales and dimensionless groups
$\varPi_Q$
,
$\textit{Pe}$
and
$\textit{Re}$
. The temperature-dependent electrical conductivity and viscosity are approximated by
Table 1. Characteristic scales used in the non-dimensionalisation of the equations.
$\rho$
,
$\gamma$
and
$c$
stand for the density, surface tension and heat capacity of the liquid;
$\mu _0$
and
$K_0$
are the viscosity and electrical conductivity at 294.15 K.
Table 2. Dimensionless numbers used in the model.
$k$
stands for the thermal conductivity.
\begin{align} K(\tilde {T})&=Y_K \exp \left (\frac {B_K}{\tilde {T}-T_K}\right)\!, \end{align}
\begin{align} \mu (\tilde {T})&=Y_\mu \exp \left (\frac {B_\mu }{\tilde {T}-T_\mu }\right)\!, \end{align}
where the constants
$Y_K$
,
$B_K$
,
$T_K$
,
$Y_\mu$
,
$B_\mu$
and
$T_\mu$
are specific to each liquid. The values of these constants are given in table 3, while Appendix A explains how we obtain them from measurements of the electrical conductivity and the viscosity. These exponential laws are accurate for most liquids (Stokes & Mills Reference Stokes and Mills1966; FogelSon & Likhachev Reference FogelSon and Likhachev2001), including ionic liquids (Okoturo & VanderNoot Reference Okoturo and VanderNoot2004; Leys et al. Reference Leys, Wübbenhorst, Preethy Menon, Rajesh, Thoen, Glorieux, Nockemann, Thijs, Binnemans and Longuemart2008).
Table 3. Physical properties, characteristic scales and dimensionless numbers of the ionic liquids studied at 294.15 K. 1-ethyl-3-methylimidazolium bis((trifluoromethyl)sulfonyl)imide (EMI-Im, CAS 174899-82-2), 1-butyl-3-methylimidazolium tricyanomethane (BMI-TCM, CAS 878027-73-7), 1-ethyl-3-methylimidazolium trifluoroacetate (EMI-TFA, CAS 174899-65-1), ethylammonium nitrate (EAN, CAS 22113-86-6). Viscosity and conductivity are reported as functions of temperature in Kelvin.
On the free surface the model includes conservation of surface charge
$\sigma$
, Gauss’ law, the surface kinematic condition, zero heat flux and the balance of stresses
\begin{align} \frac {\rm d}{{\rm d}x}\left (R\sigma v_s\right )&=\frac {2K}{K_0}RE_n^i\sqrt {1+{R'}^2}, \end{align}
\begin{align} E_n^o-\varepsilon E_n^i&=\frac {\pi \sqrt {\varPi _Q}}{2}\sigma , \end{align}
\begin{align} \boldsymbol{v} \boldsymbol{\cdot }\boldsymbol{n}&=0, \end{align}
\begin{align} \boldsymbol{\nabla }T \boldsymbol{\cdot }\boldsymbol{n}&=0, \end{align}
\begin{align} \frac {\mu }{\mu _0}\boldsymbol{t} \boldsymbol{\cdot }\big ( \boldsymbol{\nabla }\boldsymbol{v}+{\boldsymbol{\nabla }\boldsymbol{v}}^T\big )\boldsymbol{n}&=\frac { Re}{2}E_t\sigma , \end{align}
\begin{align} \boldsymbol{\nabla } \boldsymbol{\cdot }\boldsymbol{n}-p+\frac {2\mu }{\pi \mu _0 Re}\boldsymbol{n} \boldsymbol{\cdot } \boldsymbol{\nabla }\boldsymbol{v}\boldsymbol{n}&=\frac {{E_n^o}^2-\varepsilon {E_n^i}^2-{E_t}^2\left (\varepsilon -1\right )}{2\pi ^2\sqrt {\varPi _Q}}, \end{align}
where
$E_n^i$
and
$E_n^o$
are the normal components of the electric field inside and outside the surface,
$E_t$
is the tangential component, and
$\boldsymbol{n}$
and
$\boldsymbol{t}$
are the normal and tangential vectors. In summary, the model solves for
$R$
,
$\boldsymbol{v}$
,
$p$
,
$T$
,
$\varPhi _i$
,
$\varPhi _o$
and
$\sigma$
, and the solutions are parametric functions of the dimensionless groups
$\varPi _Q$
,
${\textit{Re}}$
,
$\textit{Pe}$
and
$\varepsilon$
, and of the constants parametrising the
$K(\tilde {T})$
and
$\mu (\tilde {T})$
functions.
Section 3 compares the numerical solution with the so-called ‘isothermal solution’. The viscosity and electrical conductivity are kept constant for the latter (we use their values at the upstream temperature 294.15 K). The isothermal model solves the same equations except for conservation of energy (2.1c ) and associated boundary conditions, which are not needed since the temperature is fixed. The solution including self-heating will be referred to as the ‘dissipative solution’, or simply as the ‘solution’, depending on whether we need to differentiate with the isothermal solution.
Although significant amounts of molecular ions are observed in the beams of the simulated cone jets (Gamero-Castaño & Cisquella-Serra Reference Gamero-Castaño and Cisquella-Serra2021; Caballero-Pérez, Galobardes-Esteban & Gamero-Castaño Reference Caballero-Pérez, Galobardes-Esteban and Gamero-Castaño2025), we do not include ion field emission (Magnani & Gamero-Castaño Reference Magnani and Gamero-Castaño2023) in the model, specifically in (2.3a
). The available experimental evidence indicates that, in all simulated flow rates of EMI-Im, ion emission does not take place from the surface of the jet, but from evolving droplets in the jet’s breakup region. This conclusion is based on both published results (Gamero-Castaño & Cisquella-Serra Reference Gamero-Castaño and Cisquella-Serra2021) and additional unpublished measurements of the retarding potential of ions. Perez-Lorenzo & Fernandez de la Mora (Reference Perez-Lorenzo and Fernandez de la Mora2022) report the same behaviour for cone jets of the ionic liquid 1-Ethyl-3-methylimidazolium tris(pentafluoroethyl)trifluorophosphate. In addition, the experimental characterisation by Caballero-Pérez et al. (Reference Caballero-Pérez, Galobardes-Esteban and Gamero-Castaño2025) suggests that ion emission only takes place from the breakup region in all simulated flow rates of BMI-TCM, EMI-TFA and EAN, and that ions start being emitted from the surface of the cone and the steady jet in the case of BMI-TCM at flow rates lower than presently simulated. This behaviour can be explained by the variation of
$E_n^o$
along the cone jet, and the exponential dependence of the ion emission intensity on the electric field (Iribarne & Thomson Reference Iribarne and Thomson1976). The cone jet has two regions where
$E_n^o$
exhibits local maxima: a zone near the current cross-over at the base of the jet (Gamero-Castaño & Fernandez de la Mora Reference Gamero-Castaño and Fernandez de la Mora2000; Gamero-Castaño & Magnani Reference Gamero-Castaño and Magnani2019), and the droplets and filaments in the breakup region. Because
$E_n^o$
is substantially larger in the smallest droplets than in the surface of the jet (Misra & Gamero-Castaño Reference Misra and Gamero-Castaño2022), ions are emitted from the breakup region at flow rates for which
$E_n^o$
at the base of the jet is too low to promote significant ion emission. This is also supported by Gamero-Castaño & Cisquella-Serra (Reference Gamero-Castaño and Cisquella-Serra2021), who report the existence of small EMI-Im droplets whose original net charge is depleted by ion emission at flow rates in which ion emission from the base of the jet is not observed. However, as the flow rate decreases, the electric field at the base of the jet increases (primarily due to self-heating), and we expect that ions will be emitted from both the steady jet and the breakup region at sufficiently low flow rates. We include a preliminary study of this regime in Appendix C.
Figure 1. Model domain and regions. The external boundary is a semicircle of radius
$1000 \times l_c$
centred at the origin. The regions are not to scale to better highlight their features.
Figure 1 shows the computational domain, which is divided into four regions to improve the accuracy of the solution. The cone jet (regions 1, 2 and 3) is separated from the surrounding vacuum (region 4) by the free surface
$R(x)$
. The solution is regarded hydrostatic and isothermal in region 1 due to the large cross-section available for electrical conduction and the flow, and the solution of (2.1a
)–(2.1d
) becomes trivial; in the transition region 2 the full set of equations is numerically solved in two-dimensional space; and in the jet region 3, a slender approximation yields an analytic solution. The boundaries
$\varSigma _{12}$
between regions 1 and 2, and
$\varSigma _{23}$
between regions 2 and 3, are placed at
$-200\, l_c$
and
$200\, l_c$
from the origin of coordinates, respectively. This ensures the validity of the hydrostatic and isothermal assumptions in region 1 and the slender approximation in region 3. The flow field is solved in terms of the streamfunction
$\varPsi$
and the vorticity
$\varOmega$
instead of the velocity components, which in two-dimensional problems allows decoupling of the calculation of the pressure (Higuera Reference Higuera2003; Gamero-Castaño & Magnani Reference Gamero-Castaño and Magnani2019). The equations in region 2 are discretised using finite differences in an orthogonal set of coordinates
$\{\xi ,\eta \}$
where the surface and the axis are two coordinate curves. This frame of reference is constructed using a variation of the method given by Srinivas & Fletcher (Reference Srinivas and Fletcher2002). The coordinate
$\xi$
goes from 0 on
$\varSigma _{12}$
to 1 on
$\varSigma _{23}$
;
$\eta$
, defined as
$r/R(x)$
, goes from 0 on the axis to 1 on
$\varSigma _{24}$
. The coordinate
$\eta$
is also used in region 3. Additionally, on the liquid surface we define the coordinates
$\{t,n\}$
, corresponding to the local tangential and normal directions to the surface, for variables such as
$E_n^i$
,
$E_n^o$
and
$E_t$
along
$\varSigma _{14}$
,
$\varSigma _{24}$
and
$\varSigma _{34}$
. The slender approximation in the jet, region 3, transforms (2.1a
)–(2.1d
) into
\begin{align} \varPsi _{\textit{jet}}&=\frac {\eta ^2}{2}+\frac {\eta ^4-\eta ^2}{8}\left (\frac {Re}{2\mu }R^3E_t\sigma -2{R'}^2-RR''\right)\!, \end{align}
\begin{align} \varOmega _{\textit{jet}}&=\frac {\eta }{2}\left [\frac {R''}{R^2\big (1+{R'}^2\big )}\left (4+\frac {Re}{\mu }R^3E_t\sigma -2{R'}^2-RR''\right )-\frac {Re}{2\mu }E_t\sigma \right ]\!, \end{align}
\begin{align} \varPhi _{i,\textit{jet}}&=\varPhi _s+\frac {RR'\left (\eta ^2-1\right )}{2}E_t, \end{align}
\begin{align} T_{\textit{jet}}&=f_0+\frac {1}{\sqrt {\varPi _Q}}\int _{x_{23}}^x{\left (\frac {12\mu R\prime ^2}{{Re} R^4}+KR^2{E_t}^2\right )}{\rm d}x+\sum _{n=1}^\infty {g_nJ_0 (\eta j_{1,n})e^{\tfrac {{j_{1,n}}^2(x_{23}-x)}{\textit{Pe}}}}, \end{align}
where
$\varPhi _s$
is the electric potential on the surface,
$J_0$
and
$J_1$
are the Bessel functions of the first kind of orders zero and one,
$j_{1,n}$
is the nth zero of
$J_1$
and
$x_{23}$
is the axial position that divides the surface between
$\varSigma _{24}$
and
$\varSigma _{34}$
; the coefficients
$f_0$
and
$g_n$
are given by
\begin{align} f_0&=2\int _0^1{\eta T(\varSigma _{23})}{\rm d}\eta , \end{align}
\begin{align} g_n&=\frac {2}{{J_0(j_{1,n})}^2}\int _0^1{\eta J_0(\eta j_{1,n})\left [T(\varSigma _{23})-f_0\right ]}{\rm d}\eta . \end{align}
Magnani & Gamero-Castaño (Reference Magnani and Gamero-Castaño2024) describe in detail the construction of the orthogonal grid and the jet solution.
As boundary conditions, for the electrostatic problem we impose the Taylor potential (Taylor Reference Taylor1964) on
$\varSigma _{04}$
\begin{align} \varPhi _T(x,r)&=\frac {2\pi \sqrt {\sin {2\theta _T}}}{\mathcal{K}\left (\dfrac {1+\cos \theta _T}{2}\right )}\big [\big (x^2+r^2\big )\varPi _Q\big ]^{\frac {1}{4}}\left [\mathcal{K}\left (\frac {1}{2}-\frac {x}{2\sqrt {x^2+r^2}}\right )\right .\nonumber\\&\quad \left .-2\mathcal{E}\left (\frac {1}{2} -\frac {x}{2\sqrt {x^2+r^2}}\right )\right]\!, \end{align}
where
$\mathcal{K}$
and
$\mathcal{E}$
are the complete elliptic integrals of the first and second kinds and
$\theta _T=49.29^\circ$
is the Taylor angle. Here,
$\varSigma _{01}$
is the equipotential at the value of its intersection with
$\varSigma _{04}$
. In the jet, (2.4c
) yields the potential inside the liquid from the potential on the surface, so only the potential at the intersection of
$\varSigma _{03}$
and
$\varSigma _{04}$
needs to be imposed. The temperature, streamfunction and vorticity are solved numerically only in region 2: we impose
$T=T_0$
on
$\varSigma _{12}$
, zero heat flux on
$\varSigma _{24}$
, equation (2.3d
), and
${\partial T}/{\partial \xi }= {\partial T_{\textit{jet}}}/{\partial \xi }$
on
$\varSigma _{23}$
; for the streamfunction and vorticity we impose,
\begin{align} \frac {\partial ^2\varPsi }{\partial n^2}-\frac {r'r''}{\big (1+{r'}^2\big )^{3/2}}\frac {\partial \varPsi }{\partial n}&=r\varOmega , \end{align}
\begin{align} \frac {\partial ^2\varOmega }{\partial n^2}+\left (\frac {2}{r}-\frac {r''}{1+{r'}^2}\right )\frac {r'}{\sqrt {1+{r'}^2}}\frac {\partial \varOmega }{\partial n}&=0 , \\[6pt] \nonumber \end{align}
on
$\varSigma _{12}$
; on
$\varSigma _{24}$
the kinematic condition (2.3c
) requires
$\varPsi (\xi ,1)=1/2$
, while the equilibrium of tangential surface stresses (2.3e
) provides a condition for the vorticity
\begin{equation} \varOmega (\xi ,1)=\frac {2R''}{R^2\big (1+{R'}^2\big )}\frac {\partial \varPsi }{\partial \eta }-\frac {Re}{2\mu }E_t\sigma ; \end{equation}
finally, we impose continuity,
$\varPsi = \varPsi _{\textit{jet}}$
and
$\varOmega = \varOmega _{\textit{jet}}$
, on
$\varSigma _{23}$
. Note that conditions (2.7a
) and (2.7b
) are the field equations for the streamfunction and the vorticity for
$R\gg 1$
. This facilitates the convergence of the numerical solution, as it eliminates a small jump in stress across
$\varSigma _{14}$
and
$\varSigma _{24}$
produced by the simpler Neumann boundary conditions.
Although the cone-jet model does not include the beam of droplets and its space charge, its effect on the electric field of the transition region is partially reproduced by the long jet considered in the model, which extends well beyond the breakup region and into the beam (the axial distributions of charge in this jet and in the beam are similar). Furthermore, the space charge has a small effect in the transition region, as cone jets operating in vacuum and atmospheric conditions have identical current versus flow rate characteristics despite having very different beam structures and space charge (Gamero-Castaño et al. Reference Gamero-Castaño, Aguirre-de Carcer, de Juan and Fernández de la Mora1998). Ultimately, the omission of the beam in the model is supported by the good agreement between the numerical solution and experiments for isothermal cone jets, in particular the numerical solution reproduces extremely well the scaling laws for both the current and the radius of the jet (Gañán-Calvo et al. Reference Gañán-Calvo, López-Herrera, Herrada, Ramos and Montanero2018; Gamero-Castaño & Magnani Reference Gamero-Castaño and Magnani2019).
Figure 2. Solution scheme.
Figure 2 illustrates the solution scheme. The initial guesses for streamfunction and vorticity are obtained from (2.4a
) and (2.4b
), which work well everywhere in the liquid. After discretising the domain, the thermal (2.1c
, 2.2a
, 2.2b
, 2.3d
, 2.4d
), electric (2.1d
, 2.1e
, 2.3a
, 2.3b
, 2.4c
) and fluid dynamic problems (2.1a
, 2.1b
, 2.3c
, 2.3e
, 2.4a
, 2.4b
) are solved using an iterative procedure. This partial solution (we are solving all model equations except for the balance of normal stresses (2.3f
)) is considered converged once the difference between consecutive solution vectors
$\chi$
is such that
$\sum _{j=1}^{j=N}{|\chi ^{new}_j/\chi ^{old}_j-1|}/N \lt 10^{-6}$
. We then compute the normalised residual of (2.3f
)
\begin{equation} \mathcal{R}=\frac {\int _{x_0}^{x_f}{\left |S_c-S_p+S_\mu -S_E\right |}\,{\rm d}x}{\frac {1}{4}\int _{x_0}^{x_f}{\left |S_c\right |+\left |S_p\right |+\left |S_\mu \right |+\left |S_E\right |}\,{\rm d}x}, \end{equation}
where
\begin{align} S_c&=\frac {1}{R\sqrt {1+{R'}^2}}-\frac {R''}{\big (1+{R'}^2\big )^{3/2}}, \end{align}
\begin{align} S_p&=\frac {\sqrt {\varPi _Q}}{2\pi ^2}\left [\frac {1+{R'}^2}{R^4}\frac {\partial \varPsi }{\partial \eta }^2\right ]_{x_0}^x-\frac {1}{\pi {Re}}\int _{x_0}^x{\mu \left (\frac {\varOmega }{R}+\frac {1+{R'}^2}{R}\frac {\partial \varOmega }{\partial \eta }\right )}\nonumber\\&\quad-\frac {2\sqrt {1+{R'}^2}}{R^2}\frac {\partial \mu }{\partial t}\left [\left (\frac {R'R''}{1+{R'}^2}-\frac {2R'}{R}\right )\frac {\partial \varPsi }{\partial \eta }+\sqrt {1+{R'}^2}\frac {\partial ^2\varPsi }{\partial t\partial \eta }\right ]{\rm d}x, \end{align}
\begin{align} S_{\mu }&=\frac {2\mu }{\pi {Re} R^2}\left [\left (\frac {R'}{R}-\frac {R'R''}{1+{R'}^2}\right )\frac {\partial \varPsi }{\partial \eta }-\sqrt {1+{R'}^2}\frac {\partial ^2\varPsi }{\partial t\partial \eta }\right ], \end{align}
\begin{align} S_E&=\frac {1}{2\pi ^2\sqrt {\varPi _Q}}\big [{E_n^o}^2-\varepsilon {E_n^i}^2-{E_t}^2\left (\varepsilon -1\right )\big], \end{align}
are the capillary stress, surface pressure, viscous stress and electric stress, respectively. If
$\mathcal{R} \lt 10^{-3}$
the full solution is considered converged, otherwise we improve the position of the surface
$R(x)$
by minimising the residual of (2.3f
) using the false transient method (Mallinson & de Vahl Davis Reference Mallinson and de Vahl Davis1973; Northrop et al. Reference Northrop, Ramachandran, Schiesser and Subramanian2013), which adds a pseudo-time derivative to (2.3f
) to update the position of the surface. To avoid potential oscillations during integration we add an artificial dissipation term (VonNeumann & Richtmyer Reference VonNeumann and Richtmyer1950; Mattsson & Rider Reference Mattsson and Rider2015) and use a downwind discretisation scheme for some derivatives. The resulting equation is
\begin{equation} \frac {{\rm d}R}{{\rm d}\tau }=S_c-S_p(R'_d)+S_{\mu }(R'_d)-S_E + \alpha \frac {\max \left |S_c-S_p+S_{\mu }-S_E\right |_{\tau =0}}{\max |R''|_{\tau =0}}R'', \end{equation}
where
$\tau$
stands for the pseudo-time. The last term adds artificial dissipation with the ad hoc factor
$\alpha$
controlling its magnitude. As the solution asymptotes to steady state both the time derivative and the artificial dissipation term become zero, recovering the original equation (2.3f
). We highlight the dependence of the pressure and the viscous stress terms on
$R'_d$
because this first derivative is discretised using a 3-point downwind scheme, while all other derivatives are discretised using a central scheme. In each iteration step, (2.11) is integrated from
$\tau =0$
to
$\tau _I=10^{-4}$
using a regular ordinary differential equation (ODE) solver. The values of
${\partial \varPsi }/{\partial \eta }$
,
${\partial ^2\varPsi }/{\partial t\partial \eta }$
,
$\varOmega$
,
$ {\partial \varOmega }/{\partial \eta }$
,
$\mu$
,
$ {\partial \mu }/{\partial t}$
,
$E_n^o$
,
$E_n^i$
and
$E_t$
are not updated during integration, as the position of the surface is expected to change slightly. This depends on the value of the residual, and at the start of the optimisation, when
$\mathcal{R}$
is the largest, we reduce the integration time
$\tau _I$
.
We simulate cone jets of four ionic liquids, EMI-Im, BMI-TCM, EMI-TFA and EAN. Table 3 lists their physical properties and characteristic scales. The properties of the liquids are obtained from the National Institute of Standards and Technology (NIST) ionic liquid database (Kazakov et al. Reference Kazakov, Magee, Chirico, Paulechka, Diky, Muzny, Kroenlein and Frenkel2022), with the exception of the surface tension of EMI-TFA (Fang et al. Reference Fang, Zhang, Jia, Shan, Xia and Yang2017). The dielectric constants of BMI-TCM and EMI-TFA are obtained from the list provided by Rybinska-Fryca, Sosnowska & Puzyn (Reference Rybinska-Fryca, Sosnowska and Puzyn2018), which use the quantitative structure property relationship method to estimate
$\varepsilon$
for many different liquids.
3. Results and discussion
Figure 3 illustrates a typical solution for EAN and
$\varPi _Q=50$
, comparing several features of the dissipative and isothermal solutions. The radius of the cone jet exhibits a sharp transition between the conical region and the jet, more accentuated when self-heating is taken into account. The bulk conduction current,
$\tilde {I}_b=2\pi \int _0^R{K\tilde {E}_x}\tilde {r}{\rm d}\tilde {r}$
, is dominant upstream and transforms into surface current,
$\tilde {I}_s=2\pi \tilde {R} \tilde {\sigma } \tilde {v}_s$
, along the jet. Both transport mechanisms become equal at the so-called current cross-over point. The dissipative solution yields a significantly larger total current,
$\tilde {I}=\tilde {I}_b+\tilde {I}_s=237.2$
nA, than the isothermal solution, 77.5 nA. Figure 3(b) shows the dissipation linear densities
\begin{align} P'_\varOmega (x)&=2\int _0^{\tilde {R}(x)}{r P^{\prime\prime\prime}_\varOmega (x,r)}\,{\rm d}r, \end{align}
\begin{align} P'_\mu (x)&=2\int _0^{\tilde {R}(x)}{r P^{\prime\prime\prime}_\mu (x,r)}\,{\rm d}r. \\[6pt] \nonumber \end{align}
Figure 3. Numerical solution for EAN,
$\varPi _Q=50$
. Red lines depict features of the dissipative solution, while blue lines are used for the isothermal solution: (a) cone-jet radius, bulk conduction current and surface current; (b) ohmic and viscous dissipation linear densities; (c) temperature increase, electrical conductivity and viscosity along the axis (dissipative solution); and (d) two-dimensional map of the temperature increase (dissipative solution).
When self-heating is accounted for the ohmic and viscous dissipations are respectively much larger and much smaller than in the isothermal solution. This is due to the temperature increase along the cone jet and the associated enhancement of the electrical conductivity and the reduction of the viscosity. The effect of self-heating is strong enough to reverse their relative importance, ohmic dissipation being much more intense than viscous dissipation in the dissipative solution while the latter is dominant in the isothermal solution. Figure 3(c) shows the variation of the temperature, the electrical conductivity and the viscosity along the axis. The total temperature increase is 301 K, the electrical conductivity increases by a factor of 15.4 and the viscosity decreases by a factor of 122. Such large changes in physical properties, especially in the electrical conductivity, make the isothermal assumption grossly inaccurate. The density plot in figure 3(d) shows that temperature variations along the radial coordinate are negligible in the slender jet.
Figure 4. Total current as a function of the dimensionless flow rate for EMI-Im (a), BMI-TCM (b), EMI-TFA (c) and EAN (d). The characteristic scale
$I_o$
for each liquid is
$2.77 \times 10^{-9}$
A for EMI-Im,
$4.58 \times 10^{-9}$
A for BMI-TCM,
$4.07 \times 10^{-9}$
A for EMI-TFA and
$4.26 \times 10^{-9}$
A for EAN. We also plot the best fitting function for the combined isothermal data,
$\tilde {I}/I_o= 2.52+2.32{\varPi _Q}^{1/2}-0.0049\varPi _Q$
.
Figure 4 shows the total electrospray current as a function of the dimensionless flow rate for all ionic liquids investigated. The charts include the dissipative and isothermal solutions and experimental data (Caballero-Pérez et al. Reference Caballero-Pérez, Galobardes-Esteban and Gamero-Castaño2025). The current is made dimensionless with the scale
$I_o=\sqrt {\varepsilon _0\gamma ^2/\rho }$
, to facilitate comparison with the traditional scaling law
$\tilde {I}/I_o\propto {\varPi _Q}^{1/2}$
(Fernandez de la Mora & Loscertales Reference Fernandez de la Mora and Loscertales1994; Gañán-Calvo et al. Reference Gañán-Calvo, Lasheras, Dávila and Barrero1994). The experimental values compare much better with the dissipative solution than with the isothermal solution, confirming the importance of self-heating in these cone jets. These thermal effects generate a strong positive feedback loop between the electric current and the electrical conductivity (Magnani & Gamero-Castaño Reference Magnani and Gamero-Castaño2024). The dissipative solution reproduces well the experimental values for EMI-Im and EMI-TFA, and is 10 %–20 % lower for BMI-TCM and EAN. These differences could be explained by several factors, including: (i) the need to extrapolate when estimating the values of the electrical conductivity and viscosity, since in general the experimental data used in the calculation of
$K(\tilde {T})$
and
$\mu (\tilde {T})$
extend over temperature ranges narrower than those obtained with the model. For example, we have found electrical conductivity data for EAN only in the range from 288 to 391 K (Kazakov et al. Reference Kazakov, Magee, Chirico, Paulechka, Diky, Muzny, Kroenlein and Frenkel2022), while the temperature range of the numerical solution is between 397 and 716.3 K; (ii) the uncertainty in the experimental values of the physical properties, as illustrated in Appendix A describing the calculation of the temperature-dependent functions for the physical properties; (iii) neglecting the variation with temperature of physical properties other than viscosity and electrical conductivity. Appendix B extends the present model to include the variation with temperature of all physical properties; (iv) the improved agreement for EAN and BMI-TCM at the lowest flow rates is likely an artefact of how Caballero-Pérez et al. (Reference Caballero-Pérez, Galobardes-Esteban and Gamero-Castaño2025) executed the experiments. An important goal of this study was to lower the minimum stable flow rate of these cone jets. This was achieved by using emitters with reduced inner diameters supporting smaller Taylor cones. The emitter potential was also decreased with flow rate to maintain the volume of the Taylor cone approximately constant. For example, the lowest flow rates for EAN were recorded with a Taylor cone with a base diameter of 15
$\unicode{x03BC}$
m and an emitter potential of 1150 V, while the medium and large flow rates were recorded with a 50
$\unicode{x03BC}$
m Taylor cone at 1576 V. As a result, it is likely that, although the EAN and BMI-TCM electrosprays at the lowest flow rates were stable, their electrification levels were too low, causing the electrospray current to fall below that of a cone jet. The isothermal solutions for all ionic liquids are nearly indistinguishable, and the combined data are well fitted by a single function,
$\tilde {I}/I_o= 2.52+2.32{\varPi _Q}^{1/2}-0.0049\varPi _Q$
. At flow rates
$\varPi _Q \lesssim 250$
, the data are well fitted by
$\tilde {I}/I_o = 2.45{\varPi _Q}^{1/2} + 1.01$
, in exceptional agreement with the best fit reported by Gañán-Calvo et al. (Reference Gañán-Calvo, López-Herrera, Herrada, Ramos and Montanero2018),
$\tilde {I}/I_o=2.5{\varPi _Q}^{1/2}$
, for experimental data including 54 different electrolytes with much lower electrical conductivities than EMI-IM, BMI-TCM, EMI-TFA or EAN, and therefore having negligible self-heating. The 2.45 factor also compares well with the values of 2.33 and 2.45 reported by Fernandez de la Mora & Loscertales (Reference Fernandez de la Mora and Loscertales1994) for two liquids, benzyl-OH (
$\varepsilon = 13.1$
) and 3-ethylene glycol (
$\varepsilon = 23.7$
), with similar dielectric constants to those presently simulated. In contrast,
$\tilde {I}/I_o$
for the dissipative solution depends not only on the dimensionless flow rate, but also on the ionic liquid itself. This behaviour, noted in previous studies (Gamero-Castaño Reference Gamero-Castaño2019; Magnani & Gamero-Castaño Reference Magnani and Gamero-Castaño2024), cannot be attributed to a dependency on the Peclet number, which is large in all simulations making thermal conduction negligible. Instead, the breakdown of the simple
$I(\varPi _Q)$
dependency in ionic and other highly conducting liquids is due to the spatial variation of the electrical conductivity and viscosity, the unique
$K(\tilde {T})$
and
$\mu (\tilde {T})$
relationships for each liquid and the importance of self-heating.
Although the radius of the jet varies along the axis as the electric field accelerates the liquid, a characteristic jet radius
$R_c$
can be defined due to its slow variation,
$R\propto x^{-1/8}$
(Gañán-Calvo Reference Gañán-Calvo1997). Measurements of the diameters of droplets (Gañán-Calvo et al. Reference Gañán-Calvo, López-Herrera, Herrada, Ramos and Montanero2018) and of the jet at the breakup (Gamero-Castaño & Hruby Reference Gamero-Castaño and Hruby2002; Gamero-Castaño Reference Gamero-Castaño2010) support the validity of the characteristic radius
$R_c = (\varepsilon _0\rho Q^3/(\gamma K) )^{1/6}$
, first proposed by Gañán-Calvo et al. (Reference Gañán-Calvo, Lasheras, Dávila and Barrero1994). Figure 5 plots the radius of the jet at the current cross-over point as a function of the dimensionless flow rate, for both the dissipative and the isothermal solutions. The lines
$\tilde {R}/d_0=\sqrt {\varPi _Q}$
and
$\tilde {R}/d_0=0.5\sqrt {\varPi _Q}$
bound the data, showing that
$R_c$
describes well the jet radius
\begin{equation} \tilde {R}(x)=\alpha \left (\frac {\varepsilon _0\rho Q^3}{\gamma K} \right )^{1/6} = \alpha d_0 {\varPi _Q}^{1/2}, \end{equation}
where
$\alpha$
is a factor of order unity and
$d_0$
a length defined by physical properties only,
$d_0 = ({\varepsilon _0}^2\gamma /(\rho K^2) )^{1/3}$
. However, there is a significant difference between the dissipative and the isothermal solutions: while the radii at the current cross-over collapse into a single-valued function for the latter, the dissipative radii exhibit considerable spread. As in the case of the total current, this spread is due to the importance of self-heating, the spatial variation of the electrical conductivity and viscosity and the unique
$K(\tilde {T})$
and
$\mu (\tilde {T})$
functions for each liquid.
Figure 5. Jet radius at the current cross-over point as a function of the dimensionless flow rate. The radius is normalised with
$d_0 = ({\varepsilon _0}^2\gamma /(\rho K^2))^{1/3}$
.
Figure 6. Temperature increase (red) and electrical conductivity (black) at the current cross-over and at
$I_s/I = 0.95$
, for EMI-Im (a), BMI-TCM (b) EMI-TFA (c) and EAN (d). The dashed lines show the regression law for the temperature found by Magnani & Gamero-Castaño (Reference Magnani and Gamero-Castaño2024), equation (3.3).
Figure 6 shows the increase in temperature and the electrical conductivity at the current cross-over, and at the axial position where the surface current reaches 95 % of the total current. Most dissipation has taken place before the latter point, and the temperature and conductivity have nearly reached their maximum values. The current cross-over is important because it is representative of the centre of the transition region, where variations in physical properties have the strongest effect in the total current and the jet diameter. The upstream temperature is 294.15 K in all calculations. The temperature increase is substantial for all liquids and intensifies with decreasing flow rate: at
$I_{s}/I = 0.95$
, the temperature increase ranges from 28.3 to 279.2 K, 58.7–240.4 K, 38.5–445.9 K and 124.4–443.1 K for EMI-Im, BMI-TCM, EMI-TFA and EAN, respectively; at the current cross-over point, the temperature increase ranges from 14.1 to 189.8 K, 31.4–157.1 K, 18.8–291.4 K and 67.8–292.4 K. The electrical conductivity at the current cross-over varies from 1.21 to 10.2 S m−1, 2.20–6.19 S m−1, 1.59–23.1 S m−1 and 8.46–30.1 S m−1 for EMI-Im, BMI-TCM, EMI-TFA and EAN, respectively. EAN and EMI-TFA are the two liquids which experience the highest temperature increase, over 440 K. The first one was expected, given that EAN is the liquid with the highest conductivity, while the latter is a result of both a relatively high initial conductivity (the second highest) and the steepest increase in conductivity of all liquids: an increase of almost 40 times compared with 20 for EAN. It is apparent that dissipation and self-heating significantly alter the physical properties of the liquid and strongly modify these cone jets. Figure 6 also shows a regression law for the temperature increase at
$I_{s}/I =0.95$
, obtained by Magnani & Gamero-Castaño (Reference Magnani and Gamero-Castaño2024):
\begin{equation} \Delta T\approx \frac {7.74}{ {Re}^{0.11}{\varPi _Q}^{0.11}}+\frac {12.8+4.53/ {Re}}{{\varPi _Q}^{0.73}}. \end{equation}
In the present calculations with ionic liquids, this simple expression works well at medium and high flow rates, but overestimates the temperature increase at the lowest flow rates. The disagreement is caused by the second term on the right-hand side of (3.3), associated with viscous dissipation and which becomes dominant at low flow rates. Two shortcomings in (3.3) are responsible for this: first, it is just a function of
$\varPi _Q$
and
${\textit{Re}}$
and does not account for a dependence on the specific liquid; second, it was obtained for cone jets with larger Reynolds numbers and therefore lower temperature increases than in the present study,
${Re} \geqslant 0.095$
and
$\Delta \tilde {T}\leqslant 15.7$
K. However, the very high temperatures at the lowest flow rates characteristic of ionic liquids greatly reduce the viscosity and the viscous dissipation. Under these conditions the second term in (3.3) overpredicts the temperature increase caused by viscous dissipation. Note also that the temperature increase diminishes slowly with increasing flow rate,
$\Delta T\rightarrow ( {7.74}/{{Re}^{0.11}{\varPi _Q}^{0.11}})$
. This explains the slow (or absence of) convergence at large flow rates of the ‘dissipative’ and ‘isothermal’ currents in figure 4, since their ratio for a given dimensionless flow rate approximately is
$\sqrt {K(T_o+\Delta T)/K(T_o)}$
, which is nearly constant at large flow rates.
The high increase in temperature experienced by these liquids could potentially result in significant radiation heat losses and lower temperature increase. To estimate this effect, we use the Stefan–Boltzmann law to evaluate the radiative heat flux from the surface and compute the power radiated from the cone jet as
\begin{equation} \tilde {P}_{\textit{rad}}(\tilde {x}) =\int _{-\infty }^x k_{SB}\epsilon F \big (\tilde {T}^4-{\tilde {T}_0}^4\big )2\pi \tilde {R}{\rm d}\tilde {x}, \end{equation}
where
$k_{SB}$
is the Stefan–Boltzmann constant,
$\epsilon$
the emissivity of the surface and
$F$
the view factor between the liquid and the ambient. To obtain an upper bound, we consider both
$\epsilon$
and
$F$
equal to 1. The parameter
$\tilde {P}_{\textit{rad}}$
can then be compared with the ohmic and viscous dissipation to determine whether radiation can significantly alter the temperature increase along the cone jet. Table 4 shows the maxima of
$\tilde {P}_{\textit{rad}}(\tilde {x})/\tilde {P}_\varOmega (\tilde {x})$
and
$\tilde {P}_{\textit{rad}}(\tilde {x})/\tilde {P}_\mu (\tilde {x})$
for all ionic liquids at two disparate dimensionless flow rates. The power radiated is several orders of magnitude smaller than the power dissipated, and therefore it does not have a significant effect on the temperature increase and can be neglected.
Table 4. Maximum ratio between radiated power and ohmic and viscous dissipation for dimensionless flow rates of 10 and 1000.
The high temperatures attained along the cone jet may lead to thermal decomposition of the ionic liquid. Although there is variation in the reported decomposition temperatures of ionic liquids (especially as a function of the heating rate), typical minimum values for EAN are from 196 to 278 K (Salgado et al. Reference Salgado, Parajó, Villanueva, Rodríguez, Cabeza and Varela2019), while for EMI-Im minimum decomposition temperatures range from 600 to 750 K (Beigi et al. Reference Miran Beigi, Abdouss, Yousefi, Pourmortazavi and Vahid2013; Zaitsau & Abdelaziz Reference Zaitsau and Abdelaziz2020; Yu et al. Reference Yu, Liu, Xia and Wu2022; Hung & Pan Reference Hung and Pan2023). However, thermal decomposition is a slow process. For example, the maximum decomposition rate for EAN has an associated time constant of approximately 0.01 s (Salgado et al. Reference Salgado, Parajó, Villanueva, Rodríguez, Cabeza and Varela2019), while the time constant for EMI-Im is 3.1 s (Zaitsau & Abdelaziz Reference Zaitsau and Abdelaziz2020). These times are many orders of magnitude larger than the residence time of the liquid in the jet (defined as the time it takes for the liquid to move between the points where the surface current is 10 % and 90 % of the total current), which in our simulations range between 0.23 and 224 ns depending on the flow rate. Therefore, thermal decomposition of the liquid is likely to be insignificant.
Figure 7 shows the total ohmic
$P_\varOmega$
and viscous
$P_\mu$
dissipations as functions of the dimensionless flow rate and the Reynolds number, for both the dissipative and the isothermal solutions. Here,
$P_\varOmega$
and
$P_\mu$
are the integrals of the dissipation densities in the cone jet, down to the
$I_{s}/I = 0.95$
point. Self-heating has opposite effects on the ohmic and viscous dissipations, as a higher temperature increases the electrical conductivity and decreases the viscosity (
$P^{\prime\prime\prime}_\varOmega \propto K$
and
$P^{\prime\prime\prime}_\mu \propto {\mu }$
). The ohmic and viscous dissipations for the isothermal solution are similar to those identified previously,
$P_\varOmega \approx 7.4{\varPi _Q}^{0.41}$
and
$P_\mu \propto 1/ {Re}$
, both in isothermal simulations (Gamero-Castaño Reference Gamero-Castaño2019) and in simulations including self-heating but with reduced temperature variations (Magnani & Gamero-Castaño Reference Magnani and Gamero-Castaño2024); the Reynolds numbers considered in these two studies,
$0.1\lesssim 1/ {Re}\lesssim 11$
, were much larger than in the current simulations. On the other hand, the strong self-heating occurring in ionic liquids significantly increases the ohmic dissipation due to the positive feedback with the electrical conductivity, while it reduces the viscous dissipation and even reverses the increasing
$P_\mu \propto 1/ {Re}$
trend for the smallest Reynolds number, i.e. for EAN. The EMI-TFA, which has a maximum temperature increase comparable to the one of EAN, also exhibits a trend reversal at the lowest flow rates. Furthermore, the liquid-specific
$K(\tilde {T})$
and
$\mu (\tilde {T})$
functions combined with the strong self-heating prevent the collapse of the data into functions of only the dimensionless flow rate and the Reynolds number.
Figure 7. Total ohmic (a) and viscous (b) dissipations in the cone jet.
Our cone-jet model is steady state and therefore it does not incorporate the breakup region. Cone jets are steady-state and stable systems upstream of the breakup region. The breakup instability does not propagate upstream and therefore it is appropriate to replace the time-dependent breakup region with a fictitious steady-state jet when modelling the transition region and the Taylor cone upstream. In our model, this fictitious jet extends far downstream, well beyond the breakup region. This reduces the sensitivity of the solution to the size of the computational domain. The validity of replacing the breakup region with a fictitious jet is supported by the agreement between experiments and the numerical solution of isothermal jets. For example, the numerical solution reproduces extremely well the scaling laws for both the current and the radius of the jet (Gamero-Castaño & Magnani Reference Gamero-Castaño and Magnani2019; Gañán-Calvo et al. Reference Gañán-Calvo, López-Herrera, Herrada, Ramos and Montanero2018). Although our model does not yield the axial location of the region where the jet starts to become unstable, if this axial location is known by alternative means the numerical solution can be used to analyse properties of this region such as the electric field and the radius of the jet. For example, Gamero-Castaño & Cisquella-Serra (Reference Gamero-Castaño and Cisquella-Serra2021) provide measurements and estimates for several features of the breakup region of cone jets of EMI-Im in the range
$612\lesssim \varPi _Q \lesssim 3630$
. These include the voltage difference between the emitter electrode and the breakup point
$\Delta \varPhi _{J}$
, as well as the velocity, the radius and the electric field at the position
$x_b$
of the breakup. Figure 8(a) shows the experimental voltage difference as a function of the dimensionless flow rate, which follows the power law
\begin{equation} \Delta \varPhi _{J} = 12.91+9.787\sqrt {\varPi _Q}. \end{equation}
When combined with the numerical solution, this relationship yields the axial location of the breakup point, depicted in figure 8(b). These data are approximated by the straight line
\begin{equation} \Delta \varPhi _{J} = 1.225 x_b - 51.30 ,\end{equation}
or equivalently
\begin{equation} x_b=8.133\sqrt {\varPi _Q}+48.03. \end{equation}
This linear relationship is anticipated from (3.5) and the dimensionless Taylor potential,
$\varPhi _T\propto {\varPi _Q}^{1/4} x^{1/2}$
. We do not have experimental data to validate these functions at the lowest flow rates,
$\varPi _Q \lt 612$
. However, we observe that both fittings (3.5) and (3.7) exhibit positive
$y$
-intercepts, which is consistent with a shortening transition region at decreasing flow rate, yet maintaining a finite length and potential drop at the lowest flow rates.
Figure 8. (a) Emitter potential minus potential at jet breakup for EMI-Im cone jets, experimental data (Gamero-Castaño & Cisquella-Serra Reference Gamero-Castaño and Cisquella-Serra2021); (b) emitter potential minus
$\varPhi (x)$
, and breakup position (circles). Solid circles indicate a flow rate within the range measured by Gamero-Castaño & Cisquella-Serra (Reference Gamero-Castaño and Cisquella-Serra2021). The characteristic potential for EMI-Im is
$\varPhi _c = 0.772$
V.
Figure 9 shows the velocity, the radius and the normal component of the outer electric field of EMI-Im jets at the breakup position, estimated with (3.7). It also includes values for the isothermal solution and measurements by Gamero-Castaño & Cisquella-Serra (Reference Gamero-Castaño and Cisquella-Serra2021). The dissipative solution reproduces well the experimental results, whereas the isothermal solution significantly underestimates both the jet velocity and the electric field. The differences between the dissipative and isothermal solutions increase at decreasing flow rate.
Figure 9. Jet features at the breakup for EMI-Im, dissipative and isothermal solutions and experimental data (Gamero-Castaño & Cisquella-Serra Reference Gamero-Castaño and Cisquella-Serra2021): (a) jet velocity; (b) jet radius; and (c) outward normal electric field.
Ionic liquids have extremely low vapour pressures (Horike et al. Reference Horike, Ayano, Tsuno, Fukushima, Koshiba, Misaki and Ishida2018). However, the evaporation rate increases exponentially with temperature and is further enhanced by the small radius of curvature of the surface. These combined effects may lead to significant evaporation in the cone jet, potentially invalidating the numerical solution. The total evaporation rate can be estimated by the integral
\begin{equation} \dot {m}_{\textit{ev}}=\int _S{\tilde {p}_{v\textit{ap}}\sqrt {\frac {M_m}{2\pi R_g\tilde {T}}}}{\rm d}S ,\end{equation}
over the surface of the cone jet, with the vapour pressure
\begin{equation} \tilde {p}_{v\textit{ap}}=\exp \left (\beta -\frac {\Delta H_{v\textit{ap}}}{R_g\tilde {T}}\right )\exp \left (\frac {2\gamma M_m}{\rho R_g\tilde {R}_c\tilde {T}}\right )\!. \end{equation}
Here,
$R_g$
,
$M_m$
and
$H_{v\textit{ap}}$
are the gas constant, the molecular mass and the enthalpy of vaporisation respectively;
$\beta$
is a constant, estimated from experimental data. The first exponential term in (3.9) is the Clausius–Clapeyron equation (Brown Reference Brown1951), while the second exponential term is a correction based on the curvature of the surface,
$\tilde {R}_c$
(Mitropoulos Reference Mitropoulos2008). We integrate (3.8) down to the point where the surface current is 95 % of the total current. Table 5 lists the evaporation losses for EMI-Im and EAN at the lowest flow rates studied, i.e. for the highest temperatures along the jet. Data for
$H_{v\textit{ap}}$
and
$\beta$
were obtained from Heym et al. (Reference Heym, Korth, Etzold, Kern and Jess2015), Ahrenberg et al. (Reference Ahrenberg, Beck, Neise, Keßler, Kragl, Verevkin and Schick2016) and Santos et al. (Reference Santos, Ferreira, Štejfa, Rodrigues, Rocha, Torres, Tavares and Carpinteiro2018) for EMI-Im, and from Salo et al. (Reference Salo, Westerlund, Andersson, Nielsen, D’Anna and Hallquist2011) for EAN. Despite the very high temperatures, thermal evaporation is a negligible fraction of the flow rate fed to the cone jet.
Table 5. Molecular mass, enthalpy of vaporisation, constant
$\beta$
in (3.9), evaporated mass flow rate and ratio between evaporated mass flow rate and input mass flow rate, for EMI-Im and for EAN. The evaporated mass flow rate is computed at the lowest flow rate investigated,
$\varPi _Q=10$
.
Figure 10. Cone-jet radius, normalised surface current
$I_s/I$
and terms in power balance (3.10), as functions of the axial position for EMI-Im cone jets,
$\varPi _Q =$
600 (a) and 3500 (b). The dashed section of the
$P_p$
line corresponds to negative values. The position of the breakup is indicated by a circle.
The electric power
$\tilde {P}_E$
injected in the cone jet is converted into jet kinetic power
$\tilde {P}_k$
, ohmic
$\tilde {P}_\varOmega$
and viscous
$\tilde {P}_\mu$
dissipations, and flow work associated with the pressure
$\tilde {P}_p$
and viscous
$\tilde {P}_v$
stresses. This balance, evaluated along the cone jet down to the axial position
$\tilde {x}$
is given by (Gamero-Castaño Reference Gamero-Castaño2010; Gamero-Castaño & Magnani Reference Gamero-Castaño and Magnani2019)
\begin{equation} -\underbrace {\tilde {I}\tilde {\varPhi }_s(\tilde {x})}_{\tilde {P}_E}= \underbrace {\frac {\dot {m}\tilde {v}^2(\tilde {x})}{2}}_{\tilde {P}_k} +\underbrace {Q\tilde {p}(\tilde {x})}_{\tilde {P}_p} +\underbrace {2Q\mu \left .\frac {\partial \tilde {v}}{\partial \tilde {x}}\right |_{\tilde {x}}}_{\tilde {P}_v} +\tilde {P}_\varOmega (\tilde {x}) +\tilde {P}_\mu (\tilde {x}) .\end{equation}
Figure 10 plots this balance for dimensionless flow rates of 600 and 3500. These flow rates are within the range studied by Gamero-Castaño & Cisquella-Serra (Reference Gamero-Castaño and Cisquella-Serra2021), making it possible to place the position of the jet breakup; values plotted downstream of this point correspond to the fictitious section of the jet used in the model for computational purposes. The charts also show the radius of the cone jet and the ratio between the surface and total currents. At these high flow rates, the transition region in which the conduction current becomes surface current is long. The flow work,
$\tilde {P}_p +\tilde {P}_v$
, is always negligible in the balance. The electric power injected before the current cross-over point,
$I_s/I = 0.5$
, is largely dissipated, mostly into ohmic heating. Further ohmic dissipation continues downstream, but the electric energy is increasingly converted in kinetic energy. The surface current plateaus before the kinetic energy equals dissipation (
$I_s/I \gtrsim 0.87$
and
$0.81$
), and the additional work by the electric field mostly leads to an increase in kinetic energy. Although the conduction current continues to decay asymptotically, the tangential electric field decays faster and additional ohmic dissipation is negligible. Totals of 91.7 % and 92.2 % of the total current have become convected surface charge before the breakup, for
$\varPi _Q=$
600 and 3500, respectively. Note that the pressure is negative along most of the jet, due to the high values of the surface charge and its effect on the normal electric stress,
$\tilde {E}^o_n \cong \tilde {\sigma }/\varepsilon _0$
and
${E_n^o}^2 \gt \varepsilon {E_t}^2 \gg \varepsilon {E_n^i}^2$
.
Figure 11. Values of
$\tilde {P}_E/\tilde {I}$
,
$\tilde {P}_k/\tilde {I}$
,
$\tilde {P}_{\varOmega }/\tilde {I}$
and
$\tilde {P}_{\mu }/\tilde {I}$
at the point where the surface current is 95 % of the total current, as a function of the dimensionless flow rate for EMI-Im (a), BMI-TCM, (b) EMI-TFA (c) and EAN (d).
Figure 11 shows the main terms in (3.10),
$\tilde {P}_E$
,
$\tilde {P}_k$
,
$\tilde {P}_\varOmega$
and
$\tilde {P}_\mu$
, computed at the axial position where
$I_s/I = 0.95$
and divided by the total current, as a function of the dimensionless flow rate. By dividing by the current each term has volts units, illustrating how the voltage difference between the emitter and the
$I_s/I = 0.95$
point is distributed into kinetic energy, ohmic and viscous dissipation. The sum of the latter two can be regarded as a voltage loss if the main goal is to convert electric into kinetic energy (as is the case, for example, in electrospray propulsion), and it is a characteristic of the cone-jet independent of the emitter potential and the geometry of the electrodes (Gamero-Castaño Reference Gamero-Castaño2008). In these ionic liquids the ohmic voltage losses dominate over the viscous voltage losses, and both decrease with the flow rate. The total voltage loss for the flow rates investigated ranges from 20.7 to 140 V for EMI-Im, 23.3–117 V for BMI-TCM, 23.9–204 V for EMI-TFA and 15.2–101 V for EAN. The balance of the emitter potential is converted into kinetic energy, both as the liquid gains velocity along the jet and during the flight of charged droplets.
4. Conclusions
This study has extended a previous numerical model (Magnani & Gamero-Castaño Reference Magnani and Gamero-Castaño2024) for cone jets including dissipation and self-heating, in four key areas:
-
(i) we have extensively analysed liquids in which self-heating largely drives the experimental behaviour and the numerical solution;
-
(ii) the numerical model has been improved by implementing a new optimisation method and two strategies to eliminate unphysical oscillations;
-
(iii) we have evaluated the importance of including other temperature-dependent physical properties in addition to the electrical conductivity and viscosity
-
(iv) the model has been extended to include ion field emission from the surface of the steady cone jet, enabling a preliminary study of ion field emission.
Our numerical solutions demonstrate that self-heating due to ohmic and viscous dissipation is an essential phenomenon in cone jets of ionic liquids, and more generally of liquids with high electrical conductivities. In the ionic liquids that we have studied, self-heating leads to temperature increases in the range from 28 to 446 K, at flow rates typically used in experiments. Self-heating creates a strong positive feedback between the electric current and the electrical conductivity, resulting in higher electrospray currents and smaller jet and droplet radii compared with the isothermal solution. As expected from previous studies (Gamero-Castaño Reference Gamero-Castaño2019; Magnani & Gamero-Castaño Reference Magnani and Gamero-Castaño2024), the two cone jets of the ionic liquid with highest electrical conductivity, EAN and EMI-TFA, operated at the lowest flow rates, produce the highest temperatures. Accounting for self-heating in the modelling and in the interpretation of experiments is critical for these cone jets. As a rule of thumb, temperature increases of a few centigrade degrees already occur in liquids with electrical conductivities as low as 0.01 S m−1(Magnani & Gamero-Castaño Reference Magnani and Gamero-Castaño2024). For this electrical conductivity and typical values of the dimensionless flow rate of 10, a surface tension of 0.04 kg s−2 and a density of 1100 kg m−3, the characteristic radius of the jet and droplets is
$R_c \approx 96$
nm. Clearly, self-heating is a key factor in all electrospray applications centred in the generation of nanometric-size jets and droplets.
Magnani & Gamero-Castaño (Reference Magnani and Gamero-Castaño2024) had found that ohmic dissipation is the main contributor to self-heating at most flow rates, and that only for sufficiently small flow rates,
$\varPi _Q \lesssim 1/{Re}$
, viscous dissipation becomes dominant. However, when the temperature increase is large, which is the case for the four ionic liquids investigated at low flow rates, the strong decrease of the viscosity with temperature greatly reduces viscous dissipation, making ohmic dissipation dominant at all flow rates.
The success of simple scaling laws for the current and the diameters of jets and droplets (Fernandez de la Mora & Loscertales Reference Fernandez de la Mora and Loscertales1994; Gañán-Calvo et al. Reference Gañán-Calvo, Lasheras, Dávila and Barrero1994) has played an important role in both the physical understanding and the development of applications of cone jets. Scaling laws such as
$\tilde {I} \cong 2.5 (\gamma K Q )^{1/2}$
and
$\tilde {R} \approx (\varepsilon _0\rho Q^3 / (\gamma K) )^{1/6}$
have been proven valid for values of the flow rate and the physical properties spanning many orders of magnitude (Gañán-Calvo et al. Reference Gañán-Calvo, López-Herrera, Herrada, Ramos and Montanero2018). Unfortunately, these scaling laws are only applicable under isothermal conditions, and break down when self-heating is significant. The need to account for thermal phenomena, the variation of physical properties with temperature (the electrical conductivity and the viscosity being the most influential), and the fact that these dependencies are liquid specific, increase the number of dimensionless groups needed to describe the state of a cone jet, and therefore the complexity of potential scaling laws. Thus, we anticipate that first-principles modelling and case-by-case experimental characterisation will be the main tools for investigating cone jets of liquids with high conductivities. We note that modelling is handicapped by the accuracy and range of the physical properties reported in the literature, especially in the case of ionic liquids: for example the highest temperatures for which we have found electrical conductivity data (e.g. 401 for EMI-Im, 379 for BMI-TCM, 313 for EMI-TFA and 391 K for EAN) fall short of the maximum temperatures in the numerical solutions; and there is significant variance between reported values, especially in the case of the electrical conductivity.
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. Determination of the temperature dependence of liquid properties
The cone-jet model includes several physical properties, namely the density
$\rho$
, surface tension
$\gamma$
, dielectric constant
$\varepsilon$
, thermal conductivity
$k$
, heat capacity
$c$
, viscosity
$\mu$
and the electrical conductivity
$K$
. Most of their values for the four ionic liquids studied are obtained from the experimental data in the NIST ionic liquid database (Kazakov et al. Reference Kazakov, Magee, Chirico, Paulechka, Diky, Muzny, Kroenlein and Frenkel2022). This database includes measurements of several physical properties for more than 2000 ionic liquids, listing them as functions of temperature and pressure. The fitting functions capturing the temperature dependence of the physical properties used in this article were calculated with data collected before the last update of the database, which occurred on 4 June 2024. We limited the data to measurements for the liquid phase at temperatures higher than 283 K and ambient pressure. For the dielectric constant we only use measurements at zero frequency, since the cone-jet model is steady state. The data are collected automatically using the Python library pyILT2. The NIST database does not include information for a few properties, and in these cases we obtain the data from other available publications. For all properties, an experimental datum includes a measured value
$y_i$
, the temperature
$T_i$
associated with the measurement and the measurement uncertainty
$\pm \delta _i$
. The data are separated into subsets, each one associated with a different source.
We have improved the procedure for calculating the fitting functions with respect to that used by Magnani & Gamero-Castaño (Reference Magnani and Gamero-Castaño2024), in order to reduce problems associated with the relatively large scatter of the experimental data. As a result some of the coefficients in the fittings reported by Magnani & Gamero-Castaño (Reference Magnani and Gamero-Castaño2024) for EMI-Im are slightly different from those used presently. The procedure for obtaining the temperature-dependent regression laws includes the following steps:
-
(i) We remove all data not in the liquid phase, or with temperature lower than 283 K, or with a pressure outside the range from 99 to 102 kPa.
-
(ii) If all data have a temperature range of less than 1 K we compute the average values of the property and temperature, and assume that the property does not depend on temperature.
-
(iii) If the data points are at least as many as the fitting parameters, and if they are at sufficiently different temperatures, we compute the fitting by minimising
(A1)where
\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}
$N$
is the total number of data including all subsets, and
$f(\boldsymbol{p},T_i)$
is the fitting function selected. The vector
$\boldsymbol{p}$
is the set of parameters of the fitting function.
-
(iv) Once the fitting is computed, we check if any of the subsets departs significantly from the rest. To do so, we recompute the fitting by omitting one subset at a time, and among these new solutions keep the one with the lowest residual if it is lower than the one including all subsets.
Note that (A1) computes the residual by weighing each data point with its uncertainty. The weighing factor
$1/(1+\delta _i/y_i)$
decreases as the relative uncertainty
$\delta _i/y_i$
increases. The fitting function is selected depending on the liquid property we are studying. We use
\begin{equation} f(\boldsymbol{p},T)=p_a+p_bT ,\end{equation}
for the density, surface tension, dielectric constant, thermal conductivity and heat capacity, and
\begin{equation} f(\boldsymbol{p},T)=p_a\exp \left (\frac {p_b}{T-p_c}\right)\! ,\end{equation}
for the viscosity and electrical conductivity (Okoturo & VanderNoot Reference Okoturo and VanderNoot2004; Leys et al. Reference Leys, Wübbenhorst, Preethy Menon, Rajesh, Thoen, Glorieux, Nockemann, Thijs, Binnemans and Longuemart2008). We have sufficient data distributed at different temperatures to compute the temperature dependence of all physical properties except for the dielectric constant and the heat capacity of EAN. Figures 12–15 show the different data subsets and temperature-dependent fittings of the physical properties for all four ionic liquids.
Figure 12. Experimental data and fitting functions for EMI-Im. The dashed lines represent the standard deviation associated with the fitting function.
Figure 13. Experimental measures and fitting functions for BMI-TCM. The dashed lines represent the standard deviation associated with the fitting function.
Figure 14. Experimental measures and fitting functions for EMI-TFA. The dashed lines represent the standard deviation associated with the fitting function.
Figure 15. Experimental measures and fitting functions for EAN. The dashed lines represent the standard deviation associated with the fitting function.
Figure 16. Comparison between the terms including a gradient of a physical property and the average term of the associated model equation. Both are computed using the numerical solution for EMI-Im,
$\varPi _Q=10$
and assuming constant values of all physical properties, equations (B1a
)–(B2f
) (the property gradients are evaluated using the temperature field of this numerical solution). For plotting purposes, we use the position
$x_0$
of the current cross-over point as the origin of the axial coordinate. The black line in each panel is the average term of the equation, defined by (B5), while the dashed black line corresponds to 1 % of the average term. (a) Equation of mass conservation; (b) equation of momentum conservation; (c) equation of energy conservation; (d) equation of conservation of charge in the liquid bulk; (e) equilibrium of normal stresses on the surface; (f) equilibrium of tangential stresses on the surface.
Figure 17. Relative variation of the liquid properties, computed with the temperature field from the numerical solution for EMI-Im,
$\varPi _Q=10$
and assuming constant values of all physical properties, equations (B1a
)–(B2f
). For plotting purposes, we use the position
$x_0$
of the current cross-over point as the origin of the axial coordinate. Solid lines are used for physical properties that increase with temperature, while dashed lines are used for physical properties that decrease with increasing temperature.
Appendix B. Extending temperature dependence to other physical properties
In previous studies we considered the viscosity and the electrical conductivity as the only temperature-dependent physical properties of the liquid (Magnani & Gamero-Castaño Reference Magnani and Gamero-Castaño2023, Reference Magnani and Gamero-Castaño2024). However, for large temperature variations other liquid properties can also change significantly, and it is important to evaluate the significance of temperature-dependent terms neglected in the current model. When all physical properties are considered temperature-dependent, the original model equations (2.1a )–(2.1e ) and (2.3a )–(2.3f ) become
\begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{v}&=-\frac {\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{\nabla }\rho }{\rho }, \end{align}
\begin{align} \frac {\sqrt {\varPi _Q} Re\rho }{\pi \rho _0}\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{v}+\pi Re\boldsymbol{\nabla }p&=\frac {\mu {\nabla} ^2\boldsymbol{v}}{\mu _0}+\frac {\big (\boldsymbol{\nabla }\boldsymbol{v}+{\boldsymbol{\nabla }\boldsymbol{v}}^T\big )\boldsymbol{\cdot }\boldsymbol{\nabla }\mu }{\mu _0}\nonumber \\ &\quad+\frac {2\boldsymbol{\nabla }\mu \frac {\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{\nabla }\rho }{\rho }}{3\mu _0}+\frac {\mu \boldsymbol{\nabla }\left (\frac {\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{\nabla }\rho }{\rho }\right )}{3\mu _0}, \end{align}
\begin{align} \frac {\rho c}{\rho _0 c_0}\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{\nabla }T-\frac {k{\nabla} ^2T}{k_0 \textit{Pe}}&= \frac {\mu \big (\boldsymbol{\nabla }\boldsymbol{v}+{\boldsymbol{\nabla }\boldsymbol{v}}^T\big )\boldsymbol{:}\boldsymbol{\nabla }\boldsymbol{v}}{\mu _0 Re\sqrt {\varPi _Q}}+\frac {K{\boldsymbol{\nabla }\varPhi _i}\boldsymbol{\cdot }{\boldsymbol{\nabla }\varPhi _i}}{K_0\sqrt {\varPi _Q}}+\frac {{\boldsymbol{\nabla }T}\boldsymbol{\cdot }{\boldsymbol{\nabla }k}}{k_0 \textit{Pe}}\nonumber \\ &\quad-\frac {\rho T\boldsymbol{v}\boldsymbol{\cdot }{\boldsymbol{\nabla }c}}{\rho _0 c_0}+\frac {\boldsymbol{v}\boldsymbol{\cdot }{\boldsymbol{\nabla }\rho }}{\rho \sqrt {\varPi _Q}}\left (\pi p-\frac {2\mu }{3\mu _0 Re}\frac {\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{\nabla }\rho }{\rho }\right)\!, \end{align}
\begin{align} {\nabla} ^2\varPhi _i&=-\frac {\boldsymbol{\nabla }\varPhi _i\boldsymbol{\cdot }\boldsymbol{\nabla }K}{K}, \end{align}
\begin{align} {\nabla} ^2\varPhi _o&=0, \end{align}
\begin{align} \frac {\rm d}{{\rm d}x}\left (R\sigma v_s\right )&=\frac {2K}{K_0}RE_n^i\sqrt {1+{R'}^2}, \end{align}
\begin{align} E_n^o-\varepsilon E_n^i&=\frac {\pi \sqrt {\varPi _Q}}{2}\sigma , \end{align}
\begin{align} \boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{n}&=0, \end{align}
\begin{align} T\boldsymbol{\cdot }\boldsymbol{n}&=0, \end{align}
\begin{align} \frac {\mu }{\mu _0}\boldsymbol{t}\boldsymbol{\cdot }\big(\boldsymbol{\nabla }\boldsymbol{v}+{\boldsymbol{\nabla }\boldsymbol{v}}^T\big)\boldsymbol{n}&=\frac { Re}{2}E_t\sigma -\frac {\pi Re}{\gamma _0}\frac {d\gamma }{dt}, \end{align}
\begin{align} \frac {\gamma \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{n}}{\gamma _0} +\frac {2\mu }{\pi \mu _0 Re}\left (\boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{v}\boldsymbol{n} +\frac {\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{\nabla }\rho }{3\rho }\right ) &=p+\frac {{E_n^o}^2-\varepsilon {E_n^i}^2-{E_t}^2\left (\varepsilon -1\right )}{2\pi ^2\sqrt {\varPi _Q}}. \end{align}
These equations include new terms proportional to the gradient of physical properties (e.g.
$-({(\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{\nabla }\rho)}/{\rho})$
), or to a temperature-dependent physical property (e.g.
$({\rho }{\rho _0})({\sqrt {\varPi _Q} Re}/{\pi})\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{v}$
). In order to evaluate the importance of the terms including the gradient of a physical property we solve the original equations (2.1a
)–(2.1e
) and (2.3a
)–(2.3f
) using constant values of the electrical conductivity and viscosity
\begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{v}&=0, \end{align}
\begin{align} \frac {\sqrt {\varPi _Q} Re}{\pi }\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{v}&={\nabla} ^2\boldsymbol{v}-\pi Re\boldsymbol{\nabla }p, \end{align}
\begin{align} \boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{\nabla }T-\frac {{\nabla} ^2T}{\textit{Pe}}&= \frac {\left (\boldsymbol{\nabla }\boldsymbol{v}+{\boldsymbol{\nabla }\boldsymbol{v}}^T\right )\boldsymbol{:}\boldsymbol{\nabla }\boldsymbol{v}}{ Re\sqrt {\varPi _Q}}+\frac {{\boldsymbol{\nabla }\varPhi _i}\boldsymbol{\cdot }{\boldsymbol{\nabla }\varPhi _i}}{\sqrt {\varPi _Q}}, \end{align}
\begin{align} {\nabla} ^2\varPhi _i&=0, \end{align}
\begin{align} {\nabla} ^2\varPhi _o&=0, \end{align}
\begin{align} \frac {\rm d}{{\rm d}x}\left (R\sigma v_s\right )&=2RE_n^i\sqrt {1+{R'}^2}, \end{align}
\begin{align} E_n^o-\varepsilon E_n^i&=\frac {\pi \sqrt {\varPi _Q}}{2}\sigma , \end{align}
\begin{align} \boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{n}&=0, \end{align}
\begin{align} T\boldsymbol{\cdot }\boldsymbol{n}&=0, \end{align}
\begin{align} \boldsymbol{t}\boldsymbol{\cdot }\big (\boldsymbol{\nabla }\boldsymbol{v}+{\boldsymbol{\nabla }\boldsymbol{v}}^T\big )\boldsymbol{n}&=\frac { Re}{2}E_t\sigma , \end{align}
\begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{n} +\frac {2\boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{v}\boldsymbol{n}}{\pi Re} &=p+\frac {{E_n^o}^2-\varepsilon {E_n^i}^2-{E_t}^2\left (\varepsilon -1\right )}{2\pi ^2\sqrt {\varPi _Q}}, \\[9pt] \nonumber \end{align}
evaluate each new term proportional to the gradient of a physical property with the resulting temperature field, and compare its magnitude with the dominant term
$S$
in the constant-property equation, estimated as one half of the sum of the magnitudes of all terms
\begin{align} S_{\textit{mass}}&=\frac {1}{2}\left (\left |\frac {\partial v_x}{\partial x}\right |+\left |\frac {\partial v_r}{\partial r}\right |+\left |\frac {v_r}{r}\right |\right)\!, \end{align}
\begin{align} S_{\textit{momentum}}&=\frac {1}{2}\left (\left |\left |\frac {\sqrt {\varPi _Q} Re}{\pi }\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{v}\right |\right |+\big |\big |{\nabla} ^2\boldsymbol{v}\big |\big |+\left |\left |\pi Re\boldsymbol{\nabla }p\right |\right |\right)\!, \end{align}
\begin{align} S_{\textit{energy}}&=\frac {1}{2}\left (\left |\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{\nabla }T\right |+\left |\frac {{\nabla} ^2T}{\textit{Pe}}\right |+\left |\frac {\big (\boldsymbol{\nabla }\boldsymbol{v}+{\boldsymbol{\nabla }\boldsymbol{v}}^T\big )\boldsymbol{:}\boldsymbol{\nabla }\boldsymbol{v}}{ Re\sqrt {\varPi _Q}}\right |+\left |\frac {{\boldsymbol{\nabla }\varPhi _i}\boldsymbol{\cdot }{\boldsymbol{\nabla }\varPhi _i}}{\sqrt {\varPi _Q}}\right |\right)\!, \end{align}
\begin{align} S_{\textit{Laplace}}&=\frac {1}{2}\left (\left |\frac {\partial ^2\varPhi _i}{\partial x^2}\right |+\left |\frac {\partial ^2\varPhi _i}{\partial r^2}\right |+\left |\frac {1}{r}\frac {\partial \varPhi _i}{\partial r}\right |\right)\!, \end{align}
\begin{align} S_{\textit{normal}}&=\frac {1}{2}\left (\left |\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{n}\right |+\left |\frac {2\boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{v}\boldsymbol{n}}{\pi Re}\right |+\left |p\right |+\left |\frac {{E_n^o}^2-\varepsilon {E_n^i}^2-{E_t}^2\left (\varepsilon -1\right )}{2\pi ^2\sqrt {\varPi _Q}}\right |\right)\!, \end{align}
\begin{align} S_{\textit{tangent}}&=\frac {1}{2}\left (\left |\boldsymbol{t}\boldsymbol{\cdot }\big (\boldsymbol{\nabla }\boldsymbol{v}+{\boldsymbol{\nabla }\boldsymbol{v}}^T\big )\boldsymbol{n}\right |+\left |\frac { Re}{2}E_t\sigma \right |\right)\!. \end{align}
On the other hand, the importance of including the dependence on temperature in a term proportional to a physical property
$\beta$
is gauged by the function
$|({\beta (T(x)) - \beta _0})/{\beta _0}|$
.
Figure 16(a–f) shows the importance of the terms including a property gradient in the equations of conservation of mass, momentum and energy, in the equation of conservation of charge in the bulk of the liquid, and in the balances of normal and tangential stresses. In each figure the solid black line is the dominant term in the equation (when
$S$
is a function of both the radial and axial coordinate we plot its maximum value at a given axial position), while the dashed black line represents
$S/100$
. For plotting purposes, we use the position
$x_0$
of the current cross-over point as the origin of the axial coordinate. Only the viscous
$({(\boldsymbol{\nabla }\boldsymbol{v}+{\boldsymbol{\nabla}\boldsymbol{v}}^T )\boldsymbol{\cdot }\boldsymbol{\nabla }\mu })/({\mu _0})$
term in the equation of conservation of momentum, the
$-({\boldsymbol{\nabla }\varPhi _i\boldsymbol{\cdot}\boldsymbol{\nabla }K})/({K})$
term in the equation of conservation of charge in the bulk and the heat capacitance
$({\rho T\boldsymbol{v}\boldsymbol{\cdot }{\boldsymbol{\nabla}c}})/({\rho _0 c_0})$
term in the equation of conservation of energy are significant. The
$-({\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{\nabla}\rho})/({\rho })$
term in the equation of conservation of mass is relatively significant, but given that we use the streamfunction–vorticity method to solve the equations of conservation of mass and momentum it is preferable to assume a constant density in the continuity equation.
Figure 17 shows the relative variation of the liquid properties. Terms proportional to the viscosity and the electrical conductivity must retain temperature dependency. The numerical solution of cone jets with high temperature increases would be improved by including a temperature-dependent surface tension in the balance of normal stresses, a temperature-dependent heat capacitance in the equation of conservation of energy and probably a temperature-dependent density in the equation of conservation of momentum. Including a temperature-dependent thermal conductivity is irrelevant due to the large Peclet number, which makes thermal conduction negligible compared with convection.
Based on the previous analysis, we next solve an upgraded model that includes temperature dependent viscosity, electrical conductivity, heat capacity and surface tension. When using a variable heat capacity the one-dimensional approximation (2.4d ) for the temperature along the jet needs to be modified
\begin{equation} T_{\textit{jet}}=f_0+\frac {c_0}{c\sqrt {\varPi _Q}}\int _{x_{23}}^x{\left (\frac {12\mu R'^2}{{Re} R^4}+KR^2{E_t}^2\right )}{\rm d}x+\sum _{n=1}^\infty {g_nJ_0 (\eta j_{1,n})e^{\tfrac {c_0{j_{1,n}}^2(x_{23}-x)}{cPe}}}. \end{equation}
All other jet approximations remain unchanged, since the two-dimensional equations from which they are derived do not include new terms. The same applies to
$f_0$
and
$g_n$
. Figure 18 shows the total current as a function of the dimensionless flow rate, comparing the new solutions (
$K(T)$
,
$\mu (T)$
,
$c(T)$
and
$\gamma (T)$
) with the original dissipative model (
$K(T)$
and
$\mu (T)$
only) and the isothermal solution. A temperature-dependent heat capacity is the main cause of the relatively small difference between the two dissipative solutions. Since the heat capacity of liquids typically increases with temperature, the liquid stores more thermal energy for a given temperature increase, the temperature decreases compared with a model with constant heat capacity, and so does the total current. Note that EAN, for which we lack enough heat capacity data to obtain its dependence on temperature, exhibits the smallest current difference, thus highlighting the importance of considering a temperature-dependent specific heat capacity (versus surface tension).
Figure 18. Effect on the total current of including temperature-dependent heat capacity and surface tension, in addition to temperature-dependent viscosity and electrical conductivity (the only temperature-dependent properties in the original dissipative model).
Appendix C. Ion Field Emission from the surface of steady cone-jets
As discussed in § 2, although cations constitute a significant fraction of the experimental beam current in the simulated conditions, they are emitted from droplets in the jet’s breakup region and therefore ion field emission from the steady cone jet can be neglected for the calculations in the present article. We do expect ion emission from the steady cone jet, but only at flow rates lower than simulated. We next use an extension of the present model to investigate this regime.
Ion evaporation is often modelled with the kinetic law proposed by Iribarne & Thomson (Reference Iribarne and Thomson1976). This equation relates the ion current density with the energy barrier
${\Delta G}_o-G_E$
that emitted ions must overcome
\begin{equation} \tilde {J}_\omega =\frac {k_B\tilde {T}}{h}\tilde {\sigma }\exp \left (-\frac {{\Delta G}_o-G_E}{k_B\tilde {T}}\right )\! ,\end{equation}
where
$k_B$
and
$h$
are the Boltzmann and the Planck constants,
${\Delta G}_o$
is the ion solvation energy and
$G_E$
is the electrostatic reduction of the energy barrier. To compute
$G_E$
we assume that the emitting medium is a planar dielectric (Magnani & Gamero-Castaño Reference Magnani and Gamero-Castaño2022)
\begin{equation} G_E=\left (\frac {q^{3}E_{n}^{o}\left (\varepsilon -1\right )}{4\pi \varepsilon _{0}\left (\varepsilon +1\right )}\right )^{1/2} \!, \end{equation}
where
$q$
is the net charge of the ion and
$\varepsilon$
is the dielectric constant of the liquid. We are not aware of a direct measurement of the ion solvation energy for any ionic liquid. Measuring
${\Delta G}_o$
is difficult, as it requires a precise determination not only of the ion flux, but also of the electric field and the temperature of a steady surface with a radius of curvature of tens of nanometres or smaller (these small radii of curvature are needed to sustain the high electric fields associated with significant ion emission). Ion solvation energies for a variety of ion–solvent combinations have been determined indirectly, by investigating the residues left by droplets evaporating both ions and solvent molecules at atmospheric conditions. For example, Loscertales & Fernández de la Mora (Reference Loscertales and Fernández de la Mora1995) report values between 2.05 and 2.22 eV for Li
$^+$
-based cations evaporated from ethylene glycol and formamide droplets, while Gamero-Castaño & Fernández de la Mora (2000) report a value of 1.85 eV for (C
$_7$
H
$_15$
)
$_4$
N
$^+$
evaporated from formamide. In addition,
${\Delta G}_o$
can be estimated with Born’s continuum model for emission from a planar dielectric
\begin{equation} \Delta G_{o,{B}} = \frac {3}{4} \left (\frac {\gamma q^4 (1-1/\varepsilon )^2}{\pi \varepsilon _0^2 }\right )^{1/3}\!. \end{equation}
Here,
$\gamma$
is the surface energy of the liquid. Table 6 lists
$\Delta G_{o,{B}}$
for the ionic liquids investigated. We compute values at two temperatures typical of these cone jets (dissipation significantly increases the temperature of the liquid along the jet, see figure 6). The temperature dependence of the surface tension is computed using the fittings in figures 12–15, while
$\varepsilon$
is regarded constant due to the lack of temperature-dependent data. The Born model yields ion solvation energies between 2.33 and 1.85 eV, with values decreasing at increasing temperature due to the temperature dependence of the surface tension.
Table 6. Surface tension and Born-estimated ion solvation energy for EMI-Im, BMI-TCM, EMI-TFA and EAN at 294 K and 494 K.
Figure 19. Ion current emitted up to a given axial position for EMI-Im,
$\varPi _Q = 1000$
, and different values of the ion solvation energy: (a) as a function of the axial coordinate; (b) as a function of the voltage drop from the emitter. The ion current is computed with the numerical solution a posteriori, i.e. the model equations do not include ion emission.
Figure 19(a) shows the ion current emitted from the surface of the cone jet up to a given axial position
\begin{equation} \tilde {I}_i(\tilde {x}) = \int _{-\infty }^{\tilde {x}} 2\pi \tilde {R} \sqrt {1+{\tilde {R}{'}}^2} \tilde {J}_{\omega } {\rm d}{\tilde {x}} ,\end{equation}
for EMI-Im and a dimensionless flow rate of 1000. In these calculations the model equations do not include the ion emission law (C1), instead
$\tilde {I}_i(\tilde {x})$
is estimated a posteriori using the electric field, temperature and radius of the cone jet from a numerical solution in which the viscosity and the electrical conductivity are the only physical properties that vary with temperature. The total current (without ion emission) yielded by the numerical solution is 251 nA. We consider several values of
$\Delta G_o$
between 1.3 and 1.8 eV. It is known from experiments that the ion current is near 18 % of the total current,
$I_i/I \approx 0.18$
, for this flow rate (Gamero-Castaño & Cisquella-Serra Reference Gamero-Castaño and Cisquella-Serra2021; Caballero-Pérez et al. Reference Caballero-Pérez, Galobardes-Esteban and Gamero-Castaño2025). This simple estimate conveys two important points if the ions were to be emitted from the surface of the steady cone jet: an unphysically low ion solvation energy, between 1.3 and 1.4 eV, would be needed to reproduce the experimental ion current; and ion emission would be concentrated in a section of the jet near its base, unlike in experiments where ions are emitted further downstream at a lower potential, more specifically from the jet’s breakup region (Gamero-Castaño & Cisquella-Serra Reference Gamero-Castaño and Cisquella-Serra2021; Perez-Lorenzo & Fernandez de la Mora Reference Perez-Lorenzo and Fernandez de la Mora2022). Figure 19(b) better illustrates this point by plotting the ion current as a function of the voltage drop from the emitter,
$\tilde {\phi }_e-\tilde {\phi }(\tilde {x})$
. For example, half of the ion current is emitted by the location where the voltage drop is 118 V for
$\Delta G_o = 1.3$
eV, and 121 V for
$\Delta G_o = 1.4$
eV. Meanwhile, retarding potential measurements of these ions (Gamero-Castaño & Cisquella-Serra Reference Gamero-Castaño and Cisquella-Serra2021) indicate that the emission is narrowly distributed around the voltage of the breakup region, which is 322 V lower than the voltage of the emitter (see figure 8
a). In summary, these numerical results indicate that, at the flow rates simulated, the ion current observed in experiments does not originate from the surface of the steady jet.
Figure 20. Ion currents resulting from a model that includes ion field emission, equation (C5), for EMI-Im cone jets and several values of
$\Delta G_o$
. Open symbols correspond to solutions in which only
$K$
and
$\mu$
are temperature-dependent, while solid symbols include the temperature dependencies of
$K$
,
$\mu$
,
$\gamma$
and
$c$
: (a) ratio between ion current and total current; (b) total current.
Next, we compute ion field emission self-consistently, by including the emission term in the equation of conservation of charge
\begin{equation} \frac {\rm d}{{\rm d}x}\left (R\sigma v_s\right )=2R\sqrt {1+{R'}^2} \left (\frac {K}{K_0}E_n^i - J_\omega \right)\!. \end{equation}
For simplicity we neglect the associated emission of mass, i.e. the boundary condition (2.3c
) remains unchanged. Figure 20(a) shows the ratio between the ion and total currents for EMI-Im cone jets and three values of the ion solvation energy. To gauge the importance of including the dependence on temperature of the surface tension and the specific heat, in addition to the electrical conductivity and viscosity (see Appendix B), we solve the model for both cases. The current ratio in experiments is approximately constant,
$I_i/I \approx 0.18$
, for flow rates
$\varPi _Q \gtrsim 290$
, and increases monotonically at lower flow rate (Gamero-Castaño & Cisquella-Serra Reference Gamero-Castaño and Cisquella-Serra2021; Caballero-Pérez et al. Reference Caballero-Pérez, Galobardes-Esteban and Gamero-Castaño2025). Similarly to the a posteriori estimate, the self-consistent numerical solution also requires an unphysically low ion solvation energy,
$\Delta G_o \approx 1.3$
eV, to produce the experimental ion currents. Furthermore
$I_i/I$
is not constant for
$\varPi _Q \gtrsim 290$
. Note also that the intensity of ion emission decreases significantly when retaining the temperature dependence of the specific heat and the surface tension, since this more accurate formulation yields lower temperatures. Finally, figure 20(b), which compares the total current with experimental values, shows that at sufficiently low flow rate ion emission from the base of the jet changes the monotonic trend between current and flow rate, making the total current increase at decreasing flow rate. This trend reversal is amplified when using constant specific heat and surface tension, which supports the need to use the more accurate formulation (i.e.
$K(T)$
,
$\mu (T)$
,
$\gamma (T)$
and
$c(T)$
) when ion field emission is significant.
Figure 21. Numerical solution for an EMI-Im cone jet,
$\varPi _Q = 80$
and
$\Delta G_o = 1.5$
eV. Dashed lines are for a solution that does not include ion field emission in the model equations (ion current is computed a posteriori), while solid lines are for a self-consistent solution including
$K(T), \mu (T), \gamma (T),c(T)$
.
Although ions are mostly emitted from the jet’s breakup region in the flow rates simulated, we do expect that ions will also be emitted from the surface of the steady jet at sufficiently low flow rates and it is worth investigating how ion emission will modify the cone jet. Figure 21 shows the effect of significant ion emission on the temperature and
$E_n^o$
profiles for EMI-Im,
$\varPi _Q = 80$
and
$\Delta G_o = 1.5$
eV. Under these conditions the current of ions emitted from the steady jet is 13.4 nA, representing 11.8 % of the total current. We compare these profiles with those obtained without the inclusion of ion emission in the model. The ion current estimated a posteriori is significantly larger than the ion current computed self-consistently. This is due to the lower normal component of the electric field in the self-consistent solution, a negative feedback mechanism that may help stabilising the cone jet. On the other hand, a larger conduction current develops to support ion emission, leading to higher ohmic dissipation in the emission region and an increase of the temperature. However, the larger denominator inside the exponential of (C1) associated with the higher temperature does not compensate for the increase of the energy barrier due to the lower
$E_n^o$
, and the net effect is a reduction of ion emission compared with the a posteriori estimate.
Finally, a corollary of this preliminary study is that the extension of the present model combined with accurate measurements of ion emission from cone jets of ionic liquids operating at very low flow rates will make it possible to quantify their ion solvation energies, and study ion field emission more generally. It is apparent that including dissipation effects will be required to tackle this study.
Open access
