Self-consistent modeling of laminar electrohydrodynamic plumes from ultrasharp needles in cyclohexane

1 Self-Consistent Modeling of Laminar Electrohydrodynamic Plumes from Ultrasharp Needles in Cyclohexane Marley Becerra,1,a) Henrik Frid,1 Pedro A. Vázquez2 1KTH Royal Institute of Technology, Osquldas väg 6, SE-10044 Stockholm, Sweden. 2University of Seville, Camino de los Descubrimientos s/n, 41092 Sevilla, Spain This paper presents a self-consistent model of electrohydrodynamic (EHD) laminar plumes produced by electron injection from ultra-sharp needle tips in cyclohexane. Since the density of electrons injected into the liquid is well described by the Fowler-Nordheim field emission theory, the injection law is not assumed. Furthermore, the generation of electrons in cyclohexane and their conversion into negative ions is included in the analysis. Detailed steady-state characteristics of EHD plumes under weak injection and space-charge limited injection are studied. It is found that the plume characteristics far from both electrodes and under weak injection can be accurately described with an asymptotic simplified solution proposed by Vazquez et al. Physics of Fluids 12, 2809 (2000) when the correct longitudinal electric field distribution and liquid velocity radial profile are used as input. However, this asymptotic solution deviates from the self-consistently calculated plume parameters under space-charge limited injection since it neglects the radial variations of the electric field produced by a highdensity charged core. In addition, no significant differences in the model estimates of the plume are found when the simulations are obtained either with the Finite Element Method or with a diffusion-free particle method. It is shown that the model also enables the calculation of the current-voltage (IV) characteristic of EHD laminar plumes produced by electron field emission, with good agreement with measured values reported in the literature. Keywords: Electrohydrodynamics (EHD), EHD plumes, liquid high-voltage insulation, self-consistent modeling, Fowler-Nordheim field emission, Hybrid FEM-particle method I. INTRODUCTION The injection of charges in a liquid under high voltage results in an electrohydrodynamic (EHD) flow due to motion of the fluid induced by Coulomb forces.1 Specifically, the injection of charges from a sharp electrode tip results in a charged jet propagating at high velocity, referred to as an EHD plume. Since the resulting liquid velocity is typically larger than the ion drift velocity, liquid EHD plumes are typically convection-dominated with a strong coupling between the velocity field and the charge density (see e.g. Ref 2). The first theoretical efforts to model laminar EHD plumes described the flow in terms of self-similar solutions by analogy with thermal plumes.2-3 Integral methods have later also been used in Ref 46. The initial studies in Ref 26 provided a basic understanding of EHD plumes, with equations describing the flow in the asymptotic a) Author to whom correspondence should be addressed. Electronic mail: [email protected]. 2 region far from the electrodes. These studies have shown that EHD plumes have a double-layer structure with a thin charged region surrounded by a moving layer of liquid.6,7 Analytical solutions describing EHD plumes cannot generally be obtained without significant simplifications due to the complexity of the charge distribution along the plume.8 Instead, numerical simulations of the EHD plume in the complete domain (including the asymptotic region) have been obtained by numerically solving a set of partial differential equations (PDEs).9,10 However, these previous models cannot predict both the electric and the hydrodynamic characteristics of the plume since the fully self-consistent numerical modeling of EHD plumes is an extremely difficult task.6,11 The challenges in the self-consistent modeling of EHD plumes arise due to the strong coupling between the electric field, the charge distribution and the liquid velocity field, in addition to the different length scales and the large gradients involved.6,8,11 Most importantly, the injection of charge in the liquid is difficult to predict since the functional dependence of injection with electric field is generally unknown.12 The stability of EHD flows has been studied in Refs. 8,10,1316. Based on numerical models, it has been found that EHD plumes in the blade-plane configuration have three different flow regimes.10,14 The first regime is laminar and steady while the second and third regimes are periodic and fully turbulent. The critical Reynolds number for which the transition between steady and unsteady regimes occurs has been studied for the blade-plane configuration in Ref. 15. This transition has been evaluated under different combinations of the dimensionless injection level parameter and the dimensionless ionic mobility number.15 Even though the mentioned numerical studies further improved our understanding of EHD plumes, there are still open questions regarding the modeling of the charge injection.12 The charge injection has previously been modeled based on some assumed injection law (see e.g. Refs 10,14,15,17,18). However, since the charge injection is assumed rather than calculated, such models cannot be used for calculating the current-voltage (IV) characteristic of the EHD plume. Furthermore, the mutual coupling between the charge injection and the EHD flow is not generally taken into account, in spite of this being critical to simulate space-charge-limited regimes. In order to overcome the limitations created by the lack of necessary information on actual injection laws, recent studies have attempted to study EHD phenomena in liquids based on processes with well-known theoretical description or with injection functions obtained by iterative fitting to experimental data.12, 19, 20 In such a way, current voltage characteristics are possible to be predicted in agreement with measurements for blunt needle tips.20 This paper presents a self-consistent model to simulate laminar EHD plumes generated by electron injection from ultrasharp negative needlepoints (���0.5��) in cyclohexane. For such sharp needle tips, electrons are injected into the liquid by cold field emission.21,22 Thus, the injected charge is calculated according to the Fowler-Nordheim tunneling theory instead of assuming an injection law. In addition, electron multiplication and attachment processes generating positive holes and negative ions are also considered. Thus, the model solves the coupled system of PDEs governing the EHD plume, which includes the continuity equations of three charged particle types coupled with Poisson and Navier-Stokes equations. This self-consistent model enables the calculation of the current-voltage (IV) characteristic, in good agreement with measurements reported for cyclohexane.22 Furthermore, the calculated characteristics of the hydrodynamic plume are also compared with the asymptotic solution presented in Ref 6. 3 The paper is organized as follows. The mathematical formulation of the model with the governing equations is described in Section II. Section III describes the computer implementation required for obtaining reliable numerical results. The laminar EHD plume resulting from Fowler-Nordheim electron injection is studied under steady-state in Section IV. The EHD flow close to the needlepoint electrode is studied in Section IV.A and the properties of the EHD plume in the asymptotic region far from both electrodes are studied in Section IV.B. The comparison of the estimates obtained using a FEM solver and a hybrid FEM-particle solver is also discussed. Section IV-D presents the calculated IV-characteristic, in good agreement with the experimental results in Ref 22. Finally, conclusions are drawn in Section V. II. MATHEMATICAL FORMULATION The type of EHD plume here considered is produced by non-autonomous injection of electrons from ultra-sharp needle tips (with tip radius �� smaller than 0.5 m) by cold field electron emission. According to the Fowler-Nordheim tunneling theory describing electric field-dependent electron injection23, the electron current density at the needle surface can be expressed as: �� ∗����|��|� �����������������.� |��|�� (1a) where  is the work function of the metal in the liquid, �� is the electric field of the needle tip surface and a is the boundary unit normal vector. The correction factor  is generally defined by: �������1���� �0.166�������� (1b) with ��� ��� ��� |��| �� (1c) where e is the elementary charge and � is the permittivity. The parameters aFN and bFN are the first and second Fowler- Nordheim constants, as specified in Table 1. Once the electrons are injected into the liquid, avalanches are created in cyclohexane due electron multiplication by impact ionization.24 For each ionizing event, an electron and a positive carrier is formed. Measurements indicate that the positive carriers in cyclohexane correspond to fast positive holes, with mobility about one order of magnitude higher than for positive ions.25 In addition, highly mobile electrons can be converted into slower negative ions by attachment processes.26 In order to simulate the dynamics of all the charged species created by these processes, the hydrodynamic approximation is considered. Thus, the drift, diffusion, generation and loss of electrons, positive holes and negative ions in the discharge formed in the liquid is described by a set of continuity equations: ��� �� ��∙������������ ���|�|�������|�|����|��|������ (2)  ��� �� ��∙��������������|�|���|��|������������ (3) 4  ��� �� ��∙������������ �����|�|���|��|������ (4) where n is the number density of electrons, positive holes and negative ions indicated with the subscripts e, p and n respectively. The drift velocity vector W is given by the product of the particle mobility  and the electric field � plus the liquid velocity �. The terms  (E) and  att(E) are the field-dependent ionization and attachment coefficients.  is the recombination coefficient while D is the diffusion coefficient calculated using the Einstein relation.6 The continuity equations (2-4) are coupled with Poissons equation for calculating the electric field E and the electrostatic potential V: �∙������ � � ������������ (5) and with Navier-Stokes equations for an incompressible viscous liquid under a Coulomb body-force: �∙��0 (6) ����� �� ���∙������������������������������ (7) where  L, p and  are the mass density, pressure and dynamic viscosity. Since variations in the electric permittivity due to temperature gradients are neglected, dielectrophoretic and electrostrictive forces (see e.g. Ref 27) are omitted. Since the effects of the gravitational force ��� are much smaller than the effects of the Coulomb force, its effect has been here neglected. The steady state current in the configuration is calculated considering the contributions of electrons, positive holes and negative ions given by:28 � � � ������������������ ��∙���� (8) where ��is the applied voltage and �� is the Laplacian electric field and �� is the volume of the computational domain shown in Figure 1. III. COMPUTATIONAL METHOD The numerical solution of the system of PDEs given by equations (2)-(7) is a challenging task. Due to the presence of large gradients appearing in the distribution of charged particles, the continuity equations are very stiff and therefore difficult to solve due to the strong coupling with Poisson equation. Recently, the finite element method (FEM) has been shown to provide accurate results for resolving those large gradients in the charged particle density when proper meshing and stabilization techniques are used.29 However, it is known that FEM introduces numerical diffusion that artificially widens the thin plume of negative ions created by the EHD motion of the liquid.9 5 Particle methods as described in Refs. 8, 30 are better suited to describe the asymptotic region far from both electrodes where the charged plume is thin. These particle-type methods replace the continuum charge flow assumed by the hydrodynamic approximation by a large number of discrete superparticles. Nevertheless, particle methods have the drawback of introducing artificial numerical oscillations (discrete particle noise) in the simulation of the charged particle drift.31 In order to accurately simulate both the discharge produced by the injection of electrons and the formation of the double layer of the EHD plume, a hybrid FEM-particle numerical model is used. This strategy is implemented in a commercial FEM software32 as the sequential steady-state solution of the FEM and particle transient calculations until convergence is reached. A. FEM model The FEM solution of equations (2)-(7) is performed considering a point-plane configuration using the axisymmetric computational domain illustrated in Figure 1. The radius of the needle electrode tip (���0.2μm� and the distance between the electrodes (d = 0.5 mm) are the same as used for the experiments reported in Ref. 22. The computational domain is bounded by four boundaries referred to as the needle electrode, the ground electrode, the axis of symmetry and the surrounding boundary. The distance from the axis of symmetry and the ground plate to the surrounding open boundary is set to 1 mm. Numerical experiments have shown that the distance to the surrounding open boundary does not affect the solution provided that it is larger than the gap distance. Since effective meshing is key for an accurate and efficient solution with FEM,29 the entire computational domain  3 has been divided into four sub-domains (illustrated in Figure 1) with different mesh settings. The first subdomain  1 corresponds to the region where the mutual coupling between all charged particles and the liquid is evaluated. Since electrons and positive holes are only present in the area close to the needle, equations (2)-(3) are only solved in  1. Outside this domain, the electron and positive hole density is set to zero. The second subdomain  2 defines the narrow region far from both electrodes where the fine structure of the charged plume along the axis of symmetry needs to be resolved. The third subdomain  3 corresponds to the region where negative ions spread over the ground plate. Equations (4) to (7) are solved in all the domains, including the remaining subdomain  4. Refined meshes are used to accurately resolve the high gradients in the charge particle densities under high electric fields such as in the subdomains  1 and  2. Coarser meshes are used in other subdomains. Sequential numerical experiments improving the mesh quality gradually are performed until a stable, mesh-independent solution is achieved as it will be later described in Section IV.C. An extra fine quadrilateral mesh is used in the subdomain  1 in order to resolve the electric discharge in front of the needle. For this, a structured mesh distributed exponentially along the longitudinal edges is used. The mesh element size increases from values of the order of 10-10 m close to the needle surface up to a maximum 10-6 m at the subdomain boundary. A uniform linear mapping with 130 edge elements is used for the transverse edges of this subdomain. Structured quadrilateral meshing is also used for the asymptotic region  2 and the zone close to the plate surface  3. In order to resolve the thin ionic plume, an exponentially distributed mesh is set in the radial direction, with a minimum element size of less than 510-8 m close to the axis of symmetry. About 170.000 mesh elements are used to mesh the domains  1,  2 and  3. Unstructured triangular meshing is used for the remaining domain  4 where the mesh requirements are less demanding, with 110.000 additional elements. 6 FIG 1. Illustration of the computational domains and the boundary conditions used for the simulations. The boundary conditions for the solution of the continuity equations are set as follows. For the solution of (2), Dirichlet boundary conditions are applied to all the boundaries in  1 except at the axis of symmetry. Thus, the needle surface boundary is set as the injection current density (1) divided by the product of the elementary charge and the electron drift velocity. Zerodensity is used for the other two boundaries. Observe that no additional secondary electron emission process is considered at the needle surface. The needle surface is set with a Neumann boundary defining an outflow of positive holes �∙����0, where a is the boundary unit normal vector. Zero-density for positive holes is set in the other boundaries of  1 except at the axis of symmetry. The solution of (4) considers a zero-density boundary at the needle surface, an outflow boundary on the plane surface and a zero-flux on the surrounding boundary. Dirichlet conditions are applied to the boundaries of the needle and the ground plane when solving Poissons equation, where the potential is fixed to the applied voltage �� and zero respectively. The Neumann condition �∙��� � 0is applied to the surrounding boundary. All the boundaries for the solution of Navier-stokes equations are set to the no-slip condition, except for the axis of symmetry. Numerical stabilization of the solver is crucial to solve convection-dominated, stiff equations with FEM (such as in (2)-(4)) or when solving Navier-stokes equations. For this, consistent Streamline Upwind/PetrovGalerkin (SUPG) stabilization is generally used.33 This method consistently introduces numerical diffusion only in the flow direction (e.g along the streamlines) and is less diffusive than isotropic diffusion. Numerical experiments have shown that the SUPG stabilization together with the fine mesh described above provides same results compared with more elaborated stabilization techniques such as the do Carmo and Galeao crosswind method. 7 The simulation starts with initial zero values for the density of all the charged particles and the liquid velocity. Then transient calculations are performed for each applied voltage �� until the simulated current reaches steady state. The numerical values for the material coefficients for cyclohexane used in the calculations are shown in Table I. TABLE I. Parameters used in the simulation ParameterValueSource First Fowler-Nordheim constant aFN1.5410-6 Aev/V^2 Ref.23 Second Fowler- Nordheim constant bFN6.83109 ev- 1.5V/m Ref.23 Ionization coefficient  f1(E)Ref.24 Attachment coefficient ���� f2(E)Ref.24 Electron mobility  e 410-5 m2/Vs Ref.25 Positive hole mobility  p 1.510-8 m2/Vs Ref.25 Negative ion mobility  n 10-8 m2/Vs Ref.25 Recombination coefficient  4.5 10-16 Ref.22 Liquid density  L 780kg/m3Ref.22 Dynamic viscosity  8.810-4 Pas Ref.34 Work function  4.51 eVRef.22 B. Hybrid FEM-particle model A particle method is used to improve the FEM solution of the negative charge density in the asymptotic area  2. Thus, a diffusion-free lagrangian description of the thin charged plume is used. For this, a large number of discrete particles (N=10000) are initially released from the boundary  1. The distribution of particles released is set proportional to the negative ion density from the steady-state FEM solution. In a similar manner, the initial velocity of the released particles is given by the negative ion velocity Wn calculated with FEM at the boundary. In order to track the position xj = (rj, zj) for each j-th particle released, Newtons law of motion is solved with FEM along their characteristic path.35 Thus, a second-order ordinary differential equation is solved for each component of the position: ����� ���� ������������������������ (9) where the right-hand side is the total force acting on the particle, calculated following a penalty approach.36 This force is defined by an error term to approximate the steady-state particle velocity to the negative ion drift velocity Wn calculated with FEM in the domain  2. The ratio k/m defines the penalty constant chosen large enough to reach fast convergence while being small enough to avoid divergence of the particle velocity solution from Wn. Sensitivity analysis is performed for each simulated case in order to obtain suitable penalty constants such that the error of the particle velocity estimated with (9) is 8 minimum. Different penalty constants are used in the radial and axial directions in our implementation. It is found that penalty constants k/m equal to 10��� kg s-1 and 10���kg s-1 for the radial and axial particle drift respectively are suitable values for the calculations. If each particle is assumed to represent a large number of negative ions (i.e. a superparticle approach) as in Ref. 8, a full time-domain simulation would require the injection of the particles at very small time steps until steady state is reached. This would result in a memory intensive and time consuming calculation. Since the steady-state current of the jet of negative ions injected into the domain  2 is constant, the alternative approach proposed in Ref. 35 is here followed. It considers that the particles released at the start of the calculation represent instead a continuous stream of negative ions per unit time. The number of ions per unit time represented by each particle released under the current In of the negative charged plume is given by the effective frequency frel: ���� ��� �∙� (10) Then, the derivative of the equivalent continuum number density of negative ions �� � at a mesh point xp in  2 is given by: ��� ����� �� �������� � ��� �������11� where δ is the Dirac delta function.35 The number density �� � at each mesh element is then computed by integrating equation (11) over time until the released particles can be completely traced out of the domain  2. Once the corrected number density �� � of negative ions in the plume is calculated with the particle solution, the FEM solution is solved again using the previous steady-state solution as initial value. However, the number density of negative ions in the domain  2 used for the solution of equations (57) is taken instead from the solution of �� �. Then, the updated FEM solution is used to solve the particle problem again. The calculation continues following the sequencial segregated solution of the FEM and particle problems until convergence is reached. Since the drift of electrons, positive holes and the liquid motion is weakly affected by the thickness of the ionic plume in the asymptotic zone, the solution rapidly converges. The calculation typically requires few iterations (between 3 to 4 iterations in the considered cases) until the distribution of the charged plume in the particle solution converges. IV. RESULTS A. Charge distribution close to the needlepoint Figure 2 shows representative examples of the steady-state density distribution of electrons and positive holes in the close proximity of the needle point for the charged plume at two different voltage levels of 2.2 and 2.6 kV. At these voltages the plume is generated under two different regimes: weak injection and space-charge limited injection respectively. As it can be seen, the injected electrons in both cases flow into the liquid following the electric field lines, increasing the radius of the charged plume. In contrast to EHD plumes under autonomous ionic injection,6,7 the electronic charged plume does not become constricted by the inward velocity of the liquid. Since electrons have a drift velocity significantly larger than the liquid velocity, their motion is therefore drift-dominated rather than convection-dominated. For this reason, the electronic 9 plume extends radially following the electric field lines under both weak or space charge limited injection. Therefore, the electronic charged plume expands beyond the radius of the needle electrode, contrary to what is generally assumed for plumes generated by ionic injection. 6,7 However, the highest density of electrons is mainly concentrated in front of the needle tip. As electrons drift in cyclohexane under high electric field, positive holes are created by ionizing collisions in the region close to the tip. 25 Since these charge carriers have lower velocity than electrons, they accumulate in front of the needle, reaching number densities as high as 10 21 holes/m 3 . However, the charge of the positive hole cloud does not become significantly large as to distort the electric field in front of the needle. For this reason, the positive holes do not significantly affect the rate of electron injection at the surface even under space-charge limited injection. FIG 2. Surface plot of the logarithmic concentration (in m -3 ) of positive holes n p (left) and electrons n e (right) in the close proximity of the needlepoint under the a) weak injection and b) space-charge limited injection cases. Streamlines for the liquid velocity (left) and the electric field (right) are also included. With the drift of electrons into areas of lower electric field, attachment processes become more efficient. Consequently, the density of electrons starts to decrease as they attach to neutral molecules forming negative ions. As it is shown in Figure 3, the charged layer consists mainly of negative ions at larger distances from the needlepoint electrode. Due to the lower mobility of negative ions, they start accumulating directly below the needlepoint electrode. For the case of weak injection of electrons, the maximum number density of the charge plumed is below about 10 21 ions/m 3 . Under those conditions, the electric field in front of the needle is not distorted by the charged plume. However, as the density of the charge plume core increases above 10 21 ions/m 3 , the electric field in entire gap becomes significantly distorted by the plume, defining the spacecharge limited injection regime. 16 parameter. Typical values of �� obtained from the simulation for the weak and the space-charge-limited injection cases are around 1017 and 31018 electrons/m3 respectively. Then, the parameter C for the two simulated cases is equal to 0.3 and 5.5 respectively. The dimensionless instability parameter �����������⁄ describing the ratio of the driving Coulomb force and the damping viscous forces, is in the order of about 7000. The calculated values of the mobility parameter M = �� � ⁄�� and the electric Reynolds number R= ��/M2 are 14 and 35 respectively. E. Current-voltage characteristic of the plume Once the properties of the plume are self-consistently obtained in steady state, the characteristic current-voltage (IV) curve of the configuration can be calculated. Figure 10 shows a comparison between the self-consistently estimated IV- characteristic and the measured values reported in Ref 22. Three different conduction regimes can be observed from the experimental values as the voltage increases. The first zone is defined by the ohmic regime of conduction in the liquid, dominant at low voltages (below about 1900 V). This regime is not simulated here since the low-field conductivity due to ionic-pairs in the liquid has not been considered. At higher voltages above the onset of field emission starts the second zone of weak injection of electrons (between about 1900 to 2400 V). The third zone defines the space-charge limited injection zone as the voltage further increases. FIG. 10. Simulated current compared with experimental results from Ref 22. It is found that there is an excellent agreement in the weak injection regime between the current predicted by the model and the experimental results in Figure 10. Such a good agreement is reached only when considering the contribution of collisional ionization in the simulation. Notice that calculations obtained neglecting electron multiplication in the liquid (not shown here) lead to underestimation of the currents for up to one order of magnitude. Good agreement between the model estimates and the experimental values is also found for the low-voltage end of the space-charge limited regime. However, the differences between the self-consistently estimated currents and the measured values start to increase above 2600 V. These differences may be explained by the overestimation of the liquid velocity as the laminar hydrodynamic plume in the experiment in Ref. 22 may have turned unsteady into a periodic or turbulent regime as the voltage increased. In such a case, the model cannot reproduce any 3D instability in the flow due to the assumed axial symmetry. 17 In order to estimate the current under a purely turbulent regime of the plume, an additional simulation is made. Thus, the calculations are repeated forcing both the radial and axial components of the negative ion velocity to be defined by the electrohydrodynamic turbulent mobility ����� �� ⁄,5 discarding the liquid laminar motion. This causes the increase of the radial flow of negative ions, widening of the charged plume. It also leads to the reduction in the negative ion longitudinal velocity compared with the laminar case. Interestingly, the calculations under those conditions underestimate the current in the space-charge limited regime as shown in Figure 10. This result indicates that the plume generated under the space-charge limited regime above 2.6 kV in Ref. 22 deviates from the purely laminar flow here described. However, the liquid flow in the experiment does not seem to become fully turbulent regime either otherwise the currents measured in Ref. 22 for the spacecharge limited regime had been lower. Unfortunately there is no information in the literature about the conditions for the transition of axial-symmetric EHD laminar plumes into an unsteady periodic or a fully turbulent regime. However, the stability analysis in two-dimensional plumes (e.g. in Refs. 10, 38) can provide qualitative clues about such a transition in the simulated axial-symmetric flow. Observe that the transition into an unsteady periodic regime in a blade injector has been numerically predicted to occur at dimensionless applied electric potentials Ti between 2000 and 10000.10 This range is in the same order of magnitude as the Ti values here obtained. On the other hand, a detailed theoretical analysis in Ref. 33 has shown that the transition of the plume is governed by the effective Prandlt number ����� and the stability factor G. The effective Prandlt number, defined as ����� � ����� ⁄��. The stability factor is defined as � �4���/4��/�, where ��� ����2����� ⁄ is the electrical Grashof number in a liquid with kinematic viscosity �.38 The maximum Preff for the case here simulated at 2.6 kV is 8, located at about z/d=0.05. The stability number at the same distance from the tip is about 28. Based on the results in Ref. 33 it has been shown that 2D plumes with a �����value of about 8 become unstable under G values larger than about 2. This result qualitatively supports the argument that an unsteady flow in the EHD plume in the experiment in Ref. 18 develops above about 2.6 kV. As a last remark, it is important to point out that the presented analysis is aimed to improve the understanding of EHD motion of dielectric liquids prior to breakdown under very high electric fields. It can also be relevant for the analysis of micro/nano structured electrodes producing single and two-phase fluid flow and heat transfer.39 Possible applications of injection from such electrodes can therefore be heat transfer enhancement and flow control, in macro and micro-systems.40, 41 V. CONCLUSIONS The self-consistent numerical modeling of both the charge conduction and the motion of the liquid in EHD laminar plumes produced by electron field injection from ultra-sharp needle tips in cyclohexane is presented. The model also considers electron multiplication by impact ionization and negative ion formation by attachment in the liquid. EHD plumes generated under weak injection or space-charge limited injection are analyzed. The main conclusions are summarized as follows: - The radius of the charged plume in front of the needle increases as the emitted electrons follow the electric field lines by their drift-dominated motion. Thus, the charged plume in front of the needle does not become constricted by the inward velocity of the liquid as it happens in the case of EHD plumes generated by autonomous ionic injection. Then, as electrons move away from the needle, they become efficiently attached into negative ions at lower electric fields. Consequently, negative ions start accumulating at some distance from the needle. Significant distortion of the 18 electric field by the charged ionic plume is observed under space-charge limited injection when densities larger than 1021 ions/m3 are reached. - Since the longitudinal electric field distribution and the liquid velocity radial profiles of the plume are selfconsistently estimated with the model, other simulated plume characteristics are compared with simpler asymptotic solutions. It is shown that the simulated liquid velocity and radius of the charged and hydrodynamic plumes in the asymptotic region are in excellent agreement with the integro-differential solution proposed in Ref. 6 under weak injection. However, the asymptotic solution proposed in Ref. 6 does not accurately estimate the plume characteristics under space-charge limited injection since it neglects the radial variations of the electric field in the charged core. - In agreement with previous studies, the charged core of the plume under weak injection is shown to have a rather constant density with a well-defined radius. However, the distortion of the electric field produced under spacecharge limited injection spreads the plume radially, smoothing out the radial density profile of the charged core. - Even though using the Finite Element Method (FEM) introduces numerical diffusion in the solution of the model, the obtained estimates are similar (with differences lower than 5 %) to those obtained with a diffusion-free particle method. Thus, it is shown that a FEM solution with consistent stabilization and an adequate mesh can provide reliable results of EHD plumes, even in the asymptotic region. - The self-consistent estimation of the current-voltage (IV) characteristic of EHD laminar plumes is reported for the first time. It is found that the estimated IV characteristic of laminar plumes produced by electron field emission in cyclohexane is in good agreement with experiments. Differences in the simulated and measured IV characteristics have been found only at high voltages when such plumes are expected to turn into an unstable motion regime. These differences are expected since instabilities in the plume cannot be reproduced under the axial-symmetry assumed in the model. ACKNOWLEDGMENTS MB would like to acknowledge the financial support of the Swedish strategic research program StandUp for Energy. PAV thanks the financial support from the Spanish Ministerio de Econom�a y Competitividad under Research Project No. FIS2014-54539-P. REFERENCES 1A. Castellanos, Electrohydrodynamics, New York, Springer, (1998). 2P. A. Vázquez, A. T. Pérez and A. Castellanos, "Thermal and electrohydrodynamic plumes. A comparative study," Physics of Fluids 8, 2091 (1996). 3A. Zhakin, "Electroconvective jets in liquid dielectrics," Mekh Zhidik. Gaza 6, 13 (1984). 19 4F. M. J. McCluskey and A. T. Pérez, "The electrohydrodynamic plume between a line source of ions and a flat plate," IEEE Transactions on Electrical Insulation 27, 334 (1992). 5P. Atten, B. Malraison and M. Zahn, "Electrohydrodynamic Plumes in Point-plane Geometry," IEEE Transactions on Dielectrics and Electrical Insulation 4, 710 (1997). 6P. A. Vázquez, A. T. Pérez, A. Castellanos and P. Atten, "Dynamics of electrohydrodynamic laminar plumes: Scaling analysis and integral model," Physics of Fluids 12, 2809 (2000). 7P. A. Vázquez, A. T. Pérez, A. Castellanos and P.Atten, "Origin of the double boundary layer structure of EHD plumes," Annual Report of the Conference on Electrical Insulation and Dielectric Phenomena, 211 (1998). 8R. Chicón, A. Castellanos and E. Martin, "Numerical modelling of Coulomb-driven convection in insulating liquids," J. Fluid Mech 334, 43 (1997). 9P. A. Vazquez, E. C. Vera, A. Castellanos and T. C. Rebollo, "Numerical calculations of two-dimensional EHD plumes with finite element and particle methods," Annual Report Conference on Electrical Insulation and Dielectric Phenomena, 706 (2003). 10A. T. Pérez, Ph. Traoré, D. Koulova-Nenova and H. Romat, "Numerical Study of an Electrohydrodynamic Plume between a Blade Injector and a Flat Plate," IEEE Transactions on Dielectrics and Electrical Insulation 16, 448 (2009). 11P. Atten, Electrohydrodynamic Instability and Motion Induced by Injected Space Charge in Insulating Liquids, IEEE Transactions on Dielectrics and Electrical Insulation 3, 1 (1996). 12A. Gazaryan, A. Sitnikov, V. Chirkov, and Y. Stishkov, A Method for Estimation of Functional Dependence of Injection Charge Formation on Electric Field Strength, IEEE Transactions on Industry Applications 53, 3977 (2017) 13A. T. Pérez, P. A. Vázquez and A. Castellanos, "Dynamics and Linear Stability of Charged Jets in Dielectric Liquids," IEEE Transactions Industrial Applications 31, 261 (1995). 14P. Traoré, M. Daaboul and C. Louste, "Numerical simulation and PIV experimental analysis of electrohydrodynamic plumes induced by a blade electrode," J. Phys. D: Appl. Phys. 43, 225502 (2010). 15P. Traoré, J. Wu, C Louste, Q. Pelletier and L. Dascalescu, "Electrohydrodynamic Plumes due to Autonomous and Nonautonomous Charge Injection by a Sharp Blade Electrode in a Dielectric Liquid," IEEE Transactions on Industry Applications 51, 2503 (2015). 16P. Atten and J. C. Lacroix, "Non-linear hydrodynamic stability of liquids subjected to unipolar injection," J. Mec 18, 469 (1979). 17F. J. Higuera, "Electrohydrodynamic flow of a dielectric liquid around a blade electrode," Physics of Fluids 12, 2732 (2000). 18J.Wu, P.Traore, C.Louste, D.Koulovab and H.Romata, "Direct numerical simulation of electrohydrodynamic plumes generated by a hyperbolic blade electrode," Journal of Electrostatics 71, 326 (2013). 19S. A. Vasilkova, V. A. Chirkov, and Yu. K. Stishkov, Electrohydrodynamic flow caused by field-enhanced dissociation solely, Physics of Fluids 29, 063601 (2017). 20 20Yu. K. Stishkov, V. A. Chirkov, and A. A. Sitnikov, Dynamic CurrentVoltage Characteristics of Weakly Conducting Liquids in Highly NonUniform Electric Fields, Surface Engineering and Applied Electrochemistry 50, 135 (2014). 21R. Gomer, Field emission and field ionization, New York, (1993). 22A. Denat, J. P. Gosse and B. Gosse, "Electrical Conduction of Purified Cyclohexane in a Divergent Electric Field," IEEE Transactions on Electrical Insulation 23, 545 (1988). 23R. G. Forbes and J. H.B. Deane, R. G. Forbes and J. H. B. Deane, "Reformulation of the standard theory of Fowler-Nordheim tunnelling and cold field electron emission," Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 463, 29072927, 2007. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 463, 2907 (2007). 24M. Haidara and A. Denat, "Electron multiplication in liquid cyclohexane and propane," IEEE Transactions on Electrical Insulation 26, 592 (1991). 25W. F. Schmidt, "Electronic Conduction Processes in Dielectric Liquids," IEEE Transactions on Electrical Insulation 19, 389 (1984). 26W. Schmidt, Hot electron mobility and electron attachment in non-polar fluids, in The Liquid State and Its Electrical Properties, Edited by E. Kundhardt, L. Christophoou and L. Luessen, Series B: Physics 193, New York, Plenum Press, 273 (1987). 27Y. K. Suh, "Modeling and Simulation of Ion Transport in Dielectric Liquids  Fundamentals and Review," IEEE Transactions on Dielectrics andElectrical Insulation 19, 831 (2012). 28R Morrow and N Sato, The discharge current induced by the motion of charged particles in time-dependent electric fields; Sato's equation extended, Journal of Physics D: Applied Physics 32, L20 (1999). 29S. Kumara, Y. V. Serdyuk and S. M. Gubanski, "Charging of Polymeric Surfaces by Positive Impulse Corona," IEEE Transactions on Dielectrics and Electrical Insulation 16, 726 (2009). 30R. W. Hockney and J. W. Eastwood, Computer Simulations using Particles, New York, McGraw-Hill (1981). 31O. Chanrion and T. Neubert, A PIC-MCC code for simulation of streamer propagation in air, Journal of Computational Physics 227, 7222 (2008). 32COMSOL Multiphysics 5.2, Application Library Manual (2015). 33A. N. Brooks, and T. J. R. Hughes, Streamline upwind/PetrovGalerkin formulations for convection dominated flows with particular emphasis on the incompressible NavierStokes equations, Comput. Methods Appl. Mech. Eng. 32, 199 (1982). 34P. Petrino, Y. H. Gaston-Bonhomme, J. L. E, Chevalier, Viscosity and Density of Binary Liquid Mixtures of Hydrocarbons, Esters, Ketones, and Normal Chloroalkanes, J. Chem. Eng. Data 40, 136, (1995). 35COMSOL Multiphysics 5.2, Particle Tracing Module Users Guide (2015). 36T. Chandrupatla, Finite Element Analysis, Universities Press, India, 2004. 37K.J. Bathe, D. Hendriana, F. Brezzi and G. Sangalli, Inf-Sup Testing of Upwind Methods, International Journal for Numerical Methods in Engineering 48, 745 (2000). 21 38P. Vazquez, Dynamics and Stability of Electrohydrodynamic Plumes, (In Spanish) PhD Thesis, Universidad de Sevilla, 1998. 39M.K.Russel, P.R.Selvaganapathy and C.Y.Ching, Effect of electrode surface topology on charge injection characteristics in dielectric liquids: An experimental study, Journal of Electrostatics 72, 487 (2014). 40J. Seyed-Yagoobi, Electrohydrodynamic pumping of dielectric liquids, J. Electrostatics 63, 861 (2005). 41J. Darabi and H. Wang, Development of an electrohydrodynamic injection micropump and its potential application in pumping fluids in cryogenic cooling systems, J. Microelectromech. Systems 14, 747 (2005).