Analysis of self-heating in electrosprays operating in the cone-jet mode

Abstract The electrohydrodynamic processes taking place in a cone jet cause ohmic and viscous dissipation, and ultimately self-heating of the liquid. Despite this, previous analyses have modelled cone jets as isothermal systems. To investigate the validity of this assumption, this work applies the leaky-dielectric model to cone jets, while also requiring conservation of energy to reproduce the variation of temperature caused by dissipation and temperature-dependent liquid properties. The main goals are to determine whether there exist electrospraying conditions for which the isothermal assumption is inaccurate, and quantify the temperature field under such conditions. The work confirms that self-heating and thermal effects are important in liquids with sufficiently high conductivities, which is a significant limit because these electrical conductivities are needed to produce jets and droplets with radii of tens of nanometres or smaller. The numerical solution provides accurate expressions for evaluating the dissipation and the temperature increase in cone jets, and confirms that thermal effects cause the apparent breakdown of the traditional scaling law for the current of cone jets of highly conducting liquids.


Introduction
Electrospraying is an atomization technique based on the use of electrostatic stresses to break a liquid into charged droplets.Electrosprays can be operated in several regimes (Cloupeau & Prunet-Foch 1989), among which the cone-jet mode has received significant attention due to its ability to produce droplets with narrow and controllable size distributions (Cloupeau & Prunet-Foch 1990;De La Mora & Loscertales 1994;Rosell-Llompart & De La Mora 1994;Chen, Pui & Kaufman 1995;Gañán-Calvo, Davila & Barrero 1997).A cone jet features a conical meniscus (Taylor 1964) with a long and steady jet emerging from its tip.The inherent Rayleigh instability of the jet leads to the formation of charged droplets.The electrospray current I and the average diameter of the droplets D depend on various physical properties of the liquid (surface tension γ , electrical conductivity K, density ρ, viscosity μ and dielectric constant ε), as well as the flow rate Q.These dependencies can be approximated by scaling laws derived from approximate balances and experimental data (De La Mora & Loscertales 1994;Gañán-Calvo et al. 2018), (1.1) where α and β are coefficients of order one, and ε 0 is the vacuum permittivity.The electrical conductivity is the only physical property in these scaling laws that can vary by many orders of magnitude (its value is typically fixed by dissolving the required amount of a salt), and therefore, this property plays an essential role in electrospraying.Since the minimum flow rate at which a cone jet can operate approximately scales with K −2/3 (Gañán-Calvo, Rebollo-Muñoz & Montanero 2013; Gamero-Castaño & Magnani 2019a), the diameters of the smallest primary droplets and jets scale with K −1/2 .Thus, liquids with high electrical conductivities are used to generate small droplets.For reference, conductivities near 1 S m −1 are needed to produce droplets with diameters of 10-30 nanometres.
Numerical solutions of first-principles models describe in detail the physics of cone jets.Although these models and other studies of electrospraying phenomena assume isothermal conditions (Taylor 1966;Guerrero et al. 2007;Collins et al. 2008;Herrada et al. 2012;Collins et al. 2013;Gañán-Calvo et al. 2018;Ponce-Torres et al. 2018;Gamero-Castaño & Magnani 2019b), energy analysis of droplets electrosprayed in vacuum indicates that a significant portion of the work done by the electric field on the fluid is dissipated (Gamero-Castaño 2010).Furthermore, Gamero-Castaño (2019), by extrapolating the solution of an isothermal calculation, has proposed that the temperature increase along the cone jet caused by dissipation is approximately given by where c is the specific heat capacity.This result indicates that typical liquids such as tributyl phosphate, propylene carbonate, ethylene glycol or formamide, experience temperature increases of a few degrees Celsius at conductivities near 0.05 S m −1 , while more substantial increases are expected at the electrical conductivities required for the generation of nanodroplets.A significant temperature increase impacts the operation of cone jets in multiple ways: e.g. by modifying the values of physical properties (especially the electrical conductivity and the viscosity), by increasing ion emission from the surface of the cone jet (Gallud & Lozano 2022;Magnani & Gamero-Castaño 2023) and by enhancing liquid evaporation.It also leads to the misinterpretation of experimental data based on isothermal scaling laws (Gamero-Castaño 2019; Perez-Lorenzo 2022).Consequently, it is important to consider energy dissipation and the associated self-heating to accurately capture the behaviour of cone jets, especially when operating in the nanometric regime.In this paper we describe an extension of the leaky-dielectric model (Melcher & Taylor 1969;Saville 1997) for cone jets that retains the self-heating of the liquid caused by both ohmic and viscous dissipation.To this effect, the model incorporates the equation of conservation of energy together with temperature-dependent functions for the viscosity and the electrical conductivity, which exhibit exponential behaviours in most liquids (Stokes & Mills 1966;Fogelson & Likhachev 2001;Okoturo & VanderNoot 2004;Leys et al. 2008).The numerical solution yields expressions for the dissipation and the temperature increase in the liquid, and the model is validated with existing measurements of the current and flow rate of solutions of ethylene glycol, propylene carbonate and tributyl phosphate, as well as the ionic liquid 1-ethyl-3-methylimidazolium bis((trifluoromethyl)sulfonyl)imide (EMI-Im).

