A numerical study of electro-flow-focused micro-jets

The generation of liquid micro-jets is central to a wide range of technological applications, including sample delivery systems for serial femtosecond crystallography (SFX) (Chapman et al., 2011; Knoška et al., 2020), mass spectrometry (Fenn et al., 1989), inkjet printing (Yudistira et al., 2010), drug delivery (Zhao, 2013), electrospinning of nanofibers (Loscertales et al., 2004) and micro-propulsion systems (Huh and Wirz, 2022; Thuppul et al., 2020). In these processes, achieving precise control over jet diameter, stability and breakup dynamics is essential for performance and reliability. Among the available strategies, flow focusing (FF) (Gañán-Calvo, 1998) and electrostatic focusing – electrosprays (ES) (Gañán-Calvo et al., 2018) are two of the most widely studied and applied mechanisms. While each method can independently generate thin, stable jets, both have inherent limitations that limit their performance.

FF relies on a pressurised gas stream that mechanically squeezes the liquid meniscus through a small orifice, producing jets whose diameters scale with the gas-to-liquid momentum ratio. This technique is robust and mechanically straightforward, but requires high gas flow rates, which lead to jet breakup and whipping at high gas velocities. Electrostatic focusing, in contrast, applies a strong electric field to deform and accelerate the liquid meniscus, forming jets that are orders of magnitude smaller than the capillary diameter. However, ES performance is limited by charge accumulation in the surrounding medium, which distorts the electric field and destabilises the jet at high voltages. It also exhibits a very narrow parameter space for jet stability. Adding a more viscous additive to the sample liquid (e.g. glycerol) can produce a more stable, extended ES jet. This regime, known as electrospinning, has already been successfully implemented in SFX experiments (Sierra et al., 2012).

The concept of electro-FF (EFF) (Gañán-Calvo et al., 2006; Zupan et al., 2023) combines these two mechanisms, leveraging the mechanical confinement of gas pressure and the interfacial stresses induced by an electric field. Experiments have shown that EFF can produce jets that are thinner, longer and more stable than those generated by FF or ES alone. Despite a large body of experimental results, a systematic numerical analysis of the influence of combined electrostatic and hydrodynamic forces on a micro-jet remains lacking. In particular, the roles of liquid properties (viscosity, surface tension, electrical conductivity) and process parameters (gas pressure, applied voltage, flow rate) in shaping the jet remain difficult to disentangle experimentally because of the small length scales and transient dynamics involved.

Numerical modelling provides a powerful tool to address this gap, enabling detailed resolution of coupled fluid flow, charge transport and interfacial dynamics under well-controlled conditions. It has already been successfully used for axisymmetric simulations of FF (Zahoor et al., 2018) and ES (Herrada et al., 2012), three-dimensionally resolved FF (Kovačič et al., 2024b) and ES (Cândido and Páscoa, 2023; Huh, 2023) and other types of sample delivery systems for SFX, such as multi-component FF jets (Zahoor et al., 2024a), non-Newtonian FF jets (Zahoor et al., 2024b) and liquid sheets (Belšak et al., 2021; Kovačič et al., 2026, 2024a). In a prior study, we developed, verified and experimentally validated an electrohydrodynamic (EHD) numerical model capable of capturing the gas–liquid interface, with a focus on applications to EFF (Zupan et al., 2024b, 2025). The verification comprised three no-flow tests, all of which showed excellent agreement with the analytical solutions of the electrical equations and demonstrated that the solver is charge-conservative. Experimental validation was conducted using a free-boundary problem in which the interface deformation is governed by the applied electric field. The experimental measurements were in good agreement with the numerical predictions. Following that, a systematic computational approach is therefore needed to quantify the interplay between electrostatic and hydrodynamic forces in EFF and to establish scaling relations that can guide experimental design.

In this study, we present a fully coupled multiphase numerical model of electro-flow-focused microjets, formulated within the finite volume framework of OpenFOAM. The model solves the incompressible two-phase Navier–Stokes equations, augmented with electrostatic body forces and charge conservation and utilises the volume-of-fluid (VOF) method to capture the liquid–gas interface. Using this framework, we conduct a comprehensive parametric analysis of the influence of gas pressure, electric potential, liquid flow rate, electrical conductivity, viscosity and surface tension on jet diameter and length. The results provide new insights into the governing physical mechanisms of EFF and offer design guidelines for optimising micro-jet generation in various applications.

