HIGH-FIDELITY SIMULATIONS OF A ROTARY BELL ATOMIZER
WITH ELECTROHYDRODYNAMIC EFFECTS
by
Venkata Krisshna Pydakula Narayanan
A dissertation submitted in partial fulfillment
of the requirements for the degree
of
Doctor of Philosophy
in
Mechanical Engineering
MONTANA STATE UNIVERSITY
Bozeman, Montana
July 2023
©c COPYRIGHT
by
Venkata Krisshna Pydakula Narayanan
2023
All Rights Reserved
ii
DEDICATION
Knowledge cultivates humility; humility instills worthiness,
From worthiness one earns wealth; with wealth one can perform deeds,
and so, one finds happiness.
I dedicate this dissertation to
my father, P T Narayanan, who always supports me to pursue my dreams,
my mother, K Kalpana, for all her effort and care to see me succeed,
all the teachers who inspire me to learn and remind me to stay curious,
and all the young minds who teach me to stay spirited, passionate and creative.
iii
ACKNOWLEDGEMENTS
I am immensely grateful to Dr. Mark Owkes for being an amazing mentor. Mark
has taken the time to explain complex concepts and guide me through challenging research
problems. I have hardly felt that Mark is unreachable even through the pandemic and his
year-long sabbatical. In addition to his academic guidance, Mark has cared for my well-
being in both personal and professional aspects. I especially thank Mark for providing
me with opportunities to attend conferences and present my work. I am grateful to Ford
Motor Company for their generous funding, which has provided me with the opportunity
me to undertake in an industrially driven research project. I would also like to express
my gratitude to the National Science Foundation for funding my research activity.
The material presented in this document is based upon work supported by the under NSF
Grant No. 1749779. Dr. Wanjiao Liu has been instrumental during my learning phases
and had offered me the opportunity to intern at Ford - an experience paramount to my
professional growth. I also thank Dr. Kevin Ellwood and Brandon Petrouskie for
their mentorship and guidance. Dr. Erick Johnson’s vast knowledge has encouraged me
to continue learning. Erick’s questions about my research progress are always thought-
provoking and have improved the accuracy and applicability of my research activities. Erick
does not stop with questions either - he proceeds to provide at least two paths to overcome
my challenges. Dr. Joseph Seymour’s passion for our understanding of fluids is infectious.
Joe has been appreciative and supportive of my work at all stages of research. Dr. Yaofa
Li has provided valuable feedback on my presentations. He has been instrumental in making
iv
ACKNOWLEDGEMENTS – CONTINUED
sure that my project serves as a valuable and comprehensive contribution to the scientific
community. Dr. Alan George taught me fluids in a perspective I had never seen before.
I lack sufficient words to thank her Dr. Sarah Codd who has genuinely cared for my
(and all other students’) welfare and academic success. Dr. Alan George taught me fluid
dynamics in a perspective I had never seen before. I would like to thank my labmates -
Kris, Brian, Ryan, Erin, Clark, Adam, Julian, and several others. I am grateful
for the countless discussions on matters ranging from the best paper clip shape to the rise
of artificial intelligence. Brendan has been a keystone in my growth as a scientist. Our
shared adventures both inside (various technical discussions and outreach events) and outside
the lab (rock hunts, fishing, cooking, etc.) have created unforgettable memories that I will
cherish for a long time. Coltran and Kenny from UIT have supported my high-performance
computing needs. Dr. Owkes, Dr. Johnson, George Schaible, Dr. Loribeth Evertz,
and Dr. Codd have helped me grow as a science communicator. Madhu and Toral have
joined me on many adventures. Toral has made me feel safe and loved through many seasons.
My two best friends and allies, Akshay and Aravindh, have been my beacons since we
met. My parents and family have supported me and have dedicated so much of their life to
me. I have learned to be righteous and sincere in my efforts from my father. I have learned
that it is important to be caring and perseverant from my mother. My parents have given
nothing short of their best to see me succeed and continue to guide me in making my journey
through life a wholesome experience.
v
TABLE OF CONTENTS
1. INTRODUCTION .................................................................................................. 1
1.1 Fluid Dynamics and the Use of Computer Simulations ......................................... 1
1.2 Gas-Liquid Multiphase Flows ............................................................................. 2
1.3 Atomization and Breakup .................................................................................. 4
1.4 Rotary Bell Atomizers........................................................................................ 4
1.5 Electrohydrodynamic Assisted Atomization ......................................................... 6
1.6 Literature Review: Experimental Studies ............................................................ 7
1.7 Literature Review: Numerical Studies ................................................................. 8
1.8 Motivation to Study Electrified Atomization in RBAs ........................................ 10
1.9 ERBA Problem Setup...................................................................................... 12
1.10 Contributions .................................................................................................. 12
1.11 Dissertation Layout ......................................................................................... 15
2. MATHEMATICAL FORMULATION .................................................................... 16
2.1 Governing Equations........................................................................................ 16
2.1.1 Navier Stokes Equations for Multiphase Flows........................................... 16
2.1.2 Electrohydrodynamics.............................................................................. 17
2.1.3 Rotating Reference Frame ........................................................................ 20
2.1.4 Non-Newtonian Viscosity Relationship ...................................................... 20
2.2 Validation of Physics Modules .......................................................................... 22
3. NUMERICAL METHODS .................................................................................... 25
3.1 Domain Setup ................................................................................................. 26
3.2 The Need for e-Mesh........................................................................................ 30
3.2.1 e-Mesh Framework................................................................................... 30
3.2.2 e-Mesh Domain Initialization.................................................................... 31
3.2.3 e-Mesh Boundary Conditions.................................................................... 33
3.2.4 Migration Charges to e-Mesh (NS→e Migration) ....................................... 35
3.2.5 Solving the Electric Poisson Equation on e-Mesh ....................................... 39
3.2.6 Migrating Electric Potential to NS-Mesh (e→NS Migration)....................... 40
3.2.7 Challenges and Cost With Using E-mesh................................................... 43
3.3 Mesh Sensitivity .............................................................................................. 45
3.3.1 Sensitivity of Droplet Sizes to Mesh Resolution ......................................... 46
3.3.2 Sensitivity of Droplet Charge Densities to Mesh Resolution ........................ 48
3.4 Droplet Ancestry Extraction ............................................................................ 49
3.4.1 Insights From Droplet Statistics................................................................ 49
3.4.2 Method of Extraction............................................................................... 49
vi
TABLE OF CONTENTS – CONTINUED
3.4.3 Storing Statistics in a Database ................................................................ 50
3.4.4 Post-processing the Statistics.................................................................... 51
3.5 Computational Cost......................................................................................... 53
4. DEMONSTRATION OF NUMERICAL FRAMEWORK......................................... 55
4.1 Liquid Interface ............................................................................................... 55
4.2 Electric Potential Field .................................................................................... 58
4.3 Parameter Study ............................................................................................. 58
4.4 Examining Non-Dimensional Quantities ............................................................ 61
5. EFFECT OF ELECTRIFYING THE SETUP ........................................................ 63
5.1 Liquid Interface ............................................................................................... 63
5.2 Spray Characteristics ....................................................................................... 67
5.3 Atomization Characteristics ............................................................................. 68
5.4 Droplet Size and Shape Evolution..................................................................... 71
5.5 Droplet Velocity Statistics ................................................................................ 74
5.6 Local Weber Number Analysis.......................................................................... 75
5.7 Charge Distribution Statistics........................................................................... 77
6. EFFECT OF SHEAR-THINNING FLOWS............................................................ 80
6.1 Liquid Interface ............................................................................................... 80
6.2 Atomization Characteristics ............................................................................. 82
6.3 Droplet Size Evolution ..................................................................................... 84
6.4 Viscosity Characteristics .................................................................................. 86
7. CONCLUSIONS AND OUTLOOK ........................................................................ 90
7.1 Atomization Characteristics in Rotary Atomization ........................................... 91
7.2 Impact of Research Work ................................................................................. 92
7.3 Future Work.................................................................................................... 92
REFERENCES CITED.............................................................................................. 96
APPENDIX .............................................................................................................112
APPENDIX A : Appendix ....................................................................................113
A.1 Liquid Surface Profile in a Rotating Tank ........................................................113
A.2 Trajectory of an Object in a Rotating Frame....................................................113
vii
LIST OF TABLES
Table Page
3.1 Standard values of parameters in the ERBA simulations.............................. 29
3.2 Simplified commlist corresponding to the setup shown in Fig. 3.5................. 41
3.3 Mesh sensitivity simulation domain parameters ........................................... 46
4.1 Non-dimensional numbers used in the ERBA simulations ............................ 62
viii
LIST OF FIGURES
Figure Page
1.1 A rotary bell atomizer in operation. In this image, the
shape of the spray cloud is clearly visible above the target
surface. Photograph courtesy of RISE Research Institutes
of Sweden and Fraunhofer-Chalmers Centre. ................................................. 6
1.2 A schematic representation of the domain geometry show-
ing the trajectory of the liquid jet (green arrow) and the
external body forces acting upon it (not to scale). The
device is typically operated at a distance of about 0.25 m
from the target surface. ............................................................................. 14
2.1 The viscosity-shear rate relationship of automotive paint is
approximated well using a modified power law fit. Note the
logarithmic scale on the axes...................................................................... 22
2.2 Validation of centrifugal force: Numerical results of a
rotating tank experiment where cases 1, 2, and 3 correspond
to rotation rates of 4, 10, and 14 rad/s respectively. The
analytical solution of the interface (Eq. 2.18) at various
locations is plotted with red circles. The numerically
computed liquid interface and liquid bulk are shown as black
lines and blue areas respectively. The dashed line represents
the axis of rotation for each case. ............................................................... 23
2.3 Validation of Coriolis force: Numerical results of the liquid
jet interfaces in a reference frame rotating counter-clockwise
about an axis that points out of the 2D plane. The
analytical solution (Eq. 2.19) is shown as a superimposed
solid line for each case. .............................................................................. 24
3.1 A schematic representation of the angle associated with the
ejection of liquid from the bell edge, calculated between the
axis of rotation and the edge surface. ......................................................... 28
3.2 A schematic representation of e-Mesh and available φ BCs
shown in context with NS-Mesh (not to scale)............................................. 31
3.3 A schematic representation of the domain geometry of the
numerical problem showing positions and dimensions (not
to scale) of NS-Mesh and e-Mesh................................................................ 32
ix
LIST OF FIGURES – CONTINUED
Figure Page
3.4 Schematic representation of the grids on NS-Mesh (yellow)
and e-Mesh (blue) overlapping. e-Mesh is formulated with
a stretched grid to facilitate a higher concentration of cells
near the bell edge. The nozzle is shown for reference and is
not to scale. .............................................................................................. 33
3.5 Schematic representation of the overlap of processors be-
tween NS-Mesh (green) and e-Mesh (pink). For example,
overlapping information exists on processor 3 on the NS-
Mesh (NS processor 3, highlighted in bright green) and pro-
cessor 1 on the e-Mesh (e processor 1, highlighted in bright
pink). The nozzle is shown for reference and is not to scale.......................... 34
3.6 Locations of NS-Mesh boundaries (yellow) with respect to
a cell on e-Mesh e cell (blue) ....................................................................... 36
3.7 Locations of the boundaries of a cell on NS-Mesh NS cell
(yellow) and a cell on e-Mesh e cell (blue) .................................................... 37
3.8 Electric potential field (φe−Mesh) on e-Mesh. The nozzle is
shown for reference and is not to scale. ....................................................... 39
3.9 Boundaries of regions occupied by NS processor 3 (yellow) and
e processor 1 (blue). Here, e-NS processor 1−3
overlap = true. ......................................... 40
3.10 Probability distribution of droplet sizes for varying mesh
resolutions. The lines show a log-normal fit of the data. .............................. 47
3.11 Volume weighted charge density probability distributions
for varying mesh resolutions. Note the logarithmic scale on
the y-axis. ................................................................................................ 48
3.12 Droplet ancestries in the jet shown in Fig. 4.1 represented
using the Neo4j graphical database. The colors represent
breakup levels 1 (gray), 2 (yellow), 3 (green), 4 (red), 5
(pink) and 6 (blue) breakup events. Lines connect related
droplets, i.e., droplets which split from each other. The
ligament core (black) is present at the center of the image.
Each node contains relevant statistics from the breakup event...................... 52
x
LIST OF FIGURES – CONTINUED
Figure Page
3.13 Breakdown of the computational cost of notable modules
in the numerical solver. The electric module which handles
solutions to the EHD equations presented in Section 2.1.2
is the most expensive part of the solver taking up to 77%
of the cost; the transport module which handles multiphase
flows and interface tracking accounts for 10% of the cost;
the pressure solver accounts for 12% of the cost. ......................................... 54
4.1 Liquid interface positions of a single jet exiting the edge of
the bell at different times as viewed from an x-y plane. The
nozzle is shown for reference at the top left corner of each
image and is not to scale. .......................................................................... 56
4.2 Liquid interface positions of six individual jets exiting a
serrated bell edge at different times as viewed from an x-z
plane. The nozzle is rotating clockwise and is shown for
reference at the left edge of each image and is not to scale.
The simulation domain is outlined by a dashed box. Owing
to the periodic boundary conditions along the z-axis, an arc
along the bell edge spanning six serrations is depicted in
this view by stacking the domains in an attempt to vaguely
recreate some experimental images [42, 118, 153, 176].................................. 57
4.3 The electric potential field and its corresponding contour
lines on an x-y plane in the domain. The trajectory of
the liquid jet is shown graphically as a dashed white line.
The electric force vector (blue arrow) that acts on a control
volume (blue circle) in this regime is also schematically
represented. The bell is shown at the top left corner and is
not to scale. .............................................................................................. 58
4.4 A comparison of snapshots of the liquid interface positions
after 200 µs in the parameter study simulations. The nozzle
is shown for reference at the top left corner of each image
and is not to scale. .................................................................................... 60
5.1 Liquid interface positions in a non-electrified (yellow) and
an electrified (violet) jet after exiting the edge of the bell
on an x-y plane at different times. The nozzle is shown for
reference at the top left corner and is not to scale. ...................................... 65
xi
LIST OF FIGURES – CONTINUED
Figure Page
5.2 Liquid interface positions of six individual non-electrified
(yellow) and electrified (violet) jets exiting a clockwise-
rotating serrated bell edge (left, not to scale) as viewed at
different times from an x-z plane. The simulation domain is
outlined by a dashed box. Owing to the periodic boundary
conditions along the z-axis, an arc along the bell edge
spanning six serrations is depicted in this view by stacking
the domains in an attempt to recreate experimental images
available in Shirota et al. [153], Oswald et al. [118], Domnick
[42] and Wilson et al. [176]. ....................................................................... 66
5.3 Log-normal fit of the near-bell droplet spacing probabilty density. ................ 67
5.4 The total number of droplets in the domain over time. Note
the logarithmic scale on the y-axis.............................................................. 69
5.5 The number distribution of droplets that have undergone
different breakups levels............................................................................. 70
5.6 Droplet diameter distribution. Note the logarithmic scale
on the y-axis. ............................................................................................ 71
5.7 Mean droplet diameter for each breakup level. ............................................ 72
5.8 Mean change in droplet diameters for successive breakup
levels. Here, breakup level 1 is omitted since it corresponds
to breakup from the main ligament core. .................................................... 73
5.9 Droplet sphericity distribution. .................................................................. 73
5.10 Comparing the mean velocity for different droplet sizes. .............................. 74
5.11 Comparing the mean slip velocity for different droplet sizes. ........................ 75
5.12 Local Weber number distribution associated with all struc-
tures. Note the logarithmic scale on the y-axis. ........................................... 76
5.13 Number distribution of q in liquid structres in the domain.
The inital bulk value of q is shown as a dashed line. Note
the logarithmic scale on the y-axis.............................................................. 77
5.14 Mean q in liquid structures of different sizes. The inital
bulk value of q is shown as a dashed line..................................................... 78
xii
LIST OF FIGURES – CONTINUED
Figure Page
5.15 The mean change in q between successive breakup levels. ............................ 79
5.16 The mean change in q corresponding to the change in
diameter during a breakup event. ............................................................... 79
6.1 Liquid interface positions of a Newtonian (pink) and a non-
Newtonian (green) jet at different times as viewed from an
x-y plane. The nozzle is shown for reference at the top left
corner and is not to scale. .......................................................................... 81
6.2 The total number of droplets in the domain over time in
Newtonian (pink) and non-Newtonian (green) operating
conditions in the first 300 µs of the near-bell atomization
simulation. Note the logarithmic scale on the y-axis. ................................... 82
6.3 The number distribution of droplets that have undergone
different breakups levels in Newtonian (pink) and non-
Newtonian (green) operating conditions in the first 300 µs
of the near-bell atomization simulation. ...................................................... 83
6.4 Droplet diameter distribution in the simulation. .......................................... 84
6.5 Mean droplet diameter for each breakup level in the simulation. .................. 85
6.6 Mean change in droplet diameters for successive breakup
levels. Here, breakup level 1 is omitted since it corresponds
to breakup from the main ligament core. .................................................... 85
6.7 Liquid viscosity distribution during breakup events in the
first 300 µs of the simulation...................................................................... 87
6.8 Mean viscosity of the parent droplets for each breakup level. ....................... 87
6.9 Mean viscosity of the parent droplets as a function of
droplet sizes.............................................................................................. 88
6.10 Comparing liquid viscosity for a range of local Weber
numbers during breakup events in the first 300 µs of the
simulation................................................................................................. 88
6.11 The mean change in liquid viscosity between successive
breakup levels. .......................................................................................... 89
xiii
LIST OF FIGURES – CONTINUED
Figure Page
A.1 Analytical solution of the liquid surface (red) in a tank
rotating about the y axis with an angular velocity ω
described mathematically by the function y ...............................................113
xiv
LIST OF ALGORITHMS
Algorithm Page
3.1 NS→e charge density migration algorithm .................................................. 38
3.2 Calculating index extents for overlap between NS-Mesh and
every e processor n ........................................................................................ 41
3.3 Calculating index extents for overlap between e-Mesh and
every NS processor n...................................................................................... 42
3.4 e→NS electric potential field migration algorithm ....................................... 44
xv
NOMENCLATURE
t Time
u Velocity vector
p Hydrodynamic pressure
g Gravitational acceleration
ρ Density
γ Surface tension coefficient
κ Local interface curvature
µ Dynamic viscosity
ε Electric permittivity
ω Angular velocity
E Electric field vector
f Force vector
n Interfacial normal vector
r Radial distance from rotation axis
φ Electric potential
q Volumetric charge density
J Current density
D Diffusion constant
ψ Charge mobility
I Identity tensor
Υ Viscosity coefficients
η Local shear rate
Djet Jet diameter
Ree Electric Reynolds number
Pee Electric Péclet number
σf
i Viscous stress tensor
σe
i Maxwell stress tensor
δs Interfacial Dirac-delta function
τq Charge relaxation timescale
τf Fluid advection timescale
VF Liquid volume fraction
S Structure identification number
L Liquid identification number
We local Local Weber number
xvi
NOMENCLATURE – CONTINUED
BC Boundary condition
CFD Computational fluid dynamics
CPD Cells per diameter
EHD Electrohydrodynamics
ERBA Electrostatic rotary bell atomizer
MPI Message passing interface
PDE Partial differential equation
RBA Rotary bell atomizer
TE Transfer efficiency
VOF Volume of fluid
ACES Adjustable Curvature Evaluation Scale
RPM Rotations per minute
xvii
ABSTRACT
Atomizing flows involve the breakup of a liquid into a spray of droplets. These flows
play a vital role in various industrial applications such as spray painting and fuel injection.
In particular, these processes can have severe impacts especially in automotive paint shops -
which can account for up to 50% of the total costs and 80% of the environmental concerns in
an automobile manufacturing facility. A device commonly used for painting vehicles is called
an electrostatic rotary bell atomizer (ERBA). ERBAs rotate at high speeds while electrically
charging the liquid and operating in a background electric field to direct atomized charged
droplets towards the target surface. The atomization process directly influences the transfer
efficiency (TE) and surface finish quality. Optimal spray parameters used in industry are
often obtained from expensive trial-and-error methods.
To overcome these limitations, a computational tool has been developed to simulate
three-dimensional near-bell ERBA atomization using a high-fidelity volume-of-fluid transport
scheme. Additionally, the solver is equipped with physics modules including centrifugal,
Coriolis, electrohydrodynamic (EHD), and shear-thinning viscous force models. The primary
objective of this research is to investigate the influence of EHD parameters on near-bell
atomization of paint and subsequently improve TE in ERBAs in a cost-effective manner.
Using the tools developed, numerical simulations are performed to understand the physics
of electrically assisted atomization.
The influence of various operating parameters, such as liquid flow rate, bell rotation
rate, liquid charge density, and bell electric potential, on atomization is examined. Results
from a comparative study indicate that the electric field accelerates breakup processes and
enhances secondary atomization. The droplet velocity, local Weber number and charge
density statistics are also analyzed to understand the complex physics in electrically assisted
breakup. Additionally, the effect of shear-thinning behavior of the liquid on atomization is
also explored. High-fidelity simulations allow for the extraction of breakup statistics which
are otherwise challenging to obtain experimentally.
These findings have the potential to drive improvements in the design and operation
of ERBAs, leading to enhanced TE and surface finish quality while reducing costs and
environmental concerns in automotive paint shops.
1
CHAPTER ONE
INTRODUCTION
1.1 Fluid Dynamics and the Use of Computer Simulations
Fluid dynamics is the science describing liquids and gases in motion. This field of
science helps us understand a wide range of common and complex scientific phenomena.
Some of the pioneering descriptions of how liquids behave can be traced back to ancient
Greece around the 3rd century BCE when Archimedes formulated principles related to fluid
buoyancy and the displacement of fluids. For centuries, we have described and continually
explored the principles governing fluid flow, including the conservation of mass, momentum,
and energy, as well as the effects of viscosity and turbulence. Over the past few decades, we
have seen a rise in the usage of computers in exploring fluid physics through the development
of computational fluid dynamics (CFD) [5, 11]. CFD utilizes the power of modern computers
to investigate fluid flow problems. CFD is a powerful tool that can accurately model physical
processes using suitable numerical methods and algorithms to solve the corresponding set of
equations governing the phenomenon. CFD has allowed us to simulate and analyze complex
fluid flows that would otherwise be challenging or impossible to study through traditional
experimental methods. In fact, the findings of this dissertation are arguably impossible
to obtain using current state-of-the-art experimental techniques. I have personally seen the
application of CFD tools in a wide range of scientific and industrial fields including aerospace,
chemical, agricultural, energy, medical, environmental, and automotive [99, 108, 142, 150,
178] sectors.
In computer simulations, we have the facility to model difficult geometries, impose
2
complex boundary conditions, test various fluid properties, and play with a wide range
of operating conditions. One of the many strengths of a CFD engineer lies in the
process of discretizing the complicated partial differential equations (PDEs) that provide
a mathematical description of the physics. Discretization of the PDEs is a procedure that
yields a set of (somewhat) simple algebraic equations consisting only of the basic math
operations (addition, subtraction, multiplication and division). Another interesting step in
CFD is the decomposition of the continuous flow domain into smaller discrete control volumes
called computational cells. All properties of the fluid, i.e., scalar quantities such as pressure,
charge, and temperature and vector quantities such as velocity, vorticity, and electric field,
and tensor quantities, such as the stress tensor are represented as discrete values at either
the center or faces of these computational cells.
The process of discretizing the governing equations and the continuous flow domain to
yield a set of algebraic equations that can be solved iteratively in each computational cell
using numerical methods is the core essence of every CFD algorithm. CFD enables obtaining
impeccable insights into the behavior of fluids while, most crucially, reducing costs and time
involved with expensive experimental approaches to studying the science of fluid dynamics.
1.2 Gas-Liquid Multiphase Flows
Any fluid flow consisting of more than one phase, such as solid-liquid, gas-liquid,
and solid-gas flows, within a system is what characterizes multiphase flows - a field of
fluid dynamics that has been well studied in the past [17, 19, 26]. Multiphase flows are
ubiquitous in our natural world with examples ranging from sediment transport (solid-
liquid) [93], bubbly flows (gas-liquid) [98], and dust suspension (solid-gas) [103]. CFD
algorithms have also seen continued development towards modeling complex multiphase
flows [29, 85, 132, 169, 177]. For example, gas-liquid multiphase CFD models, much like
the one used in this project, account for the presence of two or more immiscible fluids, such
3
as a gas and a liquid, within the computational domain [4]. Each phase can be modeled
separately, with their individual properties and interactions considered. Some features of
such multiphase solvers include the ability to track the interface to accurately simulate
interface dynamics and estimate the local interface curvature to compute the surface tension
force. Some multiphase models might also include framework for mass transfer and phase
interactions [21, 38]. Several interface tracking methods have been developed [28, 83].
The Volume of Fluid (VOF) method, which has been employed in the project detailed
in this dissertation, is based on a Marker-and-cell method, tracks interface motion by
storing the volume fraction of each fluid within computational cells [72, 113]. Although
this method is mass-conserving, the interface is not readily available and must be computed
using appropriate reconstruction methods, which can be complex. In a more fundamental
aspect, the modeling of a sharp interface can be challenging and erroneous if its geometric
and topological characteristics are not treated carefully [138].
The level set method represents the interface between fluids as an iso-surface of a scalar
function that evolves over time to track interface motion [116, 149]. This method is popular
due to its simplicity compared to other methods. However, scientists the method has accrued
complexities introduced by researchers aiming to address its problems with mass conservation
[117].
The front-tracking method explicitly represents interfaces using marker particles and
maintains the topology of the interface during the simulation [130, 170]. While the interface
normal and curvature are well predicted by the method, the process of fluid coalescence
or breakup, which are crucial to some applications of multiphase flows, are not modeled
elegantly [25].
4
1.3 Atomization and Breakup
The process of disintegration of bulk liquid into a cloud of droplets/particles, i.e., sprays,
in a gaseous medium is termed as atomization. We observe a wide range of sprays in our
daily lives including rain, garden sprinklers, shower-heads, hair sprays, and household cleaner
sprays. Atomization is a key area of research and development, aimed at improving process
efficiency, environmental impact, and enhancing product quality in various industries. Sprays
are commonly seen in agricultural [70], healthcare [168], automotive [47], painting [35], food
processing [13], cooling [10], fire suppression [136], power generation [15], consumer sprays
[173], snowmaking [95], and several other industrial applications [111].
The atomization process is influenced by a wide range of factors - fluid properties
(viscosity, surface tension, density), nozzle geometry and design, and operating conditions
(ambient temperature, pressure, electric field). Understanding and controlling these factors
is essential for achieving the desired droplet size distribution and spray characteristics. CFD
models have been developed to simulate atomizing flows [94, 147]. These models aim to
provide insights into the breakup and dispersion of liquid jets or films into droplets, which
is crucial for the various applications listed above.
1.4 Rotary Bell Atomizers
Sprays can be produced in various ways. The breakup of a liquid jet or sheet during
the atomization process is often achieved by the impartation of kinetic energy to the liquid.
An atomizer is a device used to convert a liquid into a spray of droplets. In general, most
atomizers utilize the high relative velocity between the liquid and the ambient gas to induce
breakup activity. Liquid is either
1. discharged at a high velocity into a relatively slow-moving stream of gas (pressure [91]
or rotary [41] atomizers) or
5
2. discharged at a low velocity into a relatively fast-moving stream of gas (air-blast
atomizers [1]) or
3. disrupted using bubbles (effervescent atomization [158]) or mechanical (ultrasonic
atomizers [89]) means to disintegrate the bulk into smaller fragments or
4. manipulated using electric fields in a technique that I will discuss in the next section.
Rotary atomizers are popular in industrial applications due to their simplicity and ease
of operation and the desirable atomization characteristics [36, 100, 124] . One such device
is called the rotary bell atomizer (RBA) (Fig. 1.1) which is characterized by a bell shaped
nozzle [114]. Liquid is injected into the center of a rotating bell. Centrifugal forces cause the
liquid to accelerate and form a film that flows radially outward along the inner surface of
the bell [45, 88, 164]. The liquid film then exits the edge of the nozzle at a high velocity to
form ligaments or sheets which eventually undergo atomization to form a cloud of uniformly
sized droplets [58]. Notable desirable features in RBAs are the uniformity in droplet sizes
and the ease of control over spray pattern.
6
Figure 1.1: A rotary bell atomizer in operation. In this image, the shape of the spray cloud
is clearly visible above the target surface. Photograph courtesy of RISE Research Institutes
of Sweden and Fraunhofer-Chalmers Centre.
1.5 Electrohydrodynamic Assisted Atomization
The fourth method of achieving atomization that I mentioned in the previous section
is called electrospraying, where an electric field is used to disintegrate a liquid into charged
droplets [27]. These charged droplets can be controlled and manipulated using the electric
field [57, 154].
The science behind electrosprays is called electrohydrodynamics (EHD). EHD provides
a mathematical description of systems with a significant interaction of fluid dynamics and
electrostatics. EHD-assisted atomization has been an increasingly important means of
producing liquid droplets that is well-established to offer advantages over other industrial
spraying processes [18, 62, 66]. It is commonly used in applications such as inkjet
printing [133], mass spectrometry [56], biochemistry [161], microfluidics [16], food technology
[6, 14, 61, 82], electrostatic precipitation [179], powder coating [80], fuel injection [155, 182],
and biomedicine [46, 50, 54].
7
Numerical methods that model EHD flows have been developed and improved over
the past few decades for various application by several research groups [97, 167]. In EHD
simulations, charges can be modeled in the liquid phase on the interface as a fully relaxed
surface charge [63, 73, 74, 107] or in the liquid volume as a uniform volumetric charge
[154, 172]. Charges can also be modeled in high-fidelity and allowed to relax and redistribute
in the liquid phase using appropriate models for charge the migration process [151].
In industrial applications of RBAs, it is common practice to electrically charge the
liquid and apply a background electric field to enhance the transfer efficiency of the device
[43, 77]. The transfer efficiency, which is an important metric for spraying devices, is the
ratio of the mass of the liquid injected into the nozzle to the mass of the liquid that is coated
onto the target surface. Atomized droplets, which also carry electric charge, move toward
the grounded target surface under the influence of the electric field. An RBA operating in
an electrified setup is called an electrostatic rotary bell atomizer (ERBA).
1.6 Literature Review: Experimental Studies
The process of atomization in RBAs has been well studied in the past. The following
is a summary of a few notable experimental studies that have characterized RBA operation.
A study conducted by Hinze and Milborn [71] have found that the droplet size
distribution can be controlled by the rotation rate, nozzle geometry, and the liquid properties.
Studies conducted by Corbeels et al. [35] and Rezayat and Farshchi [137] found similar
results and included the effects of flow rate, nozzle size and spray angle. Another study
presented in Sidawi et al. [156] measured the fraction of area covered by spray droplets for
various rotation rates. Mahmoud and Youssef [101] used the phase Doppler particle analyzer
method to study the influence of the spinning cup and disk configurations on atomization
characteristics. The effect of shaping air, which is the use of a high-speed air stream to shape
the spray, on RBAs and ERBAs was investigated by Stevenin et al. [160] and Balachandran
8
and Bailey [20], respectively. The effect of bell geometry on droplet sizes was examined by
Domnick [42]. Sidawi et al. [157] investigated the effect of the size and shape of serrations
at the bell edge on atomization quality. Shen et al. [152] investigated the effect of the
shear-thinning viscosity of paint on the initial wetting, film formation, and film thickness
distribution on the bell. A study conducted by Loch et al. [96] identified the importance
of rotation rate and the background electric field to achieve the desired coating properties.
A simple computer model was developed from experimental measurements and observations
by Yao et al. [181] to predict droplet movement in the spray cloud of rotary atomizers. Fan
et al. [53] describes experimental and mathematical methods to measure and calculate paint
droplet size, velocity, and charge-to-mass ratio distributions in ERBAs. Several other notable
studies have provided valuable insights into the physics of atomization in RBAs [7, 81].
In general, most prior studies have given little or no attention to electric effects on
multiphase transport phenomena driving atomization. High-speed cameras and imaging
equipment are seldom placed in ERBA setups to prevent electric arcing [65]. For these
reasons, experimental studies of the atomization process in ERBAs are limited in detail.
Additionally, there is currently no experimental data on the charge distribution in atomized
droplets in ERBAs. Although some research has been conducted on the effect of charge
density on the surface finish quality [166] and charge distribution in atomized droplets [104]
in electrostatic sprays, they have not been studied extensively for charged rotary atomizers.
1.7 Literature Review: Numerical Studies
Numerical research studies have been conducted on ERBAs to investigate and optimize
their performance in various spray coating applications [33, 102, 174]. The studies focus on
understanding the complex fluid dynamics, droplet formation, and electrostatic interactions
involved in the atomization process. In Yang [180], an Euler approach is used to simulate
the paint flow in the atomizer to determine the paint droplet size and velocities at the
9
edge. This is followed by an Euler-Lagrange approach to investigate the droplet transfer
process, which was the main focus of the work. Recent numerical efforts have characterized
the effects of various parameters on droplet trajectories, droplet size distribution, spray
pattern, evaporation, and deposition in ERBAs. Pendar and Páscoa [126] implemented a
comprehensive three-dimensional Eulerian-Lagrangian model to investigate the electrostatic
spray transfer processes in the ERBAs to provide a thorough understanding of the complex
airflows around the nozzle and the motion of droplets toward the target surface. Zhu et al.
[183] showed that the gas flow greatly influences droplet trajectories near the substrate.
While larger droplets are guided towards the target surface by their momentum and
electrostatic forces, small droplets continue to be carried away by the air stream. Naoki
et al. [110] found that the presence of serrations improve the atomization quality and also
that small droplets are likely to be generated when serrations are deeper and more in number.
Colbert [34] developed a mathematical model of droplet trajectories generated by ERBAs to
predict the deposition rate. In Ellwood et al. [49], the authors present a simplified analysis
method for correlating RBA performance on droplet size and coating appearance. In Ray
and Henshaw [134], the rate of evaporation of the solvent for a wide range of operating
conditions was examined. An interesting approach to modeling the paint application process
was conducted by Guettler et al. [64] where they proposed a framework that combines
experimental input data, an injection model, and a metamodel-based optimization. Im et al.
[78] made a crucial discovery that the spray shape is critically sensitive to the liquid charge
density and electric field strength. Pendar and Páscoa [128] observed that electrification
significantly impacts the droplet size distribution and the deposition rate.
The studies listed above have shed light on the performance of RBAs and helped
optimize the spray characteristics by examining the macroscopic processes of spray cloud
formation and paint deposition. However, very little work has been done towards
understanding the microscopic phenomena that drive the atomization process to create the
10
droplets in the spray cloud. The importance of atomization to the performance of this
device will be discussed in the next section. Moreover, since the device is often operated in
an electric field, it is crucial to investigate the effect of the electric field on the atomization
process.
Saye et al. [146] performed high-fidelity numerical simulations of near-bell atomization
in RBAs using a hybrid finite difference level-set method and adaptive mesh-refinement
while accounting for the non-Newtonian behavior of the fluid. The solver was also equipped
with droplet tracking to extract some atomization statistics. Using their tool, several critical
statistics which are inaccessible to experimental approaches have provided insight into rotary
atomization. The study however lacks the investigation of electrified operation which is
commonly seen in industrial applications.
In a different study conducted by Pendar and Páscoa [127], breakup models are
employed to predict the atomization process and the droplet size in ERBAs. The models
aim to capture the mechanisms responsible for droplet formation and predict the droplet
size distribution among other statistics. Pendar and Páscoa [127] have examined the Reitz-
Diwakar [135], the Pilch-Erdman [129], and the modified Taylor Analogy Breakup [12] models
in their study. While it is not relevant to discuss the details of each model, it is notable that
each of these models makes assumptions regarding the underlying physics to form empirical
correlations between aerodynamic drag, surface tension, and viscous dissipation on breakup
lengths and time scales. While the approximations provided by such models are informative,
they may not accurately capture all the complexities of the atomization process.
1.8 Motivation to Study Electrified Atomization in RBAs
An important application of ERBAs is in the automotive industry. Almost every
car today drives out of a manufacturing facility after being painted by an ERBA in the
paint shop. With the world becoming increasingly dependent on automotive means of
11
transportation and motor vehicle production approaching 100 million units per year [115],
automobile manufacturers are looking for ways to minimize production costs. The process
of painting vehicles is one of the most expensive aspects of automobile manufacturing,
accounting for up to 50% of its total costs [8, 32]. These costs are associated with the energy
used in operating and maintaining the HVAC equipment of the painting booth, paint drying,
and control of pollutants like volatile organic compounds generated by paint overspray [59].
The overspray leads to significant material waste causing environmental and cost concerns.
This makes the paint shop responsible for over 80% of the environmental concerns in a
manufacturing facility [60]. Moreover, the global automotive paints and coatings market
size was valued at $17.34 bn in 2019 and is expected to be about $26.8 bn by 2027 [131].
Small improvements can result in significant cost savings and waste reduction.
The optimal performance of ERBAs can be improved in the automotive painting
industry. Currently, painting with ERBAs often requires over-coating to ensure sufficient
finish quality [141]. This has an effective reduction in its overall TE. The operation of ERBAs
can be improved to prevent using excess paint to achieve desirable characteristics [141].
Advancing our understanding of the breakup mechanisms in ERBAs can improve the
financial and environmental impact of the device in the automotive, agricultural, and
pharmaceutical industries. Crucial metrics such as the droplet size uniformity, surface
finish quality, TE, deposition thickness, and environmental impact are directly dependent
on the atomization process of ERBAs [9, 35, 43]. Currently, expensive trial-and-error
experimentation is used to find optimal operating conditions in industrial applications
of ERBAs. The physical processes leading to the spray characteristics can be studied
computationally using an appropriate CFD solver.
Previous numerical efforts using CFD to simulate ERBA atomizing flows have employed
breakup models that limit the accuracy of the droplet size distributions. Moreover, it is
crucial to understand the charge distribution in the resulting droplet cloud and its effect on
12
atomization and surface finish quality. The goal of this work is to address the unanswered
questions regarding the effect of the electric field on droplet characteristics in industrial
applications of ERBAs.
1.9 ERBA Problem Setup
To provide some context, the goal of this research work is to understand the effect
of electrified operation on the near-bell atomization characteristics in ERBAs. Since the
breakup activity occurs very close to the bell edge, the region of interest in this work extends
only a few millimeters from the edge (Fig. 1.2). More details on the domain setup are
presented in Chapter 3. I present the following overview of the simulation problem with the
intention to guide readers through the layout of the following sections. The mathematical
framework presented in the next chapter is motivated by the need to accurately model the
relevant physical processes, i.e., multiphase gas-liquid atomizing flows with a focus on the
electrohydrodynamics. It is worth noting that the liquid flow is affected by centrifugal and
Coriolis (due to the rotating nozzle), electric (due to the charged liquid and the background
electric field) and gravity forces after being ejected from the rotating bell cup.
1.10 Contributions
Most prior studies have not considered the electric effects on multiphase transport
phenomena driving atomization in ERBAs. Experimental studies on ERBA atomization
are limited in detail due to the complexities in measuring the atomization characteristics
especially pertaining to charge densities in droplets. The numerical studies performed
on ERBA atomization have explored the overall spray pattern and deposition rate as a
macroscopic phenomenon [96]. Although they include microscopic droplets, they are often
inserted as droplets into the domain and do not form after the typical atomization processes.
13
In this work, high-fidelity simulations of primary and secondary atomization processes with
electrohydrodynamic effects are performed. Some of my key contributions include the
development of a novel method to compute an accurate electric potential field in the flow
domain called e-Mesh (Section 3.2), the implementation of a non-Newtonian viscosity model
in the liquid phase (Section 2.1.4) and the implementation of a rotating reference frame
into the preexisting framework (Section 2.1.3). Additionally, I have incorporated a droplet
ancestry extraction tool to obtain crucial statistics relevant to atomization characteristics
which are inaccessible to experimental techniques to provide an in-depth analysis into the
physics of atomization in ERBAs (Section 3.4). Using these tools, high-fidelity simulations
of the near-bell atomization process is conducted with a focus on understanding the effect
of electrified operation on the breakup processes and droplet formation. The details of my
work have been documented in a series of published articles [86, 87]. In addition, an article
on atomization statistics in ERBAs has been submitted and another article detailing the
effect of the shear-thinning behavior of paint on atomization is in preparation.
14
Figure 1.2: A schematic representation of the domain geometry showing the trajectory of
the liquid jet (green arrow) and the external body forces acting upon it (not to scale). The
device is typically operated at a distance of about 0.25 m from the target surface.
15
1.11 Dissertation Layout
This thesis dissertation is laid out as follows: In Chapter 2, the equations governing all
the multiphysics processes relevant to near-bell atomization in ERBAs are reviewed and a
mathematical formulation is laid out for the numerical model. In Chapter 3, the numerical
methods employed when solving the governing equations are presented. The discussion
includes an addition to the flow solver that is essential in computing an accurate electric
field in the computational domain, a mesh sensitivity study, and the application of a droplet
ancestry statistics extraction tool to obtain previously inaccessible information that will
help us understand the physics of atomization. In Chapter 4, the tool built to simulate high-
fidelity ERBA flows is demonstrated by presenting some results and conducting a parameter
study. Chapter 5 highlights the findings of a comparative study between a non-electrified and
an electrified jet and draws inferences regarding the effect of electrification on atomization
statistics. Chapter 6 contains a similar comparative study between a Newtonian and a non-
Newtonian jet to help isolate and understand the effects of atomization of a shear-thinning
fluid. The dissertation is summarized and concluded in Chapter 7 which also contains a
comprehensive description of several potential improvements to the setup and methods used
in this work.
16
CHAPTER TWO
MATHEMATICAL FORMULATION
This section contains mathematical descriptions of the relevant physical process that
are considered to simulate atomization in ERBAs. This section is laid out with particular
consideration to the setup described in Section 1.9.
2.1 Governing Equations
2.1.1 Navier Stokes Equations for Multiphase Flows
For low-Mach number, variable density, multiphase flows, the mass and momentum
conservation laws in both phases can be written as
∂ρi
+∇ · (ρiui) = 0, (2.1)
∂t
∂ρiui
+∇ · (ρ f
iui ⊗ ui) = −∇pi +∇ · (σi + σe
i ) + γκnδs + fexternal (2.2)
∂t
where ρ is the density, u is the velocity field vector, t is time, and p is the hydrodynamic
pressure. The subscript i denotes the fluid phase (gas or liquid). The term γκnδs denotes
surface tension force, where γ is the surface tension coefficient, κ is the local curvature of the
interface computed using the ACES technique [123], n is the normal vector to the interface,
and δs is a Dirac-delta function that is nonzero only on the interface [120]. The viscous stress
tensor σf
i in Eq. 2.2 is given by
σf 2
i = µi(∇ui +∇uT
i )− µi(∇ · ui)I (2.3)
3
17
where µ is the dynamic viscosity and I is the identity tensor. These equations form the basis
of the fluid dynamics module in this work and further details can be found in Owkes and
Desjardins [119, 122] and Desjardins et al. [39, 40]. σe
i is the Maxwell stress tensor which
will be described in the next section. fexternal is a combination of the forces experienced by
a fluid element in the system.
fexternal = fcentrifugal + fCoriolis + fgravity (2.4)
2.1.2 Electrohydrodynamics
The Maxwell stress tensor σe
i in Eq. 2.2 is given by
( )
σe
i = εiEi ⊗ Ei −
εi
Ei · Ei 1− ρi ∂εi
I (2.5)
2 εi ∂ρi
where E is the electric field vector and ε is the electric permittivity equal to the product
of the relative permittivity κi and the vacuum permittivity space ε0. Magnetic effects
have been ignored since the EHD time scale is several orders of magnitude larger than the
magnetic time scale [145] in the atomization process that is of interest here. This electrostatic
assumption eliminates the effect of the velocity of the charges (i.e., current) on the electric
field thus dictating that the electric field is only influenced by the instantaneous electric
charge distribution. The electric force can be expressed as the divergence of the Maxwell
stress tensor, ( )
∇ · σe
i = qiEi −
1 2∇ ∇ 1 ∂εi
Ei εi + ρi E2
2 2 ∂ρ i (2.6)
i
where q is the volumetric electric charge density. The first term on the right-hand side of Eq.
2.6 is the Coulomb (or Lorentz) force. The second and third terms denote the dielectric and
the electrostrictive forces respectively, which are only significant in a transient electric field
with time scales several orders of magnitude larger than what is encountered in atomization
18
problems [84].
As the electric field vector is irrotational it can be expressed in terms of the scalar
electric potential φ as
Ei = −∇φ. (2.7)
The electric potential Poisson equation describes the relationship between charge density qi
and φ as
−∇ · (εi∇φ) = qi. (2.8)
The electric Poisson equation is solved using the hypre library of routines for scalable parallel
solutions of linear systems [51, 52].
The accumulation of bulk volumetric charge as a surface charge in a thin electric
boundary layer much smaller than the hydrodynamic boundary layer is a commonly made
assumption in charged jet simulations [171, 175]. The electric charge has previously been
modeled in two ways: either a constant bulk volumetric charge or a fully relaxed surface
charge. According to the classic leaky dielectric model [105, 145, 165] that is commonly used
to describe the effects of electric charge in dielectric liquids, the fundamental underlying
assumption is that bulk volumetric electric charge has sufficient time to fully relax to a
surface charge. However, for atomizing flows the advection time scale is similar to the
charge relaxation time scale and a fully relaxed charge assumption is invalid. The charge
relaxation time τq represents the time required for volumetric charge ql to relax to the surface
[37, 84]. The fluid advection timescale τf is a characteristic time for a fluid element to move
across a relevant length scale L0. The ratio of the two time scales called the electric Reynolds
number Ree was formulated in Stuetzer [162] as
19
εl
τq = , (2.9)
ψl ql
L0
τf = , (2.10)
u
τq
Ree = (2.11)
τf
where values of the quantities that correspond to typical automotive paints and RBA
operations listed in Table 3.1 yield a value of Ree approximately equal to 8.93. This suggests
that the advection and charge relaxation timescales are comparable in this application. The
above time scale analysis highlights the necessity to model the process of charge relaxation
by migration and diffusion in addition to advection with the fluid velocity. The common
assumption of charges to be accumulated on the surface disallows the charge migration
dynamics and limits the accuracy of the droplet charge distribution after atomization. In
this work, we assume a constant bulk volumetric charge at the injector and allow charges
to relax as the jet propagates through the domain. Charge transport is described by the
conservation equation
∂qi
+∇ · Ji = 0 (2.12)
∂t
where the current density Ji is formulated as [106, 172]
Ji = qiui + qiψiEi −Di∇qi (2.13)
where ψi is the charge mobility coefficient and Di is the molecular diffusivity. The three terms
that contribute to the current density can be described as advection due to the velocity field,
advection due to the electrical velocity, and diffusion.
In summary, the charge density field (q) and boundary conditions (BCs) in electric
potential are used to obtain the electric potential field (φ) using Eq. 2.8. The electric field
20
(E) is then obtained using Eq. 2.7 after which Eq. 2.6 allows for the calculation of the
Coulomb/electric force (felectric) on the charged fluid. More details on the equations involved
in the formulation of the EHD module can be found in Sheehy and Owkes [151].
2.1.3 Rotating Reference Frame
In an ERBA undergoing ligament breakup mode, each ligament exiting the serrated bell
edge behaves similarly with statistically similar droplet size distribution and charge transport
behavior. Therefore, a single ligament ejected from the bell edge is the focus of the numerical
simulations. The reference frame is attached to the rotating nozzle and ligament.
To account for the rotating nozzle and the consequent forces, equations that describe
the rotating reference frame are included in the model. This subjects the liquid to centrifugal
and Coriolis forces which can be formulated as
fcentrifugal = −ρiω × (ω × r) (2.14)
fCoriolis = −2ρiω × ui (2.15)
where ω is the angular velocity and r is the radial distance from the axis of rotation.
2.1.4 Non-Newtonian Viscosity Relationship
Most modern automotive paints are non-Newtonian in nature. In particular, they
exhibit shear-thinning or pseudoplastic behavior. There are a number of well-established
shear-thinning viscosity models. The power law fluid model predicts a boundless shear-
thinning viscosity for high shear rates. The Carreau model first proposed by Pierre Carreau
and the Cross fluid model both predict the viscosity in a liquid that has both a low and
high viscosity limit at high and low shear rates, respectively. The Herschel-Bulkley fluid
model limits the lowest viscosity attainable by the liquid at infinitely high shear rates [69].
In this work, however, the dynamic viscosity of a liquid (µl) subject to a local shear rate
21
(η) is formulated using a modified power law. The modified power law model is a suitable
approximation for automotive paints as the empirical viscosity-shear rate data provided by
PPG Industries shows a similar behavior (Figure 2.1). Although a Herschel-Bulkley model
would have also been suitable, the unavailability of the yield stress of the liquid would have
required making an arbitrary assumption for the purpose of a fit. According to the modified
power law, liquid viscosity is related to the local shear rate as
µ = Υ · η (Υ2−1)
l 1 + Υ0 (2.16)
where Υ0, Υ1, and Υ2 are viscosity coefficients pertaining to the power law. The
modification lies in the addition of the term Υ0 which is absent in a traditional power law
equation. The coefficients are obtained by applying a fitting function on empirical viscosity-
shear rate data. Using empirical data, the set of values of the viscosity coefficients that
models its shear-thinning behavior appropriately have been determined and are listed in
Table 3.1.
Figure 2.1 shows the viscosity-shear rate relationship of the automotive paint modeled
in this study and compares it with the modified power law fit which is shown to be an
acceptable approximation for this work. For clarity, the value of viscosity in interfacial cells
(µavg) is harmonically averaged as
µl · µg
µavg = (2.17)
µg · VF + µl · (1− VF)
where µg is the dynamic viscosity of the gas phase and VF is the liquid volume fraction
in the cell. The harmonic averaging method was provided in Patankar [125] and has been
accepted as a reasonable approximation for several decades [4].
22
Figure 2.1: The viscosity-shear rate relationship of automotive paint is approximated well
using a modified power law fit. Note the logarithmic scale on the axes.
2.2 Validation of Physics Modules
The EHD module has seen continued development and validation efforts have been
presented in [151]. To validate the implementation of centrifugal force (described in Section
2.1.3), numerical experiments of a rotating tank configuration are conducted. The solution
for the steady-state interface profile h(x) of a liquid in a rotating tank is analytically derived
as
1
h(x) = h0 + (ωx)2 (2.18)
2g
where h0 is the height of the liquid at the axis of rotation and g is the gravitational
acceleration [90]. The derivation of this analytical solution is provided in Section A.1.
Numerical results show excellent agreement with the analytical solution for varying rotation
23
rates (Fig. 2.2).
Figure 2.2: Validation of centrifugal force: Numerical results of a rotating tank experiment
where cases 1, 2, and 3 correspond to rotation rates of 4, 10, and 14 rad/s respectively. The
analytical solution of the interface (Eq. 2.18) at various locations is plotted with red circles.
The numerically computed liquid interface and liquid bulk are shown as black lines and blue
areas respectively. The dashed line represents the axis of rotation for each case.
The numerical implementation of the Coriolis force is tested by simulating a liquid jet
in a rotating frame. The path of an object experiencing the centrifugal and Coriolis forces
follows a circle with a growing radius, i.e., a spiral. The analytical solution for the position
of the object (x, y) at a given time t is obtained as
x(t) = (x0 + vx0t) cos (ωt) + t (vy0 + ωx0) sin (ωt)
(2.19)
y(t) = − (x0 + vx0t) sin (ωt) + t (vy0 + ωx0) cos (ωt)
where vx0 and vy0 are the x and y components of the initial velocity of the jet and x0 is the
position of the jet at t = 0. The derivation of this solution is also provided in Section A.2.
This solution is derived for a rigid point-like object with its all mass concentrated at its center.
Liquid jets, however, are not point-like approximations and include effects of viscosity and
surface tension. Despite these challenges, satisfactory results are obtained for the the effect
of the Coriolis force on liquid jets in a rotating reference frame. Numerically simulated liquid
24
jets follow the analytically predicted trajectory for varying rotation rates (Fig. 2.3).
Figure 2.3: Validation of Coriolis force: Numerical results of the liquid jet interfaces in a
reference frame rotating counter-clockwise about an axis that points out of the 2D plane.
The analytical solution (Eq. 2.19) is shown as a superimposed solid line for each case.
25
CHAPTER THREE
NUMERICAL METHODS
The numerical setup is designed based on the various physics modules discussed in
Section 2.1. The modules have been implemented within a code called NGA - a high-order,
fully conservative, variable density, low Mach number Navier-Stokes solver that contains
various multi-physics modules implemented in parallel using message passing interface
(MPI). The formulation discretely conserves mass and momentum in a periodic domain.
NGA uses a conservative unsplit geometric volume-of-fluid (VOF) scheme as described in
the works of Blanquart et al. [24], Desjardins et al. [39, 40] and Owkes and Desjardins
[119, 120, 121, 122]. NGA has been developed by several groups to solve multiphase EHD
flows, details of which can be found in Van Poppel et al. [172], Sheehy and Owkes [151].
The finite difference operators arising from discretizing the governing equations
(Section 2.1) are formulated using a conservative finite difference scheme [39, 119]. Consistent
and conservative transport near the interface is performed using the unsplit geometric VOF
scheme [122]. The interface curvature (κ) used in the surface tension computation (Eqn. 2.2)
is calculated using the Adjustable Curvature Evaluation Scale (ACES) scheme described in
[123]. The pressure Poisson equation is solved using the ghost fluid method [55] and the
blackbox multigrid technique [79]. The electric Poisson equation (Eq. 2.8) is solved using
the LGMRES method [140] included in the HYPRE library of linear algebra solvers and
preconditioners [76].
In this chapter, the domain setup is described in Section ??. A method to obtain
accurate electric potential fields on the flow domain is presented in Section 3.2. Section 3.3
contains the details of a mesh sensitivity study that was performed to identify a suitable mesh
resolution to capture the relevant physical processes. The implementation of a novel tool that
26
extracts droplet ancestry statistics to gain more insights into the physics of atomization is
presented in Section 3.4. Finally, the computational cost of the simulation is briefly discussed
in Section 3.5.
3.1 Domain Setup
The domain is sized in such a way that it is large enough to capture most of the breakup
activity occurring near the bell edge. It has dimensions of 12.96 mm × 960 µm × 360 µm.
It is represented schematically in Fig. 1.2. In this work, a bell of radius 25 mm is studied.
Some ERBAs are equipped with serrations present at the edge of the bell which causes the
fluid to be ejected as thin jets. In the absence of these serrations, the fluid forms liquid sheets
on the inner surface of the bell. These liquid sheets may be subject to interfacial instability
and form periodic ligaments at the edge of the bell or break up in sheet mode depending
on the operating conditions. In this work, a single liquid jet exiting one of the serrations is
simulated. Appropriate boundary conditions are imposed to inform every jet of neighboring
flows from adjacent serrations which will be discussed in the following paragraphs.
The bell simulated in this work is equipped with 418 serrations along its edge, each of
which is 375 µm wide. While the liquid profile in a serration does not have a circular cross-
section, a 60 µm diameter jet that exits the bell edge is approximated for a bell operating at
40 kRPM [44]. For a nozzle flow rate of 200 mL/min, the flow rate through each serration is
7.96×10−9 m3/s. In the numerical setup, the jet is injected into the domain from the top wall
with initial velocity components that correspond to an edge angle of 25◦ through an elliptical
cross section. Breakup is induced by imparting velocity modulations to the inflow. The
inflow velocity (uinflow) is modulated by adding sinusoidal perturbations at a superposition
of three (40, 50, and 60 kHz) frequencies (ω) at an amplitude à which is 9% of its original
velocity (u) (Eq. 3.1). Several combinations of frequencies were tested and the particular
combination chosen produces a perturbation in the jet that is small enough in magnitude
27
that it does not largely affect the trajectory near the inflow. The goal of imparting these
velocity modulations is to replicate natural perturbations that would appear in experimental
setups due to instrument vibrations or fluctuations in film thickness that act as sources of
minor modifications to the flow [148, 159].
∑3
uinflow = u + Ã sin (2πωnt) . (3.1)
n=1
In ERBAs, the bell edge plane is angled away from the plane perpendicular to the axis
of rotation (Fig. 3.1). The axis of rotation, which is the center of the bell, is located one
radial length away (in the −x direction) from the inflow.
A magnitude of 2.879 C/m3 for liquid charging is applicable to ERBAs used in industrial
paint shops [48, 128]. While the charging mode in practice does not typically yield a uniform
bulk initial charge distribution, such a case is assumed for simplicity. Improvements towards
more accurate modeling of liquid charging are discussed in the future works section in
Section 7.3. It is worth noting that the charges are still free to redistribute by the processes
of advection, diffusion, and migration as dictated by Eqn. 2.12 once the liquid is in the flow
domain.
The top boundary is a wall that allows slip velocity and the left boundary is a no-flux
wall. The bottom and right boundaries are convective outflows. Periodic BCs are imposed
on the front and back walls. In addition to mass, momentum, and charge being transported
across the periodic z boundaries, the electric field is also computed with the effect of the
periodicity. This means that at any instance, the single jet simulated in the computational
domain is solved with full consideration of its nearest neighboring jets.
It is worth noting that a Cartesian grid is employed in this work for the following reason
- despite the serrations lying along the perimeter of a bell cup which is curved, the focus
of this simulation is a thin slice (360 µm wide) while the domain length along the radial
28
Figure 3.1: A schematic representation of the angle associated with the ejection of liquid
from the bell edge, calculated between the axis of rotation and the edge surface.
direction (12.96 mm) which is about half the diameter of the bell cup. This means that the
curvature of the bell is significant enough to consider using a cylindrical coordinate system
and grid. An implication of using a Cartesian grid and applying periodic BCs along the
z boundaries means that the jets (and the resulting liquid structures) interact more then
they should at the downstream part of the domain. This is a shortcoming in the current
formulation and has been addressed in Section 7.3.
Table 3.1 contains a list of parameter values used in simulations.
29
Table 3.1: Standard values of parameters in the ERBA simulations
Property Symbol Value Unit
Bell rotation rate ω 40×103 RPM
Bell radius Rbell 0.025 m
Edge angle - 25 deg
Inlet jet diameter Djet 60×10−6 m
Liq. flow rate - 7.96×10−9 m3/s
Liq. viscosity (Newtonian) µl 0.1 Pa.s
Liq. density ρl 1000 kg/m3
Surface tension γ 0.03 N/m
Bell electric potential V 3
bell 80×10 V
Liq. charge density ql 2.879 C/m3
Liq. rel. permittivity κl 50 -
Liq. diffusion constant D −6
l 2×10 m2/s
Liq. charge mobility ψl 1.79×10−8 m2/V.s
γ0 0.11 Pa.s
Liq. viscosity coefficients γ1 0.11 Pa.s
γ2 0.36 -
30
3.2 The Need for e-Mesh
An addition to the EHD module of flow solver was developed for this project and mainly
involves using a domain called e-Mesh. e-Mesh is much larger than the flow-solver domain,
henceforth called the NS-Mesh, where the Navier-Stokes equations described in Section 2.1.1
are solved. Since our interest lies in primary and secondary atomization, NS-Mesh lies at
the edge of the bell and extends only a few millimeters outside of it. NS-Mesh requires a
well-defined electric potential field (φ) to obtain the electric field (E) and solve the EHD
equations. φ is solved using appropriate boundary conditions (BCs) on NS-Mesh. However,
Fig. 3.2 highlights that φ is well defined (labeled Dirichlet) with values readily available only
at the nozzle (bell electric potential) and at the target surface (grounded) but unavailable at
the boundaries of NS-Mesh. Instead of assuming values of φ BCs on NS-Mesh, a new domain
called e-Mesh that spans between regions of well-defined φ is initialized and used to obtain
accurate φ BCs on NS-Mesh. e-Mesh is about 4000 times larger than NS-Mesh, i.e., the
region of interest in this project lies in about 0.0125% of the volume of e-Mesh. Making NS-
Mesh this large would require significantly more computational resources to solve the fluid
dynamics far from the bell edge. Such an approach would impractical and unfeasible. The
breakup activity and charge transport within the liquid, which are the processes of interest
in this work, are already completed within the NS-Mesh domain. In this way, e-Mesh is used
exclusively as an electric potential solver.
3.2.1 e-Mesh Framework
There are several steps involved in the employment of e-Mesh to obtain an electric
potential field. A well defined charge density field (q) is required to solve for φ on e-Mesh
(φe−Mesh). After obtaining φe−Mesh, the values need to be interpolated from e-Mesh to cells
that lie on the boundaries of NS-Mesh to be applied as BCs to compute φ on NS-Mesh. The
31
Figure 3.2: A schematic representation of e-Mesh and available φ BCs shown in context with
NS-Mesh (not to scale).
step-wise process can be listed as follows -
1. Initialize the e-Mesh domain geometry.
2. Impose boundary conditions on e-Mesh.
3. Migrate charge density q from NS-Mesh to e-Mesh.
4. Solve the electric Poisson equation on e-Mesh.
5. Migrate electric potential from e-Mesh to NS-Mesh.
Steps 1 - 2 are performed once at the beginning of the simulation and steps 3 - 5 are
performed every time the electric potential BCs are updated in NS-Mesh. Details of each
step are provided below.
3.2.2 e-Mesh Domain Initialization
When operating in industrial applications, the bell is typically about 0.25 m away from
the surface. This makes e-Mesh 0.25 m long in the y-direction. e-Mesh extends to a distance
of about 0.2 m away from the axis along the x-direction, φ contour lines are measured to
be parallel to the surface. This distance would allow for Neumann BCs on the right wall.
32
Figure 3.3: A schematic representation of the domain geometry of the numerical problem
showing positions and dimensions (not to scale) of NS-Mesh and e-Mesh.
e-Mesh and NS-Mesh are equal in length in the z-direction, allowing for periodic BCs on the
z+ and z− walls. Therefore, e-Mesh, which spans from the bell edge to the target surface,
has dimensions of 0.2 m × 0.25 m × 360 µm (Fig. 3.3).
e-Mesh has a stretched grid along the x and y directions to facilitate a higher
concentration of cells where it coincides with NS-Mesh near the bell edge. Grid points along
x and y directions on e-Mesh are placed according to a Gaussian function such that they
are stretched around the bell edge. Grid points along z are coincident in both meshes. As
a result, nodes and cells on both the meshes are not coincident. A schematic representation
of the stretched grid is shown in Fig. 3.4. NS-Mesh has a uniform grid consisting of cubic
cells and e-Mesh has a non-uniform grid made up of cuboidal cells. NGA, which is the code
that performs numerical calculations, is a code written with the ability to be run on multiple
processors using message passing interface (MPI). It is common practice to distribute the
computational load by allowing each processor to solve only a part of the entire domain. As a
33
Figure 3.4: Schematic representation of the grids on NS-Mesh (yellow) and e-Mesh (blue)
overlapping. e-Mesh is formulated with a stretched grid to facilitate a higher concentration
of cells near the bell edge. The nozzle is shown for reference and is not to scale.
result, NS-Mesh and e-Mesh are both split into smaller sub-domains, one for each processor.
The processor domains on e-Mesh are decomposed in such a way that the number of cells
handled by each processor is balanced. This yields a domain decomposition that also looks
non-uniform as shown in Fig. 3.5.
The discretization of the governing equations describing EHD flows (Section 2.1.2) yield
a set of several finite different operations. It is somewhat straight-forward to perform the
operations on the uniform grid on NS-Mesh. However, the finite difference operations on
e-Mesh, which has a non-uniform grid with non-constant intervals, are obtained as described
in Sundqvist and Veronis [163].
3.2.3 e-Mesh Boundary Conditions
A Dirichlet BC corresponding to the value of the bell electric potential is imposed on the
section of the top boundary of e-Mesh that overlaps the nozzle. A zero potential boundary
34
Figure 3.5: Schematic representation of the overlap of processors between NS-Mesh (green)
and e-Mesh (pink). For example, overlapping information exists on processor 3 on the NS-
Mesh (NS processor 3, highlighted in bright green) and processor 1 on the e-Mesh (e processor 1,
highlighted in bright pink). The nozzle is shown for reference and is not to scale.
35
condition is imposed on the bottom boundary which represents the grounded target surface.
Periodic BCs for potential are imposed on the z+ and z− boundaries of e-Mesh and Neumann
BCs for potential are imposed on all its remaining boundaries. The right boundary (x+) has
been chosen in such a way that it is far enough away from the axis of rotation that the
electric potential field lines are parallel to the target surface, allowing for a Neumann BC to
be imposed. A schematic representation of the BCs applied on each boundary is shown in
Fig. 3.2.
3.2.4 Migration Charges to e-Mesh (NS→e Migration)
As mentioned, the goal of e-Mesh is to obtain an electric potential field that can be
used to impose well-informed BCs on the boundaries of NS-Mesh. In order to obtain an
electric potential field on e-Mesh, the charge density field is required to be migrated from
NS-Mesh. There are two steps involved in the migration of charge density values from
NS-Mesh to e-Mesh (NS-e migration) - 1) conservative interpolation and 2) inter-processor
communication. Interpolation of information is necessary due to the non-coincidence of cells
between NS-Mesh and e-Mesh, outlined schematically in Fig. 3.4.
For interpolation, every cell on e-Mesh (e cell) is flagged by whether it shares volume
with NS-Mesh and requires to participate in the migration process by comparing its spatial
extents to the spatial extents of NS-Mesh -
if [(e cell
x− ≥ NS mesh and e cell mesh cell mesh
x− x− ≤ NS x+ ) or (e x+ ≤ NS x+ and e cell
x+ ≥ NS mesh
x− )]
and [(e cell ≥ NS mesh and e cell mesh cell mesh cell mesh
y− y− y− ≤ NS y+ ) or (e y+ ≤ NS y+ and e y+ ≥ NS y− )]
and [(e cell
z− ≥ NS mesh cell mesh
z− and e z− ≤ NS z+ ) or (e cell ≤ NS mesh cell mesh
z+ z+ and e z+ ≥ NS z− )]
then NS-e cell
overlap = true,
(3.2)
where each term represents a location outlined in Fig. 3.6. In a cell on e-Mesh (e cell) that
36
Figure 3.6: Locations of NS-Mesh boundaries (yellow) with respect to a cell on e-Mesh e cell
(blue)
needs to receive charge density, i.e., if NS-e cell
overlap = true, overlapping volume fractions in
the x, y and z directions (VFx, VFy and VFz respectively) are calculated in every cell on
NS-Mesh (NS cell) by comparing the locations outlined in Fig. 3.7. These quantities are used
to calculate the charge to be migrated (Qmigration) from every NS cell (Eq. group 3.3).
37
Figure 3.7: Locations of the boundaries of a cell on NS-Mesh NS cell (yellow) and a cell on
e-Mesh e cell (blue)
( ( ))
NS cell − e cell
VFx = min(1.0,max 0.0, x+ x−
( NS cell
x+ − NS cell
x− ))
NS cell − e cell
− min(1.0,max(0.0, x+ x+
NS cell − NS cell
x+ x− ))
NS cell
y+ − e cell
y−
VFy = min(1.0,max(0.0,
NS cell cell
y+ − NS y− ))
NS cell − e cell
− y+ y+
min(1.0,max(0.0, (3.3)
NS cell
y+ − NS cell
y− ))
NS cell − e cell
VFz = min(1.0,max(0.0, z+ z−
NS cell
z+ − NS cell
z− ))
NS cell − e cell
− min 1.0,max 0.0, z+ z+
∑( NS cell − NS cell
z+ z− )
Qmigration = VFx × VFy × VFz × NS cell
q × NS cell cell
volume × NS VF .
Here, the summation on the left hand size operates over all the cells in the domain.
NS cell and NS cell
q VF are the charge density and liquid fraction in the NS cell respectively. The
term NS cell ×NS cell ×NS cell
q volume VF represents the amount of charge in the NS cell and the term
38
VFx ×VFy ×VFz can be interpreted as its unit volumetric contribution towards migration.
The scalar charge field (Q ) can be divided by the corresponding cell volume (e cell
migration volume)
to obtain the charge density field in e-Mesh (qe−Mesh). This way, the scalar electric charge
density field q is conservatively interpolated from NS-Mesh to overlapping cells on e-Mesh.
It is evident from Fig. 3.5 that there are regions handled by different processors on
each mesh. A necessity for inter-processor communication of values during migration thus
arises. The figure shows that processors 2, 3 and 4 (NS processor 2, NS processor 3, NS processor 4)
will communicate the charge density field that exists on NS-Mesh to the overlapping region
on e-Mesh that is handled by processor 1 (e processor 1).
NS-e migration is outlined in Algorithm 3.1. Here, Qglobal migration is the total charge to
be migrated into an e cell from all overlapping NS cells and e cell
q is the charge density in the
e cell.
Algorithm 3.1
NS→e charge density migration algorithm testing testing testing testing testing
for every e cell
for every NS cell in each processor
if NS-e cell
overlap = true, then
compute volume fractions and Qmigration from Eq. 3.3
end for
sum Qmigration from all processors to get Qglobal migration
identify the processor, my proc, that contains current e cell
if currently at my proc then
set charge density, e cell cell
q = Qglobal migration/e volume
end if
end for
end for
39
Figure 3.8: Electric potential field (φe−Mesh) on e-Mesh. The nozzle is shown for reference
and is not to scale.
3.2.5 Solving the Electric Poisson Equation on e-Mesh
The electric Poisson equation on e-Mesh (Eq. 2.8) is solved using the LGMRES method
[140] included in the HYPRE library of linear algebra solvers and preconditioners [76]. In the
previous sections, we presented how e-Mesh is 1) informed with the appropriate boundary
conditions of potential and 2) assigned a charge density field by conservatively interpolating
and communicating values from NS-Mesh. Since the spatial distribution of liquid volume
is insignificant compared to the size of e-Mesh, the permittivity term (ε) is neglected when
solving the electric Poisson equation on e-Mesh. It is important to note that φ on NS-Mesh
is recomputed with the consideration of the liquid configuration. The electric potential field
(φe−Mesh) on e-Mesh is shown in Fig. 3.8 and spans from the bell potential (80 kV in this
case) to the grounded voltage at the target surface.
40
3.2.6 Migrating Electric Potential to NS-Mesh (e→NS Migration)
Once the potential field φe−Mesh is computed , its values are migrated to the boundaries
of NS-Mesh (e-NS migration) so that they can be used as BCs to solve for electric potential
on NS-Mesh. This process involves similar previously discussed steps of interpolation and
inter-processor communication. For a processor on NS-Mesh, every overlapping processor on
e-Mesh is identified by comparing their spatial extents -
( )
if [ e processor i ≤ NS processor j and e processor i ≥ NS processor j
( x− x+ x+ x− )]
and [ processor i processor j processor i processor j
(e y− ≤ NS y+ and e y+ ≥ NS y− )] (3.4)
and [ e processor i ≤ NS processor j and e processor i ≥ NS processor j
z− z+ z+ z− ]
then e-NS processor i−j
overlap = true.
where i and j are processor IDs. For example, Fig. 3.9 shows a schematic representation of
this comparison for NS processor 3 and e processor 1.
Figure 3.9: Boundaries of regions occupied by NS processor 3 (yellow) and e processor 1 (blue).
Here, e-NS processor 1−3
overlap = true.
At the beginning of the simulation, a comprehensive communications list called
41
Algorithm 3.2
Calculating index extents for overlap between NS-Mesh and every e processor n
for n = 1, nproc
for all grid points along x on NS-Mesh, i = imin, imax
if NS cell
center ≥ e processor n
x− then
set NS left index
overlap = i; exit
end if
end for
for all grid points along x on NS-Mesh, i = imax, imin, −1
if NS cell processor n
center ≤ e x+ then
set NS right index
overlap = i; exit
end if
end for
repeat in other two dimensions
end for
Table 3.2: Simplified commlist corresponding to the setup shown in
Fig. 3.5
sender
e processor 1 e processor 2 e processor 3 e processor 4
receiver
NS processor 1 true false false false
NS processor 2 true true false false
NS processor 3 true false false false
NS processor 4 true true false false
42
Algorithm 3.3
Calculating index extents for overlap between e-Mesh and every NS processor n
for n = 1, nproc
for all grid points along x on e-Mesh, i = e imin, e imax
if e cell ≥ NS processor n
center x− then
set e left index
overlap = i; exit
end if
end for
for all grid points along x on e-Mesh, i = e imax, e imin, −1
if e cell
center ≤ NS processor n
x+ then
set e right index
overlap = i; exit
end if
end for
repeat in other two dimensions
end for
43
commlist is created using e-NS processor i−j
overlap . The flag values in commlist corresponding to the
schematic shown in Fig. 3.5 is tabulated in Table 3.2. In addition to the overlap flags, the
index extents of e cells that send data to each NS processor and the index extents of NS cells that
receive data from each e processor are identified and stored as outlined in Algorithms 3.2 and
3.3 respectively. e-NS processor i−j communication is then performed as dictated by commlist.
For a cell lying on the boundary of NS-mesh, say NS boundary cell, and belonging to
processor NS processor p, the nearest cell on e-Mesh is identified and an appropriate 12-point
stencil from e-Mesh is communicated to NS processor p. The 12 points comprise of two upstream
and downstream points each in x, y and z directions. Using the electric potential values
(φe−Mesh) at each point in the stencil, the electric potential value (φ) at the center of every
NS boundary cell
center is computed using a tricubic interpolation scheme described in Lekien and
Marsden [92]. Algorithm 3.4 outlines the process of e→NS migration.
3.2.7 Challenges and Cost With Using E-mesh
The e-Mesh module has a significant impact on the overall cost of the simulation. The
electric Poisson solver on e-Mesh is much more expensive than that on NS-Mesh perhaps
due to complexities introduced by the non-uniform grid. The migration of charge densities
and electric potential stencils involves repeated communication between multiple sets of
processors, making the process expensive. On average, solving the electric Poisson equation
(Eq. 2.8) on e-Mesh accounts for up-to 74% of a timestep. For this reason, the potential
field on e-Mesh is solved every 50 timesteps. The interval between two consecutive electric
potential field solutions on e-Mesh has an non-significant impact on the accuracy of the
electric potential boundary conditions on NS-Mesh. This is because the timestep is extremely
small (about 0.01 µs) and the fluid displacement between two consecutive e-Mesh solutions
is a fraction of the cell width.
In this way, e-Mesh is used exclusively as an electric potential solver. e-Mesh spans
44
Algorithm 3.4
e→NS electric potential field migration algorithm testing testing testing testing testing
for n = 1, length(commlist)
read sender, receiver and e index
overlap index range to be communicated
if currently at sender then
send φe−Mesh data to receiver
send NS index
overlap range to receiver
end if
if currently at receiver then
receive φ data and NS index
e−Mesh overlap
for i = NS left index , NS right index
overlap overlap
for ii = e left index right index
overlap , e overlap
if e cell ii
center ≥ NS cell i
center exit
repeat in other two dimensions to identify nearest e cell
center
compute tricubic interpolated value
set φ cell
NS-Mesh
end for
end for
end if
end for
45
between regions of well-defined electric potential values that serve as boundary conditions.
Making the flow domain (NS-Mesh) this large would be uneconomical. The charge density
field is made available on e-Mesh by migrating (interpolating and communicating) the data
from NS-Mesh. Using the available electric potential boundary conditions and the scalar
charge density field, the electric potential field is solved on e-Mesh. The values of electric
potential are migrated (interpolated and communicated) back to the boundaries of NS-Mesh
so that the electric potential field can be recomputed accurately on the flow domain.
3.3 Mesh Sensitivity
To determine the effect of mesh resolution on atomization and to help choose an optimal
numerical setup, three simulations of an electrified RBA jet (Table 3.1) are performed with
varying mesh resolutions (Table 3.3). Fig. 3.10 shows that droplets larger than 10 µm are
resolved in a 5 cells-per-diameter (CPD) mesh while the 10 CPD mesh captures droplets
smaller than 10 µm. Additionally, the 15 CPD mesh captures the droplets even smaller than
5 µm. Larger droplets are similarly resolved for all meshes. Between the two finer meshes
the smaller droplets are increasingly well resolved.
The necessity for a mesh finer than 5 CPD is apparent since the coarseness of the mesh
does not allow for charge migration to be captured. All liquid structures in the 5 CPD mesh
contain the injected charge density value of 2.879 C/m3. Using 10 and 15 CPD meshes, we
are able to resolve the charge migration process and investigate the resulting droplet charge
distribution. Most droplets contain a charge density close to the initial bulk charge density
value. Compared to the initial charge density, it is more probable to find a droplet with a
lower charge density. Fig. 3.11 highlights that charge migration is more resolved in a 15 CPD
mesh and is consistent in its trend with the 10 CPD mesh. Although there is uncertainty
with the smallest scales of droplet sizes, the number of droplets larger than 20 µm does not
vary on increasing the mesh resolution, suggesting that the larger droplets are well resolved
46
Table 3.3: Mesh sensitivity simulation domain parameters
Mesh Cells per Cell Processor No. of
ID diameter width count structures
05 CPD 5 12 µm 16 74
10 CPD 10 6 µm 96 121
15 CPD 15 4 µm 256 705
in a 15 CPD mesh. The relevant physical processes are sufficiently captured on a 15 CPD
mesh.
3.3.1 Sensitivity of Droplet Sizes to Mesh Resolution
Fig. 3.10 shows that droplets larger than 10 µm are resolved in a 5 cells-per-diameter
(CPD) mesh while the 10 CPD mesh captures droplets smaller than 10 µm. Additionally, the
15 CPD mesh captures the droplets even smaller than 5 µm. Larger droplets are similarly
resolved for all meshes. Between the two finer meshes the smaller droplets are increasingly
well resolved. Although there is uncertainty with the smallest scales of droplet sizes, the
number of droplets larger than 20 µm does not vary on increasing the mesh resolution,
suggesting that the larger droplets are well resolved in a 15 CPD mesh. Since this project
is driven by an industrially relevant process, it was decided that a 15 CPD mesh is able to
capture the physics of atomization associated with the smallest droplets commonly seen in
the paint shop which are about 10 µm in diameter [7, 49]. Therefore, all subsequent results
presented in Chapters 4, 5 and 6 are obtained from simulations performed on a 15 CPD
mesh.
47
Figure 3.10: Probability distribution of droplet sizes for varying mesh resolutions. The lines
show a log-normal fit of the data.
48
(a) 10 CPD (b) 15 CPD
Figure 3.11: Volume weighted charge density probability distributions for varying mesh
resolutions. Note the logarithmic scale on the y-axis.
3.3.2 Sensitivity of Droplet Charge Densities to Mesh Resolution
The coarseness of a 5 CPD mesh does not allow for charge migration to be captured. All
liquid structures in the 5 CPD mesh contain the injected charge density value of 2.879 C/m3.
Using 10 and 15 CPD meshes, the charge migration process can be resolved to investigate
the resulting droplet charge distribution. Most of the droplets contain a charge density close
to the initial bulk charge density value. Compared to the initial charge density, it is more
probable to find a droplet with a lower charge density. Fig. 3.11 highlights that charge
migration is more resolved in a 15 CPD mesh and is consistent in its trend with the 10 CPD
mesh.
49
3.4 Droplet Ancestry Extraction
3.4.1 Insights From Droplet Statistics
High-fidelity simulations of complex systems can resolve the relevant time and length
scales required to extract droplet ancestry statistics as a valuable means of understanding the
physics of atomization. In this work, a droplet ancestry extraction tool first introduced by
Rubel and Owkes [139] and further developed by Christensen and Owkes [31] is employed to
extract and store droplet and local flow field statistics during each breakup and coalescence
event in simulations of near-bell atomization in ERBAs.
The information is stored in a Neo4j graphical database and analyzed to track the
evolution of the liquid as it moves from the liquid core to small droplets in an atomizing
flow. Graph databases are frequently used in the business industry to identify relationships
between products and consumers through links that can reveal patterns in consumer
behavior. A similar approach is used to create paths that connect droplets through breakup
and coalescence events in atomization simulations.
The analysis allows us to compute breakup and coalescence statistics, which characterize
the drop size distributions and help us understand the processes driving breakup. Such
insights into complex atomizing flows are challenging to obtain using experimental methods.
The tool allows for the analysis of droplet behavior in regions that are not easily observable
or measurable, such as near the bell edge or within the dense droplet cloud. It also enables
the study of the behavior of charged droplets under the influence of a background electric
field.
3.4.2 Method of Extraction
The droplet ancestry extraction algorithm was developed and demonstrated on
atomizing flows by Rubel and Owkes [139] and Christensen and Owkes [31]. It is incorporated
50
into the ERBA simulation tool to track, collect, and store statistics of droplet breakup and
coalescence events throughout the simulation. The droplet ancestry algorithm identifies
breakup and coalescence events at every time-step. The algorithm operates through the
implementation of two identification numbers, namely the structure (S) and liquid (L)
identification numbers. S is a unique integer computed for each liquid structure using a
connected component labeling algorithm [67, 68, 75] and provides information of the current
state of all the structures. L is a unique integer that is transported with the liquid providing
a time history of the fluid. Together, the numbers allow for breakup and coalescence events to
be identified [31]. Coalescence events occur when structures occupy the same computational
cell and breakup events occur when a gap of one computational cell forms between structures
containing the same S value. This formulation prevents the occurrence of fictitious, repetitive
coalescence and breakup events [31].
3.4.3 Storing Statistics in a Database
Once these key events are identified, the information listed below regarding parent
and child droplets are extracted and stored in a CSV file for convenient analysis and post-
processing.
1. Breakup/coalescence event identifier
2. Event time
3. Liquid identification numbers of the parent and child structures
4. Structure identification numbers of the parent and child structures
5. Volumes of the parent and child structures
6. Principal axes lengths of the parent and child structures
7. Position vector of the structure
8. Velocity vector of the structure
9. Velocity vector of the surrounding gas phase
51
10. Vorticity vector of the structure
11. Viscosity values in the parent and child structures
12. Charge density values in the parent and child structures
3.4.4 Post-processing the Statistics
The droplet ancestry statistics extraction tool operates throughout the simulation and
connects successive breakup/coalescence events to create an ancestry or family tree describing
the evolution of the spray cloud. Fig. 3.12 displays the ancestry as interpreted by a Neo4j
graphical database for the simulation of a non-electrified 15 CPD jet. The nodes represent
individual droplets and the lines connecting them represent the associated breakup events.
The ligament core is in the center of the image. To maintain the visibility of the nodes
and lines, only droplets which broke up up-to six times and their “ancestors” are displayed
in the image. The statistics associated with each node are also stored in the database and
can be processed to provide insights on the breakup mechanisms in ERBAs. Results and
insights from such analyses are presented in Chapters 5 and 6 where extracted statistics are
compared to quantify the effects of different flow conditions on atomization.
52
Figure 3.12: Droplet ancestries in the jet shown in Fig. 4.1 represented using the Neo4j
graphical database. The colors represent breakup levels 1 (gray), 2 (yellow), 3 (green), 4
(red), 5 (pink) and 6 (blue) breakup events. Lines connect related droplets, i.e., droplets
which split from each other. The ligament core (black) is present at the center of the image.
Each node contains relevant statistics from the breakup event.
53
3.5 Computational Cost
Fig. 3.13 shows the computational cost associated with each module in the solver during
a timestep in which e-Mesh does not participate. As mentioned, e-Mesh is used to update
the electric potential boundary conditions on the flow domain (NS-Mesh) once every 50
timesteps. In a timestep where e-Mesh is active, it takes up upto 74% of the cost, which is
significant.
A closer look at Fig. 3.13 shows that the two most expensive modules are the ones
that handle pressure and the electric field. Interestingly, both modules contain a step where
a Poisson equation is solved. The Poisson solver is often computationally expensive in
numerical methods. In this setup, the electric Poisson equation is solved using the LGMRES
iterative method which is available in the HYPRE library of linear algebra solvers and
preconditioners [76]. The pressure Poisson equation is solved using a multigrid technique
which employs a coarse grid to accelerate the convergence of the solution. By leveraging
the grid hierarchy and adapting the solution at multiple levels, the multigrid method can
converge to an accurate solution more quickly than a single-level LGMRES iterative solver.
The transport module handles multiphase flows, interface tracking and surface tension
computations. The 1% sector comprises of writing simulation metrics to log files, generating
outputs for post-processing, liquid and structure identifications associated with the droplet
ancestry tool, and rest periods. Each 15 CPD simulation contains 69,984,000 computational
cells on NS-Mesh and 5,400,000 computational cells on e-Mesh. The wall time for 300 µs
of simulation time was about 650 hours without electrification on 128 processors and about
1500 hours for a simulation that included EHD effects on 256 processors. The computational
efforts were performed on the Hyalite High Performance Computing System and the Tempest
High Performance Computing System, operated and supported by University Information
Technology Research Cyberinfrastructure at Montana State University.
54
Figure 3.13: Breakdown of the computational cost of notable modules in the numerical
solver. The electric module which handles solutions to the EHD equations presented in
Section 2.1.2 is the most expensive part of the solver taking up to 77% of the cost; the
transport module which handles multiphase flows and interface tracking accounts for 10%
of the cost; the pressure solver accounts for 12% of the cost.
55
CHAPTER FOUR
DEMONSTRATION OF NUMERICAL FRAMEWORK
In this chapter, the tool developed to simulate near-bell atomizing flows in ERBAs is
demonstrated to be well-posed. A parameter study is conducted to understand the effects
of different operating conditions on the flow profile. Finally, a non-dimensional analysis of
the setup is briefly presented to understand the importance of different physical processes.
4.1 Liquid Interface
Figs. 4.1 and 4.2 show the positions of the liquid interface of the jet at different instances
in time as viewed from different viewing planes. These simulations are performed on a 15
CPD mesh with parameters Table 3.1. Complex and chaotic breakup activity comprising
primary and secondary atomization that begins approximately 6 mm away from the bell
edge is observed. The breakup processes leading to primary atomization consist of ligament
thinning, Rayleigh-Plateau instabilities and aerodynamic breakup due to interaction with
quiescent air. A more detailed analysis of the breakup mechanisms is presented in the next
chapter.
56
Figure 4.1: Liquid interface positions of a single jet exiting the edge of the bell at different times as viewed from an x-y
plane. The nozzle is shown for reference at the top left corner of each image and is not to scale.
57
(a) t = 120 µs
(b) t = 200 µs
(c) t = 320 µs
Figure 4.2: Liquid interface positions of six individual jets exiting a serrated bell edge at
different times as viewed from an x-z plane. The nozzle is rotating clockwise and is shown
for reference at the left edge of each image and is not to scale. The simulation domain
is outlined by a dashed box. Owing to the periodic boundary conditions along the z-axis,
an arc along the bell edge spanning six serrations is depicted in this view by stacking the
domains in an attempt to vaguely recreate some experimental images [42, 118, 153, 176].
58
4.2 Electric Potential Field
The electric potential field (Fig. 4.3) in the domain (NS-Mesh) ranges from the bell
potential (80 kV) to about 20 kV. In the figure, the trajectory of the liquid jet is graphically
shown as a dashed white line. The contour lines of electric potential, shown as solid lines,
are almost parallel to the y-axis in this regime. This indicates that the electric field and the
electric force vectors are oriented perpendicular to the y-axis, i.e., in the direction of lower
potential. A schematic representation of the direction of the electric force (blue line) acting
on a control volume (blue circle) is shown for clarity. It is important to note that the electric
force acts along the radially outward direction of the liquid jet in this regime as this has
significant effects on the flow development which will be discussed in the following sections.
Figure 4.3: The electric potential field and its corresponding contour lines on an x-y plane in
the domain. The trajectory of the liquid jet is shown graphically as a dashed white line. The
electric force vector (blue arrow) that acts on a control volume (blue circle) in this regime is
also schematically represented. The bell is shown at the top left corner and is not to scale.
4.3 Parameter Study
A parameter study is performed on a 15 CPD mesh to understand the effects of changing
four parameters - nozzle rotation rate, liquid flow rate, liquid charge density, and bell electric
potential. 9 separate simulations were performed totaling to a wall time of about 14,000
hours. Results of the parameter study are shown in Fig. 4.4 as snapshots of the liquid
interface after 200 µs. The standard values of all varied properties are listed in Table 3.1
59
unless otherwise stated in the figure.
As the nozzle rotation rate increases, centrifugal forces are stronger on the jet and
stretch it out faster leading to early elongation and breakup (Fig. 4.4(a)). Slower rotation
rates do not stretch the jet out as much in the same duration. A higher flow rate through
the same jet diameter acts in the same way as increasing jet velocity, which initially pushes
the jet out further before centrifugal forces take over (Fig. 4.4(b)). These results are in
agreement with experimental observations conducted by [35].
An increase in either qin or φbell stretches the jet along its downstream direction
(Figs. 4.4(c), 4.4(d)). This is because the Coulomb force vectors point towards the direction
of propagation of the liquid jet, i.e., the electric potential field contour lines are approximately
perpendicular to the liquid velocity at that location as explained in Fig. 4.3.
60
(a) Varying rotation rate
(b) Varying flow rate
(c) Varying liquid charge density
(d) Varying bell electric potential
Figure 4.4: A comparison of snapshots of the liquid interface positions after 200 µs in the
parameter study simulations. The nozzle is shown for reference at the top left corner of each
image and is not to scale.
61
4.4 Examining Non-Dimensional Quantities
The non-dimensional numbers corresponding to the simulations are listed in Table 4.1.
In the interest of understanding the effect of EHD on the flow, the relevant non-dimensional
numbers are briefly discussed. The electric Reynolds number (Ree) [162] is the ratio of the
charge advection timescale to the charge mobility timescale while the electric Péclet number
(Pee) [151] is a measure of the charge mobility timescale to the charge diffusion timescale.
For high values of Ree, charge migration is insignificant and the initial charge distribution
will prevail in the liquid jet. For low values of Ree, charge migration acts quickly to relax the
charges to the surface allowing for surface charge models to be appropriate. For high values
of Pee, charge migration dominates the charge diffusion process and vice versa. Based on
the values of Ree and Pee for this setup, it is reasonable to model all three charge dynamics
processes, i.e., advection, diffusion, and migration, since their timescales are comparable to
one another. The electric Bond number [22] quantifies the importance of the deforming
electrical force compared to the restoring surface tension force. A low value signifies the
dominance of surface tension-driven breakup activity. The electro-inertial number [151]
denotes the importance of EHD forces compared to inertial advective forces. A low value
implies that inertial forces predominantly drive the hydrodynamics of the flow in this setup.
In order to better understand the effect of electrification on atomization, a more in
depth comparative study between an electrified and a non-electrified jet was conducted.
The findings are presented in the next chapter.
62
Table 4.1: Non-dimensional numbers used in the ERBA
simulations
Number Definition Value
Density ratio ρl/ρg 830.56
CFL |umax|∆t/∆x 0.30
Bulk Weber ρl|ujet|2Djet/γ 77.87
Permittivity ratio εl/εg 50.00
Electric Reynolds (Ree) εl|ujet|/(Lqlψl) 8.93
Electric Péclet (Pe ) q ψ L2
e l l / (Dlεl) 0.21
Electric Bond εl|E|2L/γ 0.091
Electro-inertial q2
l L
2/ (εlρl|u 2
jet| ) 0.0017
63
CHAPTER FIVE
EFFECT OF ELECTRIFYING THE SETUP
An overarching question among engineers operating ERBAs in industrial applications
(particularly in automotive paint shops) is - “What is the effect of the electrification on
atomization and how can it be used to improve the performance of the device?”. Crucial
metrics such as the droplet size uniformity, surface finish quality, TE, deposition thickness,
and environmental impact are directly dependent on the atomization process of ERBAs
[9, 35, 43].
Two simulations are performed - one in a non-electrified setup and one in an electrified
setup - in a domain as described in Chapter 3 with parameter values listed in Table 3.1. The
two simulations are identical except that there is no background electric potential and no
liquid charge density in the former. The two cases are simulated to help understand the effect
of charges and the background electric field on the breakup mechanisms driving atomization.
The rest of this chapter is organized as follows - a qualitative comparison between the liquid
interfaces of the two simulated cases is followed by a quantitative comparison between the
two simulated cases where inferences are drawn based on the overall spray characteristics
and droplet statistics including size, shape, velocity, and charge density.
5.1 Liquid Interface
Figs. 5.1 and 5.2 show the positions of the liquid interface of the jet at different
instances in time as viewed from different viewing planes. Long ligaments, irregular droplet
structures, and chaotic breakup activity comprising primary and secondary atomization that
begins approximately 6 mm away from the bell edge is observed. Primary atomization is
characterized by the processes of ligament thinning, Rayleigh-Plateau instability growth
64
and aerodynamic breakup due to interaction with quiescent air. Sinuous instabilities are
observed in both cases but seem suppressed in the electrified jet. The reason for suppressed
instability growth will be discussed in the following sections. As seen in Fig. 4.3, the electric
force, which points radially outward, elongates the electrified jet along the x−direction and
enhances breakup activity.
65
Figure 5.1: Liquid interface positions in a non-electrified (yellow) and an electrified (violet) jet after exiting the edge of
the bell on an x-y plane at different times. The nozzle is shown for reference at the top left corner and is not to scale.
66
(a) t = 160 µs
(b) t = 200 µs
(c) t = 240 µs
Figure 5.2: Liquid interface positions of six individual non-electrified (yellow) and electrified
(violet) jets exiting a clockwise-rotating serrated bell edge (left, not to scale) as viewed at
different times from an x-z plane. The simulation domain is outlined by a dashed box. Owing
to the periodic boundary conditions along the z-axis, an arc along the bell edge spanning
six serrations is depicted in this view by stacking the domains in an attempt to recreate
experimental images available in Shirota et al. [153], Oswald et al. [118], Domnick [42] and
Wilson et al. [176].
67
5.2 Spray Characteristics
The spray statistics point to the observation that electrification enhances breakup
activity. The radial liquid penetration after 300 µs is 5% more for the electrified jet. The
spray angle characterizes the dispersion of the spray and is computed as the angle about
the z axis on the x-y plane out of the bell edge such that 90% of the liquid is within the
prescribed region created by that angle. The spray angle is about 3.4 times greater for the
electrified jet. The distribution of droplet spacing (Fig. 5.3), which can be quantified as the
distance between every droplet pair in the spray, further reflects the larger spray coverage in
the electrified jet where the mean droplet spacing which is 12% more than the non-electrified
case. All these metrics suggest that atomization and breakup activity is enhanced by the
electric field. As seen in Fig. 4.3, the electric force points radially outward and increases
penetration and spray coverage in electrified flows. Droplet spacing is an important metric
in rotary sprays and the electric field has a significant impact on it. The optimal operating
conditions of the electrified atomizer have to be carefully selected.
Figure 5.3: Log-normal fit of the near-bell droplet spacing probabilty density.
68
In contrast, the breakup length, i.e., the length of the main ligament core, is 4%
longer for the electrified jet. Moreover, the first instance of breakup occurs 6% later in the
electrified jet. This can be attributed to the stabilizing force provided by the electric field
on the charged liquid core [107, 112, 144]. Saville [143] found that under certain conditions,
electrical shearing forces can stabilize the cylindrical interface to axisymmetric deformations.
In another study conducted by Bhuptani and Sathian [23], the axial electric field has been
quantitatively shown to have a stabilizing effect on nanoscale liquid threads. However, the
next section contains evidence that once breakup begins, the electric field accelerates the
processes leading to secondary atomization.
5.3 Atomization Characteristics
The droplet ancestry statistics extraction tool described in Section 3.4 operates
throughout the simulation and connects successive breakup/coalescence events to create an
ancestry or family tree describing the evolution of the spray cloud. Using this tool, several
atomization statistics related to the droplet size, shape, velocity, and charge density are
extracted and analyzed to gain insights into the atomization process in ERBAs. The results
are compiled after analyzing 837 liquid structures in the non-electrified setup and 1119 liquid
structures in the electrified setup.
According to Pendar and Páscoa [128], the presence of electric charges within liquid
structures causes a net reduction in the effective surface tension in large droplets and
ligaments and therefore an increase in their effective local Weber number. This leads
to enhanced breakup activity. After 300 µs, the electrified jet contains about 25% more
droplets than the non-electrified jet. Fig. 5.4 shows the time evolution of the number of
liquid structures in the domain. The delay in the first instance of breakup (discussed in
the previous section) and the accelerated breakup in the electrified jet are evident from the
plot. A closer look at the droplet statistics, which will be presented in the following sections,
69
provides more evidence for the enhanced atomization in an electrified jet.
Figure 5.4: The total number of droplets in the domain over time. Note the logarithmic
scale on the y-axis.
Two major types of atomization are identified in this system - primary and secondary.
Primary atomization occurs when droplets split from the main ligament core. Secondary
atomization comprises all breakup following primary atomization which includes second,
third, fourth, etc. breakup events of the same structure. These subsequent breakup events
are also called breakup levels in the following discussion. For clarity, breakup level 1
corresponds to liquid structures that are formed due to primary atomization, i.e., breakup
from the main ligament core; breakup levels 2 onward correspond to liquid structures created
from secondary atomization events. In a non-electrified jet, 1.64% of the total number of
droplets formed in the simulation are from primary atomization (broke up once) and 98.36%
are from secondary atomization (broke up two or more times). The prevalence of secondary
atomizing events indicates that the droplet cloud is heavily dependent on mechanisms within
this regime. The dominance is more pronounced in an electrified jet where the ratio of
secondary to primary droplets is 30% more than that in the non-electrified jet. This follows
70
the idea that the electric field accelerates secondary breakup activity. Various parameters
including the rotation rate, surface tension coefficient, and liquid viscosity affect the distance
and time at which primary atomization begins.
It is interesting to analyze the number of droplets for each breakup level. While the final
time chosen is arbitrary, it confirms the trend of atomization activity in the domain. Fig. 5.5
shows the number distribution of the droplets for each breakup level in the simulation.
Although the dominant breakup level for both cases is the same (level 4), the electrified jet,
which contains more droplets overall, also contains more droplets in deeper breakup levels.
Figure 5.5: The number distribution of droplets that have undergone different breakups
levels.
Droplet coalescence is a crucial process in atomizing flows. After 300 µs, the ratio
of coalescence to breakup events in electrified operation (0.27) is 41% more than that in a
non-electrified setup (0.19). It makes intuitive sense that more coalescence is observed in
electrified setups because although the droplets are more spread out, the spray cloud contains
more structures which in turn leads to more droplet collisions.
71
5.4 Droplet Size and Shape Evolution
The enhanced breakup activity in an electrified jet produces smaller droplets. The mean
diameter of droplets in the spray cloud near the nozzle is about 4.5% smaller when operating
in an electric field. Fig. 5.6 shows a comparison of the diameter distribution of droplets in
the domain. A higher number of smaller liquid structures are created in an electrified jet
due to the net reduction in effective surface tension in charged droplets.
Figure 5.6: Droplet diameter distribution. Note the logarithmic scale on the y-axis.
Fig. 5.7 shows the comparison of the mean droplet diameters for each breakup level. A
trend towards decreasing droplet sizes for deeper breakup levels is generally observed. Due
to coalescence events, the droplet sizes may not strictly decrease between successive breakup
levels.
Fig. 5.8 shows the mean change in diameter of droplets between successive breakup
levels. The change in the size of droplets, which is the difference between the new and
parent droplet diameters, decreases for deeper breakup levels due to the small sizes of the
droplets themselves. The mean change in diameters across all breakup levels is about 18%
72
Figure 5.7: Mean droplet diameter for each breakup level.
lower in electrified jets. There are some instances where the change in diameter during a
breakup event is a positive number. This is due to a coalescence event that occurred between
two breakup events leading to an increase in its size. After breakup level 9, the mean change
in diameter across breakup levels is seemingly constant for both cases. This is speculated to
be due to a high ratio of coalescence to breakup events in this regime.
One way to quantify the shape of the resulting droplet is to examine its sphericity
- a dimensionless quantity that represents the degree of departure from a perfect sphere.
Mathematically, sphericity can be calculated as A1/A3, where A1 ≥ A2 ≥ A3 are the
lengths of the three principal axes of an ellipsoidal approximation of the droplet structure.
A sphericity of 1, which is its lowest attainable value, indicates a perfect sphere, while larger
values indicate irregular shapes. Fig. 5.9 shows a comparison of the sphericity distribution of
droplets in the domain. The mean sphericity of all droplet structures is about 5% lower for
the electrified droplets after 300 µs of the simulation, suggesting that they are more spherical
in shape.
73
Figure 5.8: Mean change in droplet diameters for successive breakup levels. Here, breakup
level 1 is omitted since it corresponds to breakup from the main ligament core.
Figure 5.9: Droplet sphericity distribution.
74
5.5 Droplet Velocity Statistics
Analyzing the droplet velocities shows that during breakup, the mean absolute velocity
of droplets in the electrified jet is about 3% lower than the non-electrified jet during breakup
events. Fig. 5.10 shows that smaller droplets undergo breakup at lower velocities in the
electrified jet due to the net reduction in the effective surface tension force. However, larger
droplets undergo breakup at greater velocities in the electrified jet. This can be attributed
to the presence of charges which imparts an electric force on large droplets. A similar trend
is not evident in smaller droplets due to their small size and sensitivity to aerodynamic drag.
Figure 5.10: Comparing the mean velocity for different droplet sizes.
Additionally, the mean slip velocity magnitude, which is formulated as |uliquid−ugas|, is
16% higher in an electrified setup during breakup events. This suggests that charged droplets
undergo breakup at higher slip velocities. Fig. 5.11 highlights that this is true across most
droplets sizes. |uliquid| is computed as the liquid volume weighted average velocity of cells
corresponding to that structure while |ugas| is computed as the gas volume weighted average
velocity of cells neighboring the interfacial cells of that structure. It is important to note
here that this method of calculating the slip velocity is a severe underestimation. Ideally,
75
the slip velocity is calculated as the local gas velocity without the influence of velocity of
the structure itself. While its actual value is challenging to compute, our observation that
its magnitude is higher in an electrified jet is a useful statistic which will be discussed in the
next section.
Figure 5.11: Comparing the mean slip velocity for different droplet sizes.
5.6 Local Weber Number Analysis
Using the ancestry extraction tool, the local flow field surrounding breakup events is
analyzed and the local Weber number (We local) is computed as
ρgU
2
sLWe local = (5.1)
γ
where ρg is the gas density, Us is the slip velocity, L is the characteristic length which is the
equivalent spherical diameter of the droplet prior to breakup, and γ is the surface tension
coefficient. A more detailed explanation of this calculation is provided in Christensen and
Owkes [31].
76
The We local distribution corresponding to structures prior to breakup is presented in
Fig. 5.12. The mean Weber number for all breakup events is about 25% more in an electrified
jet, which likely indicates that aerodynamic forces are more actively driving atomization in
this system. This follows the findings of the slip velocity which was higher during breakup
in electrified droplets. To reiterate, the slip velocity, which is severely underestimated, is
squared in Eq. 5.1 to calculate the We local. While the accurate value of We local is challenging
to compute, the general trend that it is higher in electrified setups is a useful finding.
Since the We local is severely underestimated in this formulation, there is sufficient reason
to believe that atomization is driven by both aerodynamic breakup and instability growth
and is enhanced in electrified setups by the electric force and the presence of charges in the
liquid.
Figure 5.12: Local Weber number distribution associated with all structures. Note the
logarithmic scale on the y-axis.
77
5.7 Charge Distribution Statistics
Most droplets contain a charge density that is very close to the initial bulk charging
value of 2.879 C/m3 (Fig. 5.13). However, since the process of charge migration is included
in the model, we also see a distribution of charge densities in the droplets. Compared to the
initial charge density (Table 3.1), it is more probable to find a droplet with a lower value
of q. As charges relax towards the interface of the ligament, there is a localized reduction
of q in the bulk of the ligament. A majority of the liquid structures break off the ligament
bulk (since a greater liquid volume is present in the bulk than near the interface) to form
droplets that contain slightly lower values of q than the initial bulk value. The moderate
value of Ree (Table 4.1) which is the ratio of the residence time of liquid in the domain to
the charge relaxation timescale indicates that the charges have time to relax to some degree,
but not enough time to fully relax to the surface. Higher resolution simulations can help
better understand the charge relaxation process through the ligament core and the factors
affecting the charge distribution in the spray cloud.
Figure 5.13: Number distribution of q in liquid structres in the domain. The inital bulk
value of q is shown as a dashed line. Note the logarithmic scale on the y-axis.
78
Upon closer inspection of the charge distribution in droplets, it is evident that q is
generally higher in larger droplets (Fig. 5.14). This can be explained as follows - small
droplets are more likely to be formed from the volume of liquid in the bulk of a ligament
which, as discussed, contains a lower charge density. Larger droplets tend to break off as
long primary ligaments from the main ligament core and therefore will not necessarily have
a lower value of q.
Figure 5.14: Mean q in liquid structures of different sizes. The inital bulk value of q is shown
as a dashed line.
The change in q (∆q) for successive breakup levels also shows a decreasing trend and
approaches zero for the deepest breakup levels (Fig. 5.15). Naturally, the change in q during
a breakup event is inversely proportional to the change in the droplet sizes (Fig. 5.16).
79
Figure 5.15: The mean change in q between successive breakup levels.
Figure 5.16: The mean change in q corresponding to the change in diameter during a breakup
event.
80
CHAPTER SIX
EFFECT OF SHEAR-THINNING FLOWS
In the interest of understanding breakup physics, a similar droplet ancestry focused
study was conducted to isolate the effects of non-Newtonian flow behavior in the application
of rotary atomization. The shear-thinning property of a fluid can significantly affect its
atomization process. Deviations from the a uniform droplet size distribution and spray
pattern may lead to reduced coating quality and efficiency [35]. It is important to understand
the effects of non-Newtonian behavior on atomizing flows to achieve the desired coating
thickness and surface finish quality. By isolating and understanding these effects, more
reliable models can be developed to the predict optimal operation of rotary atomizers. Two
simulations are performed - one with a Newtonian liquid and one with a non-Newtonian,
particularly shear-thinning liquid - in a domain as described in Chapter 3 with parameter
values listed in Table 3.1. However, both simulations are conducted in the absence of a
background electric field and zero liquid charge density. The two simulations are identical
except that the former has a constant liquid viscosity of 0.1 Pa.s. The goal of this study is to
isolate and understand non-Newtonian flow effects on the breakup mechanisms. For brevity,
only the notable inferences are presented in this section.
6.1 Liquid Interface
Fig. 6.1 shows the positions of the liquid interface of the jet at different instances in time
as viewed on an x-y plane. The sprays look visually similar except for the jet penetration.
A detailed comparison based on the extracted statistics follows.
81
Figure 6.1: Liquid interface positions of a Newtonian (pink) and a non-Newtonian (green) jet at different times as viewed
from an x-y plane. The nozzle is shown for reference at the top left corner and is not to scale.
82
6.2 Atomization Characteristics
The droplet ancestry statistics extraction tool described in Section 3.4 is employed to
extract several atomization statistics related to the droplet size, shape, velocity, and charge
density can be extracted and analyzed to gain insights into the effect of shear-thinning
behavior of liquid on the atomization process in RBAs. The results are compiled after
analyzing 796 liquid structures in the Newtonian domain and 897 liquid structures in the
non-Newtonian setup.
Fig. 6.2 shows the time evolution of the number of liquid structures in the domain. The
non-Newtonian jet seemingly undergoes breakup sooner and contains 12% more droplets than
the Newtonian jet. While secondary atomization is still the dominant regime for breakup,
the ratio of secondary to primary droplets in the Newtonian jet is 23% more than that in
the non-Newtonian jet.
Figure 6.2: The total number of droplets in the domain over time in Newtonian (pink) and
non-Newtonian (green) operating conditions in the first 300 µs of the near-bell atomization
simulation. Note the logarithmic scale on the y-axis.
The non-Newtonian jet has an effectively larger value of viscosity since it’s the same at
83
high shear rates, but larger at low shear rates. Therefore, the higher viscosity at low shear
rates seemingly suppresses instabilities that cause secondary atomization. However, the
initial breakup occurs sooner in the non-Newtonian case indicating that the higher viscosity
has the opposite effect on primary atomization.
Fig. 6.3 shows the number distribution of the droplets for each breakup level in the
simulation. Although the dominant breakup level for both cases is the same (level 4), the
non-Newtonian jet, which contains more droplets overall, has a smaller spread in terms of
the number distribution of droplets across breakup levels. The deepest breakup level in the
Newtonian jet is 18 while it is only 14 for the non-Newtonian case. This further supports
the idea that the shear-thinning behavior of the liquid is suppressing secondary atomization.
After 300 µs, the ratio of coalescence to breakup events in the non-Newtonian setup is 8%
more than that in the Newtonian case suggesting that more coalescence events occur in the
non-Newtonian perhaps due to the increased number of droplets in the domain.
Figure 6.3: The number distribution of droplets that have undergone different breakups levels
in Newtonian (pink) and non-Newtonian (green) operating conditions in the first 300 µs of
the near-bell atomization simulation.
84
6.3 Droplet Size Evolution
Fig. 6.4 shows a comparison of the statistically similar diameter distribution of droplets
in the domain and the similarity is evident. This indicates that the viscosity does not have
a significant impact on the drop size distribution. The drop sizes are likely affected mainly
by the aerodynamic and surface tension forces, and not the flow within the liquid. Fig. 6.5
shows the comparison of the mean droplet diameters for each breakup level. Naturally,
a trend towards decreasing droplet sizes for deeper breakup levels is observed. However,
discrepancies in the behavior between the two cases is non-evident, reiterating that viscosity
is not playing a large role in the drop sizes. Fig. 6.6 shows the mean change in diameter of
droplets between successive breakup levels. The change in the size of droplets decreases for
deeper breakup levels due to the small sizes of the droplets themselves. Sometimes, positive
diameter changes are measured due to coalescence events, which increases the droplet size.
Interestingly, the non-Newtonian jet has a larger magnitude of the mean change in diameters
across all breakup levels.
Figure 6.4: Droplet diameter distribution in the simulation.
85
Figure 6.5: Mean droplet diameter for each breakup level in the simulation.
Figure 6.6: Mean change in droplet diameters for successive breakup levels. Here, breakup
level 1 is omitted since it corresponds to breakup from the main ligament core.
86
6.4 Viscosity Characteristics
Figure 6.7 highlights that the viscosity in an overwhelming majority of the droplets is
very close to the viscosity value for an infinite shear rate which is approximately 0.111 Pa.s
(Fig. 2.1). Higher resolution simulations can help better understand the viscosity distribution
within the ligament core and the individual droplet structures.
The mean viscosity also shows a decreasing trend during breakup into deeper breakup
levels (Fig. 6.8) which suggests that droplets breaking off at deeper breakup levels experience
higher shear rates. Consequentially, upon closer inspection of the viscosity in droplets, it
is evident that larger droplets are likely to have a higher viscosity value (Fig. 6.9). This is
speculated to be due to the fact that larger droplets undergo higher shear by virtue of their
higher surface area, whereas small droplets move with the gas flow. The mean viscosity in
droplets during breakup increases with the local Weber number (Fig. 6.10). The change
in viscosity during breakup for successive breakup levels also shows a decreasing trend and
approaches zero for the deepest breakup levels (Fig. 6.11).
87
Figure 6.7: Liquid viscosity distribution during breakup events in the first 300 µs of the
simulation.
Figure 6.8: Mean viscosity of the parent droplets for each breakup level.
88
Figure 6.9: Mean viscosity of the parent droplets as a function of droplet sizes.
Figure 6.10: Comparing liquid viscosity for a range of local Weber numbers during breakup
events in the first 300 µs of the simulation.
89
Figure 6.11: The mean change in liquid viscosity between successive breakup levels.
90
CHAPTER SEVEN
CONCLUSIONS AND OUTLOOK
This dissertation contains a detailed discussion regarding the use of numerical methods
to simulate and study gas-liquid multiphase flows, particularly atomizing flows in industrial
applications. A device called the electrostatic rotary bell atomizer (ERBA) is simulated
in this work. ERBAs are popular atomizing devices used in various industries, including
automotive painting. Multiple complex physical processes influence the mechanisms driving
atomization in the device. The simulations are characterized by an interesting set of physical
models -
1. The multiphase gas-liquid flow is modeled using the Navier-Stokes equations, a surface
tension force formulation and an interface tracking scheme.
2. The non-Newtonian flow is formulated using a modified power law fit on empirical
viscosity values corresponding to shear-thinning automotive paints.
3. The rotating flow is modeled using the centrifugal and Coriolis forces.
4. The electrohydrodynamic processes are formulated using a set of governing equations
describing the electric force, electric field, electric potential field, and charge transport
(advection, diffusion, and migration). An auxiliary domain called e-Mesh is used
exclusively as an electric potential field solver to obtain accurate boundary conditions
of electric potential on the flow domain.
Notably, a multi-physics based approach has been incorporated in this work to accurately
model the complex atomizing flows in ERBAs.
91
7.1 Atomization Characteristics in Rotary Atomization
In this work, high-fidelity simulations are performed using a numerical model developed
to simulate electrohydrodynamic atomization in ERBAs. The simulations are equipped with
a droplet ancestry extraction tool that stores droplet and local flow statistics. The extracted
statistics are analyzed to provide insights into the breakup mechanisms in an ERBA. This is a
cost-effective method to understand atomization and extract information that is challenging
to obtain experimentally.
It is found that the presence of the electric field increases jet penetration, spray angle,
and droplet spacing due to the effect of the electric force that acts along the radially outward
direction. However, the breakup length and the first instance of breakup in time are longer
in the electrified jet due to a stabilizing force provided by the electric field on the charged
ligament core. However, once breakup begins, the net reduction in the effective surface
tension in liquid structures due to the presence of charges accelerates the breakup process
leading enhanced secondary atomization in the electrified jet. The total number of droplets,
the ratio of secondary to primary droplets, and the ratio of coalescence to breakup activity
are all higher when operating in an electric field. Charged droplets seem to undergo breakup
at higher slip velocities than uncharged droplets. Most charged droplets contain a charge
density that is very close to the initial bulk charging value. While larger droplets generally
contain a higher value of charge density, it is more probable to find droplets with charge
density values lower than the initial bulk value.
In a separate study conducted in a non-electrified setup to identify the effect of the shear-
thinning behavior of liquid paints, it was found that viscosity does not have a significant
impact on the drop size distribution. The drop sizes are likely affected mainly by the
aerodynamic and surface tension forces. The higher viscosity at low shear rates seemingly
suppresses instabilities that cause secondary atomization. However, the initial breakup
92
occurs sooner in the non-Newtonian case indicating that the higher viscosity has the opposite
effect on primary atomization.
7.2 Impact of Research Work
Accurate simulations of multiphase flows in engineering applications involving multiple
physical phenomena can help create a significant impact on our understanding and the
optimization of the processes. Simulations will allow researchers to isolate certain aspects
of the flow and identify their effects on the relevant metrics. Simulations also help visualize
and measure flows in ways arguably impossible with current experimental methods. In this
work, information which is inaccessible to experimental techniques near the bell edge and
within the dense droplet cloud has been extracted and analyzed to study the behavior of
charged droplets under the influence of a background electric field. Small improvements to
the performance of this device can have significant financial and environmental benefits. The
findings of this work can help manufacturers choose optimal operating conditions of ERBAs
more carefully.
7.3 Future Work
The numerical formulation and methodology presented in this dissertation provide a
framework that can predict the atomization process in near-bell ERBA flows. However,
some aspects in framework and the numerical setup that could benefit from improvements
or inclusions are listed below.
1. A parameter study can be conducted to investigate the effect of various operating
parameters on the physics of the atomization process.
(a) The rotation rate strongly affects the droplet size distribution, spray cone angle,
and surface coverage.
93
(b) The flow rate affects the droplet size distribution and the coating thickness.
(c) Liquid viscosity plays an important role in the physics of atomization in how it
governs the resistance of flow to shear in atomizing flows. Higher viscosity fluids
may have reduced spray coverage and penetration compared to lower viscosity
fluids. Liquid viscosity can also impact the surface wetting and coating uniformity.
(d) As seen from the work done in this thesis, the liquid charge density and
background electric potential can have a significant impact on the atomization
quality in multiple ways.
(e) The size of the bell cup and its distance from target surface can affect the spray
pattern in certain regimes of operating conditions.
2. In this work, liquid was injected into the domain through an ellipsoidal inflow with a
velocity profile that corresponds to the edge angle of the bell. However, in practice, a
thin film of liquid flows out of a flat/serrated bell edge in a more complex way to form
ligaments.
(a) It would be more accurate to simulate a film flowing through serrations at the
bell edge instead of modeling an ellipsoidal inflow.
(b) Different serration geometries and widths can be simulated to better understand
the relationship between droplet sizes and bell edge geometries.
(c) In a bell without serrations, sheet breakup is the primary mode of atomization.
There exist formulations currently that can model sheet breakup for atomizing
flows. It would be interesting to explore this mode of breakup in ERBAs.
3. The initial charge distribution was assumed to be uniform through the liquid bulk
in this work. However, it would be worth exploring a more practical way to impart
charge into the liquid phase to improve the accuracy of the predicted atomizing flows
[3, 30, 109]. There are two ways of charging paint before it exits the nozzle - 1) direct
charging - where the liquid is charged by applying high voltage directly to the rotating
94
bell of the atomizer [44], and 2) external charging - where droplets are charged using
external high voltage emitting electrode needles due to free ions produced from corona
discharge at the electrodes [2].
4. A “future work” section in a CFD-based research document is incomplete without
the mention of potentially performing simulations on a finer mesh. Higher resolution
studies will shed more light on the process of charge redistribution through the ligament
core and within the droplet structures.
5. While it is well-known that the presence of charges in a droplet reduces its effective
surface tension, a quantitative analysis will be useful. In theory, the electric field,
which generates an electrostatic pressure at the interface, alters the distribution of the
surface charges and affects the surface tension force. An inclusion to the statistics
being extracted could be the electric field magnitude in each structure. This value
can be used to compute the electrostatic pressure which can help quantify the effective
reduction in the surface tension force due to the presence of charges in the droplets.
6. A feature commonly found in ERBA setups is the shaping air system which is used
as a secondary air stream to shape and control the spray pattern. It will be worth
performing simulations that include a shaping air setup to make the tool more robust
in simulating industrial applications of the device.
7. As an extension to this work, a Lagrangian droplet spray model can be implemented for
droplets that do not undergo any further breakup. Once such droplets are identified
using appropriate breakup models, the Lagrangian droplet spray model can predict
the droplet trajectories, interaction with airflow, and the deposition pattern. In the
proposed model, each droplet would be treated as a discrete particle with its own
position, velocity, charge and other relevant properties. A strong advantage to this
reduced physics model is the elimination of the necessity to solve expensive pressure and
electric Poisson equations and focus only on the centrifugal, Coriolis, aerodynamic drag
95
and electric forces on each discrete droplet structure. Moreover, it would enable the
tool to completely model the spraying process from high-fidelity injection to deposition
on the target surface - which would make it the first of its kind.
96
REFERENCES CITED
97
[1] The influence of air and liquid properties on airblast atomization. Journal of Fluids
Engineering-transactions of the ASME, 97(3):316–320, 1975.
[2] The simulation of electrostatic spray painting process with high-speed rotary bell
atomizers. part ii: External charging. Particle & Particle Systems Characterization,
23(5):408–416, 2006.
[3] Effects of electrode voltage, liquid flow rate, and liquid properties on spray charge-
ability of an air-assisted electrostatic-induction spray-charging system. Journal of
Electrostatics, 68(2):152–158, 2010.
[4] Direct Numerical Simulations of Gas-Liquid Multiphase Flows. 2011.
[5] Michael B Abbott and David R Basco. Computational fluid dynamics. Longman
Scientific & Technical Harlow, Essex, 1989.
[6] J. Abu-Ali and S.A. Barringer. Method for electrostatic atomization of emulsions in
an ehd system. Journal of Electrostatics, 63(5):361–369, 2005.
[7] Adnan Darwish Ahmad, Binit B. Singh, Mark Doerre, Ahmad M. Abubaker, Masoud
Arabghahestani, Ahmad A. Salaimeh, and Nelson K. Akafuah. Spatial positioning and
operating parameters of a rotary bell sprayer: 3D mapping of droplet size distributions.
Fluids, 4(3), 2019. ISSN 23115521.
[8] N. K. Akafuah, Kimio Toda, Abraham Salazar, and Kozo Saito. Automotive Paint
Spray Characterization and Visualization, pages 121–165. Springer Netherlands,
Dordrecht, 2013.
[9] Nelson Akafuah, Sadegh Poozesh, Ahmad Salaimeh, Gabriela Patrick, Kevin Lawler,
and Kozo Saito. Evolution of the Automotive Body Coating Process — A Review.
Coatings, 6(2), 2016. ISSN 2079-6412.
[10] Abdullah Alotaibi, Mohmoud Awad, and Ahmed Hamed. Performance of air
conditioning system using air cooled condenser with water atomization. International
Journal of Engineering Research and Science and Technology, 4(1):127–135, 2015.
[11] John David Anderson and John Wendt. Computational fluid dynamics, volume 206.
Springer, 1995.
[12] Björn Andersson, Valeri Golovitchev, Stefan Jakobsson, Andreas Mark, Fredrik
Edelvik, Lars Davidson, and Johan S Carlson. A modified tab model for simulation
of atomization in rotary bell spray painting. Journal of Mechanical Engineering and
Automation, 3(2):54–61, 2013.
[13] Ricardo D Andrade, Olivier Skurtys, and Fernando A Osorio. Atomizing spray systems
for application of edible coatings. Comprehensive Reviews in Food Science and Food
Safety, 11(3):323–337, 2012.
98
[14] T. Anukiruthika, J.A. Moses, and C. Anandharamakrishnan. Electrohydrodynamic
drying of foods: Principle, applications, and prospects. Journal of Food Engineering,
295:110449, 2021. ISSN 0260-8774.
[15] ZM Anwar, Ee Sann Tan, R Adnan, and Mohd Azree Idris. Study on atomization
characteristics for power generation application. In IOP Conference Series: Earth and
Environmental Science, volume 16, page 012016. IOP Publishing, 2013.
[16] Pooya Azizian, Milad Azarmanesh, Morteza Dejam, Mehdi Mohammadi, Milad
Shamsi, Amir Sanati-Nezhad, and Abdulmajeed A. Mohamad. Electrohydrodynamic
formation of single and double emulsions for low interfacial tension multiphase systems
within microfluidics. Chemical Engineering Science, 195:201–207, 2019. ISSN 0009-
2509.
[17] WD Bachalo. Experimental methods in multiphase flows. International journal of
multiphase flow, 20:261–295, 1994.
[18] Adrian G. Bailey. Electrostatic atomization of liquids. Science Progress (1933- ), 61
(244):555–581, 1974. ISSN 00368504, 20477163.
[19] S Balachandar and John K Eaton. Turbulent dispersed multiphase flow. Annual review
of fluid mechanics, 42:111–133, 2010.
[20] W. Balachandran and A. G. Bailey. The Dispersion of Liquids Using Centrifugal and
Electrostatic Forces. IEEE Transactions on Industry Applications, IA-20(3):682–686,
1984.
[21] Marcos AS Barrozo, Claudio R Duarte, Norman Epstein, John R Grace, and C Jim
Lim. Experimental and computational fluid dynamics study of dense-phase, transition
region, and dilute-phase spouting. Industrial & Engineering Chemistry Research, 49
(11):5102–5109, 2010.
[22] Jean Berthier and Kenneth A Brakke. The physics of microdroplets. John Wiley &
Sons, 2012.
[23] Darshak K. Bhuptani and Sarith P. Sathian. Effect of axial electric field on the Rayleigh
instability at small length scales. Physical review. E, 95(5-1):053115, may 2017. ISSN
24700053.
[24] G. Blanquart, P. Pepiot-Desjardins, and H. Pitsch. Chemical mechanism for high
temperature combustion of engine relevant fuels with emphasis on soot precursors.
Combustion and Flame, 156(3):588–607, 2009. ISSN 00102180.
[25] Victor Boniou, Thomas Schmitt, and Aymeric Vié. Comparison of interface capturing
methods for the simulation of two-phase flow in a unified low-mach framework.
International Journal of Multiphase Flow, 149:103957, 2022.
99
[26] JP Brill and SJ Arirachakaran. State of the art in multiphase flow. Journal of
Petroleum Technology, 44(05):538–541, 1992.
[27] Antonio Castellanos. Electrospray and atomization. In Electrohydrodynamics, pages
206–212. Springer, 1998.
[28] Fang Chen and Hans Hagen. A survey of interface tracking methods in multi-phase
fluid visualization. In Visualization of Large and Unstructured Data Sets-Applications
in Geospatial Planning, Modeling and Engineering (IRTG 1131 Workshop). Schloss
Dagstuhl-Leibniz-Zentrum fuer Informatik, 2011.
[29] Zhangxin Chen, Guanren Huan, and Yuanle Ma. Computational methods for
multiphase flows in porous media. SIAM, 2006.
[30] F Chowdhury, M Ray, A Sowinski, P Mehrani, and A Passalacqua. A review on
modeling approaches for the electrostatic charging of particles. Powder Technology,
389:104–118, 2021.
[31] Brendan Christensen and Mark Owkes. Efficient extraction of atomization processes
from high-fidelity simulations. Computers & Fluids, 254:105808, 2023. ISSN 0045-7930.
[32] Zémerli Clément, Edelvik Fredrik, Mark Andreas, and Hermanns Oliver. Efficient
Numerical Simulation of Spray Painting Processes in Automotive Manufacturing. In
23rd SAE Brasil International Congress and Display. SAE Technical Paper, 2014.
[33] S. A. Colbert and R. A. Cairncross. A computer simulation for predicting electrostatic
spray coating patterns. Powder Technology, 151:77–86, 2005. ISSN 0032-5910.
[34] Steven Anthony Colbert. Numerical simulations of droplet trajectories from an
electrostatic rotary-bell atomizer. PhD thesis, Drexel University, 2007.
[35] P. L. Corbeels, Dwight W. Senser, and Arthur H. Lefebvre. Atomization Characteristics
of a high speed rotary-bell paint applicator. Atomization and Sprays, 2(2):87–99, 1992.
ISSN 1044-5110.
[36] Ian P. Craig, Andrew Hewitt, and Howard Terry. Rotary atomiser design requirements
for optimum pesticide application efficiency. Crop Protection, 66:34–39, 2014. ISSN
0261-2194.
[37] Joseph M. Crowley. Fundamentals of applied electrostatics, volume 1. Wiley, 1st
edition, 1986.
[38] ZG Deng, Rui Xiao, BS Jin, QL Song, and He Huang. Multiphase cfd modeling for
a chemical looping combustion process (fuel reactor). Chemical Engineering & Tech-
nology: Industrial Chemistry-Plant Equipment-Process Engineering-Biotechnology, 31
(12):1754–1766, 2008.
100
[39] Olivier Desjardins, Guillaume Blanquart, Guillaume Balarac, and Heinz Pitsch. High
order conservative finite difference scheme for variable density low mach number
turbulent flows. Journal of Computational Physics, 227:7125–7159, 2008. ISSN 0021-
9991.
[40] Olivier Desjardins, Vincent Moureau, and Heinz Pitsch. An accurate conservative level
set/ghost fluid method for simulating turbulent atomization. Journal of Computational
Physics, 227:8395–8416, 2008. ISSN 0021-9991.
[41] N. Dombrowski and T.L. Lloyd. Atomisation of liquids by spinning cups. The Chemical
Engineering Journal, 8(1):63–81, 1974. ISSN 0300-9467.
[42] Joachim Domnick. Effect of Bell Geometry in High-Speed Rotary Bell Atomization.
2010.
[43] Joachim Domnick and M. Thieme. Atomization Characteristics of High-Speed Rotary
Bell Atomizers. Atomization and Sprays, 16(8):857–874, 2006. ISSN 1044-5110.
[44] Joachim Domnick, Andreas Scheibe, and Qiaoyan Ye. The Simulation of the
Electrostatic Spray Painting Process with High-Speed Rotary Bell Atomizers. Part I:
Direct Charging. Particle & Particle Systems Characterization, 22(2):141–150, 2005.
[45] Joachim Domnick, Z. Yang, and Q. Ye. Simulation of the film formation at a high-speed
rotary bell atomizer used in automotive spray painting processes. 2008.
[46] Peter A. M. Eagles, Amer N. Qureshi, and Suwan N. Jayasinghe. Electrohydrodynamic
jetting of mouse neuronal cells. Biochemical Journal, 394(2):375–378, 02 2006. ISSN
0264-6021.
[47] MM Elkotb. Fuel atomization for spray modelling. Progress in Energy and Combustion
Science, 8(1):61–91, 1982.
[48] Kevin R.J. Ellwood and J. Braslaw. A finite-element model for an electrostatic bell
sprayer. Journal of Electrostatics, 45(1):1–23, 1998. ISSN 03043886.
[49] Kevin R.J. Ellwood, Janice L. Tardiff, and Seyed M. Alaie. A simplified analysis
method for correlating rotary atomizer performance on droplet size and coating
appearance. Journal of Coatings Technology and Research, 11(3):303–309, 2014. ISSN
15470091.
[50] Marjan Enayati, Ming-Wei Chang, Felix Bragman, Mohan Edirisinghe, and Eleanor
Stride. Electrohydrodynamic preparation of particles, capsules and bubbles for
biomedical engineering applications. Colloids and Surfaces A: Physicochemical and
Engineering Aspects, 382(1-3):154–164, 2011.
[51] R D Falgout, J E Jones, and U M. Yang. Pursuing Scalability for hypre’s Conceptual
Interfaces. Lawrence Livermore National Laboratory, 2004.
101
[52] Robert D. Falgout, Jim E. Jones, and Ulrike Meier Yang. The design and
implementation of hypre, a library of parallel high performance preconditioners. In
Numerical Solution of Partial Differential Equations on Parallel Computers, pages
267–294. Springer Berlin Heidelberg, 2006.
[53] Hua-Tzu Fan, Harry Kuo, and Joseph Simmer. Measuring paint droplet size, velocity,
and charge-to-mass ratio distribution for electrostatic rotary bell spray simulation. In
ASME International Mechanical Engineering Congress and Exposition, volume 54921,
pages 703–709, 2011.
[54] U. Farook, E. Stride, M. J. Edirisinghe, and R. Moaleji. Microbubbling by co-axial
electrohydrodynamic atomization. Medical & Biological Engineering & Computing, 45
(8):781–789, 2007.
[55] Ronald P Fedkiw, Tariq Aslam, Barry Merriman, and Stanley Osher. A non-oscillatory
eulerian approach to interfaces in multimaterial flows (the ghost fluid method). Journal
of computational physics, 152(2):457–492, 1999.
[56] John B Fenn, Matthias Mann, Chin Kai Meng, Shek Fu Wong, and Craig M
Whitehouse. Electrospray ionization for mass spectrometry of large biomolecules.
Science, 246(4926):64–71, 1989.
[57] Daniel Fredrich, Erik Weiand, and Andrea Giusti. Electrostatic fields for the control
of evaporating charged fuel sprays. International Journal of Multiphase Flow, page
104312, 2022. ISSN 03019322.
[58] A. R. Frost. Rotary atomization in the ligament formation mode. Journal of
Agricultural Engineering Research, 26:63–78, 1981. ISSN 0021-8634.
[59] C. Galitsky and E. Worrel. Energy Efficiency Improvement and Cost Saving
Opportunities for the Vehicle Assembly Industry: An Energy Star Guide for Energy
and Plant Managers. Technical report, Lawrence Berkeley National Laboratory,
University of California: Berkeley, CA, USA, December 2008.
[60] CA Geffen and S Rothenberg. Suppliers and environmental innovation: the automotive
paint process. International Journal of Operations and Production Management, 2000.
[61] Anisha V Gorty and Sheryl A Barringer. Electrohydrodynamic spraying of chocolate.
Journal of Food Processing and Preservation, 35(4):542–549, 2011.
[62] J.M. Grace and J.C.M. Marijnissen. A review of liquid atomization by electrical means.
Journal of Aerosol Science, 25(6):1005–1019, 1994. ISSN 0021-8502.
[63] Holger Grosshans and Miltiadis V Papalexandris. A model for the non-uniform contact
charging of particles. Powder Technology, 305:518–527, 2017.
102
[64] Nico Guettler, Philipp Knee, Qiaoyan Ye, and Oliver Tiedje. Initial droplet conditions
in numerical spray painting by electrostatic rotary bell sprayers: A framework for
optimization of injection model coefficients. Journal of Coatings Technology and
Research, 17(5):1091–1104, 2020.
[65] Lutz Gödeke, Walter Oswald, Norbert Willenbacher, and Peter Ehrhard. Dimensional
analysis of droplet size and ligament length during high-speed rotary bell atomization.
Journal of Coatings Technology and Research, 18:75–81, 2021. ISSN 1935-3804.
[66] I. Hayati, A. I. Bailey, and Th. F. Tadros. Mechanism of stable jet formation in
electrohydrodynamic atomization. Nature, 319(6048):41–43, 1986.
[67] Kelli Hendrickson, Gabriel D Weymouth, and Dick K-P Yue. Informed component
label algorithm for robust identification of connected components with volume-of-fluid
method. Computers & Fluids, 197:104373, 2020.
[68] Marcus Herrmann. A parallel eulerian interface tracking/lagrangian point particle
multi-scale coupling procedure. Journal of computational physics, 229(3):745–759,
2010.
[69] Winslow H Herschel and Ronald Bulkley. Konsistenzmessungen von gummi-
benzollösungen. Kolloid-Zeitschrift, 39:291–300, 1926.
[70] Andrew J Hewitt. Droplet size and agricultural spraying, part i: Atomization, spray
transport, deposition, drift, and droplet size measurement techniques. Atomization
and Sprays, 7(3), 1997.
[71] J. O. Hinze and H. Milborn. Atomization of Liquids by Means of a Rotating Cup.
Journal of Applied Mechanics, 17(2):145–153, 04 1950. ISSN 0021-8936.
[72] C. W. Hirt and B. D. Nichols. Volume of fluid (VOF) method for the dynamics of free
boundaries. Journal of Computational Physics, 39(1):201–225, January 1981.
[73] Moses M. Hohman, Michael Shin, Gregory Rutledge, and Michael P. Brenner.
Electrospinning and electrically forced jets. I. Stability theory. Physics of Fluids, 13
(8):2201–2220, 2001. ISSN 10706631.
[74] Moses M. Hohman, Michael Shin, Gregory Rutledge, and Michael P. Brenner.
Electrospinning and electrically forced jets. II. Applications. Physics of Fluids, 13
(8):2221–2236, 2001. ISSN 10706631.
[75] Joseph Hoshen and Raoul Kopelman. Percolation and cluster distribution. i. cluster
multiple labeling technique and critical concentration algorithm. Physical Review B,
14(8):3438, 1976.
[76] hypre. hypre: High performance preconditioners. https://llnl.gov/casc/hypre,
https://github.com/hypre-space/hypre.
103
[77] Kyoung-Su Im, Ming-Chia Lai, Yi Liu, Nasy Sankagiri, Thomas Loch, and Hossein
Nivi. Visualization and Measurement of Automotive Electrostatic Rotary-Bell Paint
Spray Transfer Processes. Journal of Fluids Engineering, 123(2):237–245, 12 2000.
ISSN 0098-2202.
[78] Kyoung-Su Im, Ming-Chia Lai, Sheng-Tao John Yu, and Robert R. Jr. Matheson.
Simulation of Spray Transfer Processes in Electrostatic Rotary Bell Sprayer. Journal
of Fluids Engineering, 126(3):449–456, 07 2004. ISSN 0098-2202.
[79] J.E Dendy J. Black box multigrid. Journal of Computational Physics, 48:366–386,
1982. ISSN 0021-9991.
[80] A Jaworek, AT Sobczyk, and A Krupa. Electrospray application to powder production
and surface coating. Journal of Aerosol Science, 125:57–92, 2018.
[81] Bavand Keshavarz, Eric C. Houze, John R. Moore, Michael R. Koerner, and Gareth H.
McKinley. Rotary atomization of newtonian and viscoelastic liquids. Physical Review
Fluids, 5(3), 2020.
[82] Muhammad Kashif Iqbal Khan, Maarten A.I. Schutyser, Karin Schroën, and Remko
Boom. The potential of electrospraying for hydrophobic film coating on foods. Journal
of Food Engineering, 108:410–416, 2012. ISSN 0260-8774.
[83] Douglas B Kothe and William J Rider. A comparison of interface tracking methods.
Technical report, Los Alamos National Lab.(LANL), Los Alamos, NM (United States),
1995.
[84] Agisilaos Kourmatzis and John S. Shrimpton. Electrohydrodynamics and charge
injection atomizers: A review of the governing equations and turbulence. Atomization
and Sprays, 19(11):1045–1063, 2009. ISSN 1044-5110.
[85] M Krafczyk, P Lehmann, O Filippova, D Hänel, and U Lantermann. Lattice boltzmann
simulations of complex multiphase flows. Multifield Problems: State of the art, pages
50–57, 2000.
[86] Venkata Krisshna and Mark Owkes. High-fidelity simulation of a rotary bell atomizer
with electrohydrodynamic effects. In ICLASS 2021, 15th Triennial International
Conference on Liquid Atomization and Spray Systems, 2021.
[87] Venkata Krisshna, Wanjiao Liu, and Mark Owkes. High-fidelity simulations of
a rotary bell atomizer with electrohydrodynamic effects. International Journal of
Multiphase Flow, 168:104566, 2023. ISSN 0301-9322. doi: https://doi.org/10.1016/
j.ijmultiphaseflow.2023.104566.
[88] Maximilian Kuhnhenn, Tórstein V. Joensen, Mads Reck, Ilia V. Roisman, and
Cameron Tropea. Study of the internal flow in a rotary atomizer and its influence
104
on the properties of the resulting spray. International Journal of Multiphase Flow,
100:30–40, 2018. ISSN 0301-9322.
[89] Robert J Lang. Ultrasonic atomization of liquids. The journal of the acoustical society
of America, 34(1):6–8, 1962.
[90] Pierre-Simon Laplace. Mechanism of the Heavens. Cambridge Library Collection -
Physical Sciences. Cambridge University Press, 2009.
[91] AH Lefebvre, XF Wang, and CA Martin. Spray characteristics of aerated-liquid
pressure atomizers. Journal of Propulsion and Power, 4(4):293–298, 1988.
[92] F. Lekien and J. Marsden. Tricubic interpolation in three dimensions. International
Journal for Numerical Methods in Engineering, 63(3):455–471, 2005. ISSN 00295981.
[93] Liya Li and Masaki Sawamoto. Multi-phase model on sediment transport in sheet-flow
regime under oscillatory flow. Coastal Engineering in Japan, 38(2):157–178, 1995.
[94] PY Liang and MD Schuman. Atomization modeling in a multiphase flow environment
and comparison with experiments. In AIAA, Fluid Dynamics, 21st Plasma Dynamics
and Lasers Conference, 21st, Seattle, WA, June 18-20, 1990. 17 p., 1990.
[95] Bin Liu, Hengxiang Hu, Lisen Bi, and Panagiotis E Theodorakis. Analysis of the
characteristics of the gas–liquid mixed artificial snow-making. International Journal
of Refrigeration, 149:155–167, 2023.
[96] Thomas Loch, Hossein Nivi, Yi Liu, Ming-Chia Lai, and Kyoung-Su Im. An
Experimental Investigation of Spray Transfer Processes in an Electrostatic Rotating
Bell Applicator. In International Body Engineering Conference & Exposition. SAE
International, sep 1998.
[97] JM López-Herrera, Stéphane Popinet, and MA2764018 Herrada. A charge-conservative
approach for simulating electrohydrodynamic two-phase flows using volume-of-fluid.
Journal of Computational Physics, 230(5):1939–1955, 2011.
[98] Wilai Luewisutthichat, Atsushi Tsutsumi, and Kunio Yoshida. Bubble characteristics
in multi-phase flow systems: bubble sizes and size distributions. Journal of chemical
engineering of Japan, 30(3):461–466, 1997.
[99] L. Ma, D. B. Ingham, M. Pourkashanian, and E. Carcadea. Review of the
Computational Fluid Dynamics Modeling of Fuel Cells. Journal of Fuel Cell Science
and Technology, 2(4):246–257, 04 2005. ISSN 1550-624X.
[100] Michael B. Mackaplow, Isidro E. Zarraga, and Jeffrey F. Morris. Rotary spray
congealing of a suspension: Effect of disk speed and dispersed particle properties.
Journal of Microencapsulation, 23(7):793–809, 2006.
105
[101] Ahmed Mahmoud and M. S. Youssef. Influence of spinning cup and disk atomizer
configurations on droplet size and velocity characteristics. Chemical Engineering
Science, 107:149–157, 2014. ISSN 0009-2509.
[102] A. Mark, Bjorn Andersson, S. Tafuri, K. Engstrom, H. Sorod, F. Edelvik, and J. S.
Carlson. Simulation of electrostatic rotary bell spray painting in automotive paint
shops. Atomization and Sprays, 23(1):25–45, 2013. ISSN 1044-5110.
[103] Vidar Mathiesen, Tron Solberg, Hamid Arastoopour, and Bjørn H Hjertager. Experi-
mental and computational study of multiphase gas/particle flow in a cfb riser. AIChE
journal, 45(12):2503–2518, 1999.
[104] James E. McCarthy and Dwight W. Senser. Specific charge measurements in
electrostatic air sprays. Particulate Science and Technology, 23(1):21–32, 2005. ISSN
02726351.
[105] J R Melcher and G I Taylor. Electrohydrodynamics: A Review of the Role of Interfacial
Shear Stresses. Annual Review of Fluid Mechanics, 1(1):111–146, 1969.
[106] James R Melcher. Continuum electromechanics, volume 2. MIT press Cambridge, MA,
1981.
[107] A. J. Mestel. Electrohydrodynamic stability of a highly viscous jet. Journal of Fluid
Mechanics, 312:311–326, apr 1996. ISSN 00221120.
[108] Paul D Morris, Andrew Narracott, Hendrik von Tengg-Kobligk, Daniel Alejandro Silva
Soto, Sarah Hsiao, Angela Lungu, Paul Evans, Neil W Bressloff, Patricia V Lawford,
D Rodney Hose, et al. Computational fluid dynamics modelling in cardiovascular
medicine. Heart, 102(1):18–28, 2016.
[109] B. Mostafaie Maynagh, B. Ghobadian, T. Tavakkoli Hashjin, and M. R. and
Jahannama. Effect of electrostatic induction parameters on dropletscharging for
agricultural application. Journal of Agricultural Science and Technology, 11(3), 2009.
[110] Igari Naoki, Iso Takuro, Nishio Yu, Izawa Seiichiro, and Fukunishi Yu. Numerical
simulation of droplet-formation in rotary atomizer. Theoretical and Applied Mechanics
Letters, 9:202–205, 2019. ISSN 2095-0349.
[111] Ghasem G Nasr, Andrew J Yule, and Lothar Bendig. Industrial sprays and atomization:
design, analysis and applications. Springer Science & Business Media, 2002.
[112] N. K. Nayyar and G. S. Murty. The stability of a dielectric liquid jet in the presence of
a longitudinal electric field. Proceedings of the Physical Society, 75(3):369–373, 1960.
ISSN 03701328.
106
[113] William F Noh and Paul Woodward. Slic (simple line interface calculation).
In Proceedings of the fifth international conference on numerical methods in fluid
dynamics June 28–July 2, 1976 Twente University, Enschede, pages 330–340. Springer,
2005.
[114] Shinzi Ogasawara, Masatoshi Daikoku, Minori Shirota, Takao Inamura, Yasuhiro Saito,
Kotaro Yasumura, Masakazu Shoji, Hideyuki Aoki, and Takatoshi Miura. Liquid
Atomization Using a Rotary Bell Cup Atomizer. Journal of Fluid Science and
Technology, 5:464–474, 2010.
[115] OICA. World Motor Vehicle Production - 2019 Production Statistics. Technical report,
Organisation Internationale des Constructeurs d’Automobiles, 2019.
[116] Stanley Osher and Ronald P Fedkiw. Level set methods and dynamic implicit surfaces,
volume 1. Springer New York, 2005.
[117] Stanley Osher and James A Sethian. Fronts propagating with curvature-dependent
speed: Algorithms based on hamilton-jacobi formulations. Journal of computational
physics, 79(1):12–49, 1988.
[118] Walter Oswald, Jan Lauk, Lutz Gödeke, Peter Ehrhard, and Norbert Willenbacher.
Analysis of paint flow pulsations during high-speed rotary bell atomization. Coatings,
9(10):1–9, 2019. ISSN 20796412.
[119] Mark Owkes and Olivier Desjardins. A computational framework for conservative,
three-dimensional, unsplit, geometric transport with application to the volume-of-fluid
(VOF) method. Journal of Computational Physics, 270:587–612, 2014. ISSN 0021-
9991.
[120] Mark Owkes and Olivier Desjardins. A mesh-decoupled height function method for
computing interface curvature. Journal of Computational Physics, 281:285–300, 2015.
ISSN 0021-9991.
[121] Mark Owkes and Olivier Desjardins. Consistent and conservative framework for
incompressible multiphase flow simulations. In APS Division of Fluid Dynamics
Meeting Abstracts, pages H7–007, 2015.
[122] Mark Owkes and Olivier Desjardins. A mass and momentum conserving unsplit semi-
Lagrangian framework for simulating multiphase flows. Journal of Computational
Physics, 332:21–46, 2017. ISSN 0021-9991.
[123] Mark Owkes, Eric Cauble, Jacob Senecal, and Robert A. Currie. Importance of
curvature evaluation scale for predictive simulations of dynamic gas–liquid interfaces.
Journal of Computational Physics, 365:37–55, 2018. ISSN 0021-9991.
107
[124] JD Oxley. Spray cooling and spray chilling for food ingredient and nutraceutical
encapsulation. In Encapsulation technologies and delivery systems for food ingredients
and nutraceuticals, pages 110–130. Elsevier, 2012.
[125] Suhas V Patankar. A numerical method for conduction in composite materials, flow
in irregular geometries and conjugate heat transfer. In International Heat Transfer
Conference Digital Library. Begel House Inc., 1978.
[126] Mohammad Reza Pendar and José Carlos Páscoa. Numerical modeling of electrostatic
spray painting transfer processes in rotary bell cup for automotive painting. Interna-
tional Journal of Heat and Fluid Flow, 80:108499, 2019. ISSN 0142727X.
[127] Mohammad Reza Pendar and José Carlos Páscoa. Atomization and spray character-
istics around an ERBS using various operational models and conditions: numerical
investigation. International Journal of Heat and Mass Transfer, 161, 2020. ISSN
00179310.
[128] Mohammad-Reza Pendar and José Carlos Páscoa. Numerical analysis of charged
droplets size distribution in the electrostatic coating process: Effect of different
operational conditions. Physics of Fluids, 33(3):033317, 2021.
[129] M Pilch and CA Erdman. Use of breakup time data and velocity history data to
predict the maximum size of stable fragments for acceleration-induced breakup of a
liquid drop. International journal of multiphase flow, 13(6):741–757, 1987.
[130] Stéphane Popinet and Stéphane Zaleski. A front-tracking algorithm for accurate
representation of surface tension. International Journal for Numerical Methods in
Fluids, 30(6):775–793, 1999.
[131] Precedence Research. Automotive Paints & Coatings Market Growth, Report 2020-
2027 Automotive. Technical report, Precedence Research, 2019.
[132] Andrea Prosperetti and Grétar Tryggvason. Computational methods for multiphase
flow. Cambridge university press, 2009.
[133] Pratikkumar Raje and Naresh Chandra Murmu. A Review on Electrohydrodynamic-
inkjet Printing Technology. 2014.
[134] R. Ray and P. Henshaw. Evaporation of clearcoat solvents from a rotary bell atomizer
and its relationship with bell speed, flow rate, and electrostatic potential. Journal of
Coatings Technology and Research, 15:41–49, 2018.
[135] Rolf D Reitz and Rm Diwakar. Structure of high-pressure fuel sprays. SAE
transactions, pages 492–509, 1987.
[136] Ning Ren, A Blum, C Do, and Andre W Marshall. Atomization and dispersion
measurements in fire sprinkler sprays. Atomization and Sprays, 19(12), 2009.
108
[137] Sajjad Rezayat and Mohammad Farshchi. Spray formation by a rotary atomizer
operating in the coriolis-induced stream-mode injection. Atomization and Sprays, 29
(10):937–963, 2019.
[138] William J Rider and Douglas B Kothe. Reconstructing volume tracking. Journal of
computational physics, 141(2):112–152, 1998.
[139] Clark Rubel and Mark Owkes. Extraction of droplet genealogies from high-fidelity
atomization simulations. Atomization and Sprays, 29(8):709–739, 2019. ISSN
10445110.
[140] Youcef Saad and Martin H Schultz. Gmres: A generalized minimal residual algorithm
for solving nonsymmetric linear systems. SIAM Journal on scientific and statistical
computing, 7(3):856–869, 1986.
[141] Poozesh Sadegh, Akafuah Nelson, and Saito Kozo. Effects of automotive paint spray
technology on the paint transfer efficiency – a review. Proceedings of the Institution
of Mechanical Engineers, Part D: Journal of Automobile Engineering, 232(2):282–301,
2018.
[142] B. Sanderse, S.P. van der Pijl, and B. Koren. Review of computational fluid dynamics
for wind turbine wake aerodynamics. Wind Energy, 14(7):799–819, 2011.
[143] D. A. Saville. Electrohydrodynamic Stability: Fluid Cylinders in Longitudinal Electric
Fields. Physics of Fluids, 13(12):2987–2994, 1970.
[144] D. A. Saville. Stability of Electrically Charged Viscous Cylinders. Physics of Fluids,
14(6):1095–1099, 1971.
[145] D. A. Saville. Electrohydrodynamics: The Taylor-Melcher Leaky Dielectric Model.
Annual Review of Fluid Mechanics, 29(1):27–64, 1997.
[146] Robert I. Saye, James A. Sethian, Brandon Petrouskie, Aaron Zatorsky, Xinyu Lu, and
Reza Rock. Insights from high-fidelity modeling of industrial rotary bell atomization.
Proceedings of the National Academy of Sciences, 120(4):e2216709120, 2023.
[147] Roland Schmehl, Georg Maier, and Sigmar Wittig. Cfd analysis of fuel atomization,
secondary droplet breakup and spray dispersion in the premix duct of a lpp combustor.
In ICLASS 2000: 8th International Conference on Liquid Atomization and Spray
Systems, Pasadena, CA, USA, 16-20 July 2000. ILASS, 2000.
[148] S. Schmidt, O. Krüger, K. Göckeler, and C. O. Paschereit. Numerical investigation of
the breakup behavior of an oscillating two-phase jet. Physics of Fluids, 30(7), 2018.
ISSN 10897666.
109
[149] James Albert Sethian. Level set methods and fast marching methods: evolving
interfaces in computational geometry, fluid mechanics, computer vision, and materials
science, volume 3. Cambridge university press, 1999.
[150] C. Sharma, D. Malhotra, and A. S. Rathore. Review of computational fluid dynamics
applications in biotechnology processes. Biotechnology Progress, 27(6):1497–1510,
2011.
[151] Patrick Sheehy and Mark Owkes. Numerical study of electric Reynolds number on
electrohydrodynamic assisted atomization. Atomization and Sprays, 27(7):645–664,
2017. ISSN 1044-5110.
[152] Bo Shen, Qiaoyan Ye, Nico Guettler, Oliver Tiedje, and Joachim Domnick. Primary
breakup of a non-Newtonian liquid using a high-speed rotary bell atomizer for spray-
painting processes. Journal of Coatings Technology and Research, 16:1581–1596, 2019.
[153] M Shirota, Y Hatayama, T Haneda, T Inamura, M Daikoku, Y Saito, and H Aoki.
Formation and breakup of ligaments from a rotary bell cup atomizer. In 12th
International Conference on Liquid Atomization and Spray Systems, ICLASS 2012,
2012.
[154] John S Shrimpton. Atomization, combustion, and control of charged hydrocarbon
sprays. Atomization and Sprays, 11(4), 2001.
[155] John S. Shrimpton and Yossapong Laoonual. Dynamics of electrically charged transient
evaporating sprays. International Journal for Numerical Methods in Engineering, 67
(8):1063–1081, 2006.
[156] K. Sidawi, P. Moroz, and S. Chandra. On surface area coverage by an electrostatic
rotating bell atomizer. Journal of Coatings Technology and Research, 18:649–663, 2021.
[157] K. Sidawi, P. Moroz, and S. Chandra. Bell-cup serrations and their effect on
atomization in electrostatic rotating bell atomizers. Experiments in Fluids, 62(8):
1–16, 2021. ISSN 14321114.
[158] S.D Sovani, P.E Sojka, and A.H Lefebvre. Effervescent atomization. Progress in Energy
and Combustion Science, 27(4):483–521, 2001.
[159] Vedanth Srinivasan, Abraham J. Salazar, and Kozo Saito. Modeling the disintegration
of modulated liquid jets using volume-of-fluid (VOF) methodology. Applied Mathe-
matical Modelling, 35(8):3710–3730, 2011.
[160] Ch Stevenin, Y. Bereaux, J. Y. Charmeau, and J. Balcaen. Shaping Air Flow Char-
acteristics of a High-Speed Rotary-Bell Sprayer for Automotive Painting Processes.
Journal of Fluids Engineering, Transactions of the ASME, 137(11):1–8, 2015. ISSN
1528901X.
110
[161] Brian P. Stimpson and Charles A. Evans Jr. Electrohydrodynamic ionization mass
spectrometry of biochemical materials. Biomedical Mass Spectrometry, 5(1):52–63,
1978.
[162] Otmar M. Stuetzer. Magnetohydrodynamics and electrohydrodynamics. The Physics
of Fluids, 5(5):534–544, 1962.
[163] Hilding Sundqvist and George Veronis. A simple finite-difference grid with non-
constant intervals. Tellus, 22(1):26–31, 1970. ISSN 0040-2826.
[164] Soma Tatsuya, Katayama Tomoyuki, Tanimoto Junichi, Saito Yasuhiro, Matsushita
Yohsuke, Aoki Hideyuki, Nakai Daichi, Kitamura Genki, Miura Masanari, Asakawa
Takukatsu, Daikoku Masatoshi, Haneda Toshiki, Hatayama Yohsuke, Shirota Minori,
and Inamura Takao. Liquid film flow on a high speed rotary bell-cup atomizer.
International Journal of Multiphase Flow, 70:96–103, 2015. ISSN 0301-9322.
[165] Geoffrey Taylor. Studies in Electrohydrodynamics. I. The Circulation Produced in a
Drop by Electrical Field. Proceedings of the Royal Society of London Series A, 291
(1425):159–166, April 1966.
[166] Nikola Toljic, K. Adamiak, G. S.P. Castle, Hong Hsiang Kuo, and Hua Tzu Fan.
Three-dimensional numerical studies on the effect of the particle charge to mass ratio
distribution in the electrostatic coating process. Journal of Electrostatics, 69(3):189–
194, 2011. ISSN 03043886.
[167] Gaurav Tomar, Daniel Gerlach, Gautam Biswas, Norbert Alleborn, Ashutosh Sharma,
Franz Durst, Samuel WJ Welch, and Antonio Delgado. Two-phase electrohydrody-
namic simulations using a volume-of-fluid approach. Journal of Computational Physics,
227(2):1267–1285, 2007.
[168] MN Topp and P Eisenklam. Industrial and medical uses of ultrasonic atomizers.
Ultrasonics, 10(3):127–133, 1972.
[169] G Tryggvason, B Bunner, O Ebrat, and W Tauber. Computations of multiphase flows
by a finite difference/front tracking method. i. multi-fluid flows. Lecture Series-von
Karman Institute For Fluid Dynamics, pages 7–7, 1998.
[170] Grétar Tryggvason, Bernard Bunner, Asghar Esmaeeli, Damir Juric, N Al-Rawahi,
W Tauber, J Han, S Nas, and Y-J Jan. A front-tracking method for the computations
of multiphase flow. Journal of computational physics, 169(2):708–759, 2001.
[171] R.J. Turnbull. Self-acceleration of a charged jet. IEEE Transactions on Industry
Applications, 25(4):699–704, 1989.
[172] Bret Van Poppel, Olivier Desjardins, and J. W. Daily. A ghost fluid, level set
methodology for simulating multiphase electrohydrodynamic flows with application
111
to liquid fuel injection. Journal of Computational Physics, 229:7977–7996, 2010. ISSN
0021-9991.
[173] James Van Strien, Phred Petersen, Petros Lappas, Leslie Yeo, Amgad Rezk, Sara
Vahaji, and Kiao Inthavong. Spray characteristics from nasal spray atomization.
Journal of Aerosol Science, 165:106009, 2022.
[174] Valerio Viti, J. Kulkarni, and A. Watve. Computational fluid dynamics analysis of the
electrostatic spray painting process with a rotating bell cup. Atomization and Sprays,
20(1):1–17, 2010. ISSN 1044-5110.
[175] Xin Tao Wang, Zhi Ning, and Ming Lü. Linear instability of a charged non-Newtonian
liquid jet under an axial electric field. Journal of Applied Physics, 126(13):135301,
2019. ISSN 10897550.
[176] Jacob E. Wilson, Stephen W. Grib, Adnan Darwish Ahmad, Michael W. Renfro,
Scott A. Adams, and Ahmad A. Salaimeh. Study of Near-Cup Droplet Breakup
of an Automotive Electrostatic Rotary Bell (ESRB) Atomizer Using High-Speed
Shadowgraph Imaging. Coatings, 8(5), 2018.
[177] Wolfgang Wulff. Computational methods for multiphase flow. Multiphase Science and
Technology, 5(1-4), 1990.
[178] Bin Xia and Da-Wen Sun. Applications of computational fluid dynamics (cfd) in the
food industry: a review. Computers and Electronics in Agriculture, 34(1):5–24, 2002.
ISSN 0168-1699.
[179] T Yamamoto and HR Velkoff. Electrohydrodynamics in an electrostatic precipitator.
Journal of Fluid Mechanics, 108:1–18, 1981.
[180] Ligong Yang. Numerical simulation of paint spray process. 2006.
[181] Chunde Yao, Zhong He, Huaqing Zhu, Yan Li, and Jianjun Li. The measurement of
droplet size distributions from rotary atomizers. Journal of agricultural engineering
research, 82(4):379–386, 2005.
[182] AJ Yule, JS Shrimpton, AP Watkins, W Balachandran, and D Hu. Electrostatically
atomized hydrocarbon sprays. Fuel, 74(7):1094–1103, 1995.
[183] Yi Zhu, De-an Zhao, Fa-zhong Li, and De-yuan Kong. Numerical simulation of the gas
flow field of electrostatic rotary-bell spray system. In 2010 3rd International Conference
on Computer Science and Information Technology, volume 2, pages 1–4. IEEE, 2010.
112
APPENDIX
113
A.1 Liquid Surface Profile in a Rotating Tank
y
2 2
y = ω x
2g
δm,P(x, y)
∂fx
∂fy ∂fN
x
Figure A.1: Analytical solution of the liquid surface (red) in a tank rotating about the y
axis with an angular velocity ω described mathematically by the function y
Consider a liquid mass element located at P(x, y) on the surface of the liquid in a tank
rotating about its vertical (y) axis with an angular velocity ω (Fig. A.1). Centrifugal force,
which acts in the x-direction, can be expressed as ∂fx = δm · ω2x. Gravity, which acts in
the y-direction, can be expressed as ∂fy = δm · g. ∂fN represents the normal force acting
perpendicular to the fluid element, which is the resultant of the vector sum of ∂fx and ∂fy.
Therefore, it can be seen that the local slope of the surface at point P is
∂y ∂f 2
y ω x
= = (A.1)
∂x ∂fx g
On integrating, this gives
ω2x2
y = (A.2)
2g
which is the mathematical formulation of the liquid surface in a rotating tank.
A.2 Trajectory of an Object in a Rotating Frame
Consider an object of mass δm on a two-dimensional plane in a reference frame that
is rotating counterclockwise with an angular velocity ω about an axis perpendicular to
the plane. The object experiences centrifugal and Coriolis forces. Therefore, according
114
to Newton’s Second Law,
∂2r
δm =∂f
2 centrifugal + ∂fCoriolis
∂t (A.3)
=δm(ω × r)× ∂r
ω + 2δm × ω
∂t
where r(x, y) is its position vector. The equation of motion reduces to
∂2r 2 ∂r
= ω r + 2 × ω (A.4)
∂t2 ∂t
whose components can be simplified as
∂x2
2 ∂y
= ω x+ 2ω (A.5)
∂t2 ∂t
∂y2 ∂x
= ω2y − 2ω . (A.6)
∂t2 ∂t
Multiplying Eq. A.6 by i and adding it to Eq. A.5 gives
∂x2 ∂y2 ∂2R ∂R
+ i = = ω2R− 2iω (A.7)
∂t2 ∂t2 ∂t2 ∂t
to which the general solution is of the form
R(t) = e−iωt (C1 + C2t) . (A.8)
Applying initial conditions of r(0) = (x ∂r
0, y0) and |t=0 = (u0, v0) as R(0) = (x
∂t 0, y0) and
∂R |t=0 = (u0 + iv0), C1 and C2 can be solved to give
∂t
R(t) = e−iωt [x0 + u0t+ i(v0 + ωx0)t] . (A.9)
Isolating the real and imaginary parts, the equation of instantaneous position is obtained as
x(t) = x0 − (x0 + u0t) cos(ωt) + t(v0 + x0ω) sin(ωt)
(A.10)
y(t) = y0 − (x0 + u0t) sin(ωt) + t(v0 + x0ω) cos(ωt).