Model and solving scheme
The model domain, sketched in figure 1, is a spherical region centred on the cone-to-jet transition region.It includes the liquid phase (zones Σ 1 , Σ 2 and Σ 3 ) and the surrounding vacuum (Σ 4 ), separated by the interface R(x).A large domain radius makes the dependency of the solution on this parameter negligible.The Taylor cone potential (Taylor 1964) is imposed as a far-field boundary condition.This amounts to modelling the cone jet supported by a semi-infinite Taylor cone, yielding a universal solution independent of the particular geometry of electrodes present in any experimental configuration.We extend the leaky-dielectric model (Melcher & Taylor 1969;Saville 1997) to include the self-heating of the liquid induced by energy dissipation.The steady-state model solves for the position of the surface R(x), the velocity v, pressure p and temperature T in the liquid (with temperature-dependent viscosity and electrical conductivity), the surface charge σ and the electric potential inside the liquid Φ i and in the surrounding vacuum /(πc) as the independent set of characteristic scales (length, velocity, pressure, current and temperature, respectively) for defining dimensionless variables.The derived scales for the electric field, electric potential, surface charge and power are ) and P c = Φ c I c , respectively.Here K 0 and μ 0 are the electrical conductivity and viscosity at the inlet temperature T 0 .Note that the characteristic length l c is representative of the average diameter of the electrospray droplets and the base of the jet, (1.2).Henceforth, we designate dimensional variables with a tilde while dimensionless variables are uncapped.The model includes the equations of conservation of mass, momentum, energy and charge in the liquid, with temperature-dependent electrical conductivity and viscosity given by where Y K , B K , T K , Y μ , B μ and T μ are liquid-specific constants (see table 2).Equation (2.4) results from the standard leaky-dielectric assumptions of negligible volumetric charge and the use of Ohm's law to express the current density, J = KE, together with the irrotational nature of the electric field (∇ × E = 0 → E = −∇Φ).Note also the inclusion in (2.3) of the viscous and ohmic dissipation densities, (2.6) In the surrounding vacuum the electric potential fulfils Laplace's equation, while on the surface of the cone-jet conservation of charge, the balance of stresses in the tangential and normal directions, the surface kinematic condition and the standard condition for the jump of the electric field in the interface between dielectric media must be fulfilled, (2.12) Here n and t are the unit vectors normal and tangential to the surface, while E n and E t are the normal and tangential components of the electric field.The system of (2.1)-(2.12)contains four dimensionless numbers, namely the dimensionless flow rate Π Q , the Reynolds number Re, the Péclet number Pe and the dielectric constant ε.Table 1 compiles   980 A40-4 Self-heating in electrosprays operating in the cone-jet mode the definitions of all characteristic scales and dimensionless numbers in the model.Table 2 lists the constants in (2.5) needed to evaluate the electrical conductivity and viscosity of the liquids simulated in this paper.The cone jet is divided into three regions (see figure 1) to reduce computational cost: upstream, in region Σ 1 , the radius of the meniscus is large and the problem can be regarded hydrostatic and isothermal (Gamero-Castaño & Magnani 2019b); in the central region Σ 2 the full electrohydrodynamic model is solved using the stream function Ψ and the vorticity Ω formulation to decouple the calculation of the velocity and pressure fields (Higuera 2003;Gamero-Castaño & Magnani 2019b); in region Σ 3 we take advantage of the slenderness of the jet, R 1 and R 1, to solve analytically the differential equations for Ψ , Ω, Φ i and T using Poincaré's perturbation method (Paulsen 2013): a function of interest f (x) is expressed as a series using R(x) as the expansion parameter, f = j f j R j ; this definition is inserted in the differential equation for f ; and each factor f j is computed separately by isolating the terms of order R j in the equation.This method yields (2.14) where η is one of two coordinates in an orthogonal system perpendicular to the surface, and ranging between 0 at the axis and 1 at the surface (Gamero-Castaño & Magnani (2019b), see also Appendix B).Here Φ 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 Appendix A describes in detail the approximated jet solution.
We use the Taylor potential (Taylor 1964) as a far-field boundary condition where K and E are the complete elliptic integrals of the first and second kind, and θ T ≈ 49.29 • is the Taylor angle.This expression is numerically more accurate and easier to compute than the usual formula based on Legendre's polynomials.The boundary condition for the velocity field at the inlet of region Σ 2 is a refined version of the Neumann boundary conditions ∂Ψ/∂n = 0, ∂Ω/∂n = 0: these simple constraints introduce a small jump on the surface stress across Σ 1 and Σ 2 that hinders the convergence of the numerical solution.To eliminate this problem, we use instead the field equations for the stream function and the vorticity simplified for the case R 1, as boundary conditions, where r and r are derivatives along the inlet boundary of Σ 2 .Additional boundary conditions include T = T 0 at the inlet boundary of Σ 2 and zero heat flux at the surface of the cone jet.
In region Σ 2 all differential equations are solved using second-order central differences in an orthogonal grid, described in Appendix B (Srinivas & Fletcher 2002;

Domain discretization
Velocity: (4), ( 5), ( 12), ( 14), ( 16 Telles & Wrobel 2012; Bakr 2013).All dependent variables in region Σ 3 are computed using the analytic solutions (2.13)-(2.15).Figure 2 illustrates the iterative procedure for solving the system of equations.At the start of each iteration the four regions of the domain are discretized.Next, the model equations, separated into groups, are solved sequentially with inputs from partial solutions computed previously: first we compute the temperature field and the viscosity and electrical conductivity by solving (2.3), (2.5) and (2.15); next we compute the electric potential and the surface charge by solving (2.4), (2.7), (2.8), (2.12) and (2.14); next the stream function and the vorticity are obtained by solving (2.1), (2.2), (2.9), (2.11) and (2.13).These three blocks are recomputed until the difference between two consecutive steps is sufficiently small (we use as criterion a residual smaller than 10 −5 ).At this point, the shape of the surface is adjusted by minimizing the residual of the balance of normal stresses (2.10).This process is repeated until the residual of (2.10) is smaller than a desired value, typically 10 −5 .

Numerical solutions and discussion
We have computed numerical solutions for operational states, i.e. for sets of values of the dielectric constant, dimensionless flow rate, Reynolds number and Péclet number, previously characterized in experiments (Gamero-Castaño 2019).In particular we have computed solutions for dielectric constants of 8.25, 37.6 and 65.9, which correspond to the liquids tributyl phosphate, ethylene glycol and propylene carbonate, respectively.Figure 3 illustrates a typical solution, in this case calculated for Π Q = 100, Re = 0.19, Pe = 14.87 and ε = 65.9.The surface R(x) exhibits a sharp transition between the conical region and the jet, as shown by the narrowness of its second derivative R (x); we use the position of the maximum of R , x 0 , as the origin of the axial coordinate to plot the solution.The conduction current through the bulk of the liquid, ) dissipation linear densities, The ohmic dissipation peaks slightly upstream of the current crossover point.Here P Ω is nearly equal to the product I b E t , especially in the slender jet, and therefore, it peaks near the maximum of the tangential electric field since at this point the bulk conduction current is still dominant; P Ω decreases beyond its maximum primarily driven by the conversion of conduction current into surface current.Viscous dissipation is less intense, it has a narrower distribution and peaks slightly upstream, near x 0 ; here the flow undergoes a significant increase in velocity and a change in direction, aligning with the axis as the liquid accelerates into a jet.The viscous dissipation drops rapidly downstream of its maximum due to the slow variation of the jet's radius and the uniformity of the velocity profile.Figure 3(c) also shows the temperature increase along the axis.The temperature rises rapidly within −5 < x − x 0 < 20 due to concentration of dissipation in this region.Figure 3(d) shows a two-dimensional map of the temperature increase.Because most of the dissipation occurs in the slender jet, the dissipation densities are nearly uniform along the radial coordinate.This, combined with the adiabatic boundary condition on the surface, yield a temperature field that is nearly constant in the radial direction.Furthermore, axial diffusion is small due to the large value of the Péclet number, which is always larger than 10. Figure 4(a) shows the total ohmic dissipation, as a function of the dimensionless flow rate and the Reynolds number, and for three different values of the dielectric constant.The integral is evaluated using the axial coordinate at which the surface current is 95 % of the total current as the upper limit of integration.The ohmic dissipation increases with the dimensionless flow rate and the dielectric constant and decreases with increasing Reynolds number.The dimensionless flow rate is the stronger dependency.Figure 4(b) shows a fit to a power law that only retains the flow rate and Reynolds dependencies (we have simulated too few dielectric constants).The total ohmic dissipation is well approximated by Re 0.11 . (3.4) The exponent acting on Π Q is similar to that found by Gamero-Castaño (2019) in an isothermal calculation, and can be justified with a one-dimensional estimation of the ohmic dissipation, Since the geometry of the transition region remains virtually unchanged at varying Π Q , Re and ε, the ohmic dissipation primarily scales with E t 2 .The maximum value of the tangential electric field, which is located near the current crossover point, is constant, but otherwise E t scales as Π Q 1/4 in most of the region where conduction current becomes surface current (Gamero-Castaño 2019).Thus, the exponent acting on Π Q in the ohmic dissipation law must be smaller than 1/2, but close to this value.Equation (3.4) quantifies the dependency of the ohmic dissipation on the Reynolds number for the first time.The Reynolds number has a much smaller effect on the cone-jet solution than the dimensionless flow rate, this weaker dependency has not been analysed, and therefore, we cannot justify the exponent acting on Re.
Figure 5(a) shows the total viscous dissipation, for different values of the dimensionless flow rate, the Reynolds number and the dielectric constant.The dimensionless viscous dissipation is proportional to the inverse of the Reynolds number and decreases at increasing flow rate.The dependency on the dielectric constant is weak, and it is likely that the stronger effect mentioned by Gamero-Castaño ( 2019) is due to differences in the dimensionless flow rate investigated in this previous study.Figure 5(b) shows the data to be well fitted by The one-dimensional estimation of the viscous dissipation yields where we approximate the velocity by the ratio between flow rate and the cross-section of the cone jet.The viscous dissipation is proportional to the viscosity and, therefore, inversely proportional to the Reynolds number.The dependence on the dimensionless flow rate is better illustrated by integrating (3.8) by parts, which yields The total viscous dissipation is mostly given by the integral of μR /R 3 ; the integral of μ R /R 3 is a small and negative contribution (zero at constant temperature).As shown in  figure 3(a), R is a sharply peaked function only significant in the cone-to-jet transition.Thus, the shape of R together with (3.9) explain the narrow distribution of the linear viscous dissipation in figure 3(c).As the flow rate increases, the shape of the transition from cone to jet becomes less sharp, i.e. the second derivative R becomes smaller, and the viscous dissipation decreases at increasing flow rate as (3.7) reflects.
The temperature increase along the cone jet is approximated by where (3.4) and (3.7) are used to express the ohmic and viscous dissipations.Figure 6(a) shows a contour plot of (3.10).Solid circles indicate simulated states and provide a sense of the range of dimensionless flow rates and Reynolds numbers investigated, as well as a qualitative validation of (3.10).Figure 6(a) also includes two curves: the locus of states in which ohmic and viscous dissipations are equal, Π Q = (1.65Re 0.11 + 0.58Re −0.89 ) 1.62 , ohmic dissipation being larger than viscous dissipation to the right of this curve; and the minimum flow rate at which a cone jet is stable, which is approximately given by Π Q = 1/Re for liquids with moderate and large electrical conductivities such as those considered in this study (Gañán-Calvo et al. 2013;Gamero-Castaño & Magnani 2019a).Cone jets are stable to the right of this curve.It is apparent that ohmic dissipation is dominant in The analysis of the dimensional temperature is also important, for example, to determine when the thermal and fluid-dynamic problems are coupled.For a given liquid, i.e. for fixed physical properties, (3.10) indicates that the temperature increase augments with decreasing flow rate.Thus, the maximum temperature increase occurs at the minimum flow rate and, using Π Q,min = 1/Re, it is given by Tmax ≈ (7.74 + 12.8Re 0.73 + 4.53Re −0.27 ) Here Tmax is solely a function of the physical properties of the liquid and, since the electrical conductivity is the only property with a range of several orders of magnitude, a corollary of (3.11) is that only liquids with sufficiently high electrical conductivity undergo significant self-heating.Table 3 illustrates this by listing the estimated maximum temperature increase of ten liquids together with their physical properties.The table also includes the characteristic length of the cone jet at the minimum flow rate to correlate self-heating with the radii of the jet and droplets.The temperature increase for the propylene carbonate solutions is significant for K = 0.022 S m −1 , and reaches 27.3 • C for K = 0.176 S m −1 ; the characteristic length associated with the latter is 10.2 nm.All ionic liquids, with conductivities near 1 S m −1 , exhibit maximum temperature increases near or exceeding 100 • C, and characteristic lengths near 10 nm.Although the criterion Π Q,min = 1/Re slightly underestimates the minimum flow rates of propylene carbonate solutions (Gamero-Castaño & Magnani 2019a) and EMI-Im (Gamero-Castaño & Cisquella-Serra 2021), an electrical conductivity near or above 0.01 S m −1 can be used as a rule of thumb for expecting significant self-heating.For example, the maximum temperature increases of propylene carbonate and ethylene glycol solutions with K = 0.01 S m −1 are 4.4 • C and 3.0 • C, respectively.Figure 7 illustrates the importance of self-heating in the cone jet of a liquid with a high electrical conductivity such as EMI-Im.The dimensionless flow rate in this calculation has a value of 400, which is slightly smaller than the minimum flow rate reported by Gamero-Castaño & Cisquella-Serra (2021).The dielectric constant of EMI-Im is 12.8.Figures 7(a) and 7(b) show profiles of the currents and of the dissipation linear densities computed with both the present model and ignoring self-heating (i.e.isothermal model, we assume a constant temperature).The total current obtained with the isothermal model is 131.0 nA.When considering self-heating, the current is 62.4 % higher, 212.7 nA.Furthermore, the latter compares well with the value of 203 nA for Π Q = 400 from the fitting of experimental data reported by Gamero-Castaño & Cisquella-Serra (2021).When self-heating is accounted for, the ohmic and viscous dissipations are respectively much larger and much smaller than in the isothermal calculation.This is an obvious effect of the temperature increase along the cone jet and the associated enhancement of the electrical conductivity and the reduction of the viscosity, together with the proportionality between ohmic dissipation and electrical conductivity, (3.5), and between viscous dissipation and viscosity, (3.8). Figure 7(c) shows the variation of the temperature, electrical conductivity and viscosity along the axis computed with the present model, and the ratio R/R 0 between the radii of the surfaces for the present and the isothermal models.The total temperature increase along this section of the cone jet is 93.3 • C, the electrical conductivity increases by a factor of 5.7 and the viscosity is reduced by a factor of 7.6.The very significant changes in these two key physical properties show that an isothermal description of these cone jets would be grossly incorrect.The strong variation of the electrical conductivity and the scaling law (1.1)readily explain the enhanced value of the total current with respect to the isothermal solution, 212.7 nA vs 131.0 nA.The ratio R/R 0 is 1 upstream, increases up to a maximum value of 1.087 at x − x 0 = 45.6 and decreases downstream to a value of 0.9.The larger current and the smaller jet radius in the presence of self-heating follow the electrical conductivity trend in the scaling laws (1.1) and (1.2).The dependencies are quantitatively weaker, but this is to be expected from the gradual increase of the electrical Liquid  conductivity along the transition region.Finally, the contour plot in figure 7(d) shows the lack of radial variation of the temperature in the slender jet.Although not explored in this paper, the ability to compute the temperature field along the cone jet is key to study ion-field evaporation, an emission mechanism that strongly depends on temperature and that plays a main role in the electrosprays of highly conducting liquids, including EMI-Im (Gamero-Castaño & Cisquella-Serra 2021).The significant increase of temperature in liquids with high electrical conductivity explains their abnormal current versus flow rate behaviour.Figure 8 compares experimental values (Gamero-Castaño 2019) with the numerical solution, for electrosprays of propylene carbonate and ethylene glycol with several electrical conductivities (or, equivalently, Reynolds numbers).The experimental data for each conductivity exhibits the well-known Ĩ ∝ Q1/2 law.When plotted in dimensionless form, Ĩ/ ε 0 γ 2 /ρ vs Π Q , all points for a given liquid are expected to lie on a single straight line regardless of the electrical conductivity.Although the solutions with the largest Reynolds numbers behave as expected, the data for the smaller Reynolds numbers have a vertical offset that increases with decreasing Re.Gamero-Castaño (2019) explained this anomaly as an effect of energy dissipation and the consequent self-heating of the liquid, expected to increase at decreasing Re: when the temperature in the transition region increases significantly, using the emitter temperature to evaluate Π Q underestimates its value, shifting leftward the experimental data (Π Q is proportional to the electrical conductivity, which increases with temperature).The numerical results in figure 8 confirm this explanation.The numerical solution matches well the experimental data, reproducing the vertical offset of the current versus flow rate curves at low Re.The vertical offset starts to be significant at a Reynolds number of approximately 0.3, for which the computed temperature increase is approximately 2 • C. The liquid solutions used in the experiments were obtained by adding increasing quantities of EMI-Im to propylene carbonate and ethylene glycol, in order to increase the electrical conductivity.The uncertainty in the values of the physical properties of these mixtures, which increases at decreasing Re, are likely responsible for the differences between numerical and experimental results.Note that the comparison between the numerical solution and the experimental value in the case of EMI-Im is excellent, 212.7 nA vs 203 nA for Π Q = 400, even though self-heating is significantly more intense in EMI-Im than in any state shown in figure 8. Unlike the propylene carbonate and ethylene glycol solutions, EMI-Im is a pure liquid with accurate K(T) and μ(T) functions (Kazakov et al. 2022).