The numerical model is formulated for incompressible, immiscible flow of two Newtonian fluids. Additionally, a high-voltage electrostatic field is imposed and fully coupled with the fluid flow.

The governing equations are as follows:

where $u$ represents the fluid velocity, $ρ$ density, $P$pressure, $μ$dynamic viscosity, $σ$ surface tension, $κ$ surface curvature, $δs$ is the surface Dirac distribution function and $n$ is the interface normal. Density is retained in the divergence term, as it varies spatially. Electric force $fE$ is calculated as:

where $ρE$ represents the volumetric charge density, $E$ the electric field and $ε$ the electric permittivity. Terms on the r.h.s. in equation (3) follow directly from the divergence of the Maxwell stress tensor. The derivation is detailed in Zupan et al. (2025). The electric field is calculated from Gauss’s flux theorem:

which has the form of the Poisson equation. Owing to the irrotationality of electrostatic fields, $∇×E=0$⁠, the electric field can be rewritten in terms of a scalar potential field, $E=−∇ϕ$⁠, where $ϕ$ is the electric potential. The physical model is closed with the volumetric charge conservation equation, derived from the Poisson–Nernst–Planck equation (Zupan et al., 2025), assuming negligible thermal diffusion of ions and fully reacted solvents:

Here, $K$ denotes the electric conductivity. This approach is also valid for unipolar ion injection (Montanero, 2024). The interface is resolved using the VOF method (Hirt and Nichols, 1981), which prescribes a phase function $α(p,t)$⁠, such that α = 1 in the liquid phase (index $l$⁠), α = 0 in the gas phase (index $g$⁠) and 0 < α < 1 at the interface, representing a mixture of the two fluids. Rheological parameters are adjusted at the interface using weighted arithmetic mean (WAM) for the density and viscosity:

and using the geometric mean (GM) for the electric permitivity and conductivity:

The GM method was chosen for interpolating the electric properties based on a recent study (Zupan et al., 2025), showing that it yields a more physically accurate distribution of free charge at the interface compared to WAM (López-Herrera et al., 2011). The evolution of the phase fraction α is governed by an advection equation that ensures mass conservation while maintaining a sharp interface:

where $uc$ is the artificial compression velocity introduced to limit the spread of the interface:

Here, $Cα=1$ is a parameter adjusting the intensity of the compression, $ϕα$ is the velocity flux and $Sα$ is the interface surface area. The surface tension is modelled as a continuum surface force (Brackbill et al., 1992), computing the interface normal from the gas to the liquid phase $n$ in equations (2) and (11) as $n=∇α/∇α$ and estimating the interface curvature as $κ=−∇⋅n$⁠.

The EFF device is modelled after the experimental setup introduced by Gañán-Calvo et al. (Gañán-Calvo et al., 2006), which operates as follows. A grounded metallic cylinder with a circular orifice on one side is pressurised, ejecting gas from the orifice. A capillary needle, held at high voltage and coaxially aligned with the orifice, is inserted into the cylinder to inject the liquid, forming a thin liquid jet. A grounded impact plate is placed downstream of the cylinder to collect and discharge the liquid droplets. An axisymmetric computational domain was constructed to model the experimental setup and is shown in Figure 1. The dimensions are as follows: needle capillary radius R = 50 µm, orifice radius r = 25 µm, cylinder and needle thickness t = 25 µm and distance between the needle tip and the orifice H = 100 µm. The outlet chamber is a cylinder with a length of 1.575 mm and a diameter of 1 mm.

A structured hexahedral mesh with five levels of incremental refinement was used (see Figure 1), ranging from 5 µm at the coarsest level (L5) to 0.3125 µm at the finest level (L1). The cell size at the finest level was selected based on the mesh sensitivity study shown in Figure 2. The refinement region was chosen to ensure that the gas–liquid interface is always captured with the finest mesh. The total number of computational cells was 439,740. Boundary conditions are listed in Table 1.

