Electrified Cone Formation in Perfectly Conducting Viscous Liquids: Self-Similar Growth Irrespective of Reynolds Number
Theodore G. Albertson, Sandra M. Troian
TL;DR
This study demonstrates that the self-similar growth of electrified conic tips in perfectly conducting viscous liquids occurs regardless of Reynolds number, revealing universal shape scaling and force dynamics through finite element simulations.
Contribution
It shows that conic tip growth is always self-similar across all Reynolds numbers, challenging previous assumptions about the influence of viscous and inertial forces.
Findings
Self-similar tip growth is universal across Reynolds numbers.
Universal conic shape with a half-angle dependent on Maxwell stress.
Viscous stresses at finite Reynolds number cannot be neglected.
Abstract
Above a critical field strength, the free surface of an electrified, perfectly conducting viscous liquid, such as a liquid metal, is known to develop an accelerating protrusion resembling a cusp with a conic tip. Field self-enhancement from tip sharpening is reported to generate divergent power law growth in finite time of the forces acting in that region. Previous studies have established that tip sharpening proceeds via a self-similar process in two distinct limits - the Stokes regime at and the inviscid regime . Using finite element simulations to track the acceleration of an electrified protrusion in a perfectly conducting Newtonian liquid in vacuum held at constant capillary number, we demonstrate that the conic tip \textit{always} undergoes self-similar growth irrespective of Reynolds number. The computed blow up exponents at the tip for the…
| Quantity | Scaling | Rescaled |
|---|---|---|
| variable | ||
| Coordinates | ||
| Operators | ||
| Interface shape | ||
| Electric potential | ||
| Electric field | ||
| Pressure | ||
| Velocity | ||
| Time | ||
| Self-similar | Fitted exponent | |
| collapse time | ||
| Self-similar | ||
| time interval | ||
| Reynolds | ||
| number | ||
| Capillary | ||
| number | ||
| Bond | ||
| number |
| Re | Mesh | ||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| A | |||||||||||||
| A | |||||||||||||
| A | |||||||||||||
| A | |||||||||||||
| A | |||||||||||||
| A | |||||||||||||
| B | |||||||||||||
| A | |||||||||||||
| B | |||||||||||||
| A | |||||||||||||
| B | |||||||||||||
| A | |||||||||||||
| B | |||||||||||||
| A | – | – | |||||||||||
| B | – | – | |||||||||||
| B | – | – | |||||||||||
| B | – | – | |||||||||||
| B | – | – | |||||||||||
| B | – | – | |||||||||||
| B | |||||||||||||
| B | |||||||||||||
| B | |||||||||||||
| B | |||||||||||||
| B | |||||||||||||
| B | |||||||||||||
| B | |||||||||||||
| B |
| (N/m) | (Pas) | (kgm3) | Re | Ca | |||||
|---|---|---|---|---|---|---|---|---|---|
| Ga | 0.720 | 0.1 | 7.0834 | ||||||
| 50 | 7.0834 | ||||||||
| 500 | 7.0834 | ||||||||
| 5000 | 7.0834 | ||||||||
| 50 000 | 7.0834 | ||||||||
| In | 0.556 | 0.1 | 7.0834 | ||||||
| 50 | 7.0834 | ||||||||
| 500 | 7.0834 | ||||||||
| 5000 | 7.0834 | ||||||||
| 50 000 | 7.0834 | ||||||||
| Cs | 0.070 | 0.1 | 7.0834 | ||||||
| 50 | 7.0834 | ||||||||
| 500 | 7.0834 | ||||||||
| 5000 | 7.0834 | ||||||||
| 50 000 | 7.0834 |
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Corresponding author (URL/email) : ]www.troian.caltech.edu; [email protected]
Electrified Cone Formation in Perfectly Conducting Viscous Liquids:
Self-Similar Growth Irrespective of Reynolds Number
Theodore G. Albertson and Sandra M. Troian [ California Institute of Technology, T. J. Watson Sr. Laboratories of Applied Physics, 1200 E. California Blvd., MC 128-95, Pasadena, CA 91125, USA
(August 8, 2019)
Abstract
Above a critical field strength, the free surface of an electrified, perfectly conducting viscous liquid, such as a liquid metal, is known to develop an accelerating protrusion resembling a cusp with a conic tip. Field self-enhancement from tip sharpening is reported to generate divergent power law growth in finite time of the forces acting in that region. Previous studies have established that tip sharpening proceeds via a self-similar process in two distinct limits - the Stokes regime at and the inviscid regime . Using finite element simulations to track the acceleration of an electrified protrusion in a perfectly conducting Newtonian liquid in vacuum held at constant capillary number, we demonstrate that the conic tip always undergoes self-similar growth irrespective of Reynolds number. The computed blow up exponents at the tip for the terms in the Navier-Stokes equation and interface normal stress condition reveal the different forces at play as Re increases. Rescaling of the tip shape by the capillary stress exponent yields excellent collapse onto a universal conic tip shape with interior half-angle dependent on the magnitude of the Maxwell stress. The rapid acceleration of the liquid interface also generates a thin surface boundary layer with very high local strain rate. Additional details of the modeled flow, applicable to cone growth in systems such as liquid metal ion sources, help dispel prevailing misconceptions that dynamic cones resemble conventional Taylor cones or that viscous stresses at finite Re can be neglected.
Credit line: This article has been submitted to the Physics of Fluids - a link to the final version will be provided once published.
I Introduction
In a short and influential paper in 1884, Lord Rayleigh demonstrated by a simple argument based on linear stability analysis why a charged spherical mass of liquid acted upon solely by electrical and capillary forces represents a state of unstable equilibrium. Maxwell’s equations had only been developed a couple decades earlier so Rayleigh’s hydrodynamic analysis was especially novel and insightful. He predicted too that the instability would cause the spherical mass to distort into a peaked structure that would emit a fine jet due to the growing imbalance between the destabilizing electrical pressure and the stabilizing capillary pressure. Several years later in 1890, J. Larmor became interested in the electrification of capillary waves on a liquid surface Larmor (1890). Seemingly unaware of Rayleigh’s work, he examined the linear stability of an electrified liquid surface by appealing to the unsteady form of Bernoulli’s equation from which he derived the dispersion relation correlating the rate of growth to the wavelength of ripples along the interface of a charged liquid film. These two early studies thereby revealed the ease with which an electrified liquid surface can develop and sustain protrusions which narrow and sharpen in time.
I.1 Instability analyses by Tonks and Frenkel
For decades, this phenomenon remained mostly a curiosity until 1935 when L. Tonks Tonks (1935), a researcher at General Electric who was investigating ion emission from a solid metal electrode, turned his attention to sparking and charge emission in liquid metals. Tonks closely examined the liquid distortion process and electrical breakdown of a liquid mercury surface in close proximity to a counter electrode for field strengths approaching several megavolts/cm. In particular, he noted that a planar liquid metal surface rapidly generated an accelerating protrusion whose narrowing shape concentrated the local disruptive electric field so as to generate a liquid cusp capable of ion emission. His studies eventually confirmed that field induced ion emission from liquid metals occurs at lower field strengths than required for solid metals. Informed by these observations, he carried out a number of calculations which resulted in expressions and estimates relating the amplitude of the distortion of an electrified liquid surface to the minimal field strength required to induce a sustained deformation, and the typical time required for interface rupture to occur resulting in ion emission. Tonks emphasized two key questions that continued to drive his work - namely, what was the relation between the applied field strength, the distortion amplitude and the pressure at the tip of the accelerating protrusion, and how fast did this protrusion advance in response to the pressure at the tip? Realizing that a complete and rigorous solution to this problem was rather complicated, he instead developed an insightful approximate treatment that revealed critical aspects of the distortion runaway process. Shortly after Tonk’s work was published, J. (Ya. I.) Frenkel Frenkel (1935, 1936) provided a more general derivation of this dynamic process by investigating the velocity of surface waves triggered by an applied electric field. Such periodic distortion of the surface of a perfectly conducting liquid for sufficiently high electric field strengths is now known as the Larmor-Tonks-Frenkel instability. Despite the number of studies and contributions to this field during the past century, the two questions originally posed by Tonks have not yet been completely answered.
The investigations by Larmor and Frenkel relied on a stability analysis of capillary waves on a perfectly conducting liquid subject to an electric field, through which they uncovered that the wave frequency becomes complex valued once the critical field strength for sustained distortion is exceeded. By contrast, Tonks focused on the dynamic evolution of a single protrusion emanating from a flat liquid surface with an initial small bump modeled as a spherical cap. His analysis revealed that beyond a critical field strength, the pressure at the cap apex cannot support equilibrium since the electric pressure at the tip - now known as the Maxwell pressure or Maxwell stress - increases as the square of the apex height. This quadratic dependence cannot be balanced by the restoring hydrostatic or capillary pressures, which at best scale linearly or inversely with apex height, respectively. Therefore, any small increase in the apex amplitude simply increases the net local unbalanced pressure. This strong imbalance causes a runaway process by which the distortion progressively increases in amplitude. By relying on various approximations for the initial conditions and simplified functional forms for other quantities, he was successful in deriving an equation for the acceleration of an eccentric liquid shape with time. Remarkably, Tonks also seemed to realize that the growth of the liquid tip proceeded by a self-similar process, as is apparent from his sketch of the time sequence of an evolving interface shape in Fig. 8 of Ref. [Tonks, 1935].
I.2 Experimental studies by Krohn and co-workers
Interest in this phenomenon began to spread and by the early 1960’s, V. E. Krohn Krohn (1984) and co-workers at TRW recognized that the spontaneous process by which an electrified liquid can undergo fission to produce charged droplets and ions could be of enormous relevance to the burgeoning field of electrostatic propulsion Perel et al. (1971); Bailey et al. (1972); Zafran et al. (1973); Huberman and Rosen (1974) and in particular, development of so-called “space thrusters”. While other researchers focused on colloid thrusters based on electrospray from dielectric liquids, Krohn instead concentrated on cusp formation in liquid alkali metals. Within a few years, he successfully demonstrated ion emission from small liquid droplets of cesium, gallium and mercury supported on the tip of a tungsten needle maintained at a surface potential of a few kilovolts. By 1975, Krohn was able to demonstrate stable emission lasting several hours, a remarkable technological achievement which essentially paved the way for a number of technologies based on liquid metal ion sources (LMIS). By the late 1970’s, researchers Krohn and Ringo (1972); Seliger et al. (1979) recognized the potential of LMIS for high-resolution ion microprobe analysis, micro-milling and micro-implantation. Attention therefore soon turned to identification and control of those operational parameters affecting the diameter of the ion source area, the maximum brightness achievable and the energy spread of the ion beam. Significant industrial support led to rapid and far reaching developments in this field, so much so that LMIS now forms the basis for instrumentation ranging from high resolution focused ion beam systems for enhanced etching and deposition, to scanning ion microscopy, micro-milling, ion mass spectrometry, implantation and lithography Orloff et al. (2003). Equally exciting are numerous current developments aimed at space micropropulsion systems based on LMIS such as field emission electric propulsion Wright and Ferrer (2015); Bharti and Chalia (2017) (and references therein).
I.3 NASA interest in liquid metal ion sources for space micropropulsion
The most recent 2015 NASA Technology Roadmap NAS (2015) identifies in-space micropropulsion technologies based on micro-electrospray propulsion as key enabling technologies, defined as small scale systems which can deliver thrust using a conductive fluid and electrostatic fields to extract and accelerate ions, charged droplets, molecules or clusters of molecules depending on whether the propellant is a liquid metal or ionic fluid. These systems are intended to provide high specific thrust levels characterized by low power, small form factors and extended lifetimes. The ultimate goal is to extend current know-how gained from the fabrication and manufacture of microelectromechanical systems (MEMS) to micropropulsion platforms consisting of hundreds to thousands of needle microarrrays externally wetted by a film of liquid metal electrified on demand to emit charged particles. NASA envisions that such miniaturized systems will deliver thrust levels suitable for precision attitude, rendezvous and docking maneuvers without the use of reaction wheels. Such capability is anticipated to advance formation flight and attitude actuation for small distributed spacecraft You (2017).
The fluid dynamical process leading to cusp formation in free surface, perfectly conducting, viscous liquids has remained a key problem in the field of electrohydrodynamics for decades. Prior studies of cusp formation have focused almost exclusively on two asymptotic limits - the Stokes flow regime at and the inviscid regime at . Since the dynamic behavior of rapidly accelerating, electrified viscous fluids is a rather complex process, only these two asymptotic limits have so far yielded analytic insight by revealing which forces at the liquid tip are responsible for self-similar growth. Natural questions to ask at this point are whether the self-similar evolution of the electrified tip persists for all Reynolds numbers, how does the tip shape depend on the Reynolds number Re and capillary number Ca, and what are the exponents characterizing the divergent behavior at the tip for the different forces involved just prior to ion emission.
Using finite element simulations to track the acceleration of an initially tiny electrified protrusion in a perfectly conducting Newtonian liquid in vacuum held at constant capillary number, we demonstrate that the conic tip always undergoes self-similar growth irrespective of Reynolds number. These axisymmetric simulations allow high resolution inspection of the acceleration and sharpening of an initial rounded tiny protrusion on an otherwise flat liquid layer. The computed blow up exponents at the tip for the terms in the Navier-Stokes equation and interface normal stress condition reveal the different forces at play as Re increases. Rescaling of the tip shape by the capillary stress exponent yields excellent collapse onto a universal conic tip shape with interior half-angle dependent on the Maxwell stress magnitude. Additional details pertaining to the importance of viscous stresses even at higher Reynolds number as well as the formation of a surface boundary layer characterized by a very high local strain rate are also presented.
II Notable Studies of Electric Field Induced Cuspidal Formation in Liquids
II.1 Electro-stationary state known as a Taylor cone
The analyses of Rayleigh Rayleigh (1882) and Larmor Larmor (1890) made evident that above a critical surface charge or above a critical field strength, a perfectly conducting mass of liquid must incur deformation represented by a state of unsteady equilibrium. In 1964, Taylor Taylor (1964) examined what would be the equilibrium shape which would everywhere on the surface allow for a perfect balance between the local capillary pressure and local Maxwell pressure. This static shape was found to be a cone of constant slope with an interior half-angle of , now known as the Taylor angle. While the two opposing normal stresses are everywhere in mechanical balance, however, they undergo divergence at the conic apex, a feature which troubled subsequent researchers. Over the years, this result has tended to cause considerable confusion between Taylor’s stationary configuration and the more prevalent dynamic conic shapes, such as conic cusps, which represent non-equilibrium configurations arising from the increasing imbalance in forces at the liquid tip. Some authors Orloff et al. (2003) have attempted to clarify this issue by referring to dynamic protrusions as “Gilbert-Gray cone-jets”. Gilbert Gilbert first reported the attraction of water to a charged piece of amber with subsequent formation of a conic protuberance from whose vertex was emitted a very fine jet and Gray Gray (1731) reported similar jet formation not only in water but liquid mercury as well.
II.2 Zubarev analysis of self-similar conic cusp formation in the inviscid limit ()
In a seminal paper in 2001, Zubarev Zubarev (2001) analyzed the dynamic distortion accompanying the response of a perfectly conducting layer of liquid to an initially uniform electric field. His asymptotic analysis of the unsteady Bernoulli equation applied to the moving liquid interface was conducted under the assumptions of incompressible, inviscid and irrotational flow. Key to this analysis was the assumption that the electric field strength near the advancing tip rapidly and appreciably exceeds the applied external field value such that local conditions prevail and the governing electric field strength decreases with distance from the apex. The external field uniformity condition could also then be replaced with a far field boundary condition specifying vanishing electric field strength at infinity. This approach allowed him to analyze the fluid motion near the singular apex without reference to any particular geometry, ultimately revealing the existence of self-similar solutions for the velocity potential, electric potential and fluid interface shape when rescaled in time by the variable where is the dimensionless collapse time and the dimensionless time. We next outline Zubarev’s analysis leading to this prediction of self-similar growth and the exponents characterizing the divergence at the conic tip of the Maxwell pressure, capillary pressure and kinetic energy density.
We note that for the remainder of this work unless otherwise designated, dimensional quantities are designated by lower case letters and dimensionless quantities by upper case letters. An exception to this is the electric field vector represented in dimensional form by and in dimensionless form by . The vector denotes a unit normal vector at the moving interface pointing outward from the liquid domain. Bold faced variables are used to denote vector and tensor quantities.
Zubarev examined the evolution and growth of an electrified axisymmetric protrusion in a perfectly conducting inviscid fluid. Here, the interface function denotes the boundary separating the liquid and vacuum regions. The vacuum domain supports an electric field distribution , where the electric potential satisfies the Laplacian condition . The velocity field within the liquid region is described by where the velocity potential satisfies for an incompressible liquid subject to irrotational flow. Given that the liquid is also perfectly conducting, its interior domain defines a Gaussian volume devoid of electric field and therefore only the moving interface experiences an electrical normal force which pulls liquid toward the vacuum domain. Conservation of mass and momentum are enforced through the unsteady form of Bernoulli’s equation for inviscid irrotational flow given by
[TABLE]
for a liquid with constant density and constant surface tension . The underbrace terms represent the opposing local Maxwell and capillary pressures acting in the direction normal to the moving interface. Zubarev non-dimensionalized the governing variables and equations by introducing the following scalings:
[TABLE]
where . The boundary conditions used to solve Eq. (1) and the corresponding Laplacian fields were chosen to be
[TABLE]
where . A scaling analysis of the governing Bernoulli equation and boundary conditions revealed the existence of self-similar solutions, here indicated by a tilde, according to which
[TABLE]
where denotes the dimensionless asymptotic collapse time when the interface shape assumes a universal form.
Zubarev constructed self-consistent solutions to these equations by matching the self-similar solutions in the inner region at the liquid apex to the far field solutions describing a conventional stationary Taylor cone with interior half-angle . To leading order as (or equivalently in self-similar coordinates), he found that the asymptotic functions scaled as
[TABLE]
where
[TABLE]
where denotes the polar angle in spherical coordinates, denotes the Legendre polynomial of order 1/2, and is a constant such that . As evident from Eq. (17), the solution corresponds to fluid motion that is spherically symmetric such that liquid is transported to the conic apex by streamlines which are always tangent to the moving interface. The point therefore acts as a point sink, a feature Zubarev highlighted.
These self-similar solutions describe divergent behavior of the Maxwell pressure , capillary pressure and flow speed at the advancing conic apex according to which
[TABLE]
Through this analysis, Zubarev was the first to report that a perfectly conducting liquid subject to inviscid, irrotational and incompressible flow set by local conditions due to Maxwell and capillary forces undergoes local self-similar evolution that in the far field smoothly matches onto the conventional stationary Taylor cone shape. Zubarev described the resulting interface shape as a conic cusp.
II.3 Self-Similar Cusp Formation in the Viscous Regime
II.3.1 Numerical Studies by Suvorov, Litvinov and Zubarev
Around the time Zubarev reported that the inviscid limit supports self-similar growth of the conic tip, Suvorov and Litvinov Suvorov and Litvinov (2000) conducted a computational study based on the axisymmetric geometry shown in Fig. 2 to better understand the influence of viscous effects. They employed a cell and marker technique with a coordinate transformation method for adaptive grid generation Belozerkovskii (1994), which allowed smooth and stable evolution of the fluid interface despite development of significant and rapidly increasing curvature at the conic apex. The numerical technique they implemented prevented onset of spurious surface oscillations near the region of high curvature, a problem which had plagued previous researchers Zheng and Linsu (1988).
According to Fig. 2, the axisymmetric volume of interest was partitioned into two contacting domains - an upper domain comprising the vacuum region subject to an initial uniform electric field of strength V/m and a lower domain comprising a perfectly conducting, incompressible Newtonian liquid. Given the large disparity in density and viscosity, the hydrodynamics of the vacuum domain was neglected. The normal component of the electric field strength along the vacuum-liquid interface (held at constant electric potential) was updated in time according to its spatiotemporal location. The Navier-Stokes equation for the liquid domain was solved self-consistently subject to the kinematic boundary condition governing the motion of the moving interface. The initial thicknesses of the vacuum and liquid domains, along with their radial extent, were chosen to be 20 m. The material constants of the liquid were chosen match those of liquid mercury. The initial condition for the interface represented a flat film containing a tiny Gaussian bump centered at the origin of the form where and . The chosen dimensions and material constants therefore corresponded to laminar flow conditions at the value .
For an initial Gaussian bump of fixed lateral extent, the simulations showed that a thicker liquid layer tended to produce a more rounded protrusion advancing rapidly in time, while a thinner liquid layer produced a more pronounced cuspidal shape that advanced more slowly, presumably due to stronger viscous effects. Suvorov and Litvinov called these structures dynamic Taylor cones to emphasize that such shapes do not resemble conventional stationary Taylor cones. Instead, the initial protrusion was observed to deform into a cuspidal formation with a cone-like tip terminated by a rapidly shrinking spherical cap. Suvorov and Litvinov also carried out simulations in which the viscous term in the normal stress boundary condition was altogether omitted. They concluded from these additional studies that the viscous normal force at the advancing tip is always positive - that is, it too acts in the same direction as the capillary force but is much smaller in magnitude. They reported that the resulting shapes obtained with or without the viscous normal stress were very similar but did not include any comparative figures.
The numerical scheme upon which this worked was based, which allowed stable interface growth over half a decade in time, successfully capture divergent behavior at the conic apex of such quantities as the curvature radius, electric field strength and axial velocity. These divergent quantities led in turn to runaway behavior in the Maxwell pressure, capillary pressure, viscous normal stress and kinetic energy density. Their study, conducted at a fixed value , demonstrated that the (negative) Maxwell normal stress pulling upward on the liquid apex rapidly overtakes the (positive) normal stresses due to capillary and viscous forces. Suvorov and Litvanov concluded by urging that more such computational studies be conducted to assist theorists in developing analytical treatments for the dynamic behavior of electrified liquids used in liquid metal ion sources. While their work established key features characterizing the runaway cusp formation process in viscous fluids, their results were restricted to a single value of Re. They also did not report the actual values obtained for the exponents characterizing the various divergent forces at the accelerating tip. Soon thereafter, there followed similar computational studies by Suvorov Suvorov (2004) and Suvorov and Zubarev Suvorov and Zubarev (2004) simulating the behavior of liquid gallium at the melting point at flow conditions corresponding to a smaller value of the Reynolds number. In these studies, the viscous normal stress term in the interface boundary condition was altogether omitted based on earlier indications Suvorov and Litvinov (2000) that the viscous term in the normal stress boundary condition tended to be much smaller than the capillary term. In these simulations, the initial thicknesses of the vacuum and liquid domains, along with their radial extent, were reduced to 5 m while the initial applied electric field strength value was increased to V/m. The initial condition for the interface was again chosen to be a flat film with a tiny Gaussian bump centered at the origin of the form where and . These Gaussian parameters, corresponding to an initial radius of curvature of the liquid protrusion equal to , were chosen to approximate the radii of supporting needles typically used in liquid metal ion sources. The chosen dimensions and material constants in these studies therefore corresponded to a value .
The numerical results by Suvorov Suvorov (2004) are shown in Fig. 1. Only this set of figures will be discussed here since the results presented in Ref. [Suvorov and Zubarev, 2004] are similar. Snapshots of the interface shape depicted in Figs. 1 (a) and (b) illustrate both the far and near field views of an evolving protrusion. These images reveal how the rapidly accelerating interface undergoes self-similar collapse to an apical shape resembling a sharpening conic tip with a rapidly decreasing radius of curvature. The curves in Fig. 1(c) illustrate blow up behavior at the apical tip in the magnitudes of the Maxwell and capillary pressure as the system approaches the collapse time of 123.542 nsec. Fig. 1(d) represents a log-log plot of the magnitude of the apical values of the Maxwell pressure, capillary pressure and kinetic energy density upon approach to the collapse time. The solid lines shown were reported to yield exponent values corresponding to divergent behavior of the Maxwell pressure, capillary pressure and kinetic energy density of 0.67, 0.68 and 0.63, respectively. The dashed lines we have superposed represent the corresponding exponent values 0.70, 0.77 and 0.78, which provide a closer fit to Suvorov’s data.
II.3.2 More recent numerical studies of conic formation in electrified liquids
More recently, Collins et al. Collins et al. (2008) used the finite element method with elliptic mesh generation and adaptive time integration to simulate liquid cone formation in perfectly conducting, perfectly insulating and leaky dielectric liquids. A goal of their study was to understand why electrohydynamic tip streaming resulting in the emission of a thin fluid jet from the conical tip does not occur in perfectly conducting or perfectly insulating liquids. Their simulations of perfectly conducting liquids also revealed tip sharpening with corresponding divergent behavior in the Maxwell, capillary and viscous normal stress at the apex, in agreement with previous studies Suvorov and Litvinov (2000); Suvorov and Zubarev (2004); Suvorov (2004) - however, they did not extract the corresponding values of the blow up exponents. Most interestingly, this study confirmed that tip-streaming requires the existence of interfacial tangential stresses near the conic tip, which are ultimately responsible for the cone-jet transition observed. These numerical simulations successfully captured for the first time the various scaling laws governing the size of droplets formed from steady cone-jets as first predicted by de la Mora and Loscertales and Gañán-Calvo de la Mora and Loscertales (1994); Gañán-Calvo (2001).
Burton and Taborek Burton and Taborek (2011) have since conducted detailed boundary integral simulations to examine Coulombic fission of charged inviscid droplets. Their studies have revealed how small disturbances cause rapid elongation of an initial slightly prolate ellipsoidal droplet to produce polar cap regions resembling conic tips, which they have coined lemon-shaped drops. For perfectly conducting liquids, their results show how conic tips undergo self-similar sharpening due to local field self-enhancement with a corresponding blow up exponent of 2/3 for the capillary pressure and 1/3 for the charge density. Remarkably, they were able to confirm self-similarity over twelve decades in time - an extraordinarily long period - requiring stable integration able to resolve high resolution behavior at the advancing tip. Incidentally, they discovered that during conic sharpening, the interior half-angles never equal the Taylor cone angle value of , except perhaps at the very moment of blowup.
Different interface tracking schemes have since been implemented by other groups interesting in similar systems. For example, Garzon, Gray and Sethian Garzon et al. (2014) have examined the behavior of perfectly conducting, inviscid droplets situated in an initial uniform external electric field. A level set method was used to develop an Eulerian potential flow model for tracking drop evolution past breakup) while the flow velocity, Maxwell and capillary pressure along the free surface were extracted from axisymmetric boundary integral simulations. Depending on the value of the initial field strength, the resulting exponents at the apical tip for the Maxwell pressure, capillary pressure and kinetic energy density were found to range from 0.662 - 0.676, 0.57 - 0.62 (smaller precision reported due to larger fluctuations) and 0.672 - 0.736, respectively. In all cases examined, larger initial field strengths produced conic tips with larger interior half-angles, although the values were always smaller than the Taylor angle of .
II.4 Self-Similar Cusp Formation in the Stokes Regime ()
Additional studies during the past decade or so have examined the role of viscous forces in the formation of conic cusps. Vantzos and co-workers Betelú et al. (2006); Fontelos et al. (2008) employed boundary integral techniques for tracking the distortion of perfectly conducting droplets to quantify viscous effects. Their studies, focused exclusively on the Stokes flow regime (), have shown that for sufficiently large surface charge and/or sufficiently large external field, the polar caps of an initially rounded droplet deform into conic tips through a self-similar process as well. The interior half-angles of conic tips that form in charged or uncharged droplets range from about to , much smaller than the Taylor angle value of , irrespective of the values of the actual surface charge or field strength. These studies indicate that viscous effects tend to repress fluid tips from undergoing excessive sharpening.
In an initial study, Betel et al. Betelú et al. (2006) examined theoretically and numerically the deformation of an isolated charged droplet of a perfectly conducting viscous liquid surrounded by a dielectric viscous medium in which the only acting forces were electrostatic repulsion and capillarity. Since all the electric charge resides at the liquid/liquid interface, the surface charge density is given by , where is the electric potential and denotes the normal derivative. The repulsive electrical force per unit area is therefore given by . The interior () and exterior () liquid domains satisfy the Stokes equation and continuity equation given by
[TABLE]
[TABLE]
At the moving interface, the velocity and pressure fields must satisfy the normal stress boundary condition
[TABLE]
and the shear-free boundary condition
[TABLE]
where is the local mean interfacial curvature, is the interfacial tension and denotes the deviatoric stress tensor for a Newtonian liquid given by
[TABLE]
Here, indicates the matrix transpose, denotes the local unit normal vector pointing outward from medium 1 to 2, which is perpendicular to the two surface unit tangent vectors denoted by . The numerical results and supporting scaling arguments show that for isolated charged droplets, the interface curvature, charge density, electric field and velocity at the apical tip diverge according to as where represents the interval of time measured back from the collapse time. In a second study Fontelos et al. (2008), they examined charged droplets in an external electric field. For self-similarity to hold, it was found that the surface change density must scale as , the tip curvature as and the tip velocity as , where represents a given function dependent on the far field value of the cone interior half-angle. These studies revealed that the external field acting on a charged droplet represents only a vanishingly small perturbation to the dominant divergent field at the apical tip caused by the accumulation of local surface charge in that region. It was for this reason that the interior half-angle of charged droplets was found to be essentially independent of the value of the external field.
More importantly, Fontelos et al. realized from these studies that in the limit of vanishingly small Re, the dominant force balance leading to self-similar evolution of the conic tip is not due to opposing Maxwell and capillary stresses but opposing Maxwell and viscous stresses. From the leading order balance between viscous and Maxwell forces in Eq. (28), they concluded that the tip electric field strength scales as , which in turn requires that the tip curvature scale as and the tip velocity scale as , where denotes the self-similar coordinate.
In a separate calculation, they solved for the velocity and electric field distributions satisfying the Stokes and Laplace equation within an infinite cone of perfectly conducting fluid for variable interior half-angle. The velocity field yielded two solutions - one describing symmetric formation of cusps at opposite poles of an isolated charged droplet or an uncharged droplet in an external electric field and the other describing asymmetric formation of a single cusp. In both cases, the velocity field, found to scale as (derivation not shown), was required to smoothly and asymptotically match the far field velocity function obtained for an equipotential liquid cone subject to vanishing surface shear stress and normal stresses, which scales as , yielding the important relation . In the symmetric case, the scaling parameter was found to be constrained to the value , yielding . In the asymmetric case, was found to be a function of the interior cone angle, in which case neither nor yield simple integer ratios. Their simulations produced relatively constant values of the cone interior half-angle of about for a wide range of external fields, corresponding roughly to a value , in excellent agreement with the exponent value for the tip curvature extracted from the simulations.
III Computational Model for Perfectly Conducting Electrified Viscous Liquid in Vacuum
As discussed above, both the inviscid and Stokes flow regimes permit derivation of the exponents characterizing self-similar growth of the apical tip in perfectly conducting droplets and films. The exponents corresponding to the Maxwell pressure, capillary pressure, kinetic energy density and viscous stresses on approach to the collapse point indicate that for inviscid flow (), the dominant opposing stresses are given by the Maxwell pressure which promotes tip sharpening and the capillary pressure which represses formation of strong curvature. By contrast, the dominant opposing stresses in the Stokes flow regime () have been identified as the Maxwell pressure and the viscous normal stress. In our investigation which follows, we have explored whether self-similar evolution of the electrified tip persists for all Reynolds numbers, how the interior cone angle depends on the Reynolds number Re at fixed capillary number Ca, and extraction of exponent values characterizing blow up behavior for the various forces acting at the tip.
Shown in Fig. 2 is a sketch of the axisymmetric geometry modelled in this study. The system consists of a perfectly conducting body of liquid of radius and initial thickness held at constant potential . The initial liquid layer is separated from a grounded cylindrical electrode of radius by a vacuum layer also of thickness . This is the same geometry in the original study of Suvorov and Litvinov Suvorov and Litvinov (2000) and Suvorov and Zubarev Suvorov and Zubarev (2004) although those simulations were constrained to two values of the Reynolds number, namely and , respectively.
Since charge transport in perfectly conducting liquids is essentially instantaneous in comparison to the time scale for viscous flow, the electrostatic approximation can be used in the vacuum domain where where the electric potential then satisfies the equation
[TABLE]
subject to the boundary conditions:
[TABLE]
Accordingly, the electric field increases in time in response to the advance of the equipotential boundary toward the grounded electrode. Any protrusion along the liquid interface therefore undergoes local accelerated growth due to unbalanced field enhancement at the tip.
The hydrodynamic behavior of the electrified liquid is described by the Navier-Stokes equation for an incompressible, Newtonian liquid where the velocity field and pressure field satisfy the governing mass and momentum conservation equations given by
[TABLE]
For simplicity, gravitational effects are neglected by requiring that the Bond number . Within the liquid domain , the velocity fields and must satisfy the following boundary conditions prescribed along the axis of symmetry at and the bottom and side walls:
[TABLE]
The two remaining boundary conditions for are obtained from the constraint that the moving interface is shear-free, according to which
[TABLE]
where is the deviatoric stress tensor defined in Eq. (30). The boundary condition for the pressure field is enforced by the interface normal stress boundary condition
[TABLE]
where . The local mean curvature of the fluid interface is defined by the relation where in cylindrical coordinates
[TABLE]
The terms on the right hand side of Eq. (46) represent the capillary pressure, Maxwell pressure and viscous stress, respectively. The sign convention is such that the capillary pressure at the liquid apex is always positive and therefore repressive due to the concave shape while the Maxwell pressure is always negative, continually pulling the tip toward the counter electrode. The position of the advancing interface is obtained from the kinetic condition which requires
[TABLE]
where the surface gradient operator is given by where is the identity tensor.
The initial condition for the velocity field represents a quiescent state such that . The initial liquid interface configuration represents a flat liquid layer with a small Gaussian bump centered at the origin given by
[TABLE]
where the Gaussian peak amplitude and lateral extent are scale uniformly with the characteristic thickness such that and .
III.1 Scaling and Non-Dimensionalization
To probe the influence of inertial and viscous forces on the process of cusp formation, we implemented the scalings in Table 1 to convert equations and boundary conditions to dimensionless form. In addition to the parameter values specified in Fig. 2, the material constants characterizing the liquid metal, namely , and , represent the liquid density, viscosity and surface tension, respectively, at the specified temperature.
The dimensionless Navier-Stokes equation is then given by
[TABLE]
subject to the (dimensionless) normal stress boundary condition at the moving interface:
[TABLE]
where . In rescaled form, the initial interface shape is given by
[TABLE]
where the peak amplitude and extent are given by and . The remaining dimensionless equations and boundary conditions are easily obtained and not restated here.
III.2 Computational Study of Inertial and Viscous Effects at Fixed Ca
From Eqs. (51) and (52), it is evident that both the Reynolds and capillary number play a key role in the formation process. For too small a value of the capillary number, the flow will be overwhelmed by the force of surface tension and budding protrusions will be rapidly leveled. Too small values of Ca will therefore suppress formation of fluid elongations. By contrast, too large a value of Ca will generate ripples along the moving interface and complicate analysis because of interference effects involving multiple accelerating protrusions.
In this study, we therefore chose to fix the capillary number at a value which would allow the budding and advance of a single protrusion, as in the previous numerical studies of Suvorov Suvorov (2004) and Suvorov and Zubarev Suvorov and Zubarev (2004). To that end, we note that in their studies they chose to model the liquid metal gallium at the melting point with N/m and a layer thickness maintained at a voltage potential = 2.4 kV. From the scalings given in Table 1, the corresponding capillary number in their simulations was . We adopted this same value for all simulations presented in this work and indeed confirm that this value allowed stable formation of a single conic cusp unaccompanied by any interfacial oscillatory effects.
III.2.1 Consequences of holding Ca constant in simulations
In the results which follow, it is important to understand the consequences of holding the capillary number fixed when examining the flow behavior and interfacial shapes obtained at different Reynolds numbers. In particular, given a liquid with fixed material constants, the restriction of fixed capillary number imposes the key constraint that the ratio is a constant. Therefore, an increase in Re by a factor f reflects an increase in the voltage potential to and an increase in the characteristic length scale to . According to the scalings in Table 1), these increases lead to an effective decrease in the initial field strength value from to and a corresponding decrease in the characteristic flow speeds by a factor .
III.3 Details of Numerical Simulations
To solve the coupled equations and boundary conditions discussed in the previous section, we employed finite element simulations using the commercial software package COMSOL COMSOL Multiphysics using its arbitrary Lagrangian-Eulerian method for coupling the vacuum and liquid domains across the moving interface. The computational model was constructed assuming incompressible laminar flow conditions within the liquid domain, modeled as a perfect liquid conductor in contact with an inert vacuum domain. The liquid domain therefore contained no internal electric field and the electric field at the vacuum/liquid interface was always everywhere oriented in the direction normal and away from the moving interface. The electric potential distribution within the vacuum domain and on the liquid surface was updated at each time step according to the local position of the moving boundary modeled by a moving mesh. At each time step, all mesh element edges along the vacuum/liquid boundary were translated a distance given by the local speed of the liquid. All interior mesh elements were reconfigured by a uniform Laplace smoothing process so that boundary displacement and subsequent fluid motion proceeded smoothly throughout the liquid domain. The entire mesh was reconstituted whenever the mesh quality of the most distorted element fell below a specified value. Mesh generation relied on triangular quadratic Lagrangian elements for the velocity field and linear elements for the pressure field, so-called mixed order discretization. Second-order elements provided higher resolution with stability even at the highest Reynolds numbers.
Two different non-uniform meshes were employed, one optimized for low Re (mesh A) and the other optimized for high Re (mesh B). Flow at intermediate values of Re yielded similar results used either meshing scheme. For mesh A, the edge length of mesh elements adjacent to the moving interface were set equal to the value , where denotes the local interface radius of curvature, subject to an upper and lower bound. Additionally, at very late times where the curvature of the advancing cone increased rapidly, the mesh in the vicinity of the conic tip was iteratively refined by further dividing the triangular elements to achieve satisfactory resolution. To enforce mesh refinement, the position of the final location of the conic tip was first estimated from simulations conducted with a less refined mesh. The run was then repeated with significantly more mesh elements placed near the tip. The total number of mesh elements for the initial and final states ranged approximately from to . For mesh B, the number of elements positioned adjacent to the interface was directly prescribed, as was the ratio of the edge length of the smallest element adjacent to the conic tip to that of the largest element adjacent to the outermost domain boundary. The mesh was then iteratively refined by splitting elements in the vicinity of the tip, as done with mesh A. The number of elements for mesh B for the initial and final states ranged from approximately to .
Time integration was conducted using a variable first or second order backwards differentiation scheme with adaptive time stepping for a total of about 1000 time steps per simulation. Simulations were allowed to run just until interface acceleration near the tip led to irreparable mesh inversion. In all cases, the results reported in this work reflect simulations that were all terminated when the capillary pressure at the conic tip reached a fixed imposed value. This was enforced by stopping runs when the dimensionless radius of curvature at the conic tip (scaled by ) attained the value . In what follows, the final time at which this point was achieved is denoted by . This cutoff length scale typically corresponded to a length ten times smaller than the minimum mesh edge length along the moving interface. In previous work modeled on specific liquid metals Suvorov (2004); Suvorov and Zubarev (2004), researchers have instead terminated simulations once the maximum electric field strength at the conic tip exceeded the field emission value.
IV Results
The following simulations represent results obtained for Reynolds number ranging from at fixed capillary number Ca=7.0834. Details of the evolving interface shape and exponents characterizing the divergent behavior of governing forces at the conic tip are discussed next.
IV.1 Sharpening of conic cusp
Shown in Fig. 3 (a) - (d) are far field views of the dimensionless interface shape for values of Re spanning the viscous to inertial regime. For each value of Re shown, the two curves depicted in red in (a) - (d) bracket the time interval observed to undergo self-similar growth, to be discussed in detail below. The last curve for each time sequence shown represents the interface shape at the final simulation time , determined from the condition when the dimensionless radius of curvature at the conic tip (scaled by ) reached the value . Due to viscous retardation effects, it is evident that the time interval required for the apex to attain this small radius of curvature at is approximately two orders of magnitude longer than the time required at higher Reynolds numbers. The peak amplitudes in Fig. 3 also decrease in magnitude as Re is increased. This can be understood from the fact that at fixed Ca, an increase in Re by a factor f leads to a decrease in the initial electric field strength to , which weakens the strength of the Maxwell forces pulling on the liquid tip.
Shown in Fig. 3 (e) - (h) are magnified images of the evolving tip during the period of self-similar growth along with the values of the interior half-angle extracted. These angles were obtained from linear regression of the final interface shape over the range . (Since and are both scaled by , this angular value is a constant whether viewed in dimensional or non-dimensional coordinates.) These simulations confirm very rapid formation of conical tips with interior half-angles both above and below the conventional stationary Taylor cone value of . Larger values of Re are characterized by larger values of . This can be understood from the fact that at fixed capillary number, an increase in Reynolds number from Re to fRe leads to a decrease in the initial electric field strength from to . Weaker Maxwell forces lead to formation of wider cones, as expected.
IV.2 Velocity and strain rate profiles near conic apex
Shown in Fig. 4 are closeup images of the liquid tip at the final time for at . These are color coded by the magnitude of the dimensionless velocity . Scale bars are marked as shown in the lower left side of each image. The curves superimposed on the liquid domain represent the velocity streamlines. These indicate that the net flow at the liquid tip for all Re is oriented in the vertical direction, almost parallel to the vertical axis. The color maps indicate that the highest flow speeds are concentrated near the accelerating tip. However, while the maximum velocity undergoes a significant increase from to , there is relatively little change in the maximum value beyond that. This is due to the fact that the simulations were all terminated at about the same value of the (dimensionless) capillary pressure, which effectively establishes a cap on the corresponding maximum flow speed.
Shown in Fig. 5 are the same set of closeup images as Fig. 4 color coded by the magnitude of the dimensionless rate of strain tensor , where . Scale bars are marked as shown in the lower left side of each image. The region of largest strain rate occurs just below the surface of the conic tip. This region, which sustains the largest rate of deformation, undergoes the highest acceleration. The extent of this region of highest strain rate diminishes in size as the Reynolds number is increased.
The very large strain rates observed in the vicinity of the accelerating conic tip give to an interfacial boundary layer. Examination of the vertical velocity component along the axis of symmetry reveals the behavior of this layer with increasing Re. Shown in Fig. 6 are curves of the normalized dimensionless velocity versus distance for increasing values of Re where , denotes the final simulation time and the vertical coordinate of the conic apex at the final time. As Re is increased at fixed capillary number, the boundary layer becomes thinner, as expected. However, the correlation between the boundary layer thickness and Re observed in simple systems with constant surface velocity is not observed here since the interface is not only accelerating rapidly but also accelerating non-uniformally due to the spatially varying electric field strength along the moving interface. Further work is being conducted to elucidate the nature of this interfacial boundary layer.
Due to axisymmetry, the quantity reduces simply to the value . Therefore, the regios of highest strain rate in Fig. 5 correspond to the uppermost portions of the boundary layer where is largest. Inspection of Fig. 6 further reveals that the thickness of the ensuant boundary layers is far larger than the vertical displacements in the corresponding interface shapes in Fig. 5. The region of largest strain rate is therefore nestled deep within the boundary layer which extends much further into the bulk liquid.
IV.3 Self-similar evolution of conic tip
We now turn to evidence of self-similar growth irrespective of Reynolds number for dynamic regimes spanning the viscous to inviscid limit. The evolution of the magnitudes of the various forces acting at the conic apex point, as evaluated from the normal stress boundary condition given by Eq. (52) and the Navier-Stokes equation given by Eq. (51), are shown to sustain extended self-similar growth upon approach to the collapse time. Self-similar variable transformations obtained from exponents extracted from the numerical simulations conducted at various Re at fixed Ca confirm that the accelerating liquid tip undergoes collapse toward a conic shape whose interior angle increases with increasing Re.
IV.3.1 Determination of collapse time
In order to extract accurate power law exponents corresponding to the time at which the forces acting at the conic apex undergo divergence, it is first necessary to extract from the simulations the asymptotic collapse time . For each run corresponding to a given value Re, this collapse time was determined by examining the growth in the Maxwell pressure at the conic apex as a function of the interval plotted on double logarithmic axes and then choosing that value which produced the best linear fit over the final two and one half decades in time. Although this fitting procedure was limited to two and half decades in time, the results below indicate that in the majority of cases, power law growth was observed to occur over longer intervals, as evident in Fig. 7.
IV.4 Exponents for self-similar growth of normal stresses at conic apex
Shown in Fig. 7 (a)-(d) are representative double logarithmic plots showing the magnitude of the Maxwell, capillary and viscous normal stresses, reflecting the three terms on the right hand side of Eq. (52), for selected values of Re at . These terms describe the forces per unit area acting at the conic apex. The power law exponents listed alongside each term in (a) - (d) were extracted from those data points corresponding to the final two and half decades in time indicated by the dashed vertical lines, as also marked in images (e) - (h). These results indicate that the Maxwell (), capillary () and viscous () normal stresses at the conic apex diverge in time according to where . Shown in Table 2 are the complete set of exponents and fitting parameters for all simulations conducted in this study. In general, it was found that the exponent values and tend to decrease with increasing Re, indicating that the divergence of these opposing pressures at the conic tip is more rapid at higher Reynolds number.
As evident in Fig. 7(a)-(d), the Maxwell stress at the conic apex always exceeds the capillary and viscous normal stress irrespective of Reynolds number. At small values of Re, the Maxwell pressure exceeds the capillary and viscous normal stresses by at least an order of magnitude. As Re increases, the Maxwell and capillary pressures approach each other in magnitude, while the viscous normal stress decreases in comparison by several orders of magnitude. This trend is understandable since viscous effects always diminish in influence as Re is increased. Furthermore, while at small Re the difference in pressure between the Maxwell and capillary stress increases as , this difference in pressure is observed to decrease with increasing Re as the collapse time is approached.
With regard to the actual sign of these normal stresses, it is clear from the magnified images in Figs. 4 and 5 that the liquid shape near the apex is given by a convex curve. The curve representing the capillary stress in Fig. 7(a)-(d), computed from the first term on the right hand side of Eq. (52), is therefore always positive irrespective of Reynolds number. For the chosen value of Ca in these simulations, the capillary pressure therefore always represses enhancement of tip curvature but is overcome by the Maxwell pressure, computed from the second term on the right hand side of Eq. (52), which is always negative and increases in magnitude as the curvature increases. This self-sharpening of the tip goes hand in hand with the local self-enhancement of the electric field in that region. Were it the case that the electric field at the tip is not governed its local value but the global value set by the externally applied field, then self-similar growth would not be possible since the relative length and time scales could not therefore conform to yield self-similar scalings. This insight is what originally led Zubarev to seek the self-similar exponents reported in Ref. [Zubarev, 2001] for the inviscid case.
The viscous normal stress, computed from the final term in Eq. (52), can assume both positive and negative values depending on Re and . For example as shown in Fig. 7(a), the dip in the viscous normal stress at early time near is caused by a transition from positive to negative values. Beyond that time, the viscous normal stress becomes increasingly more negative as . Careful inspection of the velocity gradient near the tip along the central axis shown in Fig. 5 confirms that the maximum in the vertical velocity component does not occur at the conic apex but a very short distance below that point. This behavior therefore gives rise to a value . We are currently exploring whether vorticity effects near the interface are responsible for this small inward displacement of the vertical velocity maximum. Likewise in Fig. 7(b), the dip in viscous normal stress near signals a similar transition from positive to negative values. In Figs. 7(a) and (b), this transition occurs prior to the onset of the self-similar regime. In Fig. 7(c) however, this transition occurs within the self-similar period of growth, which unfortunately precludes extraction of a meaningful exponent - hence the designation “n/a”. For this reason, we did not report in Table 2 any exponent values for . For , the viscous normal stress remained positive at all times, allowing extraction of from the indicated curve in Fig. 7(d) for .
IV.5 Exponents for self-similar growth of Navier-Stokes terms at conic apex
We next evaluate the temporal behavior of each term in the Navier-Stokes equation at , which due to axisymmetric symmetry, simplifies to the form
[TABLE]
These four terms represent the forces per unit volume acting at that conic apex, which include (1) acceleration, (2) advection, (3) pressure gradient and (4) viscous forces. The double logarithmic plots shown in Fig. 7 (e)-(h) confirm self-similar growth of the conic tip as . As before, we do not report in Table 2 any exponent values for the dynamic range due to the presence of the sharp dip in the viscous stress within the self-similar regime (e.g. Fig. 7(g)) caused by a change in sign.
The results in Fig. 7(e) confirm that at low Reynolds number such as , the dominant competition is given by the pressure gradient and the viscous term, both roughly five orders of magnitude larger than the inertial terms given by Terms 1 and 2. This difference in magnitude between Terms 3-4 and 1-2 is observed to increase as . As Re increases, however, the inertial Terms 1 and 2 become more prominent, eventually approaching similar magnitude as the pressure gradient. By contrast, the viscous term decreases in magnitude, becoming rather negligible by about (not shown). For , it is evident from Fig. 7(h) that the viscous normal stress at the conic tip is at least two orders of magnitude smaller than the other three terms in the Navier-Stokes equation. In all the simulations conducted, the magnitude of the pressure gradient given by Term 3 exceeded that of the other terms throughout the period governed by self-similar growth. This likely reflects the fact that the rapidly increasing Maxwell pressure at the conic apex remains the dominant component in the pressure gradient irrespective of Reynolds number.
For completeness, we note here the changes in sign which occur in time for Terms 1 - 4 in Eq. (54) shown in Fig. 7(e)-(h). Term 1, which tracks the acceleration of the conic apex, is always positive irrespective of the value of Re. This must be the case since the dominant driving force at the tip is the Maxwell pressure which increases rapidly in time due to the increase in local electric field strength as the tip curvature increases in time. Term 2 tracks the inertial advective term whose sign depends on Re and . For , the actual sign of Term 2 is negative, which implies that at the apex. For Re = 50 and 500, there occurs a sharp dip in Term 2, which signals a transition from positive to negative values as decreases. For , there initially occurs significant scatter in Term 2 which eventually settles down to reveal steady power law growth. During this latter period, Term 2 is always positive in value.
Except for some scatter at early times, Term 3 given by the pressure gradient is always negative irrespective of the Reynolds number, which indicates that the driving pressure gradient . This behavior is expected since the dominant driving force due to the Maxwell pressure at the conic apex always pulls fluid toward the counter electrode. The viscous force per unit volume given by Term 4 is well behaved at small Reynolds number and always negative, counterbalancing the pressure gradient term. The data for Term 4 become significantly noisier with increasing Re while its magnitude decreases significantly, especially within the period of self-similar growth. This behavior can be traced to the rapid change in the vertical velocity component along the central axis within a distance of a few mesh elements of the accelerating interface. For and 500, this noisy regime occurs prior to the time interval used to extract the corresponding exponent. For , this noisy behavior persists well into the self-similar regime. For this reason, the exponents reported for Term 4 for were extracted using only data generated during the final decade in time where there occurs smooth power law growth.
IV.6 Self-similar collapse of conic tip shape
The time sequence images in Fig. 7 (i-l) confirm self-similar collapse of interface shapes obtained by rescaling the horizontal and vertical coordinates by the factor , where is the exponent for the capillary pressure growth plotted in Fig. 7 (a)-(d). This dynamic length was chosen to be the rescaling length since it is a direct measure of the diminishing radius of curvature at the conic apex. The curves in Fig. 7 (i-l) depict the evolving interface shapes during the period of self-similar growth identified in (a) - (h). The non-uniform spacing between curves reflects the adaptive time stepping scheme used in the simulations, as discussed in Section III.3.
IV.7 Normal stress exponents at conic apex vs Re
The simulations conducted in this study, whose range in Reynolds number spans the viscous to the inertial regime, confirm self-similar growth of the conic tip irrespective of Reynolds number. As clear from Fig. 7 (a) - (h), which pairs of forces predominantly contribute to the self-similar blow up behavior changes as Re increases. Figure 8 illustrates the variation in the magnitude of the blow up exponents , and for increasing Re at computed from the normal stress boundary condition given by Eq. (52). The dashed horizontal line represents the exponent values first predicted by Zubarev Zubarev (2001) for cone formation in the inviscid limit as . Error bars, which were always smaller than , are not shown since they would appear smaller than the icons used to designate the mean value. Included in Fig. 8 are the results for 5, 10, 25, 50 and 100 obtained from the two meshing schemes (A and B) discussed in Section III.3. These data are nearly indistinguishable, confirming that the two meshing schemes produced virtually identical results.
Fig. 8 illustrates which normal stresses at the moving interface maintain similar exponents. At the smallest values of Re, the Maxwell and viscous normal stresses (i.e. forces per unit area) assume a value close to one, although as indicated in Fig. 7(a), their actual magnitudes differ by more than an order of magnitude. As perhaps more evident from Fig. 7(e) which shows the resulting forces per unit volume acting at the tip, the Maxwell and viscous forces do provide the main competitive balance within the viscous dominated regime, while the capillary force per unit volume acts more as a smaller corrective term. By contrast, at much higher values of the Reynolds number, the dominant competition switches to that between the Maxwell and capillary stress while the viscous normal stress appears as a smaller correction. These results confirm the dominant pair forces reported by Fontelos et al. in the Stokes limit Fontelos et al. (2008) and Zubarev Zubarev (2001) in the inviscid limit.
IV.8 Comparison with exponents predicted by Zubarev and Fontelos et al.
The results in Fig. 8 for the exponents and indicate that as the Reynolds number increases toward the inviscid limit, both asymptote toward the value 2/3 first predicted by Zubarev Zubarev (2001). Zubarev’s scalings in Eq. (13)-(16) can also be used to estimate the corresponding asymptotic exponent for the viscous normal stress. Noting that the viscous normal stress in Eq. (52) at the apex is given by , it must therefore scale as . Our simulations confirm that at , the magnitude of this exponent has already reached a value , close to the asymptotic limit of one. The trend evident from Fig. 8 is that the magnitude of will further smoothly increase toward unity as .
Zubarev’s self-similar analysis requires that lengths scale as , velocities as and the capillary and Maxwell pressures as . In concluding the analysis in his original work, he applied these inviscid scalings to the Navier-Stokes equation given by Eq. (54) to determine what was the resulting scaling for the viscous stress and found that
[TABLE]
Inspection of the exponent values in Table 2 for from our simulations, namely , and , confirms close agreement with his inviscid prediction of . The value for the viscous stress is somewhat smaller than his asymptotic estimate of .
It is worth noting here that although Zubarev’s findings stemmed from a self-similar analysis of Bernoulli’s equation under the assumption of inviscid and irrotational flow, the same exponents can be derived from the Euler equation (i.e. the inviscid limit of the Navier-Stokes equation) and the corresponding inviscid limit of the normal stress boundary condition. A back of the envelope calculation leads to the prediction
[TABLE]
Apparently then, the irrotational constraint inherent in the Bernoulli analysis is not a requirement for deriving these exponents. Indeed, in all our studies we have observed the presence of vorticity just below the accelerating interface. In a separate work, we are investigating the behavior of this interface induced vorticity, whose presence nonetheless preserves his original predictions for these exponents.
The results of our simulations also show good agreement with the predictions of Fontelos et al. Fontelos et al. (2008) in the Stokes flow limit as . As discussed in Section II.4, Fontelos et al. demonstrated that at , the electric field strength at the apex of an electrified liquid drop must diverge as , or equivalently, that the Maxwell pressure must scale as , which leads to a value . Our simulations at yield exponent values , in excellent agreement with the asymptotic prediction. Their analysis also predicts that the exponent associated with the capillary pressure must vary with the cone interior half-angle such that . In Ref. [Fontelos et al. (2008)], the authors provide a numerical plot showing the variation in exponent values corresponding to the flow velocity (and electric potential) as a function of increasing cone interior half angle. From that plot and the value of the interior half angle extracted from our Fig.3 (e) for , we determine that . Our estimated value of for , which from the trend shown in Fig.8 will clearly increase as the Reynolds number is lowered still, is in good agreement with their asymptotic prediction that as .
Noting once again that the viscous normal stress in Eq. (52) at the apex is given by , it must then be the case from the scalings for and that the viscous normal stress scales as , the same scaling obtained in the limit of high Re. Our simulations at yield a value , in excellent agreement with this prediction. In fact, the dimensional analysis requires that regardless of Reynolds number. Inspection of the values extracted for given in Table 2 indeed confirm that all the reported values are close to one. The largest deviations from unity occur for those intermediate values of Reynolds number, which as we noted earlier, incur a change in sign within the fitting interval and a corresponding sharp dip, as shown in Fig. 7 (f) - (h). The occurrence of these dips tends to skew the extracted values of slightly away from one.
The analysis by Fontelos et al. Fontelos et al. (2008) for self-similarity in the Stokes regime also predicts that scales in time as (where ), as and the voltage potential as . The latter relation therefore requires that the electric field strength diverge as and therefore yields that the Maxwell pressure must diverge as . Substituting these scalings into Eq. (54) and the value we determined above for the runs conducted at yields
[TABLE]
Listed in Table 2 are the values of the corresponding exponents extracted from the simulations run at , namely , , and , which are all consistent with the asymptotic predictions by Fontelos et al. for .
V Discussion
V.1 Laminar flow conditions
The parameter range investigated for which at was confirmed to maintain laminar flow conditions throughout. While the Reynolds number defined in Table 1 is useful in cataloguing the general flow behavior discussed, it’s value is based on the initial electric field strength . In order to assess flow conditions in the region closest to the rapidly accelerating tip, it is also useful to define a local value of the Reynolds number given by
[TABLE]
where denotes the dimensional vertical component of the fluid velocity at the conic apex and the apex radius of curvature. Non-dimensionalizing these local parameters by and defined in Table 1 yields
[TABLE]
Our simulations and analysis indicate that at finite Reynolds number Re, the capillary pressure at the apex scales as . Since this pressure is inversely proportional to the (dimensionless) apex radius of curvature , this implies that . Furthermore, since the viscous normal stress at the apex scales as , then . These scalings therefore suggest that
[TABLE]
For the exponent values and listed in Table 2, it is clear that and therefore decreases as decreases, as confirmed by the curve shown in Fig. 9(a) and (b). In both plots, decreases monotonically as the collapse time is approached due to the sharpening of the conic tip while increases monotonically as the tip continues to accelerate. Although similar to the behavior in Fig. 9(a), the data in Fig. 9(b) indicate that the values for increase by approximately two orders of magnitude while Re increases by only two orders of magnitude (from 500 to 50,000). This correlation is clear from the product in Fig. (c) which remains fairly constant as is increased except at the smallest values of Re in the viscous dominated regime.
We note that in Fig. 9(a), the since the maximum value , which is achieved at (i.e. ), the flow at the tip maintains laminar conditions. In Fig.9(b), the maximum value , which is achieved at (i.e. ), also confirms laminar flow. The transition to turbulence for external flow typically sets in at Reynolds numbers of about . Interestingly, the maxima in the local Reynolds number in Fig. 9(a) and (b) both establish very early in the liquid deformation process, prior to the onset of self-similar growth. The early rapid increase in followed by a monotonic decrease is due to the fact that the decrease in with time is always more rapid than the corresponding increase in .
V.2 Correlations at conic apex with local Reynolds number
Shown in Fig. 9(c) is the maximum apex Reynolds number achieved during each run for at defined by . The quantities and denote the corresponding values of the tip curvature radius and tip speed. The function increases smoothly and monotonically with increasing Re such that for . It is evident from Fig. 9(c) that at intermediate and large values of Re, where the dominant opposing stresses are the Maxwell and capillary stress, the velocity and curvature radius track each other closely in magnitude. This is expected since the Maxwell pressure is always the dominant driving force and its magnitude is set by the apex curvature radius. We note that both and are rather independent of Re (except in the viscous dominated regime at low Re). This behavior confirms that these two values are set by local conditions at the apex, where the electric field strength and curvature radius undergo self-enhancement due to local conditions. This stands in contrast to the behavior at very small values of the Reynolds number , where the dominant opposing stresses are the Maxwell and viscous normal stress. The latter introduces longer range interactions mediated by the liquid viscosity which therefore causes a dependence on Re. Close inspection of the the curve for in Fig. 9(c) reveals a change in slope at a value of about . This break in slope reflects the viscous dominated regime to the left where and the inertia dominated regime to the right where .
V.3 Measurement restrictions involving liquid film thickness and electric field strength
V.3.1 Restrictions due to Bond number
In this study, the capillary number was held fixed at a value and gravity effects neglected such that the Bond number . The simulations were also terminated once the dimensionless apex radius of curvature (normalized by ) attained the value . These constraints enforce certain requirements on physical systems that may someday be used to verify the results obtained. Furthermore, it is assumed that electric field values at the apex should never exceed field ion emission values in liquid metals typically of order V/m Forbes and Mair (2008), though the actual values depend not only on the liquid metal but also the operating temperature and vacuum conditions.
V.3.2 Restrictions due to maximum field strength
To ascertain the parameter range of validity for the simulations in this study, we provide in Table 3 a list of derived values for three fluids at their melting temperature commonly used in liquid metal ion source devices, namely gallium (Ga), indium (In) and cesium (Cs). The entries were computed as follows. From the scalings in Table 1, a fixed value of Re for a given liquid system corresponds to a specific initial voltage potential . The constraint then establishes the specific thickness of the initial liquid layer, which then determines the initial field strength value . The material constants and the value then establish the corresponding value of the Bond number. For ground based experimental studies, it would be necessary to satisfy the inequality when making comparison to the findings of our study. Of course, one can always redo our numerical study using a larger value of Ca, which for a given value of would ensure a smaller value of and therefore a smaller Bond number. Too large a value of Ca, however, can trigger surface waves leading to multiple protrusions and interference effects.
Another important consideration for experimentalists is the small time window for measurement before onset of ion emission, a phenomenon that has not yet been incorporated into numerical models designed to probe conic formation. The difficulty in doing so stems from the very disparate length and time scales characterizing the hydrodynamic flow and field evaporation or related process. The initial field strength value in simulations should therefore not exceed field emission values. Shown in Table 3 are values for the initial electric field strengths for the specified values of Re and Ca used in our simulations as determined from the definitions in Table 1. The values of for the smallest Reynolds number simulated in our study, namely , would likely exceed the limit for ion emission and therefore preclude measurement; the corresponding values of are also unrealistically small. For the material constants listed, measurements at very low Reynolds numbers are likely not possible even at different values of Ca. According to the expressions for Re and Ca given in Table 1, a fixed value of Re corresponds to a fixed value of the voltage potential . Therefore smaller values of Ca correspond to larger values of and ultimately smaller values of the field strength . For , the value would have to be reduced by about two orders of magnitude to forestall ion emission, which could be accomplished by increasing by two orders of magnitude. However, that would still lead to unrealistically small values of , in the range of 0.1 nm. It is therefore unlikely that self-similar behavior in the Stokes regime can be tested using the system depicted in Fig. 2.
VI Recent experimental studies of dynamic cone formation in electrified fluids
Due to the very rapid flow speeds involved in dynamic cone formation in typical electrospray and liquid metal ion sources, some experimentalists have tried to slow the formation process by using partially conducting fluids beneath a layer of insulating liquid. This has allowed direct optical tracking and quantification of cuspidal acceleration undergoing self-similar growth. Oddershede and Nagel Oddershede and Nagel (2000) used high speed photography to track the evolution of a sharp liquid spout formed in response to a large normal electric field applied to a layer of distilled water containing NaCl beneath a layer of silicone oil. Above a critical field strength, the water generated a cuspidal shape whose apex curvature and overall height were shown to exhibit divergent self-similar growth extending almost three decades in time. The authors indicated that close to the collapse time, the shape of the conducting liquid tip was not exactly a cone but instead more closely approximated a square-root cusp. Simulations of the type described in this work could be expanded to examine such systems by replacing the perfectly conducting liquid with a partially conducting one and replacing the vacuum layer with a viscous liquid layer. These changes, of course, will modify the requisite boundary conditions at the moving interface.
More recently, Elele et al. Elele et al. (2015) have been successful in monitoring the field driven deformation of submillimeter size sessile droplets into conical shapes under ambient conditions by limiting the buildup of electric charge along the free surface. This allowed a sufficiently long period of distortion and growth before the onset of electrospray, which was achieved by separating the liquid droplet from the ground electrode with an electrically insulating film that prevented electric current from flowing between opposing electrodes. They examined the behavior of a number of partially conducting and leaky dielectric liquids in both ground based and space based studies. Their microgravity experiments carried out in the International Space Station allowed measurements on even larger droplet sizes in the centimeter range. Elele et al. confirmed the growth of conic tip shapes exhibiting a wide range of internal vertex angles dependent on operating conditions and self-similar behavior in the gap distance between the liquid apex and counter electrode, found to scale roughly as . Although this behavior is suggestive of Zubarev’s original prediction given by Eq. (13), it is important to note that the liquids used in the experiments by Elele et al. were not perfectly conducting nor the flow regime controlled by inviscid conditions, both assumptions critical to Zubarev’s findings.
VII Conclusion
The numerical simulations described in this work demonstrate how the free surface of an electrified perfectly conducting viscous liquid like a liquid metal undergoes deformation to a rapidly accelerating cusp with a conic tip for Reynolds numbers spanning the viscous to inertial regime. Field enhancement caused by the continuous sharpening of the conic apex causes divergent power law growth in finite time. Temporal tracking of the interface shape, normal stresses and forces per unit volume at the conic tip reveals persistent self-similar growth upon approach to the collapse time. The associated power laws exponents cover a dynamic range spanning six orders of magnitude in the Reynolds number. The asymptotic values of these exponents are in excellent agreement with previous predictions in the literature for the distinguished limits (Stokes regime) and (inviscid regime). The simulation results support the conclusion that at very small Reynolds number, the dominant balance at the tip as represented by the forces per unit volume acting on the conic apex represents the competition between the Maxwell and viscous forces while at intermediate and larger Reynolds number, the dominant balance switches to the competition between Maxwell and capillary stresses.
The temporal behavior of the terms in the normal stress boundary condition at the moving interface indicate that for Reynolds numbers higher than about 50, the viscous normal stress is rather insignificant in comparison to the Maxwell and capillary pressure. Similar tracking of the terms in the Navier-Stokes equation, however, reveals that viscous normal forces per unit volume acting at the conic apex cannot be neglected for values Re below about since they remain comparable in magnitude to the apical pressure gradient and apical inertial forces. Previous studies of dynamic cone formation based on Bernoulli’s equation have tended to neglect altogether the role of viscous stresses at the cone tip.
Zubarev’s Zubarev (2001) original prediction for self-similar growth in the inviscid limit led to power law exponents characterized by a ratio of two integers. In particular, he reported that the capillary pressure, the Maxwell pressure and the kinetic energy per unit volume all scaled as . By contrast, Fontelos et al. Fontelos et al. (2008) uncovered that in the Stokes regime the corresponding exponents are not simple integer ratios but instead depend functionally on the cone interior half-angle. Our study, which reveals the range and behavior of exponents characterizing the forces at the tip, confirm that the exponents at finite Re are never simple integer values. These non-integer ratios reflect the fact that the forces operating at the tip depend sensitively on the slope of the dynamic cone as the system approaches the collapse time.
In conclusion, the study presented here is based on examination of the forces operating at the conic tip and the resultant interface shapes which are shown to depend on the internal cone angle. We are currently expanding our study to evaluate the flow behavior throughout the tip region to better understand the characteristics of the interfacial boundary layer which accompanies the movement of the rapidly accelerating and highly curved interface. Such a boundary layer is expected to affect not only the stability of this interfacial region but ultimately the field emission process known to emanate from the conic apex.
Acknowledgements.
The authors gratefully acknowledge financial support from a 2014 NASA Space Technology Research Fellowship (TGA) and the NASA/Jet Propulsion Laboratory (JPL) President’s and Director’s Fund (SMT). The authors also wish to acknowledge Chengzhe Zhou for interesting discussions during the course of this work about Taylor cone formation in the inviscid limit.
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1Larmor (1890) J. Larmor, “On the influence of electrification on ripples,” Proc. Cambridge Phil. Soc. 7 , 69–71 (1890).
- 2Tonks (1935) L. Tonks, “A theory of liquid surface rupture by a uniform electric field,” Phys. Rev 48 , 562–568 (1935).
- 3Frenkel (1935) J. Frenkel, “On ‘Tonks’ theory of liquid surface rupture by a uniform electric field,” Phys. Z. Sowjetunion 8 , 675 – 683 (1935).
- 4Frenkel (1936) Ya. I. Frenkel, “About the Tonks theory of liquid surface rupture by a uniform electric field in vacuum,” Zh. Eksp. Teor. Fiz 6 , 347–50 (1936).
- 5Krohn (1984) V. E. Krohn, “Liquid metal droplets for heavy particle propulsion,” Prog. Astronautics Aeronautics 5 , 73–80 (1984).
- 6Perel et al. (1971) J. Perel, A. Y. Yahiku, J. F. Mahoney, H. L. Daley, and A. Sherman, “Operational characteristics of colloid thrusters,” J. Spacecraft Rockets 8 , 702–708 (1971).
- 7Bailey et al. (1972) A. G. Bailey, J. E. Bracher, and H. J. Von Rohden, “A capillary-fed annular colloid thruster,” J. Spacecraft Rockets 9 , 518–521 (1972).
- 8Zafran et al. (1973) S. Zafran, J.C. Beynon, P.W. Kidd, H. Shelton, and F. A Jackson, “One-mlb colloid thruster system development,” J. Spacecraft Rockets 10 , 531–533 (1973).