Conclusions
Ohmic and viscous dissipations in cone jets increase the temperature of the liquid, linking electrohydrodynamic and thermal phenomena.In this paper we self-consistently study this coupling for the first time, using an extension of the leaky-dielectric model that includes conservation of energy and temperature-dependent functions for the viscosity and the electrical conductivity.We find that the temperature increase is a function of the physical properties of the liquid and its flow rate, (3.10); for a given liquid, the maximum temperature increase only depends on the physical properties, (3.11), with the electrical conductivity playing a key role (the maximum temperature scales with K 2/3 ); and temperature increases of a few centigrade degrees occur in liquids with electrical conductivities as low as 0.01 S m −1 .Therefore, when modelling cone jets or analysing experimental data, it is important to account for self-heating when the electrical conductivity is of the order or larger than 0.01 S m −1 , which includes the conductivity values needed to produce droplets and jets with diameters of tens of nanometres or smaller.
Viscous dissipation concentrates in the transition from cone to jet.The generation of ohmic dissipation starts in this region and continues in a longer section of the jet where bulk conduction current transforms into surface current.In the operational range of cone jets where self-heating is significant the ohmic term is the main contributor to the total dissipation, and only near the minimum flow rate, i.e. near the condition Π Q = 1/Re, ohmic and viscous dissipations are comparable.Because of the distribution of dissipation, the temperature increases in the section of the jet where charge is injected from the bulk of the liquid onto the surface.The physical properties of the liquid, especially the viscosity and the electrical conductivity, are functions of the temperature, and their values change along this section where the total current and the diameter of the jet are fixed.Therefore, the traditional scaling laws for cone jets, which assume constant values of the physical properties, become increasingly inaccurate as self-heating intensifies.An example of this effect is the unexpected vertical offset of the dimensionless current versus flow rate function observed in liquids with high electrical conductivities.This work shows that this vertical offset is an artifact of employing the traditional scaling laws in cone jets in which self-heating induces significant temperature variations.expansion parameter: Equations (2.1), (2.2), (2.4) and (2.3) are written in terms of Ψ and Ω and in the orthogonal coordinates ξ , η (see Appendix B) as 980 A40-18 Self-heating in electrosprays operating in the cone-jet mode To apply the perturbation method, we first identify all the terms in these equations that are small in the jet region.Downstream of the liquid surface is slender (R 1, R 1) and the fluid velocity is well aligned with the axial coordinate (v η ∝ ∂Ψ /∂ξ 1).Also, since ξ and x are nearly parallel, R, R , R are considered a function of ξ only.These observations are summarized as the jet hypotheses: By applying (A6) to (A2)-(A5) we eliminate all the higher-order terms and obtain These equations are still too complex to be solved analytically.In order to obtain a solution we substitute Ψ , Ω, Φ i and T for their expansions (A1), and separate terms according to their order R i .
For the vorticity, we obtain the equation for the zeroth-order term.Higher-order terms are still too complex to be solved analytically.The solution of (A11), using (2.9) and the symmetry condition Ω(x, 0) = 0 as boundary conditions, is We apply the vorticity solution (A12)-(A7), to obtain the following equation set for the stream function: In this case we solve analytically for the first two terms of the Ψ expansion.The boundary conditions in this case are ) where the value of Ψ on the surface is obtained from the dimensionless flow rate.The solution of these equations is then For the electric potential, we apply the same procedure.Equation (A9) is expanded as with boundary conditions The final result is Obtaining a solution for the jet temperature is more complicated.In this case the ∂/∂ξ term cannot be eliminated from the zeroth-order equation, but the differential equation can be solved using separation of variables.Since in this case we only solve for T 0 , we use T in the equations for simplicity.The temperature is assumed to be of the form T(ξ, η) = f (ξ ) + g(ξ )h(η). (A19) Leading to the equation The solution for f , g and h is then where J 0 () and Y 0 () are the Bessel functions of order zero of the first and second kind and C is a constant to be found.As boundary conditions, we have From these conditions we obtain h 0 = 1, h 1 = 0 and C = j 1,n ε 1/2 , where j 1,n is the nth zero of the Bessel functions of the first kind of order one.These definitions result in Finally, the constants f 0 and g n are obtained by matching the temperature and its derivative at the base of the jet (x = x 23 ): Appendix B. Coordinate system orthogonal to the surface of the cone jet In region Σ 2 the equation set is solved using finite differences in an orthogonal grid with coordinates ξ , η (see figure 9).Region Σ 2 is bounded by the symmetry axis and by the liquid surface, allowing us to define an orthogonal frame of reference using a variation of where ξ and η are the two orthogonal variables, and x() is a generic function of ξ, η.The coordinate η is equal to 0 at the symmetry axis and 1 at the surface, while ξ goes from 0 on the upstream boundary x 12 to 1 at the downstream boundary x 23 .By applying (B2) into (B1), we obtain We enforce the orthogonality of the ξ, η coordinates by requiring dξ/dη = 0.This leads to ∂x ∂η = −ηRR This equation is integrated from η = 1 to η = 0, using the initial condition the distribution of points on the surface R(x), to obtain the x coordinate of all the points in the orthogonal grid.To obtain the r coordinate, we use the definition of r in (B2) together with the computed x coordinates.
All the equations used in the model are converted from x, r to ξ , η coordinates by using