Equations (1)–(11) are solved numerically using the finite volume method (Moukalled et al., 2016) discretisation, as used in the OpenFOAM (Weller et al., 1998) software. Specifically, the standard VOF solver interFoam was upgraded to include electric effects by calculating the electric and charge density fields – equations (4) and (5) and computing the electric force, equation (3), used in the momentum equation – equation (2).

All simulations used a combination of second-order accurate and bounded numerical schemes. Time discretisation is handled by the Crank–Nicolson scheme with an off-centring coefficient of 0.9, providing second-order temporal accuracy with slight numerical damping for stability (Ferziger et al., 2020). Spatial derivatives are computed using the Gauss linear scheme for gradients, which is second-order accurate on structured or smoothly varying meshes (Weller et al., 1998). Convective terms for the momentum and electric equations are discretised with the Gauss MUSCL scheme – a Monotonic Upstream-centred Scheme for Conservation Law – which is second-order accurate and bounded in resolving sharp gradients without introducing spurious oscillations (Van Leer, 1979). The volume fraction advection term uses the van Leer limiter, a Total Variation Diminishing scheme that preserves boundedness and interface sharpness, commonly used in algebraic VOF methods (Hirt and Nichols, 1981). Associated flux correction terms are treated with the Gauss linear scheme. The diffusive term employs Gauss linear discretisation, and all Laplacians are computed using Gauss linear orthogonal functions, which assume orthogonal face-cell alignment. Interpolation and surface-normal gradients are handled using linear and orthogonal schemes, respectively.

The volume fraction field is solved using the MULES algorithm (Courant and interface Courant numbers equalled 0.25), which incorporates interface compression and uses two correction steps. A Gauss–Seidel-preconditioned iterative solver (Jasak, 1996) is used for the velocity field, while the pressure correction is handled by the preconditioned conjugate gradient (PCG) solver with diagonal incomplete Cholesky preconditioning (Van Der Vorst, 1992). Electricity-related fields use either the Preconditioned Bi-Conjugate Gradient Solver with the Diagonal Incomplete LU decomposition preconditioning or PCG, balancing robustness and efficiency for asymmetric matrices. The simulation advances using the PIMPLE algorithm, which combines the PISO and SIMPLE algorithms for improved convergence in transient flows (Issa, 1986). All equations are solved without under-relaxation. The solution was computed in parallel on an HPC cluster (AMD EPYC 7402, 2.8 GHz) using two nodes with 48 cores each.

Heptane and air at 1 bar and 300 K were selected as the working fluids for the base case. To investigate the influence of fluid properties on the EFF jet, the heptane parameters were systematically varied. The base values and the ranges used in the parametric study are summarised in Table 2.

We used standard dimensionless numbers to characterise the governing physical phenomena, compare the relative importance of different forces and generalise the results across different fluid systems and scales. They are presented in Table 3 and expressed in terms of process parameters and material properties. $Q$⁠, $dj$⁠, $ΔPg$⁠, $E0$ and $κl$ are the liquid flow rate, jet diameter at the orifice, gas pressure drop through the orifice, electric field strength and the liquid ionic mobility. Specifically, we calculated the Reynolds (inertial/viscous forces), Weber (inertial/surface tension forces) and the Capillary (viscous/surface tension forces) numbers. To capture the ratio between the electrical and mechanical FF of the jet, we adopted the electrostatic influence number $χ$⁠, as proposed by (Gañán-Calvo et al., 2006). It should be stated that the parameter χ does not fully capture the relative influence of electric and hydrodynamic forces in the present configuration, as it does not explicitly account for the applied voltage. This limitation should be addressed in future experiments to provide a broad parameter space.

The diameter of the jet at the orifice stabilises after the initial transient period of jet formation (normally between 0.5 and 1 ms). We calculate its time-averaged value $dj$ as:

over a time period $Δt$ of 1 ms. Average jet length $lj$⁠, measured from the orifice exit to jet breakup, was calculated in the same way. The total simulation time ranged from 2 to 6 ms.