Figure 1 .
Figure 1.Model domain.The radius of the spherical region is set at one thousand length units to reduce dependencies on the particular placement of far-field boundary conditions.

Figure 2 .
Figure 2. Iterative scheme for solving the model.
dr, is dominant upstream and transforms into convected surface charge, I s = Rσ v s , along the jet.Both transport mechanisms become equal at x − x 0 = 7.30.The tangential and normal components of the electric field on the surface, shown in figure 3(b), exhibit maxima near this current crossover point.Here E i n is always much smaller than E o n , which is indicative of the strong screening of the electric field by the surface charge.Figure 3(c) shows the (a.u.

Figure 3 .
Figure 3. Example of the model solution for Π Q = 100, Re = 0.19, Pe = 14.87 and ε = 65.9:(a) position of the surface, its second derivative, bulk conduction current and surface current; (b) tangential and normal components of the electric field on the surface; (c) ohmic and viscous dissipation linear densities, fluid velocity and temperature increase along the axis; (d) two-dimensional map of the temperature increase.The origin of the axial coordinate is set at the maximum of R (x), x 0 , for display purposes.

εFigure 4 .
Figure 4. Total ohmic dissipation: (a) values as a function of Π Q , Re and ε; (b) fitting of the data to a power law.

Figure 5 .
Figure 5.Total viscous dissipation: (a) values as a function of Π Q , Re and ε; (b) fitting function approximating the data.

Figure 6 .
Figure 6.(a) Contour plot of the estimated dimensionless temperature increase as a function of Π Q and Re, (3.10).The plot displays simulated states as solid circles; (b) error between the estimated and computed temperature increase, ( T est − T)/ T, for propylene carbonate and ethylene glycol solutions.

Figure 7 .
Figure 7. Numerical solution for EMI-Im, Π Q = 400, Re = 8.51 × 10 −3 , Pe = 13.6,ε = 12.8.The red curves correspond to the isothermal solution while the blue curves consider self-heating: (a) bulk and surface current; (b) ohmic and viscous dissipation linear density; (c) temperature increase, electrical conductivity and viscosity along the axis for the self-heating case, and the ratio between the radii of the surfaces for the self-heating (R) and isothermal (R 0 ) cases; (d) 2-D map of the temperature increase in the cone jet (self-heating case).

Figure 8 .
Figure 8.Current versus flow rate, comparison between measurements (solid circles) and model results (open squares and triangles) for ethylene glycol solutions (a) and propylene carbonate solutions (b).

Figure 9 .
Figure 9. Example of the orthogonal grid defined in region Σ 2 .The inset zooms on the base of the jet.

Table 1 .
Characteristic scales and dimensionless numbers used in the model: length l c , velocity v c , pressure p c , current I c , temperature T c , electric field E c , electric potential Φ c , surface charge σ c , power P c , dimensionless flow rate Π Q , Reynolds number Re and Péclet number Pe.The dielectric constant is omitted.

Table 2 .
Coefficients in the temperature-dependent equations (2.5) for the liquids modelled in this paper.The pre-exponential factor of the electrical conductivities of the propylene carbonate (PC), ethylene glycol (EG) and tributyl phosphate (TBP) solutions depend on the solute concentration, and are proportional to Re −3 (we assume an inlet temperature