To investigate the combined effects of ES and FF, we first established a baseline by simulating each process independently. For FF, the gas pressure drop across the orifice was incrementally increased from 2 kPa to 6 kPa in steps of 2 kPa, while maintaining the gas Mach number below 0.3 to satisfy the incompressibility assumption. For ES, the applied voltage was increased from 0.75 kV to 1.5 kV in 250V increments. The liquid flow rate was 0.75 $mlh-1$ for all cases. The corresponding results are shown in Figures (3) and (4). At voltages above 1.25 kV, the ES meniscus is inverted because of the extremely high electric field (>10⁷ $Vm-1$⁠) at the capillary tip, exerting intense electrostatic pressure on the liquid interface. The hydrodynamic pressure is insufficient to counteract this force, causing the meniscus to retract inside the capillary at the given low flow rate and small capillary diameter.

The inversion is most likely not physical but an axisymmetric response to the applied conditions, although jets remain stable in the inverted configuration up to 2.5 kV. Another consideration regarding the maximum applied voltage is the corona discharge, which occurs above the breakup strength of air (3 MV/m). We observed that this limit is reached in the 0.75–1 kV applied potential range, depending on the charge density at the interface. Jets with higher electric conductivity reach this limit earlier. Based on these observations, we limited all subsequent simulations to 0.75 kV to ensure the physical system response and operation are outside the corona discharge region.

Figure 5 compares the FF case at 4 kPa, the ES case at 0.75 kV and the combined EFF case (⁠$χ=1$⁠) under the same conditions, all at a liquid flow rate of 0.75 ml/h. The combined effect results in a thinner (−28% compared to FF; −22% c.t. ES) and longer/shorter (−21% c.t. FF; +6% c.t. ES) jet. Figure 6 presents the radial velocity profiles at four axial positions. In all cases, the jet interface accelerates more rapidly than the bulk flow immediately downstream of the orifice. This radial velocity gradient is neutralised within the first 100 µm when gas is present (FF and EFF), whereas it reverses in the ES case due to skin friction. Overall, the EFF configuration increases the centreline velocity by 70% and 150% relative to FF and ES, respectively.

One of the limiting factors in ES performance is the accumulation of free charges in the gas surrounding the liquid jet, which distorts the applied electric field, accelerates the gas phase and diminishes the effectiveness of electrostatic focusing. This phenomenon is evident in Figure 6, where the gas-phase velocity increases downstream of the orifice. In addition to the FF effect, the gas flow enhances focusing in the combined EFF regime by convecting free charges away from the region where electrostatic focusing is strongest (Zupan et al., 2025). This is further illustrated in Figure 7, which compares the radial distribution of volumetric charge density at the orifice for the ES and EFF cases. Although the peak charge densities at the interface differ by only 10%, the charge concentration in the gas phase differs by several orders of magnitude.

The remainder of this section is divided into subsections, each analysing the influence of a different parameter in the study on the EFF.

Both increased gas pressure and higher electric potential contribute to a reduction in jet diameter, thereby enhancing the focusing effect (see Figures 8 and 9). The skin friction effect increases with increased gas pressure (Figure 8), introducing the radial velocity gradient inside of the liquid jet.

The electric potential plays a critical role in shaping the meniscus: at low values (e.g. < 0.25 kV) the meniscus is changed from the FF shape. At moderate values (e.g. ∼0.5 kV), the meniscus begins to adopt an ES-like shape and at higher voltages it begins to deform. Beyond the reversal threshold (not shown here because of non-physicality), further increases in voltage result in only minor reductions in diameter but significantly deform the meniscus, pushing it deeper into the capillary because of the intensified electric field. In our case, increasing the electric potential above 2.5 kV led to an unstable meniscus.

At a specific combination of process parameters, a spatially stable jetting regime can be established, characterised by a jet length standard deviation below 2 µm (compared to the typical range of 20–30 µm). In this regime, not only is the jet length constant but also the droplet size and their spatial positions remain fixed in time, forming a highly periodic pattern. In our simulations using heptane and air (see Table 2), these conditions correspond to Re = 61.8, We = 5.4, Ca = 0.009 and χ = 1. The resulting flow exhibited a repeating sequence of two primary droplets separated by a smaller satellite droplet, with a spatial wavelength of approximately 90 µm and a temporal period of 1 ms, as illustrated in the space–time diagram in Figure 10.

The liquid flow rate exerts a significant influence on both the jet diameter and the meniscus shape (see Figure 11). As the flow rate increases, the liquid supply counters the tendency of the meniscus to deform under the influence of the electric field. This effect forms thicker jets and widens the meniscus base.

At higher flow rates, the system becomes progressively less sensitive to the electric field, as also noted by (Gañán-Calvo et al., 2006). The electric field’s capacity to shape the interface and accelerate the liquid diminishes relative to the increasing inertial and pressure-driven flow components. This shift in force balance also leads to a notable reduction in the size of the meniscus recirculation zone, further promoting flow stability.

The observed diameter trend approaches the expected square-root scaling, and although further convergence could likely be achieved at higher flow rates, such conditions are not considered here as the focus is on regimes where electric effects remain comparable to hydrodynamic forcing.

The effect of liquid electric conductivity on the EFF jet characteristics is illustrated in Figures 12 and 13 at constant process parameters. As shown in Figure 12, electrical conductivity strongly influences jet diameter and length. At low conductivity values (⁠$K≈10-7$ S/m), the meniscus is not yet deformed, and the resulting jet is similar to the FF jet at the same parameters. An increase in conductivity by a factor of 100 deforms the meniscus and further focuses the jet. A further increase in conductivity increases the electric shear stress at the interface, further focusing the jet. Such behaviour indicates competition between electrostatic forces and hydrodynamic focusing, with the charge relaxation time relative to the hydrodynamic timescales governing the jetting process. At moderate to high conductivity values, the increased charge density at the interface raises the electric field strength above the gas-ionisation threshold. Furthermore, the jet breakup departs from the physical observations, as the produced droplets are too large. Thus, at moderate and high conductivity values, the applied electric potential should be reduced to ensure the model produces physically relevant solutions.

The charge–relaxation dynamics explain the reduction in recirculation zones with increasing conductivity. The relaxation time is relatively long at low conductivity, and induced charges penetrate deeper into the bulk liquid. This weakens the interfacial Maxwell stresses, allowing the gas flow to dominate and promote the extended recirculation zones near the nozzle entrance. As conductivity increases, the charge relaxation time decreases and charges accumulate almost instantly at the interface. The resulting localised electric stresses strongly accelerate the liquid at the meniscus, efficiently transferring momentum into the jet. This direct acceleration reduces the need for bulk recirculation and stabilises the flow in the feeding region, resulting in smaller and weaker vortical structures behind the meniscus.

Similar conclusions can be drawn from the charge density distribution at the orifice (Figure 13). For the lowest conductivity case, the peak volumetric charge density is relatively small, with a broader penetration into the bulk of the jet. At high conductivity, the charge density peak shifts closer to the surface and becomes sharper, with a maximum exceeding 10°C/m. The higher charge concentration consequently drives the electric stresses at the interface. Owing to the high charge density and relatively high electric stresses at the interface, a localised stream of the charge is discharged from the interface, which influences the velocity field (see Figure 12 at $5.2⋅10-5$ S/m).

When surface tension is raised, the liquid strongly resists deformation because of its natural tendency to minimise surface area. This enhanced restoring force hinders the jet from being significantly elongated by the electric field, resulting in stabilisation only when the jet remains relatively short. In other words, higher surface tension effectively “pinches” the liquid column, producing larger diameters and reducing the overall length of the stable jet, as illustrated in Figure 14. Moreover, the recirculation zone expands with increasing surface tension, as it becomes the principal energy sink of the system (Zahoor et al., 2020).

Viscosity, in contrast, serves as an internal damping mechanism. A more viscous liquid suppresses the amplification of disturbances along the jet and allows the liquid thread to be stretched thinner. Thus, higher viscosity generally promotes thinner electro-jets, contrary to the influence of increased surface tension (see Figure 15).

The original paper by Gañán-Calvo et al. (2006), which established the EFF method, also provided theoretical scaling relations that matched the experimental data well, with the bulk of the data points lying just below the theoretical line. They proposed the theoretical jet diameter:

and the characteristic jet diameter and flow rate for non-dimensionalisation:

The total available energy per unit volume $ΔΨ$ is $ρlσ2Kl2ε0−21/3+ΔPg$⁠. Figure 16 compares the non-dimensional simulated jet diameter to the non-dimensional flow rate. The solid line represents the analytical solution [equation (13)] and the simulated data points are drawn with empty squares. Agreement between simulations and theory is on the same order as the agreement between theory and experiments (Gañán-Calvo et al., 2006), with the bulk of the simulations placed just below the theoretical line. Pure electrospraying modes are not included in the plot.

A numerical model for electro-flow-focused micro-jets was developed by considering incompressible two-phase Navier–Stokes equations in the presence of an electrostatic field. The model is solved using FVM and a modified VOF solution procedure. The EHD approach of the model was previously validated against established analytical predictions and experimental data (Zupan et al., 2025), demonstrating its ability to resolve the interaction of electrostatic and hydrodynamic focusing mechanisms. Here, it is also validated using experimentally established scaling laws for EFF jets. A systematic parametric study is conducted to examine the effects of material properties and operating conditions on jet morphology and stability.

The main findings can be summarised as follows:

• The combined EFF regime produces jets that are thinner than those formed by either ES or FF alone.

• Increased gas pressure enhances jet focusing by convecting free charges away from the meniscus, thereby stabilising the electric field distribution.

• Increased electrical potential intensifies interfacial stresses and reduces jet diameter, but excessive voltages cause non-physical meniscus inversion and system destabilisation.

• Higher liquid flow rates thicken the jet, counteract meniscus inversion and reduce sensitivity to electrostatic forces while decreasing recirculation zones.

• Low liquid electrical conductivity leads to diffuse charge penetration and large recirculation zones, whereas high conductivity localises charge at the interface, generating stronger focusing and more stable jets.

• Elevated surface tension increases restorative forces, producing shorter and thicker jets while expanding the recirculation zone. Higher viscosity dampens interfacial instabilities, leading to thinner jets.

• A new spatially stable jetting regime was found, where droplet size and their spatial positions remain fixed in time, forming a highly periodic pattern. Such highly regular jetting modes are particularly attractive for SFX delivery systems, where synchronisation of droplet emission with the beam pulse train can substantially improve the effective hit rate and reduce sample consumption. Further experimental research is needed to validate the existence of the regime and to determine how the presence of protein crystals affects it.

However, the existence of other stable jetting regimes outside the parametric space reported here is likely, as has already been experimentally determined for a similar nozzle design (Cruz-Mazo et al., 2019; Zupan et al., 2024a, 2023). The parametric space used here was chosen based on the original experimental paper on EFF (Gañán-Calvo et al., 2006) and the corresponding range of dimensionless numbers. Furthermore, it should be noted that the jet diameter is likely more robustly predicted than the absolute breakup length, which is mostly useful for comparing trends within the simulated parameter space, as has also been noted by a recent FF study (Kovačič et al., 2024b).

Although thermal effects are not explicitly modelled in the present study, Joule heating in regions of high charge density and subcooling because of gas expansion may locally influence material properties such as viscosity and surface tension; a detailed, thermally coupled investigation of these effects is left for future work (Zupan et al., 2026). Similarly, the axisymmetric model is not capable of capturing the non-symmetric multi-jet instabilities, which are most likely present in the high voltage regimes, and the physical model is limited to cases where no gas ionisation and corona discharge occur.

Overall, this study demonstrates that EFF jet formation is governed by a delicate balance between electrostatic and hydrodynamic forces, modulated by fluid properties and operating conditions. The insights gained provide a framework for optimising micro-jet generation in applications that require precise jet diameter and stability control. This work should be regarded as a first, systematic parametric study aimed at isolating the fundamental effects of the electric field force on micro-jet formation using a model fluid system, rather than a direct optimisation study for SFX conditions, which will be addressed in future work using aqueous, biologically relevant solutions.

Table A1 compiles all reported cases from the paper, including process parameters, material properties, measured jet diameters and lengths and corresponding dimensionless numbers.