COMPUTATIONAL STUDY OF INTERNAL FLOW,

NEAR NOZZLE AND EXTERNAL SPRAY OF A GDI

INJECTOR UNDER FLASH-BOILING CONDITIONS

Item Type dissertation

Authors Mohapatra, Chinmoy krushna

DOI 10.7275/27766634

Download date 2025-03-03 02:08:21

Link to Item https://hdl.handle.net/20.500.14394/18753

 COMPUTATIONAL STUDY OF INTERNAL FLOW,
NEAR NOZZLE AND EXTERNAL SPRAY OF A GDI
INJECTOR UNDER FLASH-BOILING CONDITIONS

A Dissertation Presented

by

CHINMOY K. MOHAPATRA

Submitted to the Graduate School of the

University of Massachusetts Amherst in partial fulfillment
of the requirements for the degree of

DOCTOR OF PHILOSOPHY

February 2022

Mechanical and Industrial Engineering

© Copyright by Chinmoy K. Mohapatra 2022

All Rights Reserved

 COMPUTATIONAL STUDY OF INTERNAL FLOW,
NEAR NOZZLE AND EXTERNAL SPRAY OF A GDI
INJECTOR UNDER FLASH-BOILING CONDITIONS

A Dissertation Presented

by

CHINMOY K. MOHAPATRA

Approved as to style and content by:

David P. Schmidt, Chair

Stephen de Bruyn Kops, Member

Omar A. Abdelrahman, Member

Sundar Krishnamurty, Department Head

Mechanical and Industrial Engineering
 DEDICATION

This dissertation is dedicated to my beloved parents.

 ACKNOWLEDGMENTS

I would like to acknowledge my Ph.D. advisor, Prof. David P. Schmidt, for

his immense support, patience, wisdom and phenomenal contributions in shaping me
as an independent researcher. Prof. Schmidt, you are an inspiring character. Your

passion for teaching and research while maintaining a healthy environment around

you will be the guiding blocks for my future endeavors. I will fondly look at all the

memories of our collaboration while working under your supervision. Prof. Schmidt,

I would be forever thankful to you for giving me the opportunity to learn from you.

I would like to acknowledge the support I received from my doctoral committee,

Prof. Stephen de Bruyn Kops, and Prof. Omar AbdelRahman, in the successful

completion of my doctoral degree. Thank you very much for your probing questions

and thorough guidance in building a higher quality dissertation. It has been a great

learning experience for me.

During the course of my Ph.D., I have been fortunate enough to collaborate with

some of the leading scientists from Argonne National Lab, Sandia National Lab and
other institutions. I would like to acknowledge some of them during this time. Thank

you, Dr. Lyle Pickett (Sandia) for valuable insights and taking the initiatives for the
Engine Combustion Network (ECN). I deem myself lucky to be able to contribute
to the ECN workshops and monthly web meetings. ECN shows us the power of

open-source science, where none of the contributors get paid directly for their effort.
However, they come together annually to discuss the interesting physics behind the
fuel injection systems both computationally and experimentally. Auto industries have

greatly benefited from the discoveries of ECN. I would like to thank Dr. Chris Powell
(Argonne), Dr. Brandon Sforzo (Argonne), Dr. Katie Matusik (former Argonne),

v
Prof. Daniel Duke (former Argonne), Prof. Kaushik Saha (former Argonne), Dr.
Gina Magnotti (Argonne), Dr. Lorenzo Nocivelli (Argonne), Dr. Joshua Lacey (

former University of Melbourne), and Dr. Joonsik Hwang (former Sandia) for their
valuable insights in to the physics of spray atomization. I have personally learned a

lot by having numerous discussions with you all at various platforms.

I would also like to acknowledge the support of colleagues and friends from the
Multi-phase flow simulation laboratory including Sampath, Hannah, Peetak, Brandon

and Gabe.
A special mention to my M.S. advisor Prof. Milind A. Jog from the University of

Cincinnati,OH. His timely guidance helped me to choose the research career path.

Additionally, the work of this dissertation would not have been possible without

the generous financial support from various industrial partners. I am thankful to

the General Motors Research Centre for funding the early part of Ph.D. research. A

special mention to Dr. Ronald Grover Jr. for his technical acumen which helped

me in formulating the early problem statements of my Ph.D. work. I am sincerely

grateful to the TOYOTA motors North America -R&D division for their financial

support to the major part of my Ph.D. research. I will cherish the endless technical

discussions with Dr. Oana Nitulescu and Dr. Tom Shieh. Those discussions helped

me immensely in developing the core of my Ph.D. research. I would like to thank the

partners from the Spray Combustion Consortium (SCC) for extending the financial
support to the exploratory study of my doctoral research. The support of NSF-funded
XSEDE computing resources and the Leadership Computing Facility at Argonne, for

running CFD simulations, is greatly appreciated.

Last, but not the least, I would like to acknowledge the sacrifice and support of
my family and friends from India in getting me here. Firstly, I am thankful to my

partner Archana, especially for her continuous support during the final journey of my
doctoral studies. My cousin, Dr. GuruKrushna, for always showing me the way and

vi
for his constant encouragements throughout. Finally, I acknowledge the unparalleled
sacrifices of my parents and uncle (Ashrambapa). Their teaching of never give-up

attitude and staying humble in success helped me a lot in navigating the troubled
waters during my doctoral studies.

I truly could not have achieved this success without each and every one of you
and I am eternally grateful to you all.

vii
 ABSTRACT

COMPUTATIONAL STUDY OF INTERNAL FLOW,

NEAR NOZZLE AND EXTERNAL SPRAY OF A GDI
INJECTOR UNDER FLASH-BOILING CONDITIONS

FEBRUARY 2022

CHINMOY K. MOHAPATRA
B. Tech., NATIONAL INSTITUTE OF TECHNOLOGY,ROURKELA

M. S., UNIVERSITY OF CINCINNATI

Ph.D., UNIVERSITY OF MASSACHUSETTS AMHERST

Directed by: Professor David P. Schmidt

The early and late portions of transient fuel injection have proven to be a rich area

of research, especially since the end of injection can cause a disproportionate amount

of emissions in direct injection internal combustion engines. While simulating the

internal flow of fuel injectors, valve opening and closing events are the perennial
challenges. A typical adaptive-mesh CFD simulation is extremely computationally

expensive, as the small gap between the needle valve and the seat requires very
small cells to be resolved properly. Capturing complete closure usually involves a
topological change in the computational domain. Furthermore, Internal Combustion

Engines(ICE) operating with Gasoline Direct Injection(GDI) principle are susceptible

to flash boiling due to the volatile nature of the fuel.

The presented work simulates a gasoline direct injector operating under cavitating

conditions by employing a more gradual and easily implemented model of closure that

viii
avoids spurious water-hammer effects. The results show cavitation at low valve lift for
both flash boiling and non-flash boiling conditions. Further, this study reveals post-

closure dynamics that result in dribble, which is expected to contribute to unburnt

hydrocarbon emissions. Flashing versus non-flashing conditions are shown to cause

different sac and nozzle behavior after needle closure. In particular, a slowly boiling
sac causes spurious injection behavior.
Furthermore, a qualitative analysis of the injector tip-wetting phenomena under

both flash-boiling and non-flashing conditions are conducted and different wetting
mechanisms are identified. The jet expansion mechanism is observed to dominate the

wetting process during the main injection period, whereas the sac conditions drive the

post-closure wetting phenomena. Additionally, the effect of flash-boiling conditions on

the near-nozzle spray during the quasi-steady period of the injection cycle is explored.

The exploration captured hole-to-hole variations in the rate of injection (ROI), rate of

momentum (ROM) and hydraulic coefficients of injection. Moreover, it also indicates

influences of the in-nozzle variations on the near-nozzle spray behaviors.

Finally, a novel plume-based coupling approach is developed to couple the Eule-

rian near nozzle simulations with the Lagrangian spray simulations under both non-

flashing and flash-boiling conditions. Predictions from the novel coupling approach

are validated with the experimental observations. This coupling approach requires

running an Eulerian primary atomization model, i.e., the Σ − Y model, to initialize

the Lagrangian parcels for the secondary atomization process. Hence, this coupling
approach does not depend upon the linearized instability models to simulate the dense

spray region.

ix
 TABLE OF CONTENTS

Page

ACKNOWLEDGMENTS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . v

ABSTRACT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii

LIST OF TABLES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiv

LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xv

CHAPTER

1. INTRODUCTION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1

1.1 Gasoline Direct Injection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2

1.2 Overview of the current research . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3

2. FLASH BOILING: A REVIEW . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

2.1 Physics of flash-boiling and theoretical modeling . . . . . . . . . . . . . . . . . . . . . . 8

2.1.1 Nucleation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10

2.1.1.1 Equilibrium and stability of a dispersed phase . . . . . . . . . 11

2.1.1.2 Modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

2.1.2 Bubble growth . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

2.1.2.1 Modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

2.2 Modeling of flashing flows . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

2.2.1 Bubble dynamic based model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

2.2.2 Thermodynamic rate based models . . . . . . . . . . . . . . . . . . . . . . . . . . 21
2.2.3 Interfacial exchange based model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24

2.3 Experimental studies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

x
 2.3.1 Flash-boiling mechanism . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
2.3.2 Flash-boiling and fuel injectors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28

2.4 CFD study of flash-boiling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32

2.5 Inferences . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34

3. ATOMIZATION MODELS: A REVIEW . . . . . . . . . . . . . . . . . . . . . . . . . . 36

3.1 Lagrangian-Eulerian methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38

3.1.1 Statistical description of the spray equation . . . . . . . . . . . . . . . . . . . 38

3.2 Direct Numerical Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42

3.3 Eulerian methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
3.4 Inferences . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47

4. GOVERNING EQUATIONS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49

4.1 CFD Solver . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49

4.1.1 PISO algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52

4.2 Needle Closure Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54

4.3 Primary atomization model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56

5. CASE STUDIES: MULTIPLE INJECTIONS IN GDI

SYSTEMS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60

5.1 Results and discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66

5.1.1 Sealing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
5.1.2 Rate of Injection (ROI) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
5.1.3 Cavitation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
5.1.4 Residual fuel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70

5.2 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70

6. EFFECT OF FLASH-BOILING ON THE NEEDLE CLOSURE

EVENT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72

6.1 Simulation setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74

6.2 Results and discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75

6.2.1 Rate of Injection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75

6.2.2 EOI sac condition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
6.2.3 The λ2 -criterion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
6.2.4 Tip wetting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84

xi
 6.2.5 Dribble . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88

6.3 Summary and conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89

7. A QUALITATIVE ANALYSIS OF THE TIP-WETTING

PHENOMENA UNDER FLASH-BOILING
CONDITIONS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93

7.1 Results and discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96

7.1.1 Phase One: Needle opening phase . . . . . . . . . . . . . . . . . . . . . . . . . . . 96

7.1.2 Phase Two: Downward moving needle phase . . . . . . . . . . . . . . . . . . 97
7.1.3 Phase Three: Needle closing phase . . . . . . . . . . . . . . . . . . . . . . . . . . 100
7.1.4 Phase Four: Post-needle closure phase . . . . . . . . . . . . . . . . . . . . . . 100

7.2 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104

8. QUASI STEADY STATE ANALYSIS OF MOVING NEEDLE

GDI SIMULATIONS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105

8.1 Results and discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106

8.1.1 Time averaged hydraulic coefficients . . . . . . . . . . . . . . . . . . . . . . . . 106

8.1.2 Hole-to-hole variations in ROI and ROM . . . . . . . . . . . . . . . . . . . . 110
8.1.3 Time averaged analysis of the near nozzle (2mm away from
the injector tip) behaviors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114
8.1.4 Time-scale analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121

8.2 Inferences . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122

9. PLUME BASED COUPLING APPROACH FOR GDI

SPRAY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125

9.1 Model description . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129

9.1.1 Lagrangian model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129

9.1.2 Plume-based coupling approach . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129

9.2 Simulation setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133

9.2.1 Lagrangian spray simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133

9.3 Results and discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135

9.3.1 Qualitative validation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135

9.3.2 Liquid penetration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142
9.3.3 Vapor penetration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145

xii
 9.3.4 Axial gas velocity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 147
9.3.5 Momentum conservation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 157

9.4 Summary and conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 162

10.PILOT STUDY AND CONCLUSIONS . . . . . . . . . . . . . . . . . . . . . . . . . . 164

10.1 Pilot study: Dynamically coupled sharp and diffuse interface

approach . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 164

10.1.1 Governing equation: interFOAM . . . . . . . . . . . . . . . . . . . . . . . . . . . 165

10.1.2 Governing equations: HRMFoam with additional surface
tension models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 167
10.1.3 Test case: Oscillating drop . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 168

10.2 Summary and Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 170

11.FUTURE WORK . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 175

BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 178

xiii
 LIST OF TABLES

Table Page

5.1 ECN Spray G nominal nozzle geometry parameters. . . . . . . . . . . . . . . . . . . 61

5.2 Non-flashing (spray G) operating condition defined by the ECN . . . . . . . . 63

6.1 Operating condition for G (Non-Flashing) and G2 (Flashing)

condition. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75

8.1 Time averaged hydraulic coefficients (Cd , CM , Cv , and Ca ), ROI,

and ROM for all holes combined and the individual holes under
both flash-boiling and non-flashing conditions. . . . . . . . . . . . . . . . . . . . 109

8.2 Different time-scales for both flash-boiling and non-flash-boiling

conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121

8.3 Turbulent time scales for individual holes during different phases of
the needle motion. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123

xiv
 LIST OF FIGURES

Figure Page

2.1 Pressure vs. volume curve, showing the thermodynamics of flashing

boiling. Image from [201]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

2.2 Spinodal Curve and attainable super-heat limit of water. Image from
[110]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10

2.3 Jet transition from mechanical breakup to flare flashing, as identified

by Cleary et al. [43]. Image adopted from [43]. . . . . . . . . . . . . . . . . . . . 28

3.1 Four main regimes of round jet breakup, namely the Rayleigh regime,
the first wind-induced regime, the second wind-induced regime,
and the atomization regime characterized by the Oh number and
Re number [118] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38

3.2 The spray regions of an atomizing jet. Image adopted from [66]. . . . . . . . 39

5.1 ECN spray G computational domain with the outlet plenum of 9mm
diameter. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62

5.2 ECN spray G nominal geometry showing distribution of the 8 nozzle

holes and 5 dimples. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62

5.3 Cut plane of the CFD mesh in the needle seat area prior to the
needle lift . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63

5.4 Needle lift profile used in the CFD simulations . . . . . . . . . . . . . . . . . . . . . . 64

5.5 Initialization of the sealing field at the start of simulation . . . . . . . . . . . . . 65

5.6 Initialization of the pressure field at the start of simulation . . . . . . . . . . . . 65

5.7 Initialization of the non-condensible gas field at the start of

simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66

5.8 Sealing Region and corresponding velocity contour at different time

stamps during needle closure event. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67

xv
5.9 Rate of Injection curve for a multiple injection event . . . . . . . . . . . . . . . . . 68

5.10 Vapor volume fraction and low pressure region showing cavitation
during needle opening phase of the 1st injection cycle. . . . . . . . . . . . . . 69

5.11 Vapor volume fraction showing cavitation during needle closing phase
of the both 1st and second injection cycle. . . . . . . . . . . . . . . . . . . . . . . . 70

5.12 Liquid fuel mass fraction and velocity contour at 1.151 ms (EOI of
the first injection),2.171 ms (EOI of the second injection) showing
dribble and jet disintegration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71

6.1 Ensemble-averaged needle displacement profile generated by X-ray

measurement [59] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74

6.2 Simulated and experimental ROI vs time. Experimental ROI

obtained through the long-tube method [59]and simulated ROI
taken at the nozzle exit . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76

6.3 Liquid volume fraction just before needle closing, i.e., t = 1.038 ms
for non-flashing conditions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77

6.4 Vapor volume fraction before needle closing, i.e., t = 1.038 ms for
non-flashing conditions showing cavitation due to the nozzle inlet
corner. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77

6.5 NCG volume fraction before needle closing i.e. t = 1.038 ms for
non-flashing conditions showing gas entrainment in the counter
bore. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77

6.6 Liquid volume fraction just after needle closing i.e. t = 1.045 ms for
non-flashing conditions showing cavitation. . . . . . . . . . . . . . . . . . . . . . . . 78

6.7 Vapor volume fraction just after needle closing i.e. t = 1.045 ms for
non-flashing conditions showing cavitation. . . . . . . . . . . . . . . . . . . . . . . . 78

6.8 NCG volume fraction just after needle closing i.e. t = 1.045 ms for
non-flashing conditions showing no presence of gas in the sac. . . . . . . . 78

6.9 Liquid volume fraction long after needle closing i.e. t = 2.000 ms for
non-flashing conditions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78

6.10 Vapor volume fraction long after needle closing i.e. t = 2.000 ms for
non-flashing condition showing no traces of vapor inside the
sac. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78

xvi
6.11 NCG volume fraction long after needle closing i.e. t = 2.000 ms for
non-flashing condition showing presence of gas inside the sac. . . . . . . . 78

6.12 Pressure ratio just before needle closing i.e. t = 1.038 ms for
non-flashing conditions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79

6.13 Pressure ratio just after needle closing i.e. t = 1.042 ms for
non-flashing conditions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79

6.14 Pressure ratio after needle closing i.e. t = 1.045 ms for non-flashing
conditions showing receding low pressure waves. . . . . . . . . . . . . . . . . . . 79

6.15 Liquid volume fraction before needle closing i.e. t = 1.038 ms for
flash-boiling conditions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80

6.16 Vapor volume fraction before needle closing i.e. t = 1.038 ms for
flash-boiling conditions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80

6.17 NCG volume fraction before needle closing i.e. t = 1.038 ms for
flash-boiling conditions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80

6.18 Liquid volume fraction just after needle closing i.e. t = 1.046 ms for
flash-boiling conditions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80

6.19 Vapor volume fraction long after needle closing i.e. t = 1.046 ms for
flash-boiling conditions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80

6.20 NCG volume fraction just after needle closing i.e. t = 1.046 ms for
flash-boiling conditions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80

6.21 Liquid volume fraction long after needle closing i.e. t = 2.000 ms for
flash-boiling conditions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80

6.22 Vapor volume fraction long after needle closing i.e. t = 2.000 ms for
flash-boiling conditions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80

6.23 NCG volume fraction long after needle closing i.e. t = 2.000 ms for
flash-boiling conditions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80

6.24 Pressure ratio just before needle closing, i.e., t = 1.038 ms for
flash-boiling conditions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81

6.25 Pressure ratio just after needle closing, i.e., t = 1.046 ms for
flash-boiling conditions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81

xvii
6.26 Pressure ratio long after needle closing, i.e., t = 2.000 ms for
flash-boiling conditions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81

6.27 Liquid Volume Fraction, Vapor Volume Fraction and NCG Volume
Fraction for a non flashing condition for the sac and nozzle
region . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82

6.28 Liquid Volume Fraction, Vapor Volume Fraction and NCG Volume
Fraction for a flash boiling condition for the sac and nozzle
region . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82

6.29 Iso-surface for λ2 = −8e + 11 colored by the pressure ratio for the
non-flashing condition at t = 1.042ms, i.e., after the needle
closure. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83

6.30 Iso-surface for λ2 = −8e + 11 colored by the VVF for the non-flashing
condition at t = 1.042ms, i.e., after the needle closure. . . . . . . . . . . . . . 83

6.31 Iso-surface for λ2 = −8e + 11 colored by the pressure ratio for the
flashing condition at t = 1.042ms, i.e., after the needle closure. . . . . . 85

6.32 Iso-surface for λ2 = −8e + 11 colored by the vapor volume fraction

for the flashing condition at t = 1.042ms i.e. after the needle
closure. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85

6.33 (a) Volume averaged

  density (ρ) multiplied with the vorticity
magnitude |Ω| ⃗ vs time for the both non-flashing and flashing
conditions. (b) The gray colored area showing the sac and nozzle
regions which are used for the volume averaging. . . . . . . . . . . . . . . . . . . 85

6.34 The injector tip colored by the liquid mass fraction for the
non-flashing condition at t = 0.98ms. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86

6.35 The injector tip colored by the liquid mass fraction for the
flash-boiling condition at t = 0.98ms. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86

6.36 Mid plane clip colored by the liquid mass fraction for the non-flashing
condition at t = 0.98ms showing NCG entrainment and the liquid
jet bending. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87

6.37 Mid plane clip colored by the liquid mass fraction for the flash-boiling
condition at t = 0.98ms showing no NCG entrainment. . . . . . . . . . . . . 87

xviii
6.38 The injector tip colored by the liquid mass fraction for the
non-flashing condition at t = 2.0ms. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87

6.39 The injector tip colored by the liquid mass fraction for the
flash-boiling condition at t = 2.0ms. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87

6.40 Simulated ROI vs time. Truncated ROI to indicate the oscillations

post needle closure. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89

6.41 Mid plane clip colored by the Liquid Mass Fraction indicating dribble
at t = 1.038ms for the non-flashing condition. . . . . . . . . . . . . . . . . . . . . 90

6.42 Mid plane clip colored by the velocity magnitude indicating different
velocities in the ensuing liquid jet at t = 1.038ms for the
non-flashing condition. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90

7.1 Different phases of tip-wetting phenomena marked on the needle

displacement curve. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97

7.2 Tips colored by the liquid mass fraction (LMF) indicating different
wetting phenomena for the non-flashing(left) and flashing (right)
conditions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98

7.3 Mid-plane clip colored by the liquid mass fraction (LMF) indicating
different jet behaviors for the non-flashing(left) and flashing
(right) conditions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98

7.4 Tips colored by the liquid mass fraction (LMF) indicating different
wetting phenomena for the non-flashing(left) and flashing (right)
conditions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99

7.5 Mid-plane clip colored by the liquid mass fraction (LMF) with glyph
vectors scaled by the velocity magnitude indicating different
mechanism driving the wetting phenomena for the
non-flashing(left) and flashing (right) conditions. . . . . . . . . . . . . . . . . . 99

7.6 Mid-plane clip colored by the liquid mass fraction (LMF) with glyph
vectors scaled by the velocity magnitude indicating the mechanism
behind the wetting phenomena for the non-flashing(left) and
flashing (right) conditions during needle closure. . . . . . . . . . . . . . . . . 101

7.7 Tips colored by the liquid mass fraction (LMF) indicating different
wetting phenomena for the non-flashing(left) and flashing (right)
conditions during needle closure. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102

xix
7.8 Tips colored by the liquid mass fraction (LMF) indicating different
wetting phenomena for the non-flashing(left) and flashing (right)
conditions post-needle closure. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103

8.1 Bar chart showing the hole-to-hole variations in the coefficient of

velocity under non-flashing and flashing conditions . . . . . . . . . . . . . . . 110

8.2 Bar chart showing the hole-to-hole variations in the coefficient of area
under non-flashing and flashing conditions . . . . . . . . . . . . . . . . . . . . . . 111

8.3 Relative standard deviation of the rate of injection across all the
holes for both the flashing and non-flashing conditions . . . . . . . . . . . . 112

8.4 Relative standard deviation of the momentum rate across all the
holes for both the flashing and non-flashing conditions . . . . . . . . . . . . 113

8.5 String flash-boiling in a transient moving mesh GDI simulation under

flashing condition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113

8.6 Line plots showing the time averaged liquid volume fraction for the
hole1(left) and hole5(right) at 2mm downstream of the injector
tip for both non-flashing (black color) and flashing (red color)
conditions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114

8.7 Line plots showing the time averaged vapor volume fraction for the
hole1(left) and hole5(right) at 2mm downstream of the injector
tip for both non-flashing (black color) and flashing (red color)
conditions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114

8.8 Line plots showing the time averaged liquid volume fraction for the
hole2 (left) and hole6 (right) at 2mm downstream of the injector
tip for both non-flashing (black color) and flashing (red color)
conditions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115

8.9 Line plots showing the time averaged vapor volume fraction for the
hole2 (left) and hole6 (right) at 2mm downstream of the injector
tip for both non-flashing (black color) and flashing (red color)
conditions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115

8.10 Line plots showing the time averaged liquid volume fraction for the
hole3 (left) and hole7 (right) at 2mm downstream of the injector
tip for both non-flashing (black color) and flashing (red color)
conditions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115

xx
8.11 Line plots showing the time averaged vapor volume fraction for the
hole3 (left) and hole7 (right) at 2mm downstream of the injector
tip for both non-flashing (black color) and flashing (red color)
conditions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115

8.12 Line plots showing the time averaged liquid volume fraction for the
hole4 (left) and hole8 (right) at 2mm downstream of the injector
tip for both non-flashing (black color) and flashing (red color)
conditions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116

8.13 Line plots showing the time averaged vapor volume fraction for the
hole4 (left) and hole8 (right) at 2mm downstream of the injector
tip for both non-flashing (black color) and flashing (red color)
conditions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116

8.14 Axial velocity predictions at 2 mm downstream from the injector tip

for all nozzle holes under the non-flashing condition . . . . . . . . . . . . . . 117

8.15 Axial velocity prediction at 2 mm downstream from the injector tip

for all nozzle holes under flash-boiling condition . . . . . . . . . . . . . . . . . 117

8.16 Cut planes colored by the mixture density. (a) Density measured by
X-Ray for the non-flashing condition. (b) Mixture density
predicted by HRMFoam for the non-flashing condition. . . . . . . . . . . . 118

8.17 2-D line plot showing the comparison between the X-Ray measured
mixture density and the CFD predictions for the non-flashing
condition. The mixture density is averaged over all the 8 holes for
the quasi-steady phase. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118

8.18 2-D line plot showing the comparison between X-Ray measured
mixture density and the CFD predictions for the hole one (left)
and hole 5( right) under flash-boiling condition. The mixture
density is time averaged over the quasi-steady phase. . . . . . . . . . . . . . 119

9.1 Plumes indicated by the iso-surface of LVF=0.04 colored by the

velocity magnitude at t = 0.5 ms for the non-flashing condition.
The injector tip is colored by the liquid mass fraction. . . . . . . . . . . . . 131

9.2 Plumes indicated by the iso-surface of LVF=0.04 colored by the

velocity magnitude at t = 1.0 ms for the non-flashing condition.
The injector tip is colored by the liquid mass fraction. . . . . . . . . . . . . 131

xxi
9.3 Plumes indicated by the iso-surface of LVF=0.04 colored by the
velocity magnitude at t = 0.5 ms for the flashing condition. The
injector tip is colored by the liquid mass fraction. . . . . . . . . . . . . . . . . 132

9.4 Plumes indicated by the iso-surface of LVF=0.04 colored by the

velocity magnitude at t = 1.0 ms for the flashing condition. The
injector tip is colored by the liquid mass fraction. . . . . . . . . . . . . . . . . 132

9.5 Plumes colored by the velocity magnitude with the velocity vectors at
t = 0.89 ms for the non-flashing condition . . . . . . . . . . . . . . . . . . . . . . 133

9.6 The location in x,y,z space of each parcel created at a single time step
i.e. 0.55 ms after the start of injection. In this time step, 6327
parcels are created. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133

9.7 Flow chart indicating all the major steps and the
software/application used during the each step for the plume
based one-way coupling approach. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134

9.8 (a) Computational domain for the Eulerian internal flow simulation.
(b) Computational domain for the Lagrangian external spray
simulation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134

9.9 DBI images taken from the ECN data set [1] which show the liquid
penetration for the non-flashing condition at different time
instances, i.e., (a) 0.1 ms ASOI, (b) 0.3 ms ASOI, (c)
0.5 ms ASOI, (d) 0.7 ms ASOI, (e) 0.9 ms ASOI, and (f)
1.1 ms ASOI. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137

9.10 Binarized images showing the liquid spray boundary predicted by the
plume based coupling approach with W ecr = 1 for the
non-flashing condition at different time instances, i.e., (a)
0.1 ms ASOI, (b) 0.3 ms ASOI, (c) 0.5 ms ASOI, (d)
0.7 ms ASOI, (e) 0.9 ms ASOI, and (f) 1.1 ms ASOI. . . . . . . . . . . 137

9.11 Binarized images showing the liquid spray boundary predicted by the
plume based coupling approach with W ecr = 6 for the
non-flashing condition at different time instances, i.e., (a)
0.1 ms ASOI, (b) 0.3 ms ASOI, (c) 0.5 ms ASOI, (d)
0.7 ms ASOI, (e) 0.9 ms ASOI, and (f) 1.1 ms ASOI. . . . . . . . . . . 138

xxii
9.12 DBI images taken from the ECN data set [1] which show the liquid
penetration for the flashing condition at different time instances,
i.e., (a) 0.1 ms ASOI, (b) 0.3 ms ASOI, (c) 0.5 ms ASOI, (d)
0.7 ms ASOI, (e) 0.9 ms ASOI, (f) 1.1 ms ASOI, (g)
1.3 ms ASOI, (h) 1.5 ms ASOI, (i) 1.7 ms ASOI . . . . . . . . . . . . . . . 139

9.13 Binarized images showing the liquid spray boundary predicted by the
plume based coupling approach with W ecr = 1 for the flashing
condition at different time instances, i.e., (a) 0.1 ms ASOI, (b)
0.3 ms ASOI, (c) 0.5 ms ASOI, (d) 0.7 ms ASOI, (e)
0.9 ms ASOI, and (f) 1.1 ms ASOI. . . . . . . . . . . . . . . . . . . . . . . . . . . . 140

9.14 Binarized images showing the liquid spray boundary predicted by the
plume based coupling approach with W ecr = 6 for the flashing
condition at different time instances, i.e., (a) 0.1 ms ASOI, (b)
0.3 ms ASOI, (c) 0.5 ms ASOI, (d) 0.7 ms ASOI, (e)
0.9 ms ASOI, and (f) 1.1 ms ASOI. . . . . . . . . . . . . . . . . . . . . . . . . . . . 141

9.15 The figure includes predictions for the liquid penetration by the
plume-based coupling approach for two critical Weber number.(1
and 6) under the non-flashing condition. The predictions are
compared with the experimentally measured liquid penetration
based on two LVF thresholds for the same condition by the
Sandia National Lab (SNL) [81] and University of Melbourne
(UoM) [150]. The experimental data are obtained from the ECN
data base [1]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143

9.16 The figure includes predictions for the liquid penetration by the
plume-based coupling approach for two critical Weber number (1
and 6) under the flashing condition. The predictions are compared
with the experimentally measured liquid penetration based on two
LVF thresholds for the same condition by the Sandia National
Lab (SNL) [81] and University of Melbourne (UoM) [150]. The
experimental data are obtained from the ECN data base [1]. . . . . . . 144

9.17 The figure includes predictions for the vapor penetration by the
plume-based coupling approach for two critical Weber no.(1 and
6) under the non-flashing condition. The predictions are
compared with the experimentally measured vapor penetration
based for the same condition by the Sandia National Lab (SNL)
and University of Melbourne (UoM). The experimental data are
obtained from the ECN data base [1]. . . . . . . . . . . . . . . . . . . . . . . . . . 146

xxiii
9.18 The figure includes predictions for the vapor penetration by the
plume-based coupling approach for two critical Weber number (1
and 6) under the flashing condition. The predictions are
compared with the experimentally measured vapor penetration
based for the same condition by the Sandia National Lab
(SNL) [81] and University of Melbourne (UoM) [150]. The
experimental data are obtained from the ECN data base [1]. . . . . . . 148

9.19 The figure includes experimental measurements for the axial velocity
at a location 15 mm downstream of the injector tip along the
injector axis [209] and the computational predictions by the
plume-based coupling approach for the non-flashing condition. . . . . . 149

9.20 The figure includes comparison of the axial velocity predictions by

the plume-based coupling approach under flashing and
non-flashing conditions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 150

9.21 Vectors showing the direction of the ambient gas movement for the
non-flashing condition and W ecr = 1 at time = 0.1ms ASOI. The
red horizontal line is located at 15 mm downstream of the injector
tip. The black vertical line is the injector axis. . . . . . . . . . . . . . . . . . . . 152

9.22 Vectors showing the direction of the ambient gas movement for the
flashing condition and W ecr = 6 at time = 0.1ms ASOI. The red
horizontal line is located at 15 mm downstream of the injector tip.
The black vertical line is the injector axis. . . . . . . . . . . . . . . . . . . . . . . 153

9.23 Vectors showing the direction of the ambient gas movement for the
non-flashing condition and W ecr = 1 at time = 0.2ms ASOI. The
red horizontal line is located at 15 mm downstream of the injector
tip. The black vertical line is the injector axis. . . . . . . . . . . . . . . . . . . . 154

9.24 Vectors showing the direction of the ambient gas movement for the
non-flashing condition and W ecr = 1 at time = 0.4ms ASOI. The
red horizontal line is located at 15 mm downstream of the injector
tip. The black vertical line is the injector axis. . . . . . . . . . . . . . . . . . . . 155

9.25 Vectors showing the direction of the ambient gas movement for the
non-flashing condition and W ecr = 1 at time = 0.6ms ASOI. The
red horizontal line is located at 15 mm downstream of the injector
tip. The black vertical line is the injector axis. . . . . . . . . . . . . . . . . . . . 156

xxiv
9.26 Vectors showing the direction of the ambient gas movement for the
non-flashing condition and W ecr = 1 at time = 1.2ms ASOI. The
red horizontal line is located at 15 mm downstream of the injector
tip. The black vertical line is the injector axis. . . . . . . . . . . . . . . . . . . . 157

9.27 Vectors showing the direction of the ambient gas movement for the
non-flashing condition and W ecr = 6 at time = 1.2ms ASOI. The
red horizontal line is located at 15 mm downstream of the injector
tip. The black vertical line is the injector axis. . . . . . . . . . . . . . . . . . . . 158

9.28 Total momentum comparison between the Eulerian HRMFoam and

the Lagrangian plume-based simulations for the non-flashing
condition. The Eulerian momentum is evaluated at the exit of the
nozzle. The Lagrangian momentum is calculated based on the
parcel momentum . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 160

9.29 Total momentum comparison between the Eulerian HRMFoam and

the Lagrangian plume-based simulations for the flashing
condition. The Eulerian momentum is evaluated at the exit of the
nozzle. The Lagrangian momentum is calculated based on the
parcel momentum . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 161

10.1 Cutplanes colored by the NCG mass fraction showing the initial
droplet with different mesh resolutions. (a) coarse (b) Fine . . . . . . . . 168

10.2 Mid-plane clip colored by the pressure difference indicating a higher

pressure inside the droplet, i.e. around 4.5 Pa for the fine mesh.
The analytical pressure difference in 4.4 Pa for the same
droplet. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 170

10.3 Graphs indicating the oscillations of the droplet radius in z-direction

for two different meshes as predicted by the surface tension based
HRMFoam . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 171

xxv
 CHAPTER 1

INTRODUCTION

Since the age of the industrial revolution, the internal combustion engine(ICE)

has been one of the greatest inventions. Numerous engineers and researchers have
contributed to its continuous improvement with their growing knowledge of thermo-
dynamics and combustion. ICEs have seen a wide range of applications in utility

devices (pumps, mowers, generators, etc.), as well as all types of means of trans-
portation (tractors, propeller aircraft, ships, passenger cars, and other on/off road
vehicles). ICEs have revolutionized the transportation sector due to their high ther-

mal efficiency and high power output to volume ratio.

The wide use of ICEs can be linked to the economic growth and the industrial-
ization of the modern world. In the early days, people were not concerned about the

size and efficiency of the engines. However, the situation has changed drastically in
the last decades. Nowadays ICEs are fighting with electrical motors to coexist in the
automotive world.

We live in a world where the fossil fuel supply is finite and the population grows
every day resulting in a higher energy demand. According to a study by Exxon [3]
the global energy demand is expected to grow by 25% in the year 2040. One-third

of this requirement is expected to be met by fossil fuels. Increasing usage of fossil

fuel has given rise to a growing concern over pollution and climate change, which has
compelled engineers to develop cleaner and more efficient means of power generation.

Although ICEs have undergone a series of development, there is still a lot of

room for improvement. Under ideal conditions, modern gasoline engines have an

1
efficiency of 30%, which gets lower in real-world conditions. Furthermore, combustion
of fossil fuel in direct injection engines (DI) produces nitrogen oxides (NOx), unburnt

hydrocarbons (UHC), carbon monoxide(CO), and other greenhouse gases. As the

global emission and fuel economy regulations continue to tighten, automakers are

working diligently to develop cleaner and more efficient engines.

Design of a better and cleaner combustion system is essential for the development
of cleaner and efficient vehicles. An important step in this designing process is to op-

timize the fuel injection systems. They are active area of research, as they determine
the characteristics and distribution of the air-fuel mixture through the combustion
cycles. The current work addresses some of the most pertinent questions of the fuel

injection and spray atomization research, especially in the context of gasoline direct
injection(GDI) systems.

1.1 Gasoline Direct Injection

Gasoline Direct Injection is considered to be a more advanced version of the con-
ventional multi-port fuel injection(MPFI) system, where the fuel is directly injected

into the combustion chamber instead of the intake ports. In recent years, global auto
manufacturers have preferred GDI systems over the MPFI systems due to improved
combustion efficiencies, reduced emissions and increased fuel economy in ICEs appli-

cations [75], [244]. The injection pressure of GDI systems is higher as compared to
that of the MPFI systems. Consequently, it leads to better atomization of the fuel
entering into the combustion chambers, thus yielding a higher rate of vaporization.

Usually, a spontaneous phase change process known as flash-boiling is observed in

gasoline injection systems due to the volatile nature of the gasoline fuel, i.e., relatively
high vapor pressure at low temperature. Flash-boiling is encountered when a super-

heated liquid is subjected to sudden depressurization. This phase change process

is governed by a finite rate of heat transfer. Flash-boiling phenomena are usually

2
encountered in the gasoline engines under low and part load conditions. Although
flash-boiling phenomena have been widely witnessed in several chemical engineering

and nuclear engineering applications, their use as an effective atomization mechanism

in ICEs is relatively new. Hence, the effects of flash-boiling phenoemna on the gasoline

spray are not completely understood and constitute the essence of this dissertation
work.

1.2 Overview of the current research

The current report for my doctoral dissertation consists of a detailed literature
survey in the flash-boiling atomization field and descriptions of my research observa-
tions.

• In Chapter 2 of this dissertation, the underlying physics of the flash-boiling

process will be discussed. Additionally, different modeling approaches for the

flash-boiling process and their results will be reviewed. This chapter will also
include some of the prominent experimental observations of flash-boiling phe-
nomena in the spray applications.

• Chapter 3 will review the different atomization modeling approaches and their

advantages and shortcomings.

• Chapter 4 will include a detailed description of the governing equations for

the in-house developed CFD code for internal flow modeling. Additionally, the
novel sealing algorithm and the primary atomization model will be discussed.

• In Chapter 5 of this dissertation, observations for a first of kind multiple injec-

tion cycle simulation for the GDI injector have been presented.

• In Chapter 6 of this dissertation, the sealing algorithm will be applied to study

the end of injection (EOI) behaviors of the GDI injector.

3
• In chapter 7 of this dissertation, a qualitative analysis will be performed to
understand the tip-wetting behaviors under both non-flashing and flashing con-

ditions.

• Hole-to-hole variations in the mass flow rate and the momentum rate predictions
will be discussed in Chapter 8. The chapter will also include some of the quasi-
steady state analysis of the near-nozzle predictions.

• In chapter 9, a novel plume-based coupling approach will be developed. Predic-

tions from this approach will be validated with the available experimental data
for both flashing and non-flashing conditions.

• In chapter 10, a pilot study to include the surface tension based models in the
in-house HRMFoam solver will be performed. The chapter will also summarize

the current findings of this dissertation.

• Based on the summary and the current findings of this dissertation, a few future
research recommendations will be proposed in chapter 11.

4
 CHAPTER 2

FLASH BOILING: A REVIEW

The atomization process involves disintegration of bulk liquid mass into small

droplets in a gaseous environment. It has wide industrial applications ranging from

internal combustion engines, gas turbine engines, industrial furnaces, drug delivery
spray, spray coating of surface materials to food processing industries. Among these

applications, ICEs demand effective combustion of tiny droplets produced by the

atomization processes. Furthermore, effective combustion is driven by the evaporation
rate of droplets with less effective surface area. These droplets are produced by the

successful break up of liquid jets or annular liquid sheets during atomization processes.
Hence, liquid jet and sheet break up has been the soul of the atomization re-
search over the past century. All conventional atomizers (centrifugal, rotary, pressure

swirl, ultrasonic, piezoelectric atomizers, etc.) use mechanical energy to disintegrate

the continuous bulk liquid phase into a dispersed phase. However, some atomizers
ensure the disintegration of the bulk liquid by forcing a biphasic flow, i.e., bubbly

mixture through orifices. In effervescent atomizers, the bubbly mixture is ensured

by introducing external atomizing gases into the liquid. However, researchers have
faced difficulties in controlling the spray parameters of effervescent atomizers. The

other way of ensuring the biphasic flow is by super heating either the bulk liquid or a
specific liquid component of the mixture. The following discussions will explore the
flash-boiling atomization methods.

The liquid is said to be super heated when it is at a temperature above the

saturation condition corresponding to its pressure. As shown in the Figure 2.1, liquid

5
can gain its super heat in two different ways, i.e., either being heated (from O to
A) at a constant pressure or by depressurization (from O to B). Vaporization caused

by rapid depressurization can be classified as either cavitation or flash boiling based

upon the ambient pressure or injection temperature level.

Figure 2.1: Pressure vs. volume curve, showing the thermodynamics of flashing
boiling. Image from [201].

Cavitation is usually encountered in the regions of low temperature or low pres-

sure. In these regions, vapor density is low and a small amount of super heat is
sufficient enough to initiate the vaporization process. The cavitating bubble are

short lived, as they collapse within micro-seconds. On the contrary, flash-boiling is

a highly thermal non-equilibrium process and is witnessed in the regions with high

6
temperature and pressure conditions. Flash-boiling is usually encountered when the
downstream pressure is lesser than the liquid vapor pressure. Thus, the vapor bubbles

generated locally in the upstream regions do not collapse; rather, they stay intact and
expand in the downstream.

In flowing systems, depressurization might occur due to changes in the channel

flow area, for example, the diverging sections of converging diverging nozzles. Sudden
pressure release due to the valve opening/closing or some crack development in liquid

containers can also lead to the depressurization process. The phase change phenomena
can also occur because of the hydro static pressure drop in vertical flow paths.
Flash-boiling has a decisive presence in many industrial or technical applications.

For instance, desalination of sea water using multi-stage flash distillation [211], im-
pulse drying of grapes to improve wine quality [196], and paper pulp production [234].
In many chemical and industrial applications, accidental release of a liquefied haz-

ardous chemical or contaminant can potentially lead to fatal injuries and serious
environmental pollution. Under such conditions, the nature of the release i.e., either
sub-cooled or super-heated is the key to quantify the hazard due to the leakage [43].

In the nuclear industry, the depressurization problem is associated with the hypothet-
ical Loss of Coolant Accident(LOCA) of pressurized water. When uncontrolled, such
scenarios can be hazardous for human lives and the environment. Hence, a thorough

understanding of the flash-boiling process is also essential in addressing such issues.

The ability to produce rapid phase change has motivated the use of flash-boiling
sprays in the cooling of hot space shuttle parts [6], [7]. Flash-boiling of geothermal

fluids have been in use for the power generation in water-based geothermal power
plants [54], [205]. Lately, flash atomization has also witnessed applications in print-
ing industries and pharmaceutical sprays because of its ability to produce fine spray.

In pharmaceutical sprays, large droplets can lead to oral thrush, cough and bron-
chospasm [227].

7
 With growing concerns of global warming and climate change and the subsequent
tightening of the emission norms, design of efficient spray combustion systems has

become essential in automotive and aerospace industries. Sprays formed by flash-

boiling processes are usually accompanied by droplets of smaller mean diameter,

higher homogeneity, wider cone angles and shorter penetration depth compared to
the sprays formed by the mechanical means for the same range of operating pressure.
Flash-boiling sprays have the ability to provide the desired sprays for combustion

systems at low injection pressure ranges. Besides, flashing is also known to occur in
the engine start-up stages in the case of liquid propellant based rocket engines [99]
and also during the idle and low-speed conditions in gasoline engines [12], [240].

Irrespective of being subjected to a great deal of scrutiny both experimentally and

theoretically, the application of flash-boiling as a primary atomization mechanism is
still not understood completely.

In the following sections, some prominent experiments, modeling procedures, and

CFD-studies especially in the context of the application of flash-boiling as an atom-
ization mechanism in automotive fuel injectors are reviewed.

2.1 Physics of flash-boiling and theoretical modeling

From a physical point of view, the two-phase flow with mass and heat transfer,

i.e, the boiling flow, is more complicated compared to the single phase flow. One of the
many challenges is to understand the mechanism of bubble formation in pure liquid
and to describe the rapid transition from the pure liquid state to the vapor-liquid

mixture.
Hence, understanding the complete physics involved in the flash-boiling process is
exacting. The inception of flashing in an initially sub-cooled liquid occurs when its

super-heat limit is reached. The maximum degree of attainable liquid super-heat or

pressure undershoot determines the intensity of the subsequent boiling process. The

8
isothermal pressure drop, shown by the line OB in Fig. 2.1, corresponds to a super-
heated condition of the liquid. It implies that at a given temperature, the liquid has

a pressure that is lower than its saturation pressure. Whereas, line OA in Fig. 2.1
refers to a scenario where the liquid temperature is more than its boiling temperature

at a given pressure. The minima of the isotherms in the P-V diagram is given by
∂P
= 0, which equivalently means that the isothermal compressibility ,KT → ∞.
∂V
The locus of such points is called the liquid spinodal curve. The condition of stability

of pure fluids require the isothermal compressibility to be a positive number. Thus,

the spinodal curve determines the theoretical thermodynamics extent to which the
liquid can be brought without any vaporization. Beyond this threshold, the pure

liquid ceases to exist and will undergo spontaneous phase change as a reaction to any
significant perturbation, for example, rapid pressure drop. The region bounded by
the liquid saturation curve and the liquid spinodal curve is known as the metastable

liquid region.
Furthermore, the spinodal curve represents a fictional limit which can never be
verified experimentally, because the liquid always starts to flash-boil in experiments

before it reaches this limit. Therefore, it is often replaced by the empirical kinetic
homogeneous nucleation limit, which can be observed experimentally. Lienhard and
Karimi [?] have observed this limit to lie close to the spinodal curve. However, the

experiments for achieving the homogeneous nucleation limit must be performed in

carefully controlled laboratory conditions, free from impurities. That includes the
liquid is free from any impurities and the experimental vessel is perfectly clean to

avoid any wall effects. Otherwise, the meta-stable liquid undergoes phase change
before reaching the homogeneous nucleation limit. The limit of such nucleation pro-
cess is defined by the heterogeneous nucleation limit curve, which lies between the

liquid saturation curve and the homogeneous nucleation limit curve. In Fig. 2.2, the
blue dashed line represents the homogeneous nucleation limit, whereas the red dashed

9
line represents the heterogeneous nucleation limit. This process of rapid/spontaneous
phase change in a super-heated liquid is termed as flash-boiling and is divided into

three stages namely, nucleation, bubble growth and atomization. Each of these pro-
cesses are described below.

Figure 2.2: Spinodal Curve and attainable super-heat limit of water. Image from
[110].

2.1.1 Nucleation

The nucleation process involves formation of bubbles either from the bulk liq-
uid or from other sources such as dissolved gases, suspended particles and/or the

imperfections of the orifice walls. Based upon the source of bubble formation, the
nucleation process can be broadly classified into homogeneous nucleation and hetero-
geneous nucleation. Before describing the modeling of different nucleation processes,

it is important to understand the equilibrium of dispersed phase and the critical

work required to form a stable embryo during the nucleation process, which will be
explained in the following section.

10
2.1.1.1 Equilibrium and stability of a dispersed phase

In the current study, the liquid phase has been referred to as the continuous
phase, whereas the gaseous phase as the dispersed phase. The equilibrium condition
between the dispersed phase, i.e., the vapor bubbles and the the liquid phase can be

examined based upon the thermodynamic and mechanical point of views. Both the
equilibrium conditions are well represented by the Clausius-Clapeyron Equation and
the Laplace-Kelvin Equation, respectively.

1. The Clausius-Clapeyron Equation:

The Clausius-Clapeyron equation assumes thermodynamic equilibrium at the

vapor-liquid interface. Therefore, the vapor pressure is considered to be sat-

urated. As a result, the correlation shown in the Eqn. 2.1 is dependent upon
temperatures.
 
dp hLG
=   (2.1)
dT sat Tsat vG,sat − vL,sat
The Eqn. 2.1 is further simiplified into a linearized form, as shown in the
Eqn. 2.2.
 
TG − Tsat hLG
PG − PL =   (2.2)
Tsat vG,sat − vL,sat
In Eqn. 2.1 and 2.2, TG and Tsat are the vapor temperature and saturated liquid
temperature, respectively, i.e. TG = Tsat (PG ), Tsat = Tsat (PL ). Additionally,

vG,sat and vL,sat represents the specific volume of the vapor and the liquid,
respectively. The heat of vaporization is presented by hLG . Furthermore, PL
and PG represents the pressure in the liquid and the vapor phase.

The Clausius-Clapeyron equation assumes a fairly flat interface between phases,

hence neglects the surface force. In such cases, thermodynamic equilibrium can

exist between the phases , if the two phases are of equal temperature, pressure
and chemical potential.

11
2. The Laplace-Kelvin Equation:

For a dispersed phase, where the interface shape is significant, a mechanical equi-
librium needs to be considered in addition to the thermodynamic equilibrium.
The pressure difference between the two phases that are in thermodynamic

equilibrium is given by the Laplace’s equation, represented in Eqn. 2.3.

2σ
PG − PL = (2.3)
Rl
In the Eqn. 2.3, σ and Rl represents the surface tension and the radius of
curvature, respectively. Under equilibrium conditions, the chemical potential of

the two phases must be equal. Using this condition and the linearized Clausius-
Clapeyron Equation, Cole et al. [45] developed the following formulation.

 
vL,sat
PG − PL = (Psat − PL ) 1 − (2.4)
vG,sat

The set of Eqns. 2.3 and 2.4 are known as the Laplace-Kelvin equation. It states

that a vapor embryo grows only after reaching a critical bubble radius (Rcr ).
The work needed to create the vapor embryo of the critical radius is defined as
the critical work (Wcr ), which is further expressed in the Eqn. 2.5

4 2 16πσ 3
Wcr = πσRcr = vL 2 (2.5)
3 3(Psat − PL )2 (1 − vG
)

In the Eqn. 2.5, Psat represents the saturation pressure of the liquid. The critical
work can further be expressed based on the Gibbs number(Gb), Boltzmann’s
constant (k) and the liquid temperature (T ), as shown in the Eqn. 2.6

Wcr
Gb = (2.6)
kT

12
2.1.1.2 Modeling

As discussed earlier, the nucleation process can be subdivided into homogeneous

nucleation and heterogeneous nucleation.

• Homogeneous nucleation: The homogeneous nucleation, also referred to as bulk

nucleation, is usually encountered at high degrees of superheat and at ran-
dom locations in the bulk liquid. Fluctuations of the bulk liquid density at

the molecular level are usually the primary reasons behind the bulk nucleation
phenomena. Blander and Katz [32] observed bubble nucleation to occur when
the degree of superheat is as high as 90% of the liquid’s critical temperature.

Furthermore, the kinetic theory of homogeneous nucleation is based on the as-

sumption that at a given time, the number of nuclei growing to a visible size is
proportional to the number of nuclei formed. Using this theory, Skripov [206]

formulated the set of Eqns. 2.7 and 2.8 to predict the number of bubbles gen-
erated per unit volume per unit time(JCN T ).

 
Wcr
JCN T = J0 exp − (2.7)
kT
J0 = NM B ′ (2.8)

Where, NM is the number density of the liquid molecules and is calculated based
upon the molecular weight of the liquid. He also suggested that the molecular
interactions per second (B ′ ) to be a constant i.e. 1012 . On the contrary, Blander

and Katz [32] have suggested a widely accepted formulation for J0 , mentioned
in the Eqn. 2.9.
r
2σ
J0 = NM (2.9)
πmM B

In the Eqn. 2.9, B = 1 for cavitating conditions, otherwise B = 2/3 [32].

However, the above expression ignores the correlation between successive events

13
 which impacts the number of vapor embryos present in a cluster at a particular
instant [89]. Hence, a number of modifications to the classical nucleation theory

have been suggested [114], [169], and [90].

• Heterogeneous nucleation:

In real world scenarios, neither the bulk liquid is ever free of impurities, nor

the nozzle walls of the nozzles are perfectly smooth. Hence, bubble generation
in the bulk liquid is observed at a much lower degree of superheat compared
to the homogeneous nucleation. This nucleation process is known as hetero-

geneous nucleation. The heterogeneous nucleation process is characterized by

lower bubble radius compared to the homogeneous nucleation [201]. Unlike the

homogeneous nucleation, the heterogeneous nucleation limit curve is observed

to be farther from the spinodal curve and closer to the saturation curve. This
is observed because the depressurization during the heterogeneous nucleation

stops well before reaching the spinodal curve [161]. During the heterogeneous
nucleation (flashing process), rapid nucleation and bubble growth are encoun-
tered at relatively low superheat and tend to increase the pressure due to the

expansion of vapor bubbles, resulting in a halt in the depressurization process.

Heterogeneous nucleation necessitates a nucleation mechanism due to thermo-

dynamic fluctuations, however, at a much lower value of activation energy. To

bridge the gap of activation energy, Alamgir and Lienhard [107] and Skripov[206]
introduced a heterogeneity factor ϕ and modified the critical work needed for

nucleation. According to them, the nucleation flux(JHET,B ) can be computed

based upon the classical homogeneous nucleation theory, as shown below.

 
Wcr ϕ
JHET,B = J0 exp − (2.10)
kT

14
 Numerous studies have been conducted on bubble nucleation during surface
boiling. Among them, Bankoff [20] studied the theoretical thermodynamic

aspects of nucleation processes and found the energy required to initiate the
bubble nucleation is highly dependent upon the shape of the geometrical sur-

face, i.e., whether the surface has a cavity, a flat plane, or a protruding point.
Blander and Katz have proposed a model for heterogeneous nucleation flux
(JHET,smoothwall ) for a smooth surface mentioned in Eqn. 2.11 . They intro-

duced an extra geometrical factor η, which is a function of solid-liquid contact

angle.
r  
2/3 2σ Wcr ϕ
JHET,smoothwall = NM .η. .exp − (2.11)
πmM Bϕ kB Tl

Subsequently, several other models have been developed to account for the sur-
face roughness and cavity of the wall.

2.1.2 Bubble growth

Once bubbles which meet all the stability criteria are created, their growth period
initiates. One of the early and most popular theories which explains the bubble growth
process, is the mono-layer theory. Based on this theory, a bubble that is submerged

in a superheat thermal boundary layer removes that layer locally. Consequently, the
instantaneous superheat degree is diminished [56]. The vapor growth phenomena are
nonlinear problems, as they involve moving boundaries and a strong coupling between

the liquid momentum and the heat transfer through the interface. Hence, it is difficult
to get analytical expressions for the same without sufficient assumptions.

2.1.2.1 Modeling

The bubble growth period can be subdivided into three sub regimes [16]. The
initial bubble growth regime is known as the surface tension regime, where the forces

acting on bubbles are in equilibrium. Hence, the bubble growth rate is slow. Although

15
the initial period is very short, disregarding it leads to overpredictions of the bubble
growth rate for subsequent phases [21].

As the bubble grows, the effect of surface tension diminishes and the growth rate is
limited by the momentum transfer, i.e., how fast the bubble can push the surrounding

liquid. This growth regime is known as the inertia regime. During this regime, the
bubble attains its maximum velocity. Rayleigh [168] first came up with a model stated
in the Eqn. 2.12 to calculate the bubble growth rate for a spherical bubble without

accounting for the surface tension and viscosity.

3 Pv − Pamb.
RR̈ + Ṙ2 = . (2.12)
2 ρl
In Eqn. 2.12, R, Ṙ, and R̈ represent the bubble radius , its first and second

derivatives with respect to time, respectively. Further, Pv , Pamb , and ρl represent

the pressure inside the vapor bubble, ambient pressure, and the liquid density, re-
spectively. Eqn. 2.12 has since been generalized by including the effects of surface

tension and viscosity by Plasset and Zwick [162] and the resulting relation is often
called the Rayleigh-Plesset equation [37], given by

3 2 1 2σ 4µ 
RR̈ + Ṙ = Pv − Pamb. + − Ṙ . (2.13)
2 ρl R R
The final stage of the bubble growth happens to be the slowest of all stages and
is dominated and restricted by the thermal diffusion through the bubble-liquid inter-
face. Mikic et al. [125] were the first to develop a correlation by successfully coupling

the effects of inertia and heat transfer by neglecting the effect of the bubble growth
acceleration term (R̈). Their correlation suffers from the assumption of the linearized
Clausius-Clapeyron equation that describes the relationship between the vapor pres-

sure, temperature and the constant vapor density. This assumption of the linearized
relation between the vapor pressure and the constant vapor density is difficult to
justify in the case of thermal non-equilibrium. Hence, disagreements between experi-

16
ments and model predictions are observed for conditions with large initial super-heat,
or conditions where the variations in the vapor density are significant. This anomaly

was addressed by Miyatake et al. [126] who improved on the works of Theofanous
and Patel [220] by incorporating the R̈ term and also by accounting for the nonlinear
relation between the vapor pressure and the temperature. They also incorporated

the growth rate during the surface tension regime. Furthermore, Riznic et al. [179]
observed the bubble growth rate to be higher in flowing systems compared to static

systems due to the relative velocity between the phases. Their correlation based upon
the relative velocity between phases (vr ) is mentioned below.

r


2 vr t
Ṙ vr ̸=0
= Ṙ vr =0
1+ (2.14)
3 R

2.2 Modeling of flashing flows

Since the middle of the last century, a significant amount of research has been
dedicated to the modeling of two-phase flow mixtures produced by flashing flows.
Owing to the biphasic nature of flow, thermodynamic and hydrodynamic nonequilib-

rium pose the biggest challenge to the accurate modeling of the underlying physics
while representing practical scenarios. Modeling the velocity slip and heat transfer
between phases, vapor generation, and their interactions add complexity to the posed

problem. Based upon the interactions between phases, the flow regime can be divided
into bubbly flow, churn flow, and slug flow. In many industrial applications, there
exist no distinct boundaries between the different flow regimes, which is again an

extra disadvantage in the modelling of flash-boiling process. Hence, many researchers

have assumed different flow regimes as per their convenience to develop their models,
which hinders the universality of these models. The early flash-boiling models were

focused to estimate the mass flow rate and discharge coefficient without revealing
much details about the physics of the flow.

17
 Based upon the hydrodynamic aspects, the flash-boiling models can be classified
into homogeneous flow model, separated flow model, and two fluid models, which will

be explained below [201], [178], [46].

1. Homogeneous flow model: The model assumes hydrodynamic equilibrium be-

tween phases, that implies both the phases are assumed to have equal velocity
and pressure. These assumptions also satisfy the conditions of rapid momentum

and mass exchange. Consequently, the resulting governing equations resemble

those for a single fluid. The fluid properties such as density, viscosity, thermal
conductivity etc. for the mixture are obtained as a result of either mass or vol-

ume averaging. The homogeneous flow model is preferred for the cases where
one phase is finely dispersed in the other. However, the model runs into accu-
racy issues when applied to scenarios with high pressure fluctuations or rapid

acceleration.

2. Separated flow model: This model is applied to the cases with high density
ratio between the phases or low pressure systems, as the assumption of equal
velocities between the phases become erroneous. Hence, the velocities of the

phases are relaxed primarily by introducing a slip ratio (ratio of velocities).

This is accomplished by solving different momentum equations for both the
phases. It also accounts for the drag at the interface due to the relative velocity

between the phases, which is often calculated by empirically correlating the

interfacial shear stress and the slip ratio or the velocity difference. The drift
flux model [84] is a special case of such model, because it predicts void fractions

based upon given velocity differences.

3. Two fluid model: The two-fluid model is considered to be a sub-category of

multi-fluid models, also known as the 6-equation model. The governing equa-
tions consist of two sets of separate equations for continuity, momentum, and

18
 energy. Additional constitutive relations are solved to model the equilibrium
and interactions between the phases. Approaches similar to that of the sep-

arated flow models are used to model the differences in the velocities of the
phases, whereas the differences in the temperatures are solved through an en-

ergy balance equation. The differences in temperatures arise due to the time lag
between the phase change process and the flow inertia. The interfacial phenom-
ena, such as mass flux exchange due to evaporation or condensation, changes

in interfacial energy because of the surface tension dependency on temperature,

and interfacial deformations due to rapid depressurization in the nozzles, give
rise to the pressure non-equilibrium. Due to the above complexities, pressure

non-equilibrium is often ignored while modeling the two-phase flows. However,

the non-equilibrium becomes significant when the velocity either approaches or
exceeds the two-phase speed of sound. The model also suffers from empiri-

cism problems for complex flow situations where heat, mass, and momentum
exchange exist simultaneously [178]. To capture more underlying physics of
flash-boiling through this class of model, one needs to add more constitutive

relations. These constitutive relations are often non-trivial and modeling them
properly is the biggest challenge one faces.

In the majority of industrial applications, the vapor phase is finely dispersed

in the liquid phase in the case of flash-boiling flows. Thus, they are modeled as
homogeneous flows hydro-dynamically with an appropriate thermal nonequilibrium

model acting as a constitutive relation to close governing equations [201]. Hence, the
primary objective of a flash-boiling model is to account for the thermal nonequilibrium
between the phases by estimating nucleation and vapor generation. Based upon

the vapor generation model, the existing models can further be classified into three
categories, i.e., bubble dynamics-based models, thermodynamic rate-based models,
and interfacial exchange model.

19
2.2.1 Bubble dynamic based model

Bubble dynamics based models for vapor generation include submodels for nucle-
ation, bubble growth rate, and bubble number density. Models for the first two parts
have already been explained in the earlier sections 2.1.1.2 and 2.1.2.1, respectively.

The bubble number density, denoted as NB , affects the interfacial area per unit
volume in bubbly flow. Different correlations for NB have been developed by Lienhard
et al.[107], Riznic and Ishii [179], and by Shin and Jones [204]. However, these models

include several empirical constants.

The model developed by Kishii et al.[94] moved away from the empiricism of
calculating the local number density by solving a conservation equation in a control

volume approach , and is represented in the Eqn. 2.15.

∂NB ∂(NB U )
+ = Ψ, (2.15)
∂t ∂x
In the Eqn. 2.15, U is the velocity, and Ψ is the source term which is evaluated by

accounting for the change in the bubble number density due to different nucleation
processes and bubble interactions. Another approach, worthy of mention, corresponds
to the work of Elias and Chambre [64, 65] who solved an equation similar to the

Eqn. 2.15 for the concentration function NRB . It is a distribution function for the
bubble number density lying between an interval of bubble radius. The temporal
growth of the concentration function of Elias and Chambre [64, 65], is governed by

the bubble growth, given by:

∂NRB ∂  DR 
− = NRB (2.16)
∂t ∂R Dt

Furthermore, the vapor void fraction α is determined by using the bubble radius
(R) predicted by the bubble growth model and the bubble number density as shown

in the Eqn. 2.17.

20
 4π 3
α= R NB . (2.17)
3
Bubble dynamics based models suffer from several levels of empiricism, lack of

general validity, and are often tuned to replicate the available experimental data. In
the past, these models have assumed thermal equilibrium while being applied to one-
dimensional homogeneous flows. However, recently these models have been applied

to multidimensional CFD codes as submodels to describe the flash-boiling processes

[120], [87] and [109].
The assumption of thermal equilibrium can only be applied to processes with

very fast heat transfer rate and driven by the fluid temperature. In the case of a
cold fluid undergoing the phase change due to depressurization, the vapor density
in the bubble is insignificant when compared to the liquid density. Therefore, the

activation energy of vaporization is insignificant. This is usually encountered in the

case of cavitating flows [93]. However, the vapor density is much higher in the case
of superheated bubbles, thus limiting the rate of heat transfer. Hence, it renders

the assumption of thermal equilibrium imprecise for flash-boiling flows. Therefore,

thermal nonequilibrium becomes an important characteristic of flash-boiling models
to achieve better accuracy. A few notable nonequilibrium models based on the bubble

dynamics approach include the ones by Blinkov et al. [33], Elias et al. [64] and Ritcher
[178]. However, these models are limited only to one-dimensional flows, to date.

2.2.2 Thermodynamic rate based models

From the above discussions, it is evident that the majority of the bubble dynamics-
based models are semianalytic in nature and have a large degree of empiricism. De-

termining the void fraction (α) in such conditions is a cumbersome process. The
accuracy of these models also depends upon the reliable analytical derivation of the
bubble growth rate. To avoid these complexities, another class of models, which are

also empirical in nature, are often used to determine the vapor generation rate. These

21
models are termed as relaxation models, since they account for the transition of the
system from a state of thermodynamic nonequilibrium to equilibrium on the basis of

the vapor mass fraction.

Homogeneous Equilibrium Model (HEM), developed by Wallis et al. [228], models

one extreme of this relaxation process, where the phase change is instantaneous due
to an extremely fast heat transfer rate. HEM predicts the critical mass flow rate
accurately in the case of long pipes where the fluid has enough time to approach the

state of equilibrium. However, this model fails when applied to the flow in short pipes
or orifices or processes with very slow heat transfer rate. Flow through such a system
is so fast that the time required for the phase change process to reach an equilibrium

state is insufficient. On the contrary, the Homogeneous Frozen Model(HFM) proposed

by Henry and Fauske [76] is an empirical nonequilibrium model which is appropriate
for short channel flows. Whereas, the Homogeneous Relaxation Model (HRM) models

the scenarios which lie in between the two extremes modeled by the HEM and HFM.
HRM was originally developed by Bilicki et al. [30] for one-dimensional two-phase
flows to model the complex process of relaxation by determining the total derivative

of the quality (x), represented in the Eqn. 2.18.

Dx x̄ − x
= (2.18)
Dt Θ

In the Eqn. 2.18, x̄ refers to the equilibrium vapor mass fraction and is a function of
the local and saturation enthalpies of the liquid. It is again determined by using the

Eqn. 2.20. Furthermore, x is the instantaneous mass fraction of the vapor, and is a
function of the void fraction (α), vapor density (ρv ) and the mixture density (ρ), as
shown in the Eqn. 2.19. Finally, Θ is the relaxation time.

αρv
x= (2.19)
ρ

22
 h − hl
x̄ = (2.20)
hv − hl
In the Eqn. 2.20, h, hl , hv are the local enthalpies of the mixture, liquid, and vapor,
respectively. The void fraction in turn depends upon the local density of the mixture
(ρ), vapor (ρv ) and liquid (ρl ) and is given by Eqn. 2.21.

ρl − ρ
α= (2.21)
ρl − ρv
Downar-Zaploski et al. [58] developed a correlation for the relaxation time (Θ) by
performing an empirical fit of their experimental observations of flashing water flows

in long pipes as shown in the Eqn. 2.22.

Θ = Θ0 α a ψ b (2.22)

In the Eqn. 2.22, ψ is the non dimensional pressure, where Θ0 , a, and b are constants
which have different values for high pressure and low pressure applications as repre-

sented in the Eqn. 2.23 and 2.24 .

For lower pressure, i.e., P ≤ 10 bar, the recommended coefficients are:

Psat − P
Θ0 = 6.51 × 10−4 s, a = −0.257, b = −2.24, ψ = (2.23)
Psat

For higher pressure, i.e. P ≥ 10 bar, the recommended coefficients are:

Psat − P
Θ0 = 3.87 × 10−7 s, a = −0.54, b = −1.76, ψ = (2.24)
Pcr − Psat

HRM overcomes the need for intricate multistep modeling dictated by bubble dy-
namics based methods to determine the vapor generation. In HRM, α determines the
initial vapor generation due to phase change, whereas the non-dimensional pressure

23
term ψ is indicative of the amount of energy available for additional vapor gener-
ation due to superheat. The dependence of α on Θ corroborates the fact that the

initial vapor generation rate in the domain provides surface area for additional vapor
formation [135].

The one-dimensional HRM is appropriate for flow through pipes with large length-
to diameter ratios. In these pipes, the axial variation of the flow parameters determine
the thermal-hydraulic behavior of the jet, but the effects of cross-sectional parameters

are insignificant. However, the cross-sectional parameters become significant when the
length and diameters of the nozzles are of the same order. Such conditions necessitate
the need for multidimensional vapor generation models.

Schmidt et al. [191], [189] were the first one to extend the one-dimensional HRM to
two-dimensional forms and later to the multi-dimensional form in an Eulerian frame
of reference. As shown in the Eqn. 2.25, they used the chain rule to express the total

derivative of the density instead of the conventional equation of state. This expression
allows the pressure to respond to the compressibility of the individual phases and the
phase change phenomena.

Dρ ∂ρ Dp ∂ρ Dx ∂ρ Dh
= + + . (2.25)
Dt ∂p x,h Dt ∂x p,h Dt ∂h p,x Dt

2.2.3 Interfacial exchange based model

One of the popular approaches is to postulate that the phase change process is
induced by the interphase heat transfer. The vapor generation rate (Γg ) during the

evaporation process is related to the heat flux through the interface (q̇) via Eqn. 2.26.

q̇
Γg = Ai (2.26)
L

24
 In the Eqn. 2.26, Ai represents the interfacial area density, which requires further
modeling according to the flow regime .

Estimating the heat transfer coefficient on the liquid side is the biggest challenge
one faces while evaluating the vapor generation rate using the two-fluid approach.
Hence, the applicability of these sophisticated methods is restrained largely by the

reliability of closure models which represent interactions between the phases.

2.3 Experimental studies

The process of flash-boiling atomization has been studied experimentally for
several decades now. The current section reviews some of the early experimental

investigations to classify the flash-boiling process. It also discusses some of the ex-
perimental studies performed to gain insight into the effects of flash boiling on the
performance parameters of plain/straight nozzles and fuel injectors.

2.3.1 Flash-boiling mechanism

Brown and York [38] were the first to report flash-boiling atomization while
performing experiments with water and Freon-11 in straight nozzles of different length

to diameter ratios. They also concluded that flash boiling is an effective mode of
atomization without the need of high pressure and high velocities.
Wildgen and Straub [232] studied flashing water jets and identified two modes of

flashing, i.e., particle boiling and surface boiling. During particle boiling, the bubbles
in the water jets are formed due to the presence of suspended particles. However,
surface boiling is a wall phenomenon which is observed in the absence of suspended

particles. According to them, the surface boiling or wall boiling occurs only in longer
nozzles.
Oza et al. [146][147] observed different regimes of flash boiling while performing

experiments in an electromagnetic (poppet) injector with fluids of different volatility,

25
i.e. propane, methanol and indolene. They classified the regimes into internal flashing,
when flashing is initiated within the injector, and external flashing, which is external

to the nozzle with an intact liquid core surrounded by droplets.

However, Reitz et al. [171] refuted the claims of Oza et al. [146][147] about the

internal flashing. He conducted photographical studies to gain further insight into

the flash-boiling mechanism in a single-hole cylindrical nozzle. He used water as the
working fluid at a constant pressure gradient but at different liquid temperatures.

He also observed an intact liquid core upto certain distance from the exit of the
nozzle. The length of the liquid core reduced with the increasing inlet temperature.
His study was also able to elucidate a decrease in the drop size and mass flow rate

with the increase in inlet fluid temperature. Eventually, they observed a vapor-locked
phenomenon when the inlet fluid temperature approached its boiling temperature.
Park and Lee [151] also used transparent cylindrical nozzles to investigate the

flash-boiling mechanism with water as the test fluid. As they increased the superheat
defined by the Eqn. 2.27, they observed different internal flow regimes ranging from
bubbly flow, slug flow, to annular flow.

∆T = Tinj. − Tsat. (P∞ ), (2.27)

In the Eqn. 2.27, Tinj. is the injection/inlet temperature of the liquid and Tsat. (P∞ )
is the saturation temperature of the fluid corresponding to the chamber pressure. As

the degree of superheat increased, the bubbles created near the nozzle wall during the
bubbly flow regime started to interact with each other forming slugs of vapor in the
liquid core. This regime was classified as the slug flow regime, which later translated

into the annular flow regime with an increase in the superheat. In the annular flow
regime, a distinct vapor core surrounded by an annular liquid film near the wall is
observed. According to Park and Lee [151], slug and/or annular flow regimes can be

witnessed in long nozzles at low injection pressure and higher injection temperature

26
compared to that of the bubbly flow. Long nozzles provide longer residence time for
the bubbles to grow and coalesce compared to that of short nozzles. Consequently,

bubbly flow is observed in short nozzles. Furthermore, slug and annular flow regimes
are characterized by larger spray angles and finer droplets in comparison to the bubbly

flow. In bubbly flow, small droplets are observed around an intact liquid core. Peter
et al. [157] also observed similar effects in the shattering of liquid jets due to different
flow regimes.

Kitamura et al. [92] studied the break-up of flashing water and ethanol jets through
long straight nozzles. According to their hypothesis, bubble nucleation has negligible
effects on the jet break-up process at low superheat and the break-up is dominated

by mechanical (aerodynamic) effects. However, the contrary was true for the jet
break-up process under high superheats. Furthermore, they were able to arrive at a
correlation to determine this transition from a mechanical force-dominated jet break-

up to a bubble nucleation dominated jet break-up on the basis of the nondimensional

superheat. The nondimensional superheat is expressed by the Jakob number (Ja), as
shown in the Eqn. 2.28

Available sensible heat energy ρl Cp ∆T

Ja = = (2.28)
Available latent heat energy ρv hf g

Based on the superheat level, Cleary et al. [43] classified the jet break-up regimes
into subcooled, transitional, and flare-flashing regimes. In the subcooled regime,
the water jet shattered far from the nozzle due to mechanical forces, as shown in

the Fig. 2.3A. However, an intact liquid core for a short distance from the exit of the
nozzle was observed in the transition regime, as seen in the Fig. 2.3B. This liquid core
further shattered into tiny droplets forming an wide spray angle. In the third regime,

i.e., in the Fig. 2.3C, Cleary et al. [43] observed violent break-up of the liquid jet
immediately after the nozzle exit because of bubble interactions and explosion. This

27
condition was termed as the flare-flashing regime where thermal effects dominated
the jet break-up process.

Figure 2.3: Jet transition from mechanical breakup to flare flashing, as identified by
Cleary et al. [43]. Image adopted from [43].

Furthermore, from a thermodynamical point of view, Lamanna et al. [100] refuted

the earlier definitions of liquid superheat and suggested a new formulation involving
pressure ratio. They defined the pressure ratio as the ratio of the saturation pressure
of the fuel at the injection temperature to the ambient pressure.

2.3.2 Flash-boiling and fuel injectors

Many of the experiments described in the Sec. 2.3.1 were conducted on plain/straight
orifices with water as the test fluid and were successful in identifying the different

28
modes of flash-boiling sprays. However, they convey a little about the spray perfor-
mance parameters generated by fuel injectors in flash-boiling conditions. Vanderwege

et al. [226], and Schmitz et al. [195] conducted experiments in pressure swirl atom-
izers with gasoline as the test fluid. Their observations threw light on the effects

of fuel volatility and nozzle geometry on the spray structures. They concluded that
the higher volatile component of the liquid fuel is more susceptible to flash boiling
causing disintegration of the liquid jet, resulting in smaller droplets. Further studies

by Allen et al. [9] corroborated their observations and concluded that gasoline is more
susceptible to early inception of vapor generation than diesel primarily due to its high
volatility.

Since then, the focus of flash-boiling studies in the GDI systems has shifted to
the valve-covered orifice (VCO) and sac-type injectors. These injectors have been
the subject of an immense amount of scrutiny in the past, because of their continued

application in diesel injection systems [188]. Flow asymmetry is one of the important
features of these injectors. The rotational behavior of this asymmetric flow results in
strong vortex cores between the inlet, nozzle and sac area. The low-pressure region

in the vortex core results in the string-shaped cavitation. This kind of flow field was
initially identified by Arcoumanis et al. [15] in diesel injection and was also observed
by Gilles-Birth et al. [71] in their study using gasoline injection through VCO nozzles.

However, these studies do not provide any explanation on the coupling or the influence
of string-shaped cavitation on flash boiling.
Aleiferis et al. [8, 40, 198] used different hydrocarbon fuels (gasoline, iso-octane,

n-pentane, ethanol, and butanol) in an optical injector to study the in-nozzle cavita-
tion and its effects on flash boiling and subsequent spray formation. They identified
the higher vapor pressures because of the elevated fuel temperatures to be the cul-

prit behind the intense in-nozzle cavitation. Subsequently, the rigorous expansion
and interactions of bubbles led to rigorous flashing in low-pressure environments.

29
This coupling of cavitation and flash boiling with the effects of the injector geometry
resulted in an asymmetric flow with large spray angles. They also elucidated the sig-

nificant role that the ambient/gas pressure played while deciding the degree of spray
superheating. According to them, the cavitation increased the spray angle, but had

insignificant effects on improving the spray quality, which was more influenced by the
superheating of the fuel. Aleiferis and co-workers further suggested that the satu-
ration pressure of the most volatile component of the mixture should be considered

while determining the degree of superheat, a claim that was later refuted by Araneo
et al. [13, 14], who suggested that the saturation pressure of the multicomponent fuel
is influenced by the saturation pressure of individual fuel components.

Due to the limitations in spray diagnostic techniques, many experimental stud-

ies [129, 237, 240, 241] often report the quantitative aspects of the spray, such as
the drop size, spray angle, and penetration at several distances from the nozzle exit,

where the spray is either thin or dilute. However, the application of these information
in spray modeling is a futile exercise, since the effects of flash boiling begin to wane
out as we proceed downstream. Zeng et al. [240] were the first to report the effects

of flash boiling on the spray plume angle, which increases with the increase in su-
perheat level during the transition regime, but starts decreasing in hard flash-boiling
conditions. This difference in trend can be attributed to the improper measurement

of the spray width without considering the transition distance from nozzle exit [13].
Recent experimental study by Araneo et al. [13] shows that the spray angle expands
in the near field region and starts contracting in the far field for multihole injectors

under flash-boiling conditions. Therefore, they have suggested to be careful while

considering the transition region, which is a function of the superheat. This also sup-
ports some of the past ambiguities observed by Lamanna et al. [100], Wu et al. [235]

and Montanaro et al. [130]. Wu et al. [235] have also reported a wider near-nozzle
spray for the 2 hole GDI injector compared to the single hole injector. A similar

30
trend for the 3-hole and 6-hole injectors were also observed by Kramer et al. [96].
Khan et al. [91] have reported an adverse effect of the injector opening velocity on

the near-nozzle spray width in flashing conditions.

Similarly, there exist ambiguities in the literature of spray penetration measure-

ments. Montanaro et al. [130] found a lesser spray penetration in the flash-boiling
condition compared to that in the non-flashing condition. Chen et al. [41] have ob-
tained a similar trend in their novel study of heated tip GDI injectors. Khan et

al. [91] found an insignificant effect of the injector opening velocity on the spray pen-
etration under hard flashing conditions. However, Zeng et al. [240] observed favorable
effects of the superheating of the fuel on the spray penetration under hard flashing

conditions.
Furthermore, experimentalists have observed the spray collapse event for multihole
injectors under flashing conditions. It should be noted that the spray collapse event is

detrimental to the combustion process. Collapsing sprays are often characterized by

narrowing spray width, which prohibits proper air-fuel mixing, resulting in improper
combustion and increased emissions. During spray collapse events, multiple spray

jets interact with each other to form a single jet of spray. Consequently, the axial
spray penetration increases, which could potentially lead to the film formation on
the piston head. However, the process of spray collapse and its inception is less

understood among the communities and the researchers are of different opinion for
the reason of spray collapse inception. Some of these reasons include, but are not
limited to, the possibility of a low pressure region in the injector center-line [242],

smaller droplets experiencing lesser resistance from the aerodynamic forces [226, 195],
nozzle hole configuration [129, 11], and jet-to-jet interactions [96]. Furthermore, Guo
et al. [74] suggested jet-to-jet interactions during the flashing condition might induce

a local low-pressure region at the injector tip, resulting in the bending of spray plumes
towards the center-line of the injector axis. They also suggested that the condensation

31
process inside the spray plumes can also play a significant role in inducing the spray
collapse, which needs further validation. Recently, Lacey et al. [97] have found that

the nozzle diameter, drill angle, and specific volume of the fuel play an important
role in inducing jet-to-jet interaction in multihole injectors, which results in spray

collapse.

2.4 CFD study of flash-boiling

The previous section 2.3.2 shows that the majority of experimental studies talk
about the effects of flash boiling on external spray characteristics, but the studies on
the internal flow behaviors of nozzles and near nozzle characteristics are pretty limited

because of the experimental constraints. The near nozzle flow physics are important
driving parameters for the external flow characteristics. Hence, many CFD studies
try to simulate the near-nozzle behaviors coupled with the internal flow. A few of the

prominent CFD studies with a concentration on flash boiling have been discussed in
the current section.
Janet et al. [87] applied a bubble dynamics-based two-fluid model to simulate

the flow through a converging diverging nozzle in flash-boiling conditions. They also
incorporated a heterogeneous wall nucleation model and an interfacial heat transfer
model to account for the vapor generation. Their observations of the wall nucleation

and axial pressure distribution were in good agreement with the experimental inves-
tigations. However, their simulations under predicted the radial distribution of the
vapor phase owing to the low bulk interphase mass transfer rate.

Similarly, Liao et al. [109] applied a bubble dynamics-based two-fluid model with
a constant bubble number density to simulate the flashing flow through a converging
diverging nozzle. They were able to capture the flashing inception with a good agree-

ment with the experiments but predicted a lower pressure undershoot and uniform

32
radial distribution of the vapor phase because of the assumption of constant number
density.

The computationally expensive nature of the bubble dynamics models limits their
applications in complex fuel injectors. Recently, Shen et al. [200] applied a bubble

dynamics-based model to study the flash-boiling spray of a multihole injector without

simulating the complete geometry. They were able to predict a low-pressure region
near the injector tip resulting in spray collapse.

Schmidt et al. [189] have successfully validated the two-dimensional Homogeneous

Relaxation Model against the experiments conducted by Reitz [171] and were able
to accurately predict the choking flow and vapor lock phenomena. Since its multidi-

mensional extension, HRM has been successfully applied to the study of flash boiling
in pressure swirl atomizers [135], and super-heated jet fuel [105]. Since then, many
researchers have successfully simulated the cavitating flow through single-hole and

multi-hole diesel injectors [238, 245, 24, 22, 23]. Subsequently, the original capability
of HRM code has been improved to simulate the effect of ambient gas on near nozzle
flow characteristics by solving an extra transport equation for the mass fraction of

the non-condensible gas (NCG) as the third phase. The improved model has been
successfully validated against the phase-contrast X-Ray imaging (PCI) of the internal
nozzle flow [60] in both the flashing and non-flashing conditions.

Moulai et al. [132], Strek et al. [213], Saha et al. [183, 185] and Baldwin et al. [19]
have simulated the flashing and non-flashing sprays through a multihole GDI noz-
zle for the Spray G target conditions put forward by the Engine Combustion Net-

work(ECN) [1] using the HRM model. Moulai et al. [132] showed the influence of
the counterbore and the differences in the gas ingestion process between flashing and
non-flashing conditions. Strek et al. [213] compared the hole-to-hole variations in

the mass flow rate and density with the X-Ray tomography data. Saha et al. [183]
studied the effect of needle lift on the internal flow pattern in both flashing and non-

33
flashing conditions using static grids of different needle elevations. They observed
that the flashing flow pattern is more affected by the low needle lift compared to the

non-flashing one. Baldwin et al. [19] applied the HRM model to the ECN prescribed
GDI injector with a transient needle motion and predicted the string flash boiling

and swirling spray with a nonequilibrium thermal core.

2.5 Inferences
The following inferences can be drawn from the review of the various experi-

mental investigations of flash-boiling sprays:

• Flash-boiling significantly reduces the drop size while improving the spray qual-
ity.

• Internal flash boiling is primarily observed in long nozzles subjected to higher

degrees of superheat.

• The occurrence of cavitation increases the degree of flashing and near nozzle
spray width.

• In flash boiling conditions, the near nozzle spray width has a different trend

compared to the far nozzle spray width.

• In multihole injectors, an increase in superheat although helpful in reducing the

drop size, but is often characterized by spray collapse because of the higher jet-
to-jet interactions. Collapsing spray may result in increased spray penetration

and may lead to spray impingement on the piston head.

The following can be summarized from the discussion of the modeling aspects of
flash-boiling and two-phase flows in general:

• The bubble dynamics-based models suffer from severe empiricism and are often
restricted to one-dimensional flow. Mostly, two-fluid models lack the desired

34
 accuracy because of the unavailability of proper closure models. These models
are computationally expensive and are often considered for simpler geometries.

• Multidimensional modeling is necessary for fuel injectors, primarily to study

the influence of the cross-sectional flow parameters.

• Thermodynamic rate-based models can be applied to a range of thermody-

namic behaviors, i.e., ranging from equilibrium to frozen flows. These models
have been observed to perform better compared to the bubble-dynamics-based

models under flash-boiling conditions in automotive spray applications [39].

• Multidimensional nature, a wide range of applicability, and validity of the ho-

mogeneous relaxation model (HRM) can be stated as the reasons for claiming
that a level of maturity has been achieved in modeling internal nozzles both in

diesel and gasoline direct injection systems.

35
 CHAPTER 3

ATOMIZATION MODELS: A REVIEW

The spray and atomization process is a multiphase flow phenomenon, which

includes the gas phase as the continuum, whereas the liquid phase is in the form of
discrete droplets and ligaments. The kinetic energy of the spray represents the main
source for the production of turbulence which governs the microscale air-fuel mixing.

The disintegration of the continuous liquid phase occurs when the external disruptive
forces such as aerodynamic force, surface shear force, and centrifugal force acting
on the liquid surface exceed the surface tension force. This initial breakup process is

often termed as primary break-up or primary atomization. Liquid ligaments and large
droplets produced during the primary atomization event are unstable and undergo
further break-up to form smaller stable droplets. This process is identified as the

secondary breakup or secondary atomization.

In traditional fluid mechanics, the spray behaviors are conveniently character-
ized with several nondimensional parameters, such as: the Liquid Reynolds Number

(Rel ), Liquid Weber Number (W el ), Aerodynamic Weber Number (W eg ), and the

Ohnesorge Number (Oh), as mentioned below.

r
ρl ul dl ρl u2l dl 2
ρg Urel dl W el µl
Rel = , W el = , W eg = , Oh = =√ (3.1)
µl σ σ Rel ρl σdl

In Eqn. 3.1, ρ, u, and d represent the density, velocity, and diameter of each phase,

respectively. Whereas, the subscripts l and g represent the liquid and gaseous phases.
Furthermore, Urel represents the relative velocity between the phases.

36
 The primary break-up process can be classified into different regimes based on the
previously mentioned nondimensional numbers, as discussed below.

• Rayleigh regime: This regime is usually witnessed for flows with low Reynolds
number, where identical droplets are produced. In this regime, the jet disinte-

gration is driven by the surface tension force.

• Aerodynamic regime: At intermediate Reynolds number, the drop formation is

influenced by the aerodynamic forces. These forces result in symmetric (first

wind-induced mode) and asymmetric (second wind-induced mode) wave growth

on the gas-liquid interface leading to jet disintegration.

• Atomization regime: This regime is generally encountered for flows with very
high Reynolds number, where the jets get disintegrated immediately after the
nozzle exit. The sizes of the resulting droplets are much smaller than the nozzle

diameter. All practical automotive sprays fall in this regime, thus, the current
thesis explores models for the atomization regime.

Images of the different breakup regimes can be seen in Fig. 3.1.

The atomization regime can be described by two regions, i.e., the dense spray
region and the dilute spray region, as shown in Fig. 3.2. The dense spray region is
observed close to the nozzle exit, where a mixing layer surrounds the distinct liquid

core. The mixing layer consists of liquid ligaments and big unstable droplets, which
undergo further break-up processes. As we proceed farther from the exit of the nozzle,
the droplets formed from the primary and secondary atomization processes undergo

further break-up to form stable droplets. Consequently, the spray region gets wider
and sparse. This region is known as the dilute spray region. Droplet evaporation
becomes dominant in this region. Computational scientists have put a lot of effort

to model these different regions of atomization. A few of the famous approaches are
discussed below.

37
Figure 3.1: Four main regimes of round jet breakup, namely the Rayleigh regime,
the first wind-induced regime, the second wind-induced regime, and the atomization
regime characterized by the Oh number and Re number [118]

3.1 Lagrangian-Eulerian methods

As the name suggests, the Lagrangian-Eulerian (LE) approach treats the dis-
crete liquid phase as the Lagrangian particle, whereas the continuous gas phase is
treated as the Eulerian phase. Since its inception, the LE approach has become an

industrial standard to model the spray dynamics, owing to the simplicity of imple-
mentation and high computational efficiency. The LE spray models usually focus
on tracing the trajectories of individual droplets or parcels alongside the continuum.

However, the suitability of the LE method to the dense spray region is a topic of
debate, as these models are droplet oriented [18]. Hence, the formulations for these
model don’t consider the interactions between neighboring droplets.

3.1.1 Statistical description of the spray equation

Williams [233] was the first one to introduce a statistical description for the

stochastic behavior of sprays. The Williams spray equation includes formulations for

38
 Figure 3.2: The spray regions of an atomizing jet. Image adopted from [66].

the formation of new droplets, their growth rate, collision between droplets, and the
aerodynamic forces acting on them. He used a probability density approach to model
the spray droplets, represented by the PDF, f (t, x, y, z). Here, f is the probable

number of droplets per unit volume in time and space, described by its location (x),
velocity (v), radius (r), the temperature (Td ), the deformation parameter (y), and
the rate of deformation (ẏ). The Williams spray equation is mentioned below.
 
∂f ∂ (f ṙ) ∂ f T˙d ∂ (f ẏ) ∂ (f ÿ)
+ ∇x . (f v) + ∇v . (f v̇) + + + + ˙ + S˙bu (3.2)
= Scoll
∂t ∂r ∂Td ∂y ∂ ẏ

In Eqn. 3.2, ṙ, v̇, and T˙d are the time rate of change of the droplet radius, droplet

velocity and the droplet temperature, respectively. The droplet acceleration (v̇) is
further calculated by determining the drag force exerted on the liquid droplets by the
continuum. Two notable drag models, i.e., the spherical drag model or the droplet

drag model are usually used in engineering applications to calculate the drag force.
These models compute the drag force by determining the relative velocities between
the droplets and their drag coefficients. Further, ṙ is evaluated by applying the

principles of convective heat and mass transfer to the phenomena of evaporation or

condensation. Additionally, the fundamentals of energy balance help in determining

39
the T˙d , where the energy supplied to the droplet either raises the droplet temperature
or complements the latent heat of vaporization. Drop distortion parameters (y, ẏ, and

ÿ) are modeled by the Taylor’s drop oscillator model (TAB) [145]. It describes the
distortion of the droplets by a forced, damped harmonic oscillator. Here, the liquid
viscosity and surface tension act as the damping and restoring forces, respectively,

whereas the aerodynamic drag force acts as the forcing term.

The right hand side of the Williams spray equation includes two source terms,
˙ , and S˙bu . The first term, Scoll
i.e., Scoll ˙ , accounts for the collision of droplets and the

second term, S˙bu , considers the droplet break-up process. Two different approaches
have been suggested by O’Rourke et al. [145] and Schmidt et al. [193] for modeling
˙ . The former approach calculates the probability of parcels colliding with
the Scoll

other parcels based upon the size, velocity, and number of droplets present inside
them [145]. However, the latter approach is based upon the No Time Counter (NTC)
method inspired by the gas dynamics fundamentals used in the direct simulations

Monte Carlo techniques (DSMC) [31].

The source term, S˙bu , accounts for breakups due to both primary and secondary at-
omization. Several modeling approaches have been proposed to model these breakup

phenomena. A few of the popular breakup models are listed below.

• Linearized instability models: This class of models are based on the linear sta-
bility theory [170], [25], which assumes the breakup of a liquid jet is dominated
by the growth of the most unstable surface waves. Following this assumption,

Reitz and Diwakar [172] performed a first order linear analysis of a Kelvin-
Helmotz (KH) instability growing on the surface of cylindrical liquid jet pen-
etrating into a quiescent incompressible gaseous environment. Subsequently,

they postulated the KH breakup model for the primary atomization regime.
This breakup model is also known as the blob injection model as the initial

40
 liquid structures are modeled as spherical blobs which have the same size as the
diameter of the nozzle.

For the secondary breakup, Patterson and Reitz [154] proposed the Rayleigh-
Taylor(RT) breakup model inspired by the theoretical work of Taylor [219],

who studied the stability of a liquid-gas interface when subjected to accelera-

tion/deacceleration in the normal direction. According to them [154], the RT
instabilities grows when the fluid acceleration is in the opposite direction of the

density gradient. Subsequently, Baele and Reitz [25] proposed a hybrid KH-RT
model for engine applications, where both the KH and RT breakup model are

implemented in a competing manner.

Although, these models are popular among the industry users owing to their

simplicity, they don’t account for the in-nozzle turbulence or cavitation. Ad-
ditionally, these models are based on linear theory which has excellent success
in predicting the most unstable modes in various canonical two-phase flows,

such as liquid sheets, cylindrical jets, annular jets, liquid films, and liquid
threads [203], [111], which are often laminar in nature. However, application of
these break up models in automotive sprays, which are often characterized by

high Reynolds number and very small nozzle diameters is questionable [4].

• Turbulence and cavitation models: Based on the blob injection models, Huh
and Gosman [80] proposed a turbulence induced breakup model. According to
them, the initial surface perturbations are driven by the turbulent forces within

the liquid. Subsequently, Nishimura and Assanis [139] extended the work of Huh
and Gosman [80] and postulated a cavitation and turbulence-induced primary
breakup model for the diesel spray applications. To account for the cavitation,

they introduced parcels containing bubbles to the spray chamber simulations.

Furthermore, Som and Aggarwal [207] introduced the KH-Aerodynamic Cavi-

41
 tation Turbulence (KH-ACT) model for the primary atomization regime. How-
ever, all these models are still based on the linear stability theory at their core.

Additionally, performances of these models are driven by the heavy tuning of

the model constants for automotive spray applications. Models tuned for one

test condition often fail to accurately predict the spray behaviors in other con-
ditions without tuning of model constants. Hence, these models are not truly
predictive in nature and lacks general applicabilities.

• Taylor-Analogy Breakup model (TAB model): This model was first postulated

by O’Rourke and Amsden [145] based on the the analogy between and oscillat-
ing and distorting droplet and a damped spring mass system. According to this
analogy, the external aerodynamic drag forces and the surface tension forces

acting on a droplet are equivalent to the external force and the restoring force
of a spring mass system, respectively. Whereas, the viscous force of the liquid
droplet acts as a damping force. Based on this analogy, the amplitude of os-

cillations are evaluated for the forced-harmonic oscillator system. Breakup of a

parent droplet is only considered when the amplitude of this oscillation exceeds
the droplet radius. Post-breakup droplet radius is evaluated by accounting for

the energy conservation. Subsequently, Tanner [217] showed under-predictions

of droplet diameters for full cone sprays by the TAB model and proposed an
enhanced version known as the ETAB model. This model reflects a cascade of

droplet breakups where droplet size is reduced in a continuous manner until the
child droplets attain stability.

3.2 Direct Numerical Simulation

The Direct Numerical Simulation (DNS) involves solving the governing equations
for both phases, with the objective to resolve all necessary length and time scales

associated with the flow field. For high pressure injection systems, the flow velocity

42
can approach the speed of sound. Hence, the DNS approach requires to solve the
compressible formulation of the Navier-Stokes equation, as shown in the following

Eqs. 3.3 and 3.4

∇. (ρU) = 0 (3.3)
∂ρU
+ ∇. (ρU ⊗ U) = −∇p + ∇. µ ∇U + ∇T U + ρg + Tσ



(3.4)
∂t
In the above equation, Tσ is the surface tension force, which is nonzero only on

the interface. Hence, it poses a further challenge of tracking the interface in the high-
fidelity modeling of the multiphase flow. Furthermore, in spray modeling, the fluid
interface might be discontinuous, as it constantly moves, deforms, and breaks apart.

Thus, it becomes cumbersome to directly resolve the smallest length scale associated
with the droplets and the ligament breakup process. Therefore, two different ap-
proaches are generally used to describe the motion and the location of the interface.

One of these approaches involves modeling and tracking of the interface with marker
particles in a Lagrangian framework. This approach is known as the interface tracking
method [177]. However, these marker particles require constant rearranging to ensure

their conformity to the interface. As a result, the interface tracking method suffers
from mass conservation and faces difficulties in handling topological changes [187].
The other approach involves capturing the interface either by the volume of fluid

method (VOF) [223] or the level-set method [143].

The VOF methods [78], [142], [187] define an indicator function (α) for the volume
fraction, which carries a value ranging from 0 to 1. In a pure phase, i.e., in gas or

liquid, α carries a value of 0 or 1, whereas any decimal value between 0 and 1 represents
the interface. Furthermore, the VOF methods solve an advection equation for α to
track the interface. However, due to the step behavior of the indicator function,

a geometric or algebraic reconstruction of the interface is applied. It prevents the

smearing of the interface due to numerical diffusion [239], [153], [173], [176], [182].

43
 In level-set methods, a scalar field ϕ is defined in the domain. At the interface, ϕ
is assigned a fixed value of ϕ0 , and in the pure phase ϕ is either greater than ϕ0 or

lesser than ϕ0 . Away from the interface, smoothing functions for ϕ are used to avoid
numerical diffusion. However, level-set methods are not volume conserved in nature.

To avoid these drawbacks, a fine grid resolution or coupling level-set method with the
VOF method has been proposed [77], [216].
Despite the advantages of capturing the spray physics accurately using sharp

interface approaches, these methods are still not popular in engineering applications
because of the following reasons.

• It is expensive to fully resolve the interface for the flow with high Reynolds
number. Automotive sprays such as diesel and gasoline sprays operate in the

high Reynolds number and Weber number regime. hence, these flows are highly
turbulent in nature. DNS simulations of turbulent flows are expensive and
numerically challenging [73].

• Very small bubbles can be present in the spray under flash-boiling conditions,

which render further complexity in the tracking and modeling interface.

• In the near nozzle regions of the diesel and gasoline sprays, a distinct sharp
interface is not physically present [49, 48].

• During atomization, the interface undergoes severe deformations such that the
length scale reduces to zero, e.g., pinching of droplets or breakup of ligaments.

Because of the above mentioned difficulties associated with the atomization pro-

cess, 3-D DNS are not preferred for primary atomization. Rather, many of the stud-
ies pertaining to the primary atomization process are Large Eddy Simulations (LES)
for the single phase region coupled with the VOF methods for the two-phase re-

gion. These approaches are similar to the under-resolved DNS method, hence known

44
as the Quasi-DNS approach. The Quasi-DNS approaches have appropriately cap-
tured the in-nozzle turbulence and their effects on the primary atomization pro-

cess [28], [29], [50]. However, the Quasi-DNS approach is still expensive for engi-
neering applications.

3.3 Eulerian methods

One of the famous Eulerian methods is the Eulerian Liquid Eulerian Gas (ELEG) [85]

or the two-fluid approach, where both phases are treated as the continuum through-
out the domain. The ELEG method considers additional source terms for the lift,
drag, and virtual mass forces in the momentum equation. The literature shows a good

agreement of the near field spray predictions by the ELEG approach with the exper-
imental observations, but the accuracy drops significantly in the sparse regions of
the spray. The ELEG approach is also computationally expensive compared to other

approaches. To avoid these limitations, Subramaniam and O’Rourke [215] suggested

to use the EE approach in the near-field region and the LE approach in the dispersed
spray region [148]. However, coupling the Eulerian liquid phase with the Lagrangian

liquid phase is perplexing, as it requires knowledge about the relationships between

the two representations. Recently, Pai et al. [148] sampled the Eulerian liquid phase
data by probability distribution functions to achieve the coupling between the Eule-

rian and Lagrangian liquid phases. They have also provided consistency relationships
which must be satisfied for a proper EE-LE transition.
Another well-known Eulerian approach, the Σ − Y approach, models the primary

atomization using the diffused interface method. The Σ−Y model was first postulated
by Vallet and Borghi [225]. It solves a density-averaged flow field while incorporating
models for interfacial effects. The Σ−Y model revolves around four basic assumptions

as mentioned below.

45
 1. At very high Reynolds number, as seen in automotive sprays, large-scale features
are independent of the surface tension and viscous forces. These forces are

predominant at smaller length scales, where a large velocity gradient occurs or

the curvature of the interface is large.

2. The mean velocity of the RANS turbulence closure model (k − ϵ) can be used
to study the velocity field of the flow.

3. The dispersion of the liquid phase in the gaseous phase is mixing dominated and

can be computed by solving a transport equation for the liquid mass fraction.
The closure of the liquid mass fraction transport equation can be modeled by
the turbulent diffusion liquid flux term.

4. The average droplet size for the dispersed liquid phase can be predicted by

solving a transport equation for the mean surface area of the liquid/gas interface
per unit volume (Σ).

Based on the first two assumptions, the mean flow for both phases can be treated
as a mixture with variable density which shares a single velocity field. As a result,

the mass and momentum conservation equation for the mixture can be represented
in the following forms.

∂ ρ̄ ∂ ρ̄ũi
+ = 0 (3.5)
∂t ∂xi
′ ′
∂ ρ̄ũj ∂ ρ̄ũi ũj ∂ p̄ ∂ ρ̄ug
i uj
+ = − − (3.6)
∂t ∂xi ∂xj ∂xi

In Eqn. 3.5 and 3.6,˜indicates Favre/density averaging,¯indicates time averaging, and

′ ′ ′
indicates turbulent fluctuations. The Reynolds stress term (ug
i uj ) in the Eqn. 3.6 can

be closed using the Boussinesq eddy viscosity model following the second assumption.

46
 Using the third assumption, the transport equation for the mean liquid mass
¯
ρY
fraction (Ỹ = ) can be represented as seen in the Eqn. 3.7.
ρ̄

′
∂ ρ̄Ỹ ∂ ρ̄ũi Ỹ ∂ ρ̄ug
iY
+ =− (3.7)
∂t ∂xi ∂xi

′
The closure of the turbulent mixing term (ug
i Y ) in the Eqn. 3.7, is modeled by the

Fick’s diffusion law given by Eqn. 3.8, where, Sc is the turbulent Schmidt number.

′ µt ∂ Ỹ
ρ̄ug
iY = (3.8)
Sc ∂xi

Following the fourth assumption, a transport equation for Σ is solved and represented
in the Eq. 3.9.

!
∂ Σ̄ ∂ ũj Σ̄ ∂ ∂ Σ̄
+ = DΣ + (A + a)Σ̄ − Vs Σ̄2 , (3.9)
∂t ∂xj ∂xj ∂xj

The first term on the right hand side of the Eqn.3.9 represents the diffusion term,

where the second term represents the production of the interface. The interface
production term considers the contributions from the collisions of droplets and the
stretching of the interface due to the gradient of the mean velocity. The third term

on the right hand side of the equation models the destruction of the interface due to
coalescence. Since its conceptualization, the Σ − Y equation has been successfully ap-
plied to model the primary atomization in diesel spray and gasoline sprays. However,

there is still a lot of scope left to improve the model.

3.4 Inferences
• The LE approaches are popular among the industry users, as they offer solutions

with reasonable accuracy while incurring low computational cost for engineering
applications.

47
• The LE models struggle with accuracy in the dense spray region.

• It is cumbersome to initialize the LE model due to the lack of available experi-

mental data in the near nozzle region.

• The LE models also suffer from numerical stability problems and mesh depen-

dency issues.

• True DNS are costly in highly atomizing conditions due to the wide range of
scales involved in the process.

• Quasi-DNS methods can be used to validate the engineering level models in

regions where experimental results are unattainable.

• Quasi-DNS methods are computationally exorbitant, because of the stringent

requirements on the grid resolution.

• The Σ − Y model overcomes the restrictions posed by the quasi-DNS approach.

48
 CHAPTER 4

GOVERNING EQUATIONS

4.1 CFD Solver

This chapter includes a detailed description of the governing equations of the
in-house Eulerian CFD solver, HRMFoam. The solver is implemented using the li-

braries of the foam-extend [2] branch of the OpenFOAM CFD framework. HRMFoam
solver was initially developed by Schmidt et al. [189], later extended by Gopalakr-
ishnan et al. [72] and Nerooorkar et al. [136], [137], [138]. Several past publica-

tions [194], [19], [132], [213], and [44] have successfully validated the HRMFoam code
for different flow applications such as channel flow, condensing two-phase injector
flow, cavitating diesel flow, and flash-boiling GDI flow.

The HRMFoam solver implements a diffuse interface based approach to ensure

mixing between the individual phases, i.e., pure liquid phase, pure non-condensible
gas (NCG) phase and the vapor phase. To ensure mixing, the solver defines a mixture

density (ρ) which is formulated based on the density of fuel (ρf ), density of the non-
condensible gas (ρg ) and the NCG mass fraction (y), as shown in Eqn. 4.1. The fuel
density is further determined by accounting for the liquid fuel density (ρl ), vapor

density (ρv ) and the vapor mass fraction (x), as seen in Eqn. 4.2.

1 − y y −1
ρ= + (4.1)
ρf ρg
1 − x x −1
ρf = + (4.2)
ρl ρv

49
 From the above discussion, it is evident that the mixture density (ρ) is determined
by the fuel quality (x), NCG mass fraction, and the density of individual phases. Fur-

thermore, the densities of the phases depend upon the pressure and enthalpy values.
Hence, the mixture density can be denoted by ρ (x, y, p, h). The diffuse interface

approach assumes mechanical equilibrium between the phases which means all the
phases experience same velocity and pressure fields. The governing equations for
mass and momentum conservation are represented in Eqns. 4.3, and 4.4, respectively.

The transported variables in these equations are Favre averaged as the formulation
accounts for turbulence and variable density.

∂ρ
+ ∇ · (Φ) = 0 (4.3)
∂t
→
−
∂ρ U  →− →
−
+ ∇ · ΦU = −∇p + ∇→
−
τ (4.4)
∂t

→
−
In the preceding equations, Φ represents the face centered mass flux, and →
−
τ corre-
sponds to the stress tensor, which includes both the viscous and turbulent stresses.

To account for the mixing of ambient air and fuel, a transport equation for the mean
non-condensible mass fraction, y, is solved, as shown in Eqn. 4.5.

∂ρy 
′Y ′

+ ∇ · (Φy) = ∇ · ρug (4.5)
∂t
′Y ′ =
µt
ρug ∇y (4.6)
Scr

′ Y ′ , accounts for the mixing

In Eqn. 4.5, the turbulent diffusion gas flux term, ug
effect of the relative velocity between the two phases. This term is closed using Fick’s
law of diffusion [53], as seen in Eqn. 4.6, where µt and Scr are the turbulent viscosity

and realizable Schmidt number fields.

The HRMFoam solver uses the homogeneous relaxation model (HRM) to define
the amount of vapor generation. In HRM, “Homogeneous” refers to the assumption

50
of homogeneous mixing of the two phases within any computational cell. Whereas,
“Relaxation” refers to the process that local flows develop towards thermal equilib-

rium. This model was proposed by Billicki and Kestin [30] based on the modification
to the homogeneous equilibrium model to incorporate the one-dimensional variability

of the vapor generation and condensation rates. This model governs the local rate of
change of the dryness fraction (x), which tends towards its equilibrium value (x̄) over
a specific time scale (Θ), as shown in Eqn. 4.7. The local enthalpy and pressure are

inputs for a lookup table generated by the REFPROP database [108] that governs
the equilibrium quality.
Dx x̄ − x
= (4.7)
Dt Θ
Based upon Reocreux’s “Moby Dick” experiments [174], Downar-Zapolski et al. [58]

proposed the time scale for phase change to be a function of the fuel void fraction
(α) and a non-dimensional pressure (ψ). Expressions to calculate these quantities are
included in Eqs. 4.8 and 4.9.

xρl
α = (4.8)
ρv + x(ρl − ρv )
psat − p
ψ = (4.9)
pcrit − psat

In Eqn. 4.9, psat is the saturation pressure and pcrit is the critical pressure of the
fuel. Here, p refers to the local total pressure instead of partial pressure. They also

proposed the following correlations for the time scale for operating pressures exceeding
10 bar, Eqn. 4.10.

Θ = Θ0 α−0.54 ψ −1.76 (4.10)

Θ0 = 3.84 × 10−7 s (4.11)

51
 To ensure energy conservation, a transport equation for the enthalpy (h)is included
in the HRMFoam solver, refer Eqn. 4.12.

∂ρh Dp
+ ∇.(Φh) = + ∇.J + D (4.12)
∂t Dt

In Eqn. 4.12, J accounts for the diffusion of the enthalpy and D represents the
conversion of kinetic energy to sensible enthalpy due to viscosity and turbulence.

4.1.1 PISO algorithm

The coupling between the momentum and the mass balance equations is ensured
by a modified Pressure Implicit Splitting Operator(PISO) algorithm. The discretized

form of the momentum equation 4.4 in a quasi-linear formulation is mentioned in

Eqn. 4.13.
X
ap up + aN uN = −∇p (4.13)
N

In Eqn. 4.13, the subscript p is used to represent the current computational cell, while

N represents the neighboring cells. The term a accounts for the coefficient of velocity
(u) which appears in the linearized momentum equation due to the discretization of
the advection, diffusion and the source terms. Furthermore, the term accounting for

the neighboring cells is replaced by another operator (H (u)), Eqn. 4.14.

ap up = H(u) − ∇p (4.14)

To allow for the pressure to respond to compressibility, fluctuations of density due

to phase change, and density change from turbulent mixing with the NCG, HRM is
extended to the three dimensions(3-D). The 3-D formulation is achieved by connecting

52
the predicted phase change with the conservation of mass and momentum through a
simple chain rule, as shown in Eqn. 4.15

Dρ ∂ρ Dp ∂ρ Dh ∂ρ Dx ∂ρ Dy
= + + + (4.15)
Dt ∂p x,y,h Dt ∂h x,y,p Dt ∂x y,p,h Dt ∂y x,p,h Dt

An assumption of adiabatic condition leads to dh = 0. Using this assumption and

Dρ
the mass conservation formulation for the term Dt
, Eqn. 4.15 reduces to the following

form. The formulation ensures response of the pressure field to the divergence of the
velocity.
∂ρ Dp ∂ρ Dx ∂ρ Dy
−ρ∇.u = + + (4.16)
∂p x,y,h Dt ∂x y,p,h Dt ∂y x,p,h Dt

Furthermore, substituting Eqn. 4.14 results in Eqn. 4.17.

!
H(u) 1 ∂ρ Dp ∂ρ Dx ∂ρ Dy
−ρ∇. + ρ∇. ∇p = + + (4.17)
ap ap ∂p x,y,h Dt ∂x y,p,h Dt ∂y x,p,h Dt

Dy Dx
In Eqn. 4.17, the terms Dt
and Dt
are closed using the formulations presented in
Dp
Eqns. 4.5, and 4.7, respectively. Furthermore, the term Dt
in Eqn. 4.17 is expanded

resulting in Eqn. 4.18.

!
1 ∂ρ ∂p H(u) 1
+ ∇.(pu) − p∇.u + ∇. − ∇2 p+
ρ ∂p x,y,h ∂t ap ap
1 ∂ρ Dx 1 ∂ρ Dy
+ = 0 (4.18)
ρ ∂x y,p,h Dt ρ ∂y x,p,h Dt

The partial derivatives of ρ with respect to p, x and y are further written in the

following forms.
!
∂ρ 1 − y (1 − y)(1 − x) y
= ρ2 + + , (4.19)
∂p x,y,h ρv a2v ρl a2l ρv a2g
!
∂ρ 1 1
= (y − 1)ρ2 − , (4.20)
∂x y,p,h ρv ρl
!
∂ρ 1 1 − x x
= −ρ2 − − (4.21)
∂y x,p,h ρg ρl ρl

53
 In Eqn. 4.19, av , al and ag are the speed of sound in fuel vapor, liquid fuel, and
ambient gas, respectively. Based upon Brennen’s analysis [37], the partial derivative
∂ρ
is assumed to be related to the enthalpic sonic speed of the two-phase mixture
∂p x,y,h
and is given by Eqn. 4.22
∂ρ 1
= 2 (4.22)
∂p x,y,h a
However, the expression for compressibility in Eqn. 4.19 is found to be numerically un-
stable due to the very abrupt changes in two-phase regions. Consequently, a smoother
mass-weighted expression for compressibility is adopted in the present work and is

respresented in Eqn. 4.23:

!
∂ρ (1 − y)x (1 − y)(1 − x) y
=ρ + + (4.23)
∂p x,y,h ρv a2v ρl a2l p

After implicitly solving the pressure by the help of Eqn. 4.18, the flow velocity u is
updated using Eqn. 4.14.

4.2 Needle Closure Model

The challenge of modeling sealing is how to effectively cut off the communication

between two parts of the domain at the needle seat. While the obvious approach is to
topologically sever the domain, this approach would encounter a few disadvantages.
The first is the abruptness of such an operation would create massive disturbances in

the flow, triggering strong waves that could result in spurious waterhammer effects.
While such effects are expected in rapid valve closure, an instantaneous severance
would be unnaturally violent. Secondly, it is unclear what a linear solver that is

designed for solving a single domain problem would do when confronted with two
disjoint problems stored in a single linear system.
Instead, we propose to model what cannot be directly resolved in the simulation.

When the needle is nearly closed to the seat, the drag force produced by the walls

54
begins to dominate the flow. The closure model will employ a gradual increase in drag
that halts the flow through the narrow gap. To avoid degenerate cells, this sealing

model will active at low, but finite, needle valve lift. To model this event, a volumetric
source term f¯ is added to the right hand side of the momentum equation, Eqn. 4.4.

The modified momentum equation is mentioned below.

∂ρu
+ ∇.(Φu) = −∇p + ∇τ̄¯ + f¯ (4.24)
∂t
Sf
f¯ = ρu (4.25)
Sd

Eqn. 4.25 describes the artificial drag force, with Sf and Sd declared as the sealing
factor and drag constant, respectively. The drag constant represents how vigorously

the imposed force opposes the motion of the fluid. A very small value of the order
10−8 s or 10−9 s is generally assigned to Sd . The value must be sufficiently high to
bring the velocity of the fluid down to near-zero values in the sealing region.

Furthermore, the drag force is applied gradually to avoid an unrealistically violent

closure or opening. This is accomplished by modeling Sf as a continuous function
of time and seal constant S∞ , as shown in Eqn. 4.26. The discrete S∞ is unity only

when the needle-seat sealing is being applied.

 t

Sf = S∞ 1 − e−( τ ) (4.26)

In Eqn. 4.26, τ represents the time relaxation factor. Changing τ changes the rate

at which the sealing is applied or removed. The expression in Eqn. 4.25 creates an
exponential decay of velocity dependent on the strength of the drag constant, Sd . The
drag is turned on or off at a rate governed by τ . This gradual switching is controlled

by Sf , which represents a relaxation towards the state indicated by S∞ . In the seat

region, which is designated by an user input of a bounding box, the value of S∞ is

55
switched between zero and unity to represent opening or closing, respectively. For
rest of the domain, the value is always zero so that the artificial drag force is zero at

all times.
Because of the form of Eqn. 4.25, when treated implicitly, the source term appears

on the diagonal of the matrix in the linear system of equations for velocity, the ap
term of Eqn. 4.13. This formulation intrinsically benefits stability by increasing the
diagonal dominance of the matrix. If the drag term and temporal term dominate the

equation, then the decay of velocity should, at least for momentum, unconditionally
stable. Considering the coupling with the other transport equations, one cannot a
priori know the stability of the system as a whole.

4.3 Primary atomization model

Rachakonda et al. [166], [165] expanded the capabilities of the in-house HRM

solver to model the primary atomization by including a modified vesrion of the Σ − Y

solver discussed in Sec. 3.3. The formulation includes a transport equation for the
variable Ω, mass based liquid-gas interface density, defined in the Equation 4.27.

The formulation is numerically more stable compared to other formulations proposed

by several researchers [69] , [225], [10], [51], and [42] due to its volume conserving
nature. However, unphysical droplet diameters are predicted by the formulation

in the primary atomization regime. Hence, the current research uses an improved
formulation, which will be discussed in the latter part of this section.

Σ
Ω= (4.27)
ρ

In Eqn. 4.27, Σ is the interfacial area density which has a unit of m−1 . Further,

Σ and the liquid volume fraction (Y ) can be used to predict the Sauter mean diame-
ter(SMD) in the primary atomization regime, as seen in Eqn. 4.28 . The formulation

56
applies the assumption that the liquid phase in a computational cell can be modeled
with spherical droplets to determine SMD in the transition region from the dense

spray to sparse spray region, as shown in Eqn. 4.28.

6Y
SM D = (4.28)
Σ

Recently, several revisions of the initial Σ − Y formulation [225] have been pro-
posed [26], [104], [63], [79] to expand the modeling paradigm into the secondary
atomization regime. The notable formulation proposed by Lebas et al. [104] is men-

tioned in Eqn. 4.29. The equation includes an indicator function (Ψ) that acts as
 
a switch between the terms for secondary atomization, i.e., Ω̇breakup , Ω̇coal/coll and
 
the primary atomization terms Ω̇mixture , Ω̇stress . The term Ωvap accounts for the

production/destruction of the interface due to vaporization.

 
∂ρΩ µt  
+ ∇. (ρuΩ) = ∇. ∇Ω + ρΨ Ω̇mixture + Ω̇stress
∂t ScΣ
  (4.29)
+ρΩ̇vap + ρ (1 − Ψ) Ω̇breakup + Ω̇coal/coll

However, including models for secondary atomization in an Eulerian primary at-

omization model is computationally expensive. Moreover, the secondary break up

processes can be modeled by droplet oriented approaches. Hence, a simplified ver-

sion of the Σ − Y formulation, proposed by Lebas et al. [104] is used in the current
research. The formulation only accounts for the primary atomization regime in the

GDI applications, as seen in Eqn. 4.30.

 
∂ρΩ µt  
+ ∇. (ρuΩ) = ∇. ∇Ω + ρ Ω̇mixture + Ω̇stress (4.30)
∂t ScΣ

In Eqn. 4.30, Ω̇mixture ensures production of the liquid/gas interface due to mixing
of phases. This term further enables the initialization of the mean interface surface

57
density. Nevertheless, it is of importance only in the vicinity of the nozzles. The
mixture term is modeled based on the assumption that the initial length scale of the

produced liquid blob should be proportional to the turbulent length scale (lt ) [104],
[63]. The expression for the mixture term is included in Eqn. 4.31, and it depends

upon the gradient of liquid mass fraction (Ỹ ). Liquid mass fraction is further deter-
mined by Eqn. 4.32.

12νt ρ2    
Ω̇mixture = ∇Ỹ . ∇Ỹ (4.31)
ScΣ ρl ρg lt

     
ml ml mf mv mg
Ỹ = = = 1− 1− = (1 − x) (1 − y) (4.32)
mt mf mt mf mt

In Eqn. 4.30, Ω̇stress represents the source term that accounts for the produc-
tion/destruction of interface due to the turbulent stretching. It also considers the

effects of collision and coalescence in the dense spray region [104]. The stress term
attempts to attain an equilibrium interfacial surface density Ω∗eq in a duration pro-

portional to the turbulent time scale (τt ). Further, the equilibrium interfacial surface

density corresponds to the quantity of surface that would be obtained by keeping

the liquid volume fraction Ȳ , and the mixture turbulent kinetic energy (k) con-
stant at a critical Weber number (W ecr ). We hypothesize that the droplets in the

primary atomization regime attain the largest expected radius based on the critical
Weber numbers in the range of 1 − 6, as observed by Pilch and Erdman for a wide
range of Ohnesorge number [160]. The terms, Ω̇stress and Ω̇∗eq are calculated based on

Eqns. 4.34, and 4.35. In these equations, k and σ correspond to the turbulent kinetic
energy and the surface tension force.

58
  
Ω Ω
Ω̇stress = 1− ∗ (4.33)
τt Ωeq
Ȳ k
Ω∗eq = (4.34)
σW ecr
W ecr = 1 : 6 (4.35)

These auxiliary transport equations for the interfacial surface density do not feed-
back into the momentum or energy equations. Hence, to reduce computational cost,

these equations can be turned off when modeling primary atomization is not the
prime goal of the computation. The future chapter about the plume-based coupling
approach uses this Σ − Y formulation to model primary atomization proicess.

59
 CHAPTER 5

CASE STUDIES: MULTIPLE INJECTIONS IN GDI

SYSTEMS

The early and late portions of transient fuel injection have proven to be a rich
area of research, especially since the end of injection can create a disproportionate

amount of emissions in direct injection internal combustion engines. A perennial

challenge in simulating the internal flow of fuel injectors is the valve opening and
closure event. In a typical adaptive mesh CFD simulation, the small gap between

the needle valve and the seat must be resolved to very small cells, resulting in costly
computations. Capturing a complete closure usually involves a topological change
in the computational domain. In this chapter, we have applied the needle closure

algorithm, discussed in the previous section 4.2, to a gasoline direct injection system
operating under cavitating conditions. The closure algorithm is designed to avoid
spurious water-hammer effects.

Traditionally, multiple fuel injection pulses have been employed in diesel engines
to improve emissions and heat release rate [134], and [221]. Similarly, in GDI ap-
plications, multiple injection strategies have shown benefits in dramatically improv-

ing combustion stability and spark authority when compared to a single injection
event [180]. The current chapter demonstrates the first of a kind single simulation of
a multiple gasoline injection event known to the authors.

60
SIMULATION SETUP
Geometric details

The ECN Spray G nozzle consists of a multihole nozzle with eight counterbore

orifices. Its geometric details are listed in Table 5.1. For the current simulations,
the specified nominal values are used to build the computational domain (named
Generation 1 geometry), shown in Figure 5.1 and Figure 5.2. The spherical shape of

the needle is also built from the nominal manufacturing dimensions. A semispherical
volume of 9 mm diameter is added as the discharge plenum to simulate the first
millimeters of the spray development.

Table 5.1: ECN Spray G nominal nozzle geometry parameters.

Nozzle type valve covered orifice (VCO)

Bend angle 0°
Number of holes 8
Orifice shape circular
Hole shape straight
Nozzle shape step hole
L/D ratio 1.4
Orifice diameter 0.165 mm
Orifice length 0.16 − 0.18 mm
Step diameter 0.388 mm
Orifice drill angle 37° relative to nozzle axis
Full outer spray angle 80°

Mesh

A primarily hexahedral mesh for the nominal geometry with a cell count of

1.44 million has been created using the GridPro meshing tool. The grid spacing is
roughly 7µm in the sac region, as opposed to the spacing of around 10µm in the nozzle
region. To avoid higher cell counts, an anisotropic refinement in the narrow region

between the needle and the nozzle was performed. The mesh in the seat region can
be seen in Fig. 5.3. The mesh is created at an initial displacement of 5µm. A prior
work investigated mesh independence [132] and in the present work the cell count is

61
 Figure 5.2: ECN spray G nomi-
Figure 5.1: ECN spray G computational
nal geometry showing distribution
domain with the outlet plenum of 9mm di-
of the 8 nozzle holes and 5 dimples.
ameter.

comparable to or higher than that of other diffuse interface simulations reported by

Moulai et al. [132] and Saha et al. [183],[185]. In the study by Moulai et al.[132],
these spacings were 10µm and 40µm, respectively. For Saha et al., they were 15µm

and 30µm [183], and [185].

Needle motion

Needle motion in the simulation is controlled by a profile of a typical double

injection for these conditions. This profile is generated using the experimental data.
The profile is adjusted to account for the initial lift of the mesh, as shown in Fig. 5.4.
The layer addition and removal algorithm is used to accomplish the prescribed mesh

motion.
In this simulation, the sealing algorithm is activated whenever the needle lift is

less then 1µm from its initial position. The mesh is built with the needle lifted 5µm
from the seat, so the algorithm is effectively disengaged at 6µm of lift. When engaged,
the algorithm brings the fluid within the needle-seat region to rest by applying the

large artificial drag force previously described.

62
Figure 5.3: Cut plane of the CFD mesh in the needle seat area prior to the needle lift

Table 5.2: Non-flashing (spray G) operating condition defined by the ECN

Spray G operating condition

Fuel Iso-octane
Upstream pressure 20 MPa
Fuel temperature 90 0 C
Ambient temperature 300 0 C
Ambient density 3.5 kg/m3
Back pressure 600 kPa

Initial and boundary condition

The computational study is based upon the Engine Combustion Network (ECN)

Spray G operating conditions, mentioned in Table 5.2.Under these conditions, the

elevated back pressure leads to a non-flashing flow.
Zero-gradient boundary conditions are used for velocity at both the inlet and the

outlet. For the pressure at the inlet, a total pressure boundary condition is chosen. At
the outlet, a transonic total pressure boundary condition is applied. This boundary

63
 Figure 5.4: Needle lift profile used in the CFD simulations

condition switches between the zero-gradient and total pressure boundary conditions

based upon the Mach number at the exit.

During the initial condition, the needle is closed and the fluid is stationary, so all
pressure drop should occur at the needle seat. To avoid the formation of spurious

pressure waves in the flow domain during the start of the simulation, a hyperbolic
tangent drop of pressure across the needle seat area is applied, as shown in the
Figure 5.6. The interface between the fuel and the non-condensible gas is placed at

the lower bound of the sealing region, as seen in Fig. 5.7. This leaves the sac initially
filled with the non-condensible gas. The initial values of the sealing constant (S∞ )
and the sealing factor (Sf ) are both unity in the seat region and zero elsewhere. The

sealing field at the start of the simulation can be seen in the Fig. 5.5.

Turbulence closure and model assumptions

To model turbulence, the SST k − ω RANS model in conjunction with the

log-layer wall function is used. This model is well known for resolving separating
flows with adverse pressure gradients. In the present simulation, iso-octane is used

as the fuel. To account for the compressibility effects, the current work assumes the

64
Figure 5.5: Initialization of the sealing field at the start of simulation

Figure 5.6: Initialization of the pressure field at the start of simulation

65
Figure 5.7: Initialization of the non-condensible gas field at the start of simulation

liquid fuel to be compressible. Furthermore, the walls are assumed to be adiabatic,

which is consistent with the assumptions by Gavaises et al. [70] in the case of diesel
injection systems. For further simplification of the posed problem, no heat transfer

is considered between the non-condensible gas and the liquid fuel. To model these
phenomena, a careful evaluation of the interfacial heat transfer coefficient over the
complete injection cycle is required, which is beyond the scope of the current research.

5.1 Results and discussion

5.1.1 Sealing

As discussed in the previous chapter, the sealing algorithm switches the value
of S∞ to unity in the seat region whenever the gap between the needle wall and the

injector wall is less than 6 µm. The value of Sf plummets to zero everywhere, shortly
after the needle lift exceeds the needle displacement limit, i.e., 1µm. As intended,
once the needle is closed after the main injection event, the value of S∞ switches to

66
unity and the value of Sf quickly follows. The right-bottom image in Fig. 5.8 shows
the value of Sf during the dwell period between the individual injection cycles.

The effect of the sealing field on the flow domain is as intended. In Fig. 5.8,
velocity contours in the sealing region are presented for the different instances of
the dwell phase. Careful inspection shows that the velocity in the sealing region

diminishes to negligible values while the downstream jets persist for a short period.

Figure 5.8: Sealing Region and corresponding velocity contour at different time
stamps during needle closure event.

5.1.2 Rate of Injection (ROI)

The Rate of Injection (ROI) for the multiple injection event is measured at the

nozzle exit, as shown in Fig. 5.9 . The ROI is found to follow the needle displacement
curve. However, an interesting phenomenon during the second injection is observed.
Although the maximum needle displacement during the second injection is half of

that of the first injection cycle, the maximum predicted ROI for the injection cycles
are found to be of similar order.

67
Figure 5.9: Rate of Injection curve for a multiple injection event

68
5.1.3 Cavitation

The simulation captures cavitation during the low lift period of the needle’s
transient motion(during both needle opening and closing phases). The cavitation
bubbles created at the end of both the first and second injection cycles are found to

linger until the start of the dwell period. This cavitation can be a major contributor
to the injector failure in the needle-seat region [186]. Fig. 5.10 presents the fuel
vapor volume fraction and the corresponding local low pressure regions during the

needle opening phase indicating the cavitation event. Similarly, cavitation is also
demonstrated during the needle closure event for both injection shots, as seen in
Fig. 5.11. After completion of the 2nd injection, fuel vapor has been observed in the

sac in addition to the seat region. A detailed analysis of the seat and sac bubbles
during the needle closure events will be included in the future chapters.

Figure 5.10: Vapor volume fraction and low pressure region showing cavitation during
needle opening phase of the 1st injection cycle.

69
Figure 5.11: Vapor volume fraction showing cavitation during needle closing phase of
the both 1st and second injection cycle.

5.1.4 Residual fuel

Just after needle closure, the velocity at the sealing region decreases to zero
and the pressure inside the sac is reduced. However, residual fuel exists in the nozzle

holes. Although the inertia of this residual fuel causes it to exit the hole, it decelerates
because of the low pressure gradient present. As a result, the spray plume coming out
of the exit plane starts to oscillate leading to eventual jet disintegration. This slow

moving liquid fuel jet is unlikely to atomize well. Hence, it is most likely to contribute
to the hydrocarbon emissions during combustion. Figure 5.12 indicates possible liquid
fuel dribbling events during the EOI of both the first and second injection pulses.

5.2 Conclusions
A sealing algorithm has been developed and successfully applied to a non-flashing

Spray G injection with a double injection cycle. The primary results are summarized
below:

70
Figure 5.12: Liquid fuel mass fraction and velocity contour at 1.151 ms (EOI of
the first injection),2.171 ms (EOI of the second injection) showing dribble and jet
disintegration

• Results demonstrate the ability of the sealing algorithm to completely block

the needle-seat area at low needle lifts, resulting in zero flow velocity within the
sealing region.

• At low needle lift and during both SOI and EOI, fuel vapor generation is ob-
served as the liquid fuel passes through the narrow region between the needle

and seat.

• Towards the end of the second injection, fuel vapors are captured within the
sac region.

• After the end of injection, the possible liquid fuel dribble leading to eventual
jet disintegration has been demonstrated.

71
 CHAPTER 6

EFFECT OF FLASH-BOILING ON THE NEEDLE

CLOSURE EVENT

Gasoline Direct Injection (GDI) engines offer the potential of dense and effi-

cient power generation, but incur difficulties with emissions. GDI engines have made
significant progress since their inception, but more study is required to understand
the impact of the injector’s internal flow physics on the overall durability and perfor-

mance of the engine. The Start of Injection(SOI) and End of Injection(EOI) are two
crucial phases of the needle’s transient motion. Both of these injection phases gen-
erate pressure waves and flow turbulence attributable to the opening and closing of

the valve. Cavitation, a critical concern for the health of the needle seat, is often en-
countered at both the SOI and EOI [186]. Liquid fuel dribble, a major contributor to
unburned hydrocarbon emissions, is encountered at the EOI phase. Hot gas ingestion

occurs in the counterbore and nozzle regions during the needle’s dwell time between
injections. This ingestion of gas can have adverse effects on the next cycle’s SOI, as
it both delays the liquid injection and leads to more liquid vaporization[47, 159].

The computational research on SOI and EOI phases for the GDI injector is lim-
ited comapred to the diesel injection system. Ishii et al.[83] studied the flow inside
a multihole GDI injector with transient needle motion under non-flashing conditions

using a VOF model with the particle/grid hybrid method. They captured the liq-
uid fuel dribble and film formation at the EOI, but also showed a leakage velocity
through the seat area during this injection phase. Their model was computationally

expensive, as very fine grids were required to properly simulate the particles. Cav-
itation and turbulence, which are quite common during the EOI for GDI injectors,

72
were not modeled. The EOI event of a GDI injector has also been simulated by [133]
while using a VOF approach. They implemented the 1-D Sauer and Schnerr model

in conjunction with the Continuum Surface Model (CSM) to capture the sac and
tip conditions post-needle closure. However, the simulated ambient conditions were

limited to non-flashing conditions, and their simulations only included the needle lift
curve for the needle closing phase due to the computationally expensive nature of the
VOF approach. Recently, Hwang et al. [82] employed the diffuse interface approach

to predict the boil-off characteristics of a single axial-hole transparent nozzle under

non-flashing, mildly flashing, and flare-flashing conditions. They were able to get
good agreement for the boil-off time under non-flashing and flare-flashing conditions

with the experimental observations. Yet, they overpredicted the boil-off time for
mildly flashing conditions. The needle valve opening and closing represents a general
challenge in computational fluid dynamics. The valve closure represents an abrupt

change in domain topology, and the flash boiling process in highly turbulent flow
involves sudden vapor generation and collapse that is difficult to model.
In the present chapter, the needle closure model in tandem with the homogeneous

relaxation model has been applied to a single injection cycle of GDI injector operating
in different ambient conditions. This allows for a novel study of the EOI for the
multihole GDI injector under flash boiling conditions. The details about the model

implementation are discussed in the previous chapters. Moreover, the current study
includes the different sac conditions present, post-needle closure, for different ambient
conditions, and the mechanism driving these sac conditions. The studies about the

tip-wetting phenomena for both flash-boiling and non-flashing conditions are also
included. The study indicates dribbling in the non-flashing condition. Additionally,
the CFD predictions are compared with the existing experimental observations.

73
6.1 Simulation setup
The current simulations use the same geometry and mesh described in Ch. 5.
However, the needle displacement curve is different, as shown in Fig. 6.1. The needle
displacement curve is adjusted to account for the initial mesh lift. The maximum lift

for the single injection event is 55µm. The simulations are performed for ECN pre-
scribed non-flashing(G) and flash-boiling (G2) conditions. The operating conditions
are mentioned in the Table 6.1. Both conditions have the same injection pressure

and fuel temperature (20M P a and 900 C respectively), but the ambient density and
back pressure are different. The back pressure and ambient density are 600kP a and
3.5kg/m3 for the non-flashing case, whereas they are 53kP a and 0.5kg/m3 for the

flashing condition.

Figure 6.1: Ensemble-averaged needle displacement profile generated by X-ray mea-

surement [59]

In the current simulations, the sac is initialized with non-condensible gas at the
start of injection. The sealing algorithm is activated whenever the needle is displaced

74
Table 6.1: Operating condition for G (Non-Flashing) and G2 (Flashing) condition.

Parameter G-Standard G2-Flash boiling

Fuel Iso-Octane Iso-Octane
Fuel temperature 90 °C 90 °C
Ambient temperature 573 K 333 K
Fuel pressure 20 MPa 20 MPa
Ambient density 3.5 kg m−3 0.5 kg m−3
Ambient Pressure 600 kPa. 53 kPa

less then 1µm from its initial position. The lalplacian mesh motion library of FOAM

extend is used to lift the needle according to the experimentally measured needle
displacement curve. The same model assumptions discussed in the previous chapter 5
are used for the current simulations.

6.2 Results and discussion

6.2.1 Rate of Injection

For validation, Rate of Injection (ROI) data from the flashing and non-flashing
cases are compared with the experimentally measured ROI profile provided by the

ECN. These measurements were obtained by Duke et al. [59] through the Bosch long
tube method. The simulated maximum rate of injection and the ROI profile both
match well with the experimental observations, as seen Fig. 6.2. The percentage dif-

ference between the experimental observations and the computational predictions for
the maximum ROI and the ROI during the steady injection phase are 3% and 4%,
respectively. In the current simulation, the needle is assumed to be closed until it

reaches a displacement of 6µm from the seat. This assumption leads to a difference
in the observed needle opening duration between the experiment and the simula-
tions, i.e., the simulations assumed the needle opening duration to be 735µs versus

the observed needle opening duration of 760µs. Consequently, the simulated ROI
is observed to diminish to a zero value faster compared to the experimental ROI.

75
The differences in the injection duration can also be attributed to the fact that the
simulations assumed the sac to be empty at the start of the simulation, whereas the

long tube method starts with a liquid-filled tube. The simulations predict a total
injected mass of 10.01 mg for both the flashing and non-flashing conditions. The

Spray G target condition specifies a total injected mass of 10 mg, and experiments
have observed 10.1 mg being injected. Furthermore, the ROI predicted by the non-
flashing condition is similar to the ROI predicted by the flashing condition. It implies

that the different downstream pressures have no effects on the mass flow rate. This
corroborates the fact that these nozzles are operating in choking conditions [132].

Figure 6.2: Simulated and experimental ROI vs time. Experimental ROI obtained
through the long-tube method [59]and simulated ROI taken at the nozzle exit

76
6.2.2 EOI sac condition

Post injection sac and nozzle conditions can be detrimental to the emissions
behavior of GDI engines. Traces of liquid or vapor fuel in the sac will alter the
behavior of a second injection compared to the empty sac condition, and can also

affect the hydraulic injection duration.

Figure 6.4: Vapor vol- Figure 6.5: NCG vol-

Figure 6.3: Liquid vol- ume fraction before nee- ume fraction before
ume fraction just before dle closing, i.e., t = needle closing i.e. t =
needle closing, i.e., t = 1.038 ms for non-flashing 1.038 ms for non-flashing
1.038 ms for non-flashing conditions showing cavi- conditions showing gas
conditions. tation due to the nozzle entrainment in the
inlet corner. counter bore.

For reference, we define a pressure ratio as the ratio of local static pressure, i.e,

the pressure at the cell center to the saturation pressure of the fuel at its injected
temperature. For the non-flashing case, the ambient pressure ratio is around 7.7.
Figures 6.3, 6.4, and 6.5 captures the gas ingestion during the phase of injection

where the needle is open. As described by Fig. 6.4, 6.5, 6.12 a local low pressure
region is also created in the counterbore and at the nozzle inlet corner, resulting in
vapor generation. When the needle begins to close, the low pressure waves created

in the seat region and at the nozzle exit start to move into the sac, resulting in
low pressure spots inside the sac. This water hammer type transient phenomenon is
clearly visible in a simple pressure contour, as shown in figures 6.12, 6.13, and 6.14.

These low pressure spots dip below the prevailing saturation pressure, resulting in
the cavitation seen in figures 6.6 and 6.7, but these bubbles are very short lived as

77
the back pressure is greater than the saturation pressure. Post needle closure, i.e.,
during the dwell phase, the sac reaches the ambient pressure and entrains ambient

gas. This leads to a mixture of liquid fuel and ambient gas for the initial condition
of the next injection (see Fig. 6.9, 6.10, 6.11). Similar observations have been made

by Moon et al. [131] experimentally. However, they injected n-heptane from a 3-hole
GDI injector under non-flashing conditions.

Figure 6.6: Liquid vol- Figure 6.7: Vapor vol- Figure 6.8: NCG volume
ume fraction just after ume fraction just after fraction just after needle
needle closing i.e. t = needle closing i.e. t = closing i.e. t = 1.045 ms
1.045 ms for non-flashing 1.045 ms for non-flashing for non-flashing condi-
conditions showing cavi- conditions showing cavi- tions showing no pres-
tation. tation. ence of gas in the sac.

Figure 6.10: Vapor vol-

Figure 6.11: NCG vol-
Figure 6.9: Liquid vol- ume fraction long after
ume fraction long after
ume fraction long after needle closing i.e. t
needle closing i.e. t =
needle closing i.e. t = = 2.000 ms for non-
2.000 ms for non-flashing
2.000 ms for non-flashing flashing condition show-
condition showing pres-
conditions. ing no traces of vapor in-
ence of gas inside the sac.
side the sac.

Drastically different sac and nozzle conditions are observed for the flash boiling
case, where the ambient pressure ratio is 0.68. The simulation captures low pres-

78
 Figure 6.14: Pressure ra-
Figure 6.12: Pressure
Figure 6.13: Pressure ra- tio after needle closing
ratio just before needle
tio just after needle clos- i.e. t = 1.045 ms for non-
closing i.e. t = 1.038 ms
ing i.e. t = 1.042 ms for flashing conditions show-
for non-flashing condi-
non-flashing conditions. ing receding low pressure
tions.
waves.

sure regions inside the counter bore region. These regions further increase vapor
generation and create a vapor core in the jet emanating from the nozzle, as seen in
Fig. 6.15, 6.16, 6.18, 6.19, 6.24, and 6.25. Post needle closure, a local low pressure

region is observed inside the sac leading to the vapor generation seen in figures 6.19
and 6.25. Further investigation shows similar reflecting pressure waves compared to
the non-flash boiling case. However, the vapors generated during the flash boiling

case are much longer lived than the non-flashing case due to lower ambient pressure
waves. The pressure in the sac region reaches the ambient pressure gradually during
the needle’s closed dwell phase as before. Nevertheless, the residual fuel inside the

sac starts boiling off slowly as the ambient pressure is now below the saturation pres-
sure. The resultant expanding fuel vapor inside the sac and nozzle does not allow the
ambient gas to be entrained into the nozzle and sac. The post-injection boiling sac

phenomena can be viewed in figures 6.21, 6.22, 6.23, and 6.26, taken at EOI. A sim-
ilar boiling sac phenomenon was witnessed by [82] in a single axial-hole transparent
nozzle with n-pentane as the test fluid under flashing conditions.

For further comparison between the flash boiling and non-flashing cases, the total
volume of liquid and vapor fuel and non-condensible gas present in the sac and nozzle
are calculated separately for each time step using the empty volume of the sac for

79
Figure 6.15: Liquid vol- Figure 6.16: Vapor vol- Figure 6.17: NCG vol-
ume fraction before nee- ume fraction before nee- ume fraction before nee-
dle closing i.e. t = dle closing i.e. t = dle closing i.e. t =
1.038 ms for flash-boiling 1.038 ms for flash-boiling 1.038 ms for flash-boiling
conditions. conditions. conditions.

Figure 6.18: Liquid vol- Figure 6.19: Vapor vol- Figure 6.20: NCG vol-
ume fraction just after ume fraction long after ume fraction just after
needle closing i.e. t = needle closing i.e. t = needle closing i.e. t =
1.046 ms for flash-boiling 1.046 ms for flash-boiling 1.046 ms for flash-boiling
conditions. conditions. conditions.

Figure 6.21: Liquid vol- Figure 6.22: Vapor vol- Figure 6.23: NCG vol-
ume fraction long after ume fraction long after ume fraction long after
needle closing i.e. t = needle closing i.e. t = needle closing i.e. t =
2.000 ms for flash-boiling 2.000 ms for flash-boiling 2.000 ms for flash-boiling
conditions. conditions. conditions.

80
Figure 6.24: Pressure
Figure 6.25: Pressure ra- Figure 6.26: Pressure ra-
ratio just before needle
tio just after needle clos- tio long after needle clos-
closing, i.e., t = 1.038 ms
ing, i.e., t = 1.046 ms for ing, i.e., t = 2.000 ms for
for flash-boiling condi-
flash-boiling conditions. flash-boiling conditions.
tions.

reference. Fig. 6.27 and 6.28 represent the liquid volume fraction, vapor volume
fraction, and NCG volume fraction for the complete injection cycle under the non-

flashing and flashing conditions, respectively. From these figures, it can be concluded
that the sac and nozzle reach a non-flashing postinjection condition where almost
50% of the space is occupied by the liquid fuel, while the other 50% is occupied

by the non-condensible gas. Alternatively, for the flashing condition, 90% of the
liquid fuel present in the sac and nozzle vaporizes to fill the sac and keep the non-
condensible gas away from the nozzles. The total volume of the sac and nozzle for

the 8-hole GDI injector is 0.09 mm3 . Hwang et al. [82] observed a boil-off time of 1
ms for the sac volume of 0.07 mm3 with the pressure ratio of 0.69. From Fig. 6.28,
it can be observed that the sac takes roughly 0.96 ms to get filled with 90% fuel

vapor under flash-boiling condition. The boiling time is of similar order compared to
the experimentally observed sac boil-off time. Even so, it should be noted that the
experiments were performed with a different fuel and nozzle configuration.

6.2.3 The λ2 -criterion

The ensuing spray from the multihole GDI injector is highly turbulent. Thus,
there exist many vortical coherent structures in the internal nozzle region. Fur-

thermore, these vortices interact among themselves while either extending between

81
Figure 6.27: Liquid Volume Fraction, Va- Figure 6.28: Liquid Volume Fraction, Va-
por Volume Fraction and NCG Volume por Volume Fraction and NCG Volume
Fraction for a non flashing condition for Fraction for a flash boiling condition for
the sac and nozzle region the sac and nozzle region

neighboring holes, or just ending in a single nozzle hole. In this section, we intend
to capture the vortices which are present during the needle closing phase and impact

the sac boil-off process. The λ2 -criterion defines a vortex as a connected fluid region
with two eigen values of the velocity gradient tensor. It looks for a pressure mini-
mum but removes the effects of unsteady straining and viscosity by discarding them.

The definitions for the Q-criterion and the λ2 -criterion tend to be similar, however,
their vortex predictions differ significantly for the regions where vortex stretching or
compression are dominant. The λ2 -criterion is the median of the three eigen vectors

(λ1 , λ2 , λ3 ) of S 2 + Ω2 . Here, S, and Ω represent the symmetric and antisymmetric

components of the velocity gradient. Nonetheless, the Q-criterion is evaluated as,
Q = − 12 (λ1 + λ2 + λ3 ). For vortices in the variable density flows, Jeong et al. [88]

observed better identification with the λ2 -criterion compared to the Q-criterion. Ac-
cording to them, the vortices can be visualized as iso-surfaces for different negative
values of λ2 .

Fig. 6.29, and 6.30 include the iso-surfaces of the λ2 criteria colored by the pressure
ratio and the vapor volume fraction for the non-flashing condition at t = 1.042ms,

82
respectively. The time stamp corresponds to the needle closure phase. The figures
capture vortices extending from the sac region into the nozzle holes. Furthermore, the

pressure in the vortex core falls below the saturation pressure, resulting in the forma-
tion of vapor bubbles inside the sac. This observation of sac cavitation corroborates

our previous observations in Fig. 6.7.

Figure 6.29: Iso-surface for λ2 = Figure 6.30: Iso-surface for λ2 =

−8e + 11 colored by the pressure ra- −8e + 11 colored by the VVF for
tio for the non-flashing condition at the non-flashing condition at t =
t = 1.042ms, i.e., after the needle 1.042ms, i.e., after the needle clo-
closure. sure.

Similarly, vortices extending between the sac and the nozzle holes are also iden-
tified for the flash-boiling condition through the λ2 criteria, as shown in Fig. 6.31
and 6.32. Nonetheless, the strengths of vortices are found to be lower compared to

the non-flashing condition. Subsequently, these vortices tend to terminate early in

the flash-boiling condition. These contrasting observations can be explained by the

hypothesis that the flash-boiling condition has a lower angular velocity due to a higher
moment of inertia compared to the non-flashing condition. The presented graph in
Fig. 6.33 further supports our hypothesis. Fig. 6.33 represents the variations of the

volume-averaged vorticity magnitude multiplied with the mixture density over the
EOI phase. During the steady injection phase, i.e., when the needle is at the maxi-
mum lift, the ambient pressure does not impact the pressure inside the sac. However,

83
post injection, the sac region experiences a lower pressure front for the flashing con-
dition compared to the non-flashing case. Subsequently, the residual fuel inside the

sac region vaporizes and expands, resulting in a higher moment of inertia for the
flash-boiling condition. Hence, Fig. 6.33 indicates a lower magnitude of Enstrophy

for the flashing condition. The evaluation procedures for Enstrophy in a rotating flow
are mentioned below.

Z
Enstrophy = ρ (||Ω||2tr ) d∀ (6.1)
→
−
Z
Enstrophy = ρ (| Ω |2 ) d∀ (6.2)
→
−
Ω = [Ωi , Ωj , Ωk ] (6.3)
Z
Ω2i + Ω2j + Ω2k d∀


Enstrophy = ρ (6.4)
1 T

Ω = ∇U − (∇U) (6.5)
2
 
 0 Uj,i − Ui,j Uk,i − Ui,k 
 
Ω = 
 Ui,j − Uj,i 0 Uk,j − Uj,k 
 (6.6)
 
Ui,k − Uk,i Uj,k − Uk,j 0
 
 0 Ωk −Ωj 
 
Ω = 
−Ωk 0 Ωi  (6.7)
 
Ωj −Ωi 0

→
−
In Eqns. 6.1, and 6.2, ||Ω|| and Ω correspond to the spin tensor and the vorticity
vector, respectively. The calculation procedures for the spin tensor are mentioned in

Eqs. 6.5, 6.6.

6.2.4 Tip wetting

In an ideal scenario, the sac should be free of any trace of liquid fuel and should

be filled with the non-condensible gas post-needle closure. If not, this liquid leads
to tip wetting post-injection. Eventually, the unatomized liquid film on the injector

84
Figure 6.31: Iso-surface for λ2 = Figure 6.32: Iso-surface for λ2 =
−8e + 11 colored by the pressure −8e + 11 colored by the vapor vol-
ratio for the flashing condition at ume fraction for the flashing condi-
t = 1.042ms, i.e., after the needle tion at t = 1.042ms i.e. after the
closure. needle closure.

Figure 6.33:
 (a) Volume averaged density (ρ) multiplied with the vorticity magni-
⃗
tude |Ω| vs time for the both non-flashing and flashing conditions. (b) The gray
colored area showing the sac and nozzle regions which are used for the volume aver-
aging.

85
tip leads to tip coking and soot formation [95]. Adverse effects of coked injectors on
further spray and combustion cycles have also been reported [17], and [202]. In this

section, we have identified the different tip-wetting phenomena for the non-flashing
and flash-boiling conditions.

Figure 6.34: The injector tip col- Figure 6.35: The injector tip col-
ored by the liquid mass fraction for ored by the liquid mass fraction for
the non-flashing condition at t = the flash-boiling condition at t =
0.98ms. 0.98ms.

Contours of liquid mass fraction on the injector tip for both the non-flashing and
flashing conditions during the needle closing phase, i.e., at t = 0.98ms are included in
Fig. 6.34 and 6.35, respectively. Clearly, for the flash-boiling condition, the injector

tip near the counterbore is more wetted during the needle closing phase compared to
the non-flashing condition. This phenomenon can be explained by the underexpanded
jet theory, postulated by Lacey et al. [98] for flash-boiling spray. Based on this theory,

when the underexpanded jet exits the counterbore region, it expands in the nozzle
vicinity. However, for the non-flash boiling condition, the ambient gas entrains into
the counterbore area creating buffer regions. These buffer regions inhibit the ensuing

liquid jet to wet the tip regions around the counter bore exit which are away from the
injector axis. Differences in jet behavior including the gas entrainment phenomena for
both conditions can further be observed via the mid-plane clips in Fig. 6.36, and 6.37.

86
Figure 6.36: Mid plane clip colored Figure 6.37: Mid plane clip col-
by the liquid mass fraction for the ored by the liquid mass fraction for
non-flashing condition at t = 0.98ms the flash-boiling condition at t =
showing NCG entrainment and the 0.98ms showing no NCG entrain-
liquid jet bending. ment.

Figure 6.38: The injector tip colored Figure 6.39: The injector tip colored
by the liquid mass fraction for the by the liquid mass fraction for the
non-flashing condition at t = 2.0ms. flash-boiling condition at t = 2.0ms.

87
 Furthermore, the tip-wetting behaviors are compared at t = 2.0ms, i.e., long after
the needle closure. Again, more tip wetting is witnessed for the flash-boiling condition

compared to the non-flashing condition, as seen in Fig. 6.38, and 6.39. These results
complement the findings by Leick et al. [106]. Fig. 6.39 also captures the boiling film

on the injector tip for the flash-boiling condition. Our current model formulations do
not include surface tension terms. Hence, it is impossible to quantify the amount of
fuel that sticks to the injector tip wall. Nonetheless, more studies are needed to fully

understand the tip-wetting phenomena.

6.2.5 Dribble

Dribble is another EOI phenomenon which is undesirable, as it involves large

droplets and unatomized ligaments. These unatomized fluid structures can be one
of the major sources of soot emission in the GDI engines. The diffused interface
approach coupled with the RANS turbulence model is not adequate to capture these

detailed dribbling structures. Even so, an indication of dribble can be obtained by

analyzing the ROI post injection [82]. Oscillations in the ROI can be observed in
Fig. 6.40.

To obtain further indications of the dribble in the non-flash-boiling condition,

the contours for the liquid mass fraction and the velocity magnitude are compared
on a mid-plane clip, as seen in Fig. 6.41, and 6.42. In Fig 6.41, the left plume

contains regions with high liquid mass fraction (Liquid M ass F raction ≈ 0.75) away
from the injector indicating liquid fuel. Fig. 6.42 shows high velocity magnitude
in these regions. Furthermore, the local maxima of the liquid mass fraction along

the plume are separated by local minima. Velocities corresponding to these local
minima of the liquid mass fraction are also found to be low, as shown in Fig 6.42.
These intermittencies in the velocity magnitude and the liquid mass fraction profile

along the liquid plume ensuing from the injector represent the disintegration of the

88
Figure 6.40: Simulated ROI vs time. Truncated ROI to indicate the oscillations post
needle closure.

liquid plume. The current modeling approach does not track the interface, hence,
cannot capture the jet break-up phenomena. Nonetheless, the above dispersion can be

treated as indications for the dribbling phenomena. However, our current simulations
in flashing conditions do not capture similar indications for the dribbling phenomena.
A pilot study to include the surface tension models in the in-house HRMFoam solver

will be expolred in the future chapter

6.3 Summary and conclusions

The current work includes the study of internal nozzle flow dynamics in the
EOI phase for a multihole GDI injector operating at both non-flash boiling and flash-
boiling conditions. Furthermore, the numerical predictions for the rate of injection

are validated with the experimental observations. The main findings of this study are
summarized as follows:

89
 Figure 6.42: Mid plane clip colored
Figure 6.41: Mid plane clip colored
by the velocity magnitude indicating
by the Liquid Mass Fraction indicat-
different velocities in the ensuing liq-
ing dribble at t = 1.038ms for the
uid jet at t = 1.038ms for the non-
non-flashing condition.
flashing condition.

• The ratio of ambient pressure to the saturation pressure has minimal effect on
the rate of injection, corroborating the fact that these nozzles operate in choked

conditions. The numerically observed spray injection duration is smaller com-

pared to the experimental observation because of the difference in the assumed
injector opening duration. These limitations can be avoided either by modify-

ing the experimentally measured needle displacement curve or by changing the

sealing constant. The modified curve should account for the sealing algorithm
and the limitation of the initial mesh lift. However, we choose not to modify

the experimentally observed boundary condition, i.e., the needle lift curve.

• Post need closure, the downstream low pressure waves are observed propagating
into the nozzles. Additionally, the ambient gas travels into the sac for the
non-flash boiling condition. Hence, a mixed sac condition is observed, where

liquid fuel and ambient gas coexist. This observation supports the experimental
findings in literature.

• For the flash-boiling condition, a boiling sac condition is observed at the EOI.
Similar observations have been made by the experimentalists with different fuel

90
 and a different injector. Furthermore, the expanding vapors do not allow the
ambient gas to enter into the sac region.

• During the needle closing phase, cavitating bubbles are found in the needle seat

area due to the water-hammer effect. Additionally, just after needle closure,
cavitating bubbles get generated in the sac due to the interactive sac vortices.
These observations are confirmed with the analysis of λ2 criterion. These bub-

bles can be detrimental to the injector’s needle seat and sac.

• Furthermore, the amount of gas and vapor present in the sac and nozzle post-
needle closure are also quantified. Under non-flashing conditions, the nozzle and
sac are found to have almost equal amounts of liquid fuel and non-condensible

gas long after needle closure. Yet, almost 90% fuel vapor is found in the flash-
boiling condition. These results need further confirmation by experiments.

• Different tip-wetting characteristics are also observed in the flash-boiling con-

dition compared to the non-flash boiling condition. During the needle closing

phase, the tip area surrounding the nozzle outlet is found to be more wetted
in the flash boiling condition compared to the non-flashing condition. This dif-
ference in observation can be attributed to the expansion characteristics of the

underexpanded flashing jets. Furthermore, long after completion of the injec-

tion, the injector tip is found to be covered by liquid fuel and its vapor in the
flash-boiling condition. This wetted tip can lead to tip coking during the com-

bustion cycle. Nonetheless, our current simulations lack the ability to quantify
the thickness of this liquid film.

• Indications of dribbling during non-flashing conditions are also captured.

The investigation of the EOI in the microsac GDI injector for a single injection
event under both flashing and non-flashing conditions is a novel and interesting result.

91
The current investigation captures the nuances of EOI in different ambient conditions.
The tip wetting and dribbling phenomena need further exploration and validation in

different ambient conditions. In the future chapter, the current modeling capabilities
of the diffused interface approach will be expanded to include formulations to model

interfacial phenomena.

92
 CHAPTER 7

A QUALITATIVE ANALYSIS OF THE TIP-WETTING

PHENOMENA UNDER FLASH-BOILING CONDITIONS

The efficiency of an engine, as well as the composition of the exhaust emis-

sion are driven by the in-cylinder combustion process. Effective combustion further

depends upon the preceding air-fuel mixing process [167]. The macroscopic character-
istics of the the fuel spray, i.e., the shape of jet, entrainment of hot gas, and behaviors
of individual liquid structures drive the mixing process [218]. Since the introduction

of the direct injection strategy, there have been clear advantages and disadvantages.
The direct injection strategy has improved the fuel economy while avoiding abnormal
knocking [244], [75]. Hence, GDI engines emit less green house gases (CO2 ) while

achieving a specific power output in the lower engine speeds [224], [117]. Nonetheless,
gasoline engines equipped with GDI injectors are associated with higher particulate
emissions compared to the port injected engines [102].

In general, wetted walls such as spray impinged combustion chamber, wetted pis-
ton crown and wetted injector tip are considered to be the major sources of particulate
matter (PM) and particulate number (PN) emissions [208]. These wetted areas gives

rise to pool fires [243], [115]. Furthermore, in endurance tests, PN emissions are
found to increase with the operation time, although the tests are started with clean
injectors [208], [202]. This increase in PN values is accredited to the thin layer of car-

bon deposit on the injector tip [208]. Injectors with a coked tip, i.e., fouled injectors
are observed to have different macroscopic spray characteristics compared to clean
injectors [112].

93
 The fuel film, developed on the injector tip can not fully evaporate before the
onset of combustion. Eventually, the film encounters the flame front which leads

to higher PN emissions and tip coking due to lack of oxygen and high in-cylinder
temperature [212]. Additionally, this porous coking structure aids accumulation of

more liquid fuel on the tip which further slows down the evaporation process, thus
facilitating the growth of carbon deposits [212]. Consequently, researchers [158] have
found direct relationships between the engine soot emissions and the wetting degree

of the injector tip.

However, the physics of tip-wetting phenomena under different ambient condi-
tions is poorly understood. Recently, a few experimentalists have made efforts to

understand the complex wall wetting phenomena. Leick et al. [106] from their ex-
periments with the Laser-induced Fluorescence (LIF) technique concluded that the
injection pressure and the flash-boiling conditions are the biggest factors in influenc-

ing the tip-wetting phenomena. Park et al. [152] and Song et al. [208] indicate more
tip-wetting in the flash-boiling conditions. However, they don’t explain the mech-
anism of the tip-wetting process. Based on the hypothesis proposed by Medina et

al. [123], the tip-wetting mechanism includes several distinct wetting processes such
as wide plume wetting, voretx tip wetting due to recirculation/air entrainment, and
fuel dribble wetting. According to them, wide plume wetting is encountered when the

liquid fuel reaches the outlet corner edge. Consequently, a film around the outlet is
formed due to surface roughness, wide plume and adhesive forces. Furthermore, they
theorize that the droplets stripping off the spray plume move towards the tip due to

the air entrainment vortices leading to the vortex tip-wetting phenomena. Similarly,
during the end of injection, they postulate the dribbling fuel to attach itself to the
injector tip because of large size and low axial momentum. However, their engine

experiments do not reveal the macroscopic/microscopic fluid dynamics behaviors to

support their hypothesis due to the lack of physical and optical access to the injector

94
tip. They also fail to explain the effects of different ambient conditions on the pro-
posed hypothesis. Recent LIF and long microscopic experiments by Xiao et al. [236]

support the theory of expansion wetting and the dribble wetting proposed by Medina
et al. [123]. Furthermore, they also indicate that higher fuel temperatures negates

the effects of flash-boiling on the tip-wetting phenomena. However, their experiments

are performed in a two-hole GDI injector, where plume-to-plume interactions are in-
significant. They also do not capture the internal flow dynamics of the injector which

would impact the wetting process because of optical inaccessibility.

Therefore, computational studies are required to understand the complex wetting
phenomena. Recently, Mouvanal et al. [133] simulated the EOI phase of a three hole

GDI injector using the commercial CFD tool, ANSYS Fluent. Nonetheless, they only
indicate certain wetted tip area around the counter-bore exit for the non-flashing con-
dition. Other computational efforts include the work by Fischer et al. [67] who employ

a Lagrangian framework coupled with the wall film model to visualize the wetting be-
havior. Their modeling framework is informed by the in-nozzle Eulerian simulations.
However, they generate the Lagrangian parcels inside the counter-bore regions where

the dense spray exists. Validity of the Lagrangian models in the near-nozzle regions
will be reviewed in the future chapters. Furthermore, Mohapatra et al. [127] indicate
different wetting behaviors under flashing and non-flashing conditions while simulat-

ing the internal flow by a diffuse interface based approach. The tip-wetting analysis
of the current chapter will be premised on the suggestions by Mohapatra et al. [127].
The chapter includes a detailed qualitative analysis of the wetting behaviors under

different ambient conditions. Besides, the simulation setups of the single injection
event described in the previous Sec. 6.1 are used for the current analysis. The current
modeling capabilities of the in-house code, discussed in the previous chapter 4 do not

include the surface tension model. These models are essential in low Weber number
flows, hence are not preferred for the spray modeling paradigm which involves high

95
Reynolds number and Weber number. Hence, the current modeling methods can not
fully resolve the slow moving tip-wetting flow. Nevertheless, in the current chapter,

indications of the tip-wetting flow are obtained by projecting the liquid fuel mixture
of the boundary cells on to the tip-wall.

7.1 Results and discussion

To study the detailed mechanisms involved in the tip-wetting phenomena un-

der different ambient conditions, the total simulation duration is subdivided into four
distinct subphases based on the needle motion transients. The first subphase corre-
sponds to the needle opening and quasi-steady state transients, whereas, the second

phase captures the mechanisms involved during the quasi steady state in addition to
the needle downward motion period. Subsequently, the third and fourth phases in-
clude the pre and post needle closure tip-wetting phenomena, respectively. The four

distinct phases are marked on the needle displacement curve in Fig. 7.1. The differ-
ent wetting mechanisms involved during these phases are discussed in the following
subsections.

7.1.1 Phase One: Needle opening phase

In the current simulations, phase one includes the transients between 0.34 ms
to 0.69 ms of the simulation period. Hence, during this phase, the needle moves from

the seated condition and attains the maximum steady needle lift. In the flashing
condition, indications of wetted tip around the counterbore exit is captured, as seen
in Fig. 7.2. However, the wetting is not prominent in the non-flashing condition. This

difference in observation is explained by the plume-wetting mechanism postulated by

Medina et al. [123], and further supported by Fig. 7.3. Fig. 7.3 captures the jet

expansion phenomena for the flashing condition. In contrast, the ambient air enters
into the nozzles for the non-flashing condition resulting in buffer layers. These buffer

96
Figure 7.1: Different phases of tip-wetting phenomena marked on the needle displace-
ment curve.

layers do not let the non-flashing jets to expand in the near nozzle regions resulting

in different tip conditions.

7.1.2 Phase Two: Downward moving needle phase

The downward moving needle phase corresponds to the interval of 0.7 ms to

1.05ms. This phase also includes 0.1 ms of the quasi-steady state. Fig. 7.2 indicates
the tip-wetting behaviors under both non-flashing (left) and flashing (right) condi-
tions. The figure exhibits traces of fuel mass on the tip around the injector axis for

the non-flashing condition. For non-flashing, the plumes are observed to be disper-
sive in nature due to the buffer layers present in the counterbore regions. Moreover,
vortices created around the plumes impart an upward motion to the ambient gas.

This upward moving gas carries the dispersed fuel mass towards the injector tip, refer
Fig. 7.5.The described wetting mechanism is similar to the vortex tip wetting theory
proposed by Medina et al. [123]. However, similar wetting process is not observed for

the flash-boiling condition. The expansion jet wetting mechanism still dominates the
tip-wetting process in flashing.

97
Figure 7.2: Tips colored by the liquid mass fraction (LMF) indicating different wetting
phenomena for the non-flashing(left) and flashing (right) conditions.

Figure 7.3: Mid-plane clip colored by the liquid mass fraction (LMF) indicating
different jet behaviors for the non-flashing(left) and flashing (right) conditions.

98
Figure 7.4: Tips colored by the liquid mass fraction (LMF) indicating different wetting
phenomena for the non-flashing(left) and flashing (right) conditions.

Figure 7.5: Mid-plane clip colored by the liquid mass fraction (LMF) with glyph
vectors scaled by the velocity magnitude indicating different mechanism driving the
wetting phenomena for the non-flashing(left) and flashing (right) conditions.

99
7.1.3 Phase Three: Needle closing phase

As discussed in previous Ch. 6, this phase corresponds to the duration when

the gap between the needle and the seat region is small enough to trigger the sealing
algorithm, i.e., 1.06ms − 1.2ms. Eventually the needle closes, disconnecting the

upstream flow domain from the sac region. As a result, the momentum of the spray
plumes starts to decrease creating dribbling like pheonmena, refer Sec. 6.2.5. However,
the ambient gas still maintains the motion towards the injector holes. Thus, vortices

get formed in the ambient, as shown in Fig. 7.6. These vortical structures carry
the slow moving fuel mass towards the injector. Hence, a significant amount of tip-
wetting is predicted immediately after the needle closure for both non-flashing and

flashing conditions, refer Fig. 7.7. This tip-wetting mechanism further supports the
fuel dribble wetting hypothesis of Medina et al. [123].

7.1.4 Phase Four: Post-needle closure phase

Fig. 7.8 includes the injector tip colored by the liquid fuel mass fraction for the
time stamps of 1.25 ms, 1.5 ms, and 2.00ms for different ambient conditions. At

this point, the needle has already been closed. The sealing algorithm discussed in the
Sec. 4.2 ensures no flow from the upstream of the seat to the sac-nozzle regions. The
figure suggests a gradual increment of fuel mass on the tip. Even so, for the non-

flashing condition, footprints of the fuel mass are found to decrease on the injector
tip. The previous discussion in Ch. 6 shows that the sac region attains the pressure
equilibrium with respect to the ambient condition. The discussion further indicates

mixed sac and boiling sac conditions for the non-flashing and flashing conditions,
respectively. Hence, it would be safe to conclude that post needle closure, the tip-
wetting phenomena is driven by the sac conditions. This conclusion is consistent with

the hypothesis by Medina et al. [123]. Thus, in flashing, the post needle closure tip
wetting is driven by the vapor dribbling phenomena.

100
Figure 7.6: Mid-plane clip colored by the liquid mass fraction (LMF) with glyph
vectors scaled by the velocity magnitude indicating the mechanism behind the wetting
phenomena for the non-flashing(left) and flashing (right) conditions during needle
closure.

101
Figure 7.7: Tips colored by the liquid mass fraction (LMF) indicating different wetting
phenomena for the non-flashing(left) and flashing (right) conditions during needle
closure.

102
Figure 7.8: Tips colored by the liquid mass fraction (LMF) indicating different wet-
ting phenomena for the non-flashing(left) and flashing (right) conditions post-needle
closure.

103
7.2 Conclusions
In the current chapter, the in-house solver HRMFoam in conjunction with the
sealing model is applied to simulate the tip-wetting phenomena. The current solver
does not include surface tension models which are essential to simulate the low Weber

number wall bounded flows. Hence, the solver can not fully resolve and quantify the
fuel mass sticking to the injector wall. However, a qualitative study of the tip-wetting
behaviors under different ambient conditions is performed by projecting the mixture

mass of the boundary cells onto the injector wall. This study is able to corroborate
the hypothesis put forward by experimentalists to explain the complex wall wetting
phenomena. Major findings of this study are summarized below.

• The complete injection duration is subdivided into multiple subphases to ex-

plain the distinct wetting mechanisms which are prevalent in different ambient

conditions.

• The flashing spray plume expansion drives the wall wetting behaviors around the

counterbore exit. However, in the non-flashing condition, wetting is inhibited

by the buffer regions between the spray plumes and the non-condensible gas.

• During the downward motion of the needle, traces of fuel mass are observed
on the injector tip for the non-flashing condition. Upward moving ambient gas
carrying the dispersed fuel is accounted for such tip wetting.

• A significant amount of tip-wetting is observed during needle closure phase

for both flashing and non-flashing conditions. During this period, wetting is
dominated by the dribbling fuel and upward moving ambient gas.

• Post needle closure, tip wetting is driven by the sac conditions. Hence, fuel mass
on the injector tip tend to decrease because of mixed sac conditions. However,

in the flash-boiling condition, boiling sac results in vapor dribbling which drives
the increasing tip-wetting behaviors.

104
 CHAPTER 8

QUASI STEADY STATE ANALYSIS OF MOVING

NEEDLE GDI SIMULATIONS

Gasoline direct injection (GDI) offers the opportunity to carefully tailor the
delivery of fuel. A designer could hypothetically control where and when the fuel

arrives in the cylinder. However, to take advantage of this opportunity, we must

have the capability of predicting how a specific injector will deliver fuel and how this
fuel will atomize and mix with the surrounding air. This predictive ability depends

on our understanding of the internal flow intricacies. The internal details of fuel
injectors have a profound impact on the emissions from gasoline direct injection (GDI)
engines [229]. However, the impact of injector design features is not completely

understood due to the difficulty in observing and modeling internal injector flows.
GDI flows involve moving geometry, flash boiling, and high levels of turbulent two-
phase mixing.

To understand the effects of the internal flow field of the GDI injector on the
near nozzle behavior, quasi-steady state analysis of the moving mesh simulations are
included in this chapter. In moving needle simulations, quasi-steady state corresponds

to the period when the needle is at the maximum lift. Hence, the flow field variables
are time averaged over the same duration, i.e., 0.4-0.9 ms in the current simulations
to eliminate the transient variability. This allows for a better understanding of the

effects of different ambient conditions (non-flashing and flashing) on the internal and
near nozzle spray behaviors during the steady injection phase. The simulation setups
described in the previous section 6.1 are used for the current analysis.

105
8.1 Results and discussion
8.1.1 Time averaged hydraulic coefficients

Mass flux measurements in addition to the momentum flux measurements, have

proven to be essential to characterize the hydraulic behavior of injection systems [155],

and [156]. Payri et al. [155], and [156] defined dimensionless coefficients such as dis-
charge coefficient (Cd ), momentum coefficient (CM ), velocity coefficient (Cv ), and area

coefficient (Ca ) to describe the internal flow behaviors based on these measurements.

• Discharge coefficient: This coefficient is defined as the ratio of real or mea-

sured mass flux (ṁ) to the maximum theoretical mass flux (ṁth ). Additionally,
the theoretical mass flux is calculated by considering the maximum Bernoulli

velocity (Uth ) crossing the total nozzle cross-sectional area (A0 ). The correla-
tions for the computationally and theoretically predicted mass flow rates are
represented in Eqn. 8.1, and 8.2.

Z
ṁ = ρ U dA (8.1)
A0

s
2∆P
ṁth = A0 ρf (8.2)
ρf
ṁ
Cd = (8.3)
ṁth
In Eq. 8.1, and 8.2, ρ, ρf , U, and ∆P represent the mixture density, density of
the liquid fuel, velocity at the nozzle exit, and the pressure drop, respectively.

• Momentum coefficient: Similarly, the momentum coefficient is defined as

 
the ratio of the measured momentum flux Ṁ to the maximum theoretical

106
  
momentum flux Ṁth , Eqn. 8.4. Again, the theoretical momentum flux is
evaluated based on the Bernoulli assumption, as shown in Eqns. 8.5 and 8.6.

Ṁ
CM = (8.4)
Ṁth

Ṁth = 2A0 ∆P (8.5)

Z
Ṁ = ρ U 2 dA (8.6)
A0

• Velocity coefficient This coefficient correlates the effective velocity (Uef ) to

the maximum theoretical Bernoulli velocity. It is calculated with the following

equation,

Uef Uef
Cv = =p (8.7)
Uth 2∆P/ρf
• Area coefficient. This coefficient is used for evaluating the reduction of the
effective area (Aef ) with respect to the geometric one, and it is calculated as

Aef
Ca = (8.8)
Ao
Based on these definitions, Mohapatra et al. [128] show alternative approaches to

evaluate the coefficient of velocity and the coefficient of area, as shown in Eq. 8.9,
and 8.10. These coefficients help in identifying the head loss at the nozzle outlets that

would otherwise be obscured by the coefficient of discharge. Furthermore, Lagrangian

models described in section 3.1 use these dimensionless coefficients to initialize parcels
in engine simulations.
CM
Cv = (8.9)
Cd
Cd
Ca = (8.10)
Cv
In the current section, the transient flow field variables such as density and veloc-
ity are averaged over the quasi-steady phase and the corresponding rate of injection

107
(ROI), rate of momentum (ROM), and dimensionless hydraulic coefficients are eval-
uated for both the flashing and non-flashing conditions. Furthermore, the overall

computational predictions for the non-flashing conditions are validated against the
available experimental data [155], as shown in Table 8.1. The percentage difference

between the experimentally measured ROI and the computational prediction for the
non-flashing condition is one percent; as a result the corresponding difference between
the coefficient of discharge is also less than one percent. However, a higher percentage

of difference is observed for the ROM, which translates into higher differences (around
10 %) in other dimensionless coefficients (CM , Cv , Ca ). These observed differences
can be attributed to two factors. The first factor is the experimental uncertainities in

spray momentum measurements, as indicated by Payri et al. [155]. The experiments

are performed while mounting the spray G injector in a frontal configuration, i.e.,
the injector axis aligns with the axis of the piezoelectric sensor. Hence, the spray

impact angle has to be considered for the momentum measurements. For the non-
flashing condition, the spray impact angle is different from the injector drill angle
and is highly transient [132]. Payri et al. [155] also observed an uncertainty of six

percent in the momentum measurement due to an error of 50 in the spray impact

angle measurement. The second factor is the modeling errors, as indicated by Moha-
patra et al. [128]. According to them, the CFD models over-predict the downstream

fuel density. Moreover, they suggest to revisit the mixing model of non-condensible
gas with the liquid fuel. A detailed comparison of the experimentally observed fuel
density and the CFD predictions will be included in the subsequent subsection.

Furthermore, Table 8.1 captures the differences in the hydraulic coefficients for dif-
ferent ambient conditions. In comparison to the flashing condition, higher coefficients
for the momentum and velocity are observed under the non-flashing condition. The

percentage of difference is around four percent. The Table 8.1 further indicates the
hole-to-hole variations in all computationally predicted hydraulic coefficients. For the

108
 - Cd CM Cv Ca ROI (g/s) ROM (N )
Experiment Non-flashing 0.52 0.40 0.77 0.67 13.82 2.63
CFD Non-flashing 0.52 0.44 0.85 0.61 13.63 2.89
All hole Flashing 0.51 0.42 0.82 0.62 13.60 2.86
Non-flashing 0.53 0.44 0.83 0.64 1.74 0.36
H1
Flashing 0.51 0.44 0.86 0.59 1.70 0.37
Non-flashing 0.54 0.46 0.85 0.63 1.76 0.37
H2
Flashing 0.52 0.44 0.85 0.61 1.72 0.37
Non-flashing 0.50 0.42 0.84 0.60 1.63 0.34
H3
Flashing 0.51 0.41 0.80 0.64 1.69 0.34
Non-flashing 0.52 0.44 0.85 0.61 1.70 0.36
H4
Flashing 0.51 0.43 0.84 0.61 1.70 0.36
Non-flashing 0.51 0.43 0.84 0.61 1.68 0.35
H5
Flashing 0.50 0.42 0.84 0.60 1.66 0.35
Non-flashing 0.53 0.46 0.87 0.61 1.75 0.38
H6
Flashing 0.52 0.41 0.79 0.66 1.72 0.35
Non-flashing 0.52 0.45 0.86 0.60 1.70 0.37
H7
Flashing 0.50 0.41 0.82 0.61 1.68 0.34
Non-flashing 0.51 0.43 0.84 0.61 1.66 0.35
H8
Flashing 0.51 0.44 0.86 0.59 1.72 0.37

Table 8.1: Time averaged hydraulic coefficients (Cd , CM , Cv , and Ca ), ROI, and
ROM for all holes combined and the individual holes under both flash-boiling and
non-flashing conditions.

109
Figure 8.1: Bar chart showing the hole-to-hole variations in the coefficient of velocity
under non-flashing and flashing conditions

flash-boiling condition, holes numbered one and six have the highest and the lowest
coefficient of velocity, respectively. Whereas, the same holes have the lowest and the

highest Cv under the non-flashing condition, as seen in Fig. 8.1. Similarly, the hole
numbered one has the highest coefficient of area in the non-flashing condition, but
the lowest in the flashing condition, refer Fig. 8.2. Under flashing, the hole numbered

six has the highest coefficient of area.

8.1.2 Hole-to-hole variations in ROI and ROM

Hole-to-hole variations in the predicted ROI and ROM are observed during both
the flashing and non-flashing conditions. To quantify these variations, the relative
standard deviations (RSD) in the ROI across all 8 holes are calculated using the

formulation represented in Eqn. 8.11.

s
P8

100 i=1 ṁi − ṁ
RSD = ∗ (8.11)
ṁ 7

110
Figure 8.2: Bar chart showing the hole-to-hole variations in the coefficient of area
under non-flashing and flashing conditions

In Eqn. 8.11, ṁ and ṁi refer to the instantaneous mass flow rate averaged over
all eight holes and the instantaneous mass flow rate through the hole numbered i,

respectively. The hole-to-hole variations in the ROI in terms of the RSD can be seen
in Fig. 8.3
A similar analysis for the momentum rate prediction is also performed. Momen-

tum rate is a good indicator of the discharge through individual holes. The RSD of the
ROI and ROM are transient in nature and have higher values during the quasi-steady
state in comparison to the needle opening/closing transients. These hole-to-hole vari-

ations under non-flashing conditions corroborate the findings by Baldwin et al. [19].
They attribute the variations to transient interactive sac vortices, extending between
multiple nozzle holes. Similar transient interacting vortex structures are also captured

under flash boiling conditions, as shown in Fig. 8.5.

Fig. 8.5 consists of the total pressure (14 M P a) isosurface colored by the static
pressure, indicating the string flash boiling phenomenon. Furthermore, it has a

mid-plane clip colored by the mixture density. The density contour indicates the

111
Figure 8.3: Relative standard deviation of the rate of injection across all the holes for
both the flashing and non-flashing conditions

112
Figure 8.4: Relative standard deviation of the momentum rate across all the holes
for both the flashing and non-flashing conditions

Figure 8.5: String flash-boiling in a transient moving mesh GDI simulation under
flashing condition

113
Figure 8.6: Line plots showing the Figure 8.7: Line plots showing the
time averaged liquid volume fraction for time averaged vapor volume fraction for
the hole1(left) and hole5(right) at 2mm the hole1(left) and hole5(right) at 2mm
downstream of the injector tip for both downstream of the injector tip for both
non-flashing (black color) and flashing non-flashing (black color) and flashing
(red color) conditions. (red color) conditions.

hollow-cone spray due to string flash boiling. The figure includes velocity streamlines
representing the vortex generation in individual holes. It also captures the vortex

interactions between the neighboring holes. These transient vortices are the primary
reason behind the hole-to-hole variations in the eight-hole GDI injector.

8.1.3 Time averaged analysis of the near nozzle (2mm away from the

injector tip) behaviors

For a better understanding of these in-nozzle transients and their effects on the

near-nozzle spray behavior, a cut plane is chosen at the 2mm downstream location
away from the injector tip. The chosen cut plane is perpendicular to the injector axis.
The concerned flow variables (LVF, VVF, and axial velocity) are interpolated to the

chosen 2-D cut plane and averaged over the quasi-steady state phase.
Figures 8.6, 8.8, 8.10, and 8.12 represent the predicted liquid volume fraction
at the downstream location for each nozzle, whereas figures 8.7, 8.9, 8.11, and 8.13

represent the vapor volume fraction prediction at the same locations. As expected, the

114
Figure 8.8: Line plots showing the time Figure 8.9: Line plots showing the time
averaged liquid volume fraction for the averaged vapor volume fraction for the
hole2 (left) and hole6 (right) at 2mm hole2 (left) and hole6 (right) at 2mm
downstream of the injector tip for both downstream of the injector tip for both
non-flashing (black color) and flashing non-flashing (black color) and flashing
(red color) conditions. (red color) conditions.

Figure 8.10: Line plots showing the time Figure 8.11: Line plots showing the time
averaged liquid volume fraction for the averaged vapor volume fraction for the
hole3 (left) and hole7 (right) at 2mm hole3 (left) and hole7 (right) at 2mm
downstream of the injector tip for both downstream of the injector tip for both
non-flashing (black color) and flashing non-flashing (black color) and flashing
(red color) conditions. (red color) conditions.

115
Figure 8.12: Line plots showing the time Figure 8.13: Line plots showing the time
averaged liquid volume fraction for the averaged vapor volume fraction for the
hole4 (left) and hole8 (right) at 2mm hole4 (left) and hole8 (right) at 2mm
downstream of the injector tip for both downstream of the injector tip for both
non-flashing (black color) and flashing non-flashing (black color) and flashing
(red color) conditions. (red color) conditions.

non-flashing condition predicts more liquid fuel at the downstream location compared
to the flashing condition. Furthermore, different slopes of these line plots ascertain

different spray angles for flashing compared to non-flashing conditions. Similarly,

higher vapor volume fraction distributions are predicted for the flashing conditions.
It should be noted that, although both flashing and non-flashing conditions predict

similar ROIs because of the choked state, contrasting downstream flow behaviors and
spray jet compositions are observed due to the dissimilar thermodynamic states.
Additionally, the axial velocity predictions for the individual plumes emanating

from the 8-hole injector are compared for different ambient conditions. Figures 8.14
and 8.15 represent the downstream axial velocity predictions for the non-flashing
and flash-boiling conditions, respectively. The negative velocities indicate the plume

going away from the injector, whereas the positive velocities indicate the entrained
chamber gas going towards the injector tip. For the flash-boiling conditions, the

plume velocities are predicted to be higher compared to the non-flashing conditions.

116
Figure 8.14: Axial velocity predictions at Figure 8.15: Axial velocity prediction at 2
2 mm downstream from the injector tip mm downstream from the injector tip for
for all nozzle holes under the non-flashing all nozzle holes under flash-boiling condi-
condition tion

Similarly, higher velocities are observed for the entrained chamber gas in the flashing
scenario.
Argonne measures the projected mass density at a distance of 2mm downstream

away from the injector tip using an X-ray measurement techniques [122]. The mea-
sured quantity is time averaged over the quasi-steady state phase, which is similar to
the time-averaged fuel density simulated by the CFD models. For the non-flashing

condition, the projected fuel mass density [122] is compared to the mixture density
at 2mm downstream location, as shown in Fig. 8.16. The CFD predictions of plume
shapes qualitatively match the X-ray tomography measurements for the non-flashing

condition, though the magnitude of the predicted mixture density is higher. Similar
evidences are captured when the density is averaged over all eight holes, as shown
in Fig. 8.17. In a similar fashion, the predictions for the flash-boiling condition for

the hole one and five are also compared and higher differences between the density
predictions are observed, refer Fig. 8.18.
These figures clearly show the predicted plumes have a very high concentration of

fuel as consequence of the difference in the momentum predictions, though the mass

117
Figure 8.16: Cut planes colored by the mixture density. (a) Density measured by
X-Ray for the non-flashing condition. (b) Mixture density predicted by HRMFoam
for the non-flashing condition.

Figure 8.17: 2-D line plot showing the comparison between the X-Ray measured
mixture density and the CFD predictions for the non-flashing condition. The mixture
density is averaged over all the 8 holes for the quasi-steady phase.

118
Figure 8.18: 2-D line plot showing the comparison between X-Ray measured mixture
density and the CFD predictions for the hole one (left) and hole 5( right) under flash-
boiling condition. The mixture density is time averaged over the quasi-steady phase.

119
flow rate matches with the experiments for both flash-boiling and non-flash-boiling
conditions. The mass flow rate corresponds to the product of fuel density multiplied a

velocity, integrated over an area. The fact that the mass flow rates are well predicted
by both CFD simulations indicates a good level of accuracy of these quantities within

the nozzle, which is the choke point of the flow. If a quasi-stable flow is assumed, the
the CFD code should presumably be predicting the correct mass flow rate at the 2
mm observation plane.

The mass flow rate constrains the product of density and velocity. Thus, the
downstream over-prediction in density is indicative that the predicted velocity is far
too low. This discrepancy in density predictions requires a revisit of the mixing

model of the non-condensible gas with the liquid fuel and the subsequent evaporation
due to heat transfer. The current solver uses the simple Fickian diffusion law to
account for the mixing of non-condensible gas and the liquid fuel. This law might

not be applicable to the highly vaporizing and turbulent sprays as observed in the
GDI injectors. Furthermore, the experimental predictions have shoulders at the peak
density predictions, as shown in Fig. 8.18. However, the bi-modal curves predicted by

the computation are found to be smooth in nature. This discrepancy indicates that
the experiments capture the plume-to-plume interactions which are not observed in
the computations.

Additionally, the current simulations use the k − ω SST to model the turbulence
closure. In general, conventional two-equation models such as the k − ω SST model
perform poorly for flows which are not in near-equilibrium state [27]. This class of

turbulence models are based on single time-scale schemes [124]. However, the GDI
flows involve multiple time-scales for both flashing and non-flashing conditions. A
detailed time-scale analysis for both thermodynamic conditions will be included in

the following sub-section.

120
 Time scales Non-flash-boiling Flash-boiling
Bernoulli time scale 6.94 × 10−7 s 6.84 × 10−7 s
Vaporization time scale 3.84 × 10−7 s 3.84 × 10−7 s
Sealing model time scale 2 × 10−6 s 2 × 10−6 s
Time scale of needle opening phase 4 × 10−5 s 4 × 10−5 s
Time scale of quasi steady phase 5.4 × 10−4 s 5.4 × 10−4 s
Time scale of pre-needle closure phase 1.6 × 10−4 s 1.6 × 10−4 s
Time scale of needle closure phase 2.0 × 10−5 s 2.0 × 10−5 s
Time scale of post-needle closure phase 1.5 × 10−4 s 1.5 × 10−4 s
Turbulent time scale for needle opening 5.41 × 10−7 s 1.83 × 10−7 s
Turbulent time scale for quasi-steady phase 1.33 × 10−7 s 1.43 × 10−7 s
Turbulent time scale for pre-needle closure 1.09 × 10−7 s 1.09 × 10−7 s
Turbulent time scale for needle closure 1.29 × 10−7 s 3.65 × 10−7 s
Turbulent time scale for post-needle closure 1.45 × 10−6 s 1.76 × 10−6

Table 8.2: Different time-scales for both flash-boiling and non-flash-boiling conditions

8.1.4 Time-scale analysis

The complex GDI flows involve multiple time-scales due to the phase-change
process, turbulence modeling, needle closure models and the needle motion. It is

difficult to comment on the interactive nature of these time scales. Hence, Table 8.2
includes a comparison of the order of magnitude of these time-scales for both non-
flashing and flashing conditions.

In tab. 8.2, the Bernoulli time scale refers to the theoretical time taken by the fluid
to cross the defined length of the nozzle (Lnoz ) for a prescribed pressure gradient (Lnoz )
under both non-flash-boiling and flash-boiling conditions, as shown in Eqn. 8.12.

r
ρ
TBern. = Lnoz (8.12)
2∆p

In Eqn. 8.12, ρ corresponds to the saturated liquid fuel density at the inlet tempera-
ture and pressure.

The time scale of vaporization is assumed to be the same order as the time constant
of the HRM model, as the other parameters such as volume fraction and the non-

121
dimensional pressure are of unity order. Hence, flashing and non-flashing conditions
have the same time-scale for vaporization.

The complete injection duration is subdivided into distinct phases based on their
impact on the flow field. These individual phases are classified as the needle opening

phase, quasi-steady phase, pre-needle closure phase, needle closure phase and the
post-needle closure phase. Furthermore, the turbulent time scale corresponds to the
time-scale of dissipation of the turbulent kinetic energy (k) at a dissipation rate of

ϵ. This time scale is also known as the eddy life time. As the turbulent time-scale is
formulated based on the over all energy and its dissipation, it is a scale of the larger
and more energetic eddies[62]. From the comparison presented in Table 8.2, it is

evident that the turbulent time scales during the injection phase are of similar order
as the vaporization time scales and the Bernoulli time scales. However, post-needle
closure the turbulent time scales are found to be larger and of the order of the sealing

model time scale. Furthermore, during needle opening phase, the eddy life time are
found to be lager in the non-flashing condition compared to the flashing condition.
In contrast, the trend reverses during the needle closure phase. This time scale also

supports our hypothesis that the fluid in the sac and nozzle regions are less rotational
in the flashing compared to the non-flashing condition. Hole-to-hole variations in
the turbulent time scales for different thermodynamic conditions are also captured.

These variations are tabulated in Table 8.3.

8.2 Inferences
The current work includes the study of internal nozzle flow dynamics in the
quasi-steady phase for a multi-hole GDI injector operating at both non-flash boiling
and flash-boiling conditions. The numerical predictions for the time-averaged rate of

injection, rate of momentum, and dimensionless hydraulic coefficients are validated

122
 Opening Quasi-steady Pre-closure Closure Post-closure
G1 5.37 × 10−7 s 1.29 × 10−7 s 1.01 × 10−7 s 1.31 × 10−7 s 1.44 × 10−6 s
H1
G2 1.81 × 10−7 s 1.47 × 10−7 s 1.25 × 10−7 s 3.66 × 10−7 s 1.67 × 10−6 s
G1 5.35 × 10−7 s 1.19 × 10−7 s 1.13 × 10−7 s 1.30 × 10−7 s 1.45 × 10−6 s
H2
G2 1.86 × 10−7 s 1.23 × 10−7 s 1.01 × 10−7 s 3.16 × 10−7 s 1.95 × 10−6 s
G1 5.41 × 10−7 s 1.38 × 10−7 s 1.13 × 10−7 s 1.35 × 10−7 s 1.68 × 10−6 s
H3
G2 1.85 × 10−7 s 1.48 × 10−7 s 0.89 × 10−7 s 3.29 × 10−7 s 1.83 × 10−6 s
G1 5.43 × 10−7 s 1.47 × 10−7 s 1.09 × 10−7 s 1.22 × 10−7 s 1.40 × 10−6 s
H4
G2 1.88 × 10−7 s 1.79 × 10−7 s 1.05 × 10−7 s 4.57 × 10−7 s 2.04 × 10−6 s
G1 5.41 × 10−7 s 1.21 × 10−7 s 0.97 × 10−7 s 1.37 × 10−7 s 1.66 × 10−6 s
H5
G2 1.80 × 10−7 s 1.39 × 10−7 s 1.08 × 10−7 s 3.81 × 10−7 s 1.76 × 10−6 s
G1 5.44 × 10−7 s 1.38 × 10−7 s 1.12 × 10−7 s 1.39 × 10−7 s 1.56 × 10−6 s
H6
G2 1.78 × 10−7 s 1.37 × 10−7 s 1.22 × 10−7 s 4.08 × 10−7 s 1.85 × 10−6 s
G1 5.47 × 10−7 s 1.34 × 10−7 s 1.14 × 10−7 s 1.38 × 10−7 s 1.33 × 10−6 s
H7
G2 1.89 × 10−7 s 1.37 × 10−7 s 1.13 × 10−7 s 3.17 × 10−7 s 1.56 × 10−6 s
G1 5.41 × 10−7 s 1.38 × 10−7 s 1.13 × 10−7 s 0.97 × 10−7 s 1.11 × 10−6 s
H8
G2 1.76 × 10−7 s 1.33 × 10−7 s 1.15 × 10−7 s 3.55 × 10−7 s 1.45 × 10−6 s

Table 8.3: Turbulent time scales for individual holes during different phases of the
needle motion.

with the experimental observations for the non-flashing condition. The main findings

of this study are summarized as follows:

• For different ambient conditions, the percentage differences in the predicted

ROI, ROM, Cd , CM , Cv , and Ca are around four percent. However, when the
ambient condition changes, the holes which have the highest and the lowest

dimensionless coefficients also change.

• Distinct hole-to-hole variations are observed in the predicted ROI and ROM.
A maximum relative standard deviation of fifteen percent is observed in the
predicted ROI for the non-flashing condition. These hole-to-hole variations are

driven by the transient sac vortices. These interactive vortices extend between
multiple nozzle holes.

• The hole-to-hole variations in the internal nozzle flow influence the downstream
flow behavior. Downstream flow behaviors are further impacted by the different

123
 ambient conditions. Consequently, hollow cone sprays with the vapor core at the
centre are observed during the flash-boiling condition. Moreover, the flashing

jets are observed to have higher axial velocities compared to the non-flashing
jets. These high jet velocities additionally result in higher entrained chamber

gas velocities.

• Discrepancies between the momentum measurements and the CFD predictions

for both flash-boiling and non-flashing conditions can majorly be attributed

to the modeling errors the CFD solver. The large error in the predicted and

measured fuel concentration downstream of the injector likely indicates an issue

with the fuel dispersion model. The Eulerian model of the fuel/air mixing
typically rely on the simplest possible closure, Fickian diffusion. This closure

presumes that turbulent mixing is entirely responsible for the fuel-air mixing.
Further, the fuel dispersion models are not cognizant of the large density ratio
between phases.

124
 CHAPTER 9

PLUME BASED COUPLING APPROACH FOR GDI

SPRAY

Modeling the complex air-fuel mixing process of GDI engines is a perennial chal-

lenge for computational fluid dynamics (CFD) simulations, as it involves predicting

the complete evaluation of the fuel spray from the liquid phase to the vapor phase with
adequate accuracy. Furthermore, complex phenomena such as break up and coales-

cence occur during the spray injection process. When the liquid fuel exits the nozzle
and penetrates into the ambient air, it goes through the primary atomization pro-
cess. During this process, the liquid jet disintegrates into ligaments and large droplets

which represent the dense spray core regions. These large liquid structures undergo
further breakup and coalescence processes to form stable droplets which are present
in the dilute spray regions. Recently, several attempts have been made to model these

complicated processes together. Among these approaches, the Lagrangian-Eulerian

(LE) methods such as the Discrete Droplet Model (DDM) [25] have dominated the
spray modeling in CFD. Efficient samplings of the actual drops have been enabled

by the stochastic nature of the Lagrangian-Eulerian model [190] while avoiding the
annoyance of numerical diffusion [61].
The LE approaches employ parcels (stochastic representation of a group of droplets)

to represent the Lagrangian liquid fuel drops. The accuracy of the predictions by the
DDM depends upon several submodels such as the droplet break-up model, drag
model, collision model, evaporation model, and turbulent dispersion model [233].

However, the spray models are droplet oriented and take little account of the neigh-
boring drops in the dense spray region [18]. It is very difficult to obtain the desired

125
predictions without adjusting the numerous model parameters of the LE models [57].
Hence, it can be concluded that these models are not truly predictive as the model

parameters are not known a priori, but require continuous manipulation. A recent
published work by Agrawal et al. [5] questioned the validity of the linear stability the-

ory, which is the central assumption of the KH-RT model, in the primary atomization
region, i.e., the dense spray region. Besides, these models require careful definitions for
droplet diameter, droplet temperature, droplet velocity, coefficient of discharge, spray

cone angle, and turbulent intensity to initialize the parcels. Defining these physical
parameters in GDI applications is a cumbersome process considering the fact that the
ensuing spray is highly transient in the near-nozzle region and features hole-to-hole

variations due to interactive vortices, as shown in Sec. 8.1. Furthermore, Secs. 8.1.1,
and 5.1 have shown that the GDI injectors operating at different ambient conditions
predict nearly similar mass flow rate and coefficient of discharge, although they ex-

hibit completely different internal and near-nozzle flow behaviors. Most of the LE
spray calculations are grossly under-resolved, as a higher number of cell counts [197]
and parcel counts [190] compared to the existing practices are required to achieve

converged results. Recently, researchers [209] [199] have applied the LE approach
to GDI applications and are able to match the experimental predictions after heavy
tuning of the model parameters, spray cone angle, and the turbulence model con-

stants. Moreover, their simulation setups lack general applicabilities, as the setups
tuned for one ambient condition require further tuning for better predictions in dif-
ferent ambient conditions. Paredi et al. [150] have observed good agreement between

the experimental observations and the LE predictions with minimal tuning of the
model constants for different ambient conditions. However, they vary the spray cone
angles and use the Pope turbulence correction [163] to fine tune the computational

predictions. Additionally, they have observed differences between the experimentally

generated PLV (Projected Liquid Volume) maps and the computationally predicted

126
PLV maps, although the predicted 1-D spray analysis such as liquid length, vapor
penetration, and axial velocity are close to the experimental measurements.

The second approach involves treating both the liquid and gaseous phases as
the Eulerian phase for numerical calculations. The two-fluid approach treats both

phases as a continuum throughout the domain, hence, requiring interface treat-

ment. Therefore, the Eulerian-Eulerian (EE) approaches which use sharp inter-
face methods are computationally expensive. To avoid these difficulties, Eulerian

paradigms [238], [34], [69] that model the interface instead of resolving the inter-
face have gained prominence. Other advantages include a lesser mesh dependency
compared to the LE approach [238] and no requirements of initial parcels’ defini-

tion. Furthermore, the Eulerian models are volume conserved as they account for the
volume occupied by the spray, which is generally ignored by the LE simulations [149].
The literature [34], [86], [116], [103], and [53] shows good agreement between the

experimental observations and the computational predictions by the EE approach

in the dense spray region, but the accuracy of these models drops significantly in
the sparse regions of the spray. Additionally, coupling of the internal injector flow

with in-cylinder simulations is essential for predictive CFD for engine development.
Practical experience has shown that small details in the injector flow can make a large
difference in the engine emissions. Simultaneously, the internal injector flow evolves

over small time and length scales compared to the in-cylinder flow. In consequence,
the EE calculation will be prohibitively expensive, if applied to a complete engine
simulation.

To address the above limitations, Subramanian at al. [214] have suggested to use
the EE approach in the near-field region and the LE approach in the dispersed spray
region. However, coupling the Eulerian liquid phase with the Lagrangian liquid phase

is perplexing, as it requires knowledge about the overlap of the mesh and physical
quantities between the two representations. Furthermore, Saha et al. [184] have used

127
a diffused interface-based Eulerian approach to capture the flow field at the nozzle
exit of a 8-hole GDI injector and used those data to initialize the parcels for further

LE calculations. Nonetheless, the simulated discharge coefficient and the prescribed

nozzle diameter are used to determine the initial droplet diameters. Additionally,

they have used models based on the linear stability theory for primary atomization.
Recently, Nocivelli et al. [141] [140] have considered the Sauter mean diameter (SMD)
predicted by Ultra-Small-Angle X-ray Scattering (USAXS) measurements to define

the distribution of droplet diameters during the primary atomization process. Parcels
initialized with these droplets are considered for subsequent atomization. The US-
AXS measurements are taken at 1 mm downstream of the injector tip, yet the parcels

are generated at the counter-bore exit which represent the dense spray region. Sph-
icas et al. [209] have used similar locations to inject parcels while initializing the
droplet diameters by the predictions from the Eulerian-Lagrangian spray atomization

(ELSA) approach. Even so, these one-way coupling approaches [209] [184] [141] [140]
heavily depend upon the tuning of spray cone angles for better predictions in engine
simulations.

In the current chapter, a new plume-based coupling approach inspired by the one-
way coupling approach [184] is proposed. The one-way coupling approach applies
the droplet oriented LE models in the dense spray region, where they struggle the

most. However, the plume-based coupling approach runs the EE primary atomization
model in the dense spray region. Followed by the suggestions by Subramaniam and
O’Rourke [215], the plume-based coupling approach employs the LE models in the

sparse spray region, where they work the best because of their inherent formulations.
Because of the unique formulation of the plume-based coupling approach, it captures
the nuances of the in-nozzle and near nozzle behaviors and drives the LE simulations

accordingly. Thus, this approach does not depend upon the primary atomization
models which are based upon the linear stability theory. The plume-based coupling

128
approach is applied to the 8-hole GDI injector [1] operating under both flash-boiling
and non-flashing conditions. Predictions from the novel approach are also validated

against the existing experimental observations for different ambient conditions.

9.1 Model description

The current simulation approach includes two modeling frameworks. The first
modeling framework simulates the internal and near nozzle flow behaviors including

the primary atomization process using an Eulerian mixture model approach, that is
based upon the Homogeneous Relaxation Model (HRM), described in secs. 4.1, 4.2,
and 4.3. The second modeling framework incorporates the Lagrangian modeling

approach implemented in the commercial 3D CFD solver, Converge 2.4 [175], to

simulate the dilute spray regions. Furthermore, transition criteria are applied to
couple both modeling frameworks.

9.1.1 Lagrangian model

Our current approach includes the discrete droplet models [61] of CONVERGE
2.4 [175] to simulate the sparse spray region. This method couples the liquid spray

parcels to the Eulerian gas phase, which is discretized based on the finite volume
approach. The k − ϵ turbulence model is used to model the gas phase turbulence.
Subsequently, fluctuations in the parcel properties are modeled by the O’Rourke

turbulent dispersion model [144]. Besides, the dynamic droplet drag model, Tay-
lor analogy breakup (TAB) model [145], No Time Counter (NTC) droplet collision
model [193], and the Frossling evaporation model [68] are employed to model the

secondary break-up processes.

9.1.2 Plume-based coupling approach

It is computationally exorbitant to simulate the internal nozzle flow simula-

tion coupled with the in-cylinder engine simulations in a single attempt, because

129
the internal injector flow evolves over small time and length scales compared to the
in-cylinder flow. Thus, our current approach involves running and tabulating the

internal injector flow simulation apriori. Then, the approach uses a plume-based
one-way coupling to inform the subsequent Lagrangian calculations. The first step

of the coupling approach is to define a transition criteria, and the current transition
criteria is based on the inter-droplet spacing. A specific inter-droplet spacing can
be translated in to a unique liquid volume fraction (LVF) with the assumption of

uniform cubic lattice arrangement of droplets in a three-dimensional space with a

uniform droplet diameter. Direct numerical simulation (DNS) by Quan et al. [164]
reveals strong droplet-droplet interactions for separation distances of less than four

droplet diameters, which translates into a LVF of 0.01. However, we establish the
plume boundaries, i.e., the transition surface, as the isosurface of LVF=0.04, which
roughly correlates with the inter-droplet spacing of three droplet diameters. This

isosurface is tessellated with triangles and exported as a polygon file format (.PLY)
in the current workflow. While choosing a specific LVF to define the plume surfaces,
special attention need to be given to ensure that the Eulerian plumes are closed at the

tip. Closed tip plumes inherently ensure finite volume conservation while coupling
the Eulerian-Eulerian simulations with the Lagrangian-Eulerian simulations. A lower
value LVF, i.e., less than 0.04, demands a bigger Eulerian domain in comparison to

the current internal flow simulations. A bigger computational domain for the inter-
nal flow simulations significantly increases the present computational cost, hence, it
justifies the current choice of the inter-droplet spacing.

The next step of the coupling approach involves interpolating the flow field vari-
ables such as the temperature and velocity to the locations of the plume faces using
the post-processing methods. The current work uses the commercial post-processing

software Ensight to perform the flow field interpolation. Due to the inherently tran-
sient nature of the injection process, the plume data are written to the disk frequently.

130
Each of these data files includes the detailed information about the droplet size, ve-
locity, temperature, and the location of each plume face. The mass flow rate through

these triangular faces is also noted. As seen in the Figures 9.1, and 9.2, the plume
sizes and shapes are highly variable both temporally and spatially. Nonetheless, this

approach inherently captures the differences in the spray plume bending and the spray
cone angle for different ambient conditions, as shown in Figures 9.1 - 9.4.

Figure 9.1: Plumes indicated by the iso- Figure 9.2: Plumes indicated by the iso-
surface of LVF=0.04 colored by the ve- surface of LVF=0.04 colored by the ve-
locity magnitude at t = 0.5 ms for the locity magnitude at t = 1.0 ms for the
non-flashing condition. The injector tip non-flashing condition. The injector tip
is colored by the liquid mass fraction. is colored by the liquid mass fraction.

In the final step of the coupling approach, a C++ application named Ensight-
ToConverge is developed. This application reads the postprocessed plume data and

applies subsequent mass flow rate-based filtering to generate a single map to be used
for the parcel injection process in CONVERGE-2.4. The filtering process filters out
the triangular faces of the plume surface which have the velocity vector directed

into the plume. It ensures that the parcels have the initial velocity away from the
plume surface. The next challenge of the parcel injection process is to define the
parcel injection location at each time step of the simulation. The parcels are created

stochastically with probabilities proportional to the mass flow rate. For any given
point on a plume surface, the probability of a parcel appearing at that location is

131
Figure 9.3: Plumes indicated by the iso- Figure 9.4: Plumes indicated by the iso-
surface of LVF=0.04 colored by the veloc- surface of LVF=0.04 colored by the veloc-
ity magnitude at t = 0.5 ms for the flash- ity magnitude at t = 1.0 ms for the flash-
ing condition. The injector tip is colored ing condition. The injector tip is colored
by the liquid mass fraction. by the liquid mass fraction.

equal to the mass flow rate through the plume surface at that point divided by the
total mass flow rate. Hence, this algorithm creates a representative population of
parcels around the plume location. Once the algorithm defines a specific location in

the Lagrangian mesh to inject the parcel, it looks for the corresponding tabulated
Eulerian mesh data, i.e., the data interpolated to the plume surface. Furthermore,
this information is used to initialize each parcel with the corresponding droplet di-

ameter, velocity, temperature, and TKE based on its defined location. Additionally,
the droplet diameters are an output of the Σ − Y formulation, as explained in the
Sec. 4.3. Examples for the plume surface and the plume representation of the injected

parcels are shown in Figs. 9.5, and 9.6, respectively.

Mass conservation is ensured by coupling the mass flow rate predicted by the
internal flow simulations with the total mass of parcels to be injected at each instance.

The current simulation uses the computational ROI curves presented in Fig. 6.2.
Additionally, the mass of the individual parcel is kept at a fixed default value of
4.5 × 10−12 kg. This allows the injected parcel count in the Lagrangian calculation to

132
be responsive to the Eulerian ROI predictions. Figure 9.7 represents a detailed work
flow of the plume-based one-way coupling approach.

Figure 9.5: Plumes colored by the Figure 9.6: The location in x,y,z space of each
velocity magnitude with the ve- parcel created at a single time step i.e. 0.55 ms
locity vectors at t = 0.89 ms for after the start of injection. In this time step,
the non-flashing condition 6327 parcels are created.

9.2 Simulation setup

The plume-based coupling approach assumes that the injector flow is a driver

of the in-cylinder flow, but the details of the in-cylinder flow do not alter the injector
flow. Consequently, the Lagrangian calculations are performed at a much larger
time step with much coarser resolution compared to the internal injector simulations.

The simulation setup described in sec. 6.1 is used for the internal flow simulations.
Additionally, to account for the primary atomization process, the transport equation
for the Σ − Y model with two different critical Weber numbers (W ecr = 1 and 6) is

solved. The details of the Σ − Y model with the critical Weber number are described
in sec. 4.3.

9.2.1 Lagrangian spray simulations

To avoid interactions between the spray and solid walls, a cylindrical domain

with the diameter and height equal to 240 mm is used to simulate the constant volume
spray chamber, as seen in Fig. 9.8. The base mesh size is kept at a fixed value of

133
Figure 9.7: Flow chart indicating all the major steps and the software/application
used during the each step for the plume based one-way coupling approach.

Figure 9.8: (a) Computational domain for the Eulerian internal flow simulation. (b)
Computational domain for the Lagrangian external spray simulation.

134
2.5 mm. Furthermore, a fixed conical embedding region extending 50 mm from the
nozzle exit with three levels of refinement is defined. This embedding region accounts

for the strong gradient of momentum near the nozzles due to liquid spray injection.
Additionally, adaptive mesh refinement (AMR) based on the second derivative of

velocity is applied to keep the computational cost minimum while maintaining a fine
mesh resolution in the far field of the spray. The Lagrangian simulations do not
include the geometry of the injector or nozzles. However, the domain is assumed to

be filled with the non-condensible gas, which is consistent with the assumptions for
the Eulerian internal flow simulations. In addition, the boundary walls are assumed to
have the same temperature as the inlet liquid fuel temperature for both non-flashing

and flashing conditions.

9.3 Results and discussion

9.3.1 Qualitative validation

To understand the evolution of spray morphology under different ambient condi-

tions, the CFD predictions from the plume-based coupling approach are first qual-

itatively compared between the non-flashing and flashing conditions. Furthermore,

they are validated against images obtained from the Diffused Back Illumination (DBI)

technique. These images have been captured at the University of Melbourne and are
made publicly available by the Engine Combustion Network (ECN) data set [1].
In DBI experiments, a light source and detection optical instruments including

the high speed camera are placed at either side of the spray. Before injection, the
camera records an image to characterize the intensity of the reference light source, i.e.,
I0 . During the injection event, the light intensity gets attenuated by the spray. The

corresponding light intensity (I) is determined by the detection optical instruments.

I
Moreover, the ratio of the attenuated to reference light intensity, i.e., I0
is correlated

135
with the measured optical thickness (τ ) via Beer-Lambart law, as shown in Eqn. 9.1.
Details of this spray diagnostic technique are well explained by Westlye et al. [231].

I
= e−τ (9.1)
I0

To compare the Lagrangian parcel simulation results with the experimentally mea-
sured optical thickness, Magnotti et al. [118], [119] propose a postprocessing method

based on the Mie theory. The proposed correlation is presented in Eqn. 9.2. Before
applying the correlation, the 3-D parcel simulation data is projected on to a 2-D inter-
rogation plane with a finite mesh size of ∆x∆y by performing a line-of-sight integral.

Then, for each cell, the optical thickness (τpredicted ) is determined by considering the
extinction coefficient (Qext ), total no. of droplets (Nj ) in each parcel in the concerned
cell, and the corresponding droplet diameter (dj ). A value of τ > 1 reflects the liquid

region. Based on the literature [121], a value of 2.1 is selected for the extinction
coefficient, Qext .
#parcels
Qext π X
τpredicted = d2j Nj (9.2)
4∆x∆y j=1

Fig. 9.9, and 9.12 include the DBI images quoted from the ECN data set [1] for
both non-flashing and flashing conditions, respectively. These figures further capture
the time evolution of the spray morphology under different ambient conditions. In

the DBI images, the blue and red curves indicate the lower and upper bounds of the
liquid spray structure. The transient evolution of the spray morphology predicted
by the plume-based coupling approach with two different critical Weber numbers

under both non-flashing and flashing conditions are included in Figs. 9.10, 9.11, 9.13,
and 9.14. Visual inspections of these figures indicate good agreement between the
CFD predictions and the experimental observations for the temporal spray evolution

under contrasting thermodynamic states. In the non-flashing condition, the spray

morphology is found to be less sensitive to the critical Weber number. In contrast,

136
Figure 9.9: DBI images taken from the ECN data set [1] which show the liq-
uid penetration for the non-flashing condition at different time instances, i.e., (a)
0.1 ms ASOI, (b) 0.3 ms ASOI, (c) 0.5 ms ASOI, (d) 0.7 ms ASOI, (e)
0.9 ms ASOI, and (f) 1.1 ms ASOI.

Figure 9.10: Binarized images showing the liquid spray boundary predicted by the
plume based coupling approach with W ecr = 1 for the non-flashing condition at
different time instances, i.e., (a) 0.1 ms ASOI, (b) 0.3 ms ASOI, (c) 0.5 ms ASOI,
(d) 0.7 ms ASOI, (e) 0.9 ms ASOI, and (f) 1.1 ms ASOI.

137
Figure 9.11: Binarized images showing the liquid spray boundary predicted by the
plume based coupling approach with W ecr = 6 for the non-flashing condition at
different time instances, i.e., (a) 0.1 ms ASOI, (b) 0.3 ms ASOI, (c) 0.5 ms ASOI,
(d) 0.7 ms ASOI, (e) 0.9 ms ASOI, and (f) 1.1 ms ASOI.

in the flashing condition, when the critical Weber number increases, the predicted

spray plumes are found to be wider. Consequently, the liquid spray with the lower
critical Weber number penetrates the farther distance. As evident from sec. 4.3, a
higher critical Weber number leads to a lower equilibrium interfacial surface density

(Ωeq ) which in turn results in a higher SMD value on the plume surface. The change
in penetration distance due to a change in the critical Weber number indicates the
sensitivity of the plume-based coupling approach to the initial droplet diameter in the

flash-boiling condition. Furthermore, the spray plumes are found to penetrate farther
distances in flashing compared to the non-flashing condition, attributing to the fact
that the flashing condition has a lower ambient density than the non-flashing. These

visual agreements are better than other reported comparisons [140], [141], [150], [184],
and [209] for the same operating conditions.

138
Figure 9.12: DBI images taken from the ECN data set [1] which show the liquid pen-
etration for the flashing condition at different time instances, i.e., (a) 0.1 ms ASOI,
(b) 0.3 ms ASOI, (c) 0.5 ms ASOI, (d) 0.7 ms ASOI, (e) 0.9 ms ASOI, (f)
1.1 ms ASOI, (g) 1.3 ms ASOI, (h) 1.5 ms ASOI, (i) 1.7 ms ASOI

139
Figure 9.13: Binarized images showing the liquid spray boundary predicted by the
plume based coupling approach with W ecr = 1 for the flashing condition at different
time instances, i.e., (a) 0.1 ms ASOI, (b) 0.3 ms ASOI, (c) 0.5 ms ASOI, (d)
0.7 ms ASOI, (e) 0.9 ms ASOI, and (f) 1.1 ms ASOI.

140
Figure 9.14: Binarized images showing the liquid spray boundary predicted by the
plume based coupling approach with W ecr = 6 for the flashing condition at different
time instances, i.e., (a) 0.1 ms ASOI, (b) 0.3 ms ASOI, (c) 0.5 ms ASOI, (d)
0.7 ms ASOI, (e) 0.9 ms ASOI, and (f) 1.1 ms ASOI.

141
9.3.2 Liquid penetration

In an engine environment, liquid spray penetration influences the wetting of

internal components (cylinder and piston). For an effective comparison between the
experimental observations and computational predictions, the ECN community [1],

[150] has defined an approach based on the the projected liquid volume (PLV) method
for the DBI measurements. According to this approach, the measured optical thick-
ness (τ ) is correlated with the PLV, which is the integral of the liquid volume fraction

(LVF) along the cross-stream direction, y, as shown in Eqn. 9.3. The extinction
∗
coefficient (Cext ) and droplet diameter are defined based on the experimental mea-
surements [81]. In addition, ECN recommends two thresholds for the PLV to define

the liquid penetration, i.e., 2 × 10−4 mm3 /mm2 and 2 × 10−3 mm3 /mm2 .

3
πd y∞
Z
P LV = τ ∗6 = LV F dy (9.3)
Cext −y∞

To compare the predictions from the parcel simulations with experimentally pre-
dicted liquid penetration, first, an Eulerian LVF field is generated from the Lagrangian

parcel simulation. Then, the line-of-sight integrations of the 3-D LVF field along the
cross stream direction are performed to project the data onto a 2-D background
mesh. Finally, the liquid penetration is defined as the maximum axial position of any

plume with the projected value lower than the two experimentally defined thresh-
olds. Comparisons between the computationally predicted liquid penetration and
the experimentally measured penetration by two institutions, i.e., Sandia National

Lab.(SNL) and University of Melbourne (UoM) for both the LVF thresholds under
non-flashing and flashing conditions are included in Figs. 9.15, and 9.16, respectively.

For the non-flashing condition, the discrepancies in the experimental measure-

ments can be attributed to the fact that the UoM results maintain a nozzle temper-

142
Figure 9.15: The figure includes predictions for the liquid penetration by the plume-
based coupling approach for two critical Weber number.(1 and 6) under the non-
flashing condition. The predictions are compared with the experimentally measured
liquid penetration based on two LVF thresholds for the same condition by the Sandia
National Lab (SNL) [81] and University of Melbourne (UoM) [150]. The experimental
data are obtained from the ECN data base [1].

143
Figure 9.16: The figure includes predictions for the liquid penetration by the plume-
based coupling approach for two critical Weber number (1 and 6) under the flashing
condition. The predictions are compared with the experimentally measured liquid
penetration based on two LVF thresholds for the same condition by the Sandia Na-
tional Lab (SNL) [81] and University of Melbourne (UoM) [150]. The experimental
data are obtained from the ECN data base [1].

144
ature of 383K, whereas SNL maintains a nozzle temperature of 363K [1], which is
the nominal value. Again, the predicted liquid penetration and residence time match

well with the SNL’s data for the lower LVF threshold (2 × 10−4 mm3 /mm2 ). The
percentage differences in the maximum predicted liquid penetration and the predicted

liquid fuel residence time between the computational predictions for the unity crit-
ical Weber number and the SNL measurements for the lower LVF threshold under
the non-flashing condition are around one percent. Furthermore, more liquid pen-

etration and liquid fuel residence time are observed in the flashing condition both
computationally and experimentally compared to the non-flashing condition, as seen
in Fig. 9.16. However, the values predicted by the plume-based coupling approach are

higher compared to the experimental measurements. CFD simulation with W ecr = 6

has the closest penetration value to the SNL experiments for the lower LVF threshold,
and the percentage difference in the highest predicted liquid penetration is around

twenty percent. Thus, it can be inferred that the computationally predicted spray is
less evaporative than the experiment.

9.3.3 Vapor penetration

Vapor penetration is another macroscopic spray characteristic which influences

the combustion quality. Experimentally, the vapor penetration is measured by means
of a Schlieren technique. The details of the experimental measurements are explained

by Hwang et al. [81]. For an effective comparison between the computational pre-
dictions and the experimental findings, vapor penetration is defined as the maximum
axial distance from the injector tip where a mixture fraction 0.1% is found.

Figs. 9.17, and 9.18 include the vapor penetration comparison between the exper-
imental data [1] and the computational predictions for different ambient conditions.
Less difference in the vapor penetration measurements by the different institutions are

observed unlike the liquid penetration results. In non-flashing conditions, predictions

145
Figure 9.17: The figure includes predictions for the vapor penetration by the plume-
based coupling approach for two critical Weber no.(1 and 6) under the non-flashing
condition. The predictions are compared with the experimentally measured vapor
penetration based for the same condition by the Sandia National Lab (SNL) and
University of Melbourne (UoM). The experimental data are obtained from the ECN
data base [1].

146
by the plume-based coupling approach match well with the measured values. How-
ever, in flashing conditions, the plume-based coupling approach predicts higher vapor

penetration values compared to the experimental observations. The percentage dif-

ference between the measured and predicted vapor penetration values are four and six

percent in the non-flashing and flashing conditions, respectively. Similar to the liquid
penetration predictions, higher vapor penetrations are observed in the flashing con-
dition compared to the non-flashing due to lower ambient density. The non-flashing

CFD predictions are found to be insensitive to the initial droplet diameter, however,
lower vapor penetration for a higher critical Weber number is observed in flashing.

9.3.4 Axial gas velocity

Sphicas et al. [209], and [210] have performed particle image velocimetry (PIV)
measurements and reported the temporal evolution of the gas velocity on the injector
axis at 15 mm downstream of the injector tip for the non-flashing condition. For

further validation, the gas velocity predictions by the plume-based coupling approach
are compared with the experiments, as shown in Fig. 9.19. In Fig. 9.19, the negative
axial velocity refers to the condition when the ambient gas has an upward motion,

i.e., the gas is moving towards the tip of the injector. As evident from the figure,
during the injection phase (780µs ASOI), the predicted axial velocity is within the
experimental uncertainties. However, the plume-based simulations predict less posi-

tive axial velocities compared to the experiments. Additionally, in the non-flashing

condition, the axial velocity predictions are found to be sensitive to the critical Weber
number, unlike the penetration predictions. Henceforth, the initial droplet diameter

effects the momentum exchange between the ambient gas and the parcels. Similarly,
the comparisons between the predictions for the non-flashing and flashing conditions
are included in Fig. 9.20. Higher axial gas velocities are observed in the flashing con-

dition compared to the non-flashing. This difference in observations is attributable

147
Figure 9.18: The figure includes predictions for the vapor penetration by the plume-
based coupling approach for two critical Weber number (1 and 6) under the flashing
condition. The predictions are compared with the experimentally measured vapor
penetration based for the same condition by the Sandia National Lab (SNL) [81] and
University of Melbourne (UoM) [150]. The experimental data are obtained from the
ECN data base [1].

148
to the higher plume velocity due to the lower ambient gas density in flashing. In
flashing, the gas velocity is also found to be sensitive to the initial droplet diameter.

Figure 9.19: The figure includes experimental measurements for the axial velocity at
a location 15 mm downstream of the injector tip along the injector axis [209] and the
computational predictions by the plume-based coupling approach for the non-flashing
condition.

Additionally, Sphicas et al. [210] have subdivided the injection duration into sev-
eral subphases based on the observed gas velocity to determine the jet-to-jet inter-

actions under different ambient conditions. Similar phenomena are also captured
by the plume-based simulations for both flash-boiling and non-flashing conditions as
explained below.

• Pre-plume arrival period: This phase is defined as the phase before the heads of
the individual plumes reach 15 mm downstream of the injector tip [210]. During

149
Figure 9.20: The figure includes comparison of the axial velocity predictions by the
plume-based coupling approach under flashing and non-flashing conditions.

150
 this phase, although the plumes have not arrived at the predetermined locations,
the downstream gas already starts responding to the momentum exchange and

the pressure waves created in the upstream locations, as seen in Fig. 9.21. In
addition, in the flash-boiling condition, this response is quicker due to higher

plume velocities, refer Fig. 9.22. Furthermore, the direction of the gas moment
is purely vertical for the non-flashing conditions. This supports the hypothesis
proposed by Sphicas et al. [210] that the spray is creating a low pressure regime

that pulls the ambient gas up towards the injector tip, establishing a central
recirculation zone between the plumes. As time increases, the magnitude of the
upstream velocity increases.

• Time of plume arrival at the 15mm location: Similar to the experiments [210]

the plumes arrive around 200µs after the start of injection for the non-flashing
condition. As the plumes arrive at the 15mm location, an upward moving
ambient gas gets established, as seen in Fig. 9.23.

• Profile relaxation period: As seen in experiments [210], the velocity profile of

the upward moving gas becomes homogeneous with respect to the axial position

post arrival of the plumes. Moreover, the recirculation velocity plateaus after
increasing for certain durations, for the non-flashing condition. Evidence of
similar phenomena are captured in Fig. 9.24. In flashing, this state is attained

at a much faster rate.

• Uniform upward motion period: During this period, uniform gas velocity profiles
are observed on the 15 mm line both experimentally [210] and computationally.
Evidence of the same for the non-flashing condition is captured in Fig. 9.25

• Reversing period: Sphicas et al. [210] further defines this period postcomple-

tion of the injection. During this phase, the magnitude of the axial velocity
decreases and eventually it changes sign. Hence, the ambient gas no longer

151
Figure 9.21: Vectors showing the direction of the ambient gas movement for the non-
flashing condition and W ecr = 1 at time = 0.1ms ASOI. The red horizontal line
is located at 15 mm downstream of the injector tip. The black vertical line is the
injector axis.

152
Figure 9.22: Vectors showing the direction of the ambient gas movement for the
flashing condition and W ecr = 6 at time = 0.1ms ASOI. The red horizontal line
is located at 15 mm downstream of the injector tip. The black vertical line is the
injector axis.

153
Figure 9.23: Vectors showing the direction of the ambient gas movement for the non-
flashing condition and W ecr = 1 at time = 0.2ms ASOI. The red horizontal line
is located at 15 mm downstream of the injector tip. The black vertical line is the
injector axis.

154
Figure 9.24: Vectors showing the direction of the ambient gas movement for the non-
flashing condition and W ecr = 1 at time = 0.4ms ASOI. The red horizontal line
is located at 15 mm downstream of the injector tip. The black vertical line is the
injector axis.

155
Figure 9.25: Vectors showing the direction of the ambient gas movement for the non-
flashing condition and W ecr = 1 at time = 0.6ms ASOI. The red horizontal line
is located at 15 mm downstream of the injector tip. The black vertical line is the
injector axis.

156
 moves towards the injector tip, rather it moves in the direction of the plume
motion. Merged/connected plumes are accredited for this behavior of the am-

bient gas [210]. Similar evidences are captured by the plume-based coupling
approach for both flashing and non-flashing conditions, Figs. 9.26 and 9.27.

Figure 9.26: Vectors showing the direction of the ambient gas movement for the non-
flashing condition and W ecr = 1 at time = 1.2ms ASOI. The red horizontal line
is located at 15 mm downstream of the injector tip. The black vertical line is the
injector axis.

9.3.5 Momentum conservation

The plume-based coupling approach inherently accounts for the mass conserva-
tion, as the total number of parcels injected into the Lagrangian domain at a specific

time is driven by the ROI predictions from the Eulerian simulation. However, the

157
Figure 9.27: Vectors showing the direction of the ambient gas movement for the non-
flashing condition and W ecr = 6 at time = 1.2ms ASOI. The red horizontal line
is located at 15 mm downstream of the injector tip. The black vertical line is the
injector axis.

158
same cannot be said for the momentum conservation. The fuel emanating from the
nozzle mixes with the ambient gas due to turbulence and results in the momentum

exchange between the Non-condensible gas and the spray plumes. The plume-based
coupling approach does not account for this Eulerian gas momentum while initializing

the parcels in the LE simulations. Hence, there would be a definite difference between
the total momentum calculated at the nozzle exit and the total momentum of the
parcels. To account for this discrepancy, the total rate of momentum is evaluated

at each nozzle’s exit. The calculated rate of momentum is further integrated over
the whole injection duration to determine the total Eulerian momentum exiting the
individual nozzles. Similarly, for the plume-based simulations, the momentum of an

individual parcel injected during a specific time step is determined based on its fixed
mass and the initialized velocity. Then, the total momentum injected at the spe-
cific time step for the individual plume is computed by segregating the total number

of injected parcels based on their position vectors for that time instance. Finally,
the momentum of the parcels for each plume group are integrated over the injection
duration to evaluate the total Lagrangian momentum injected during a single spray

event.
Fig. 9.28 includes a comparison of the above described momentum for the non-
flashing condition. As expected, the total momentum injected into the Lagrangian

domain by the plume-based coupling approach is lesser compared to the Eulerian

momentum evaluated at the nozzles’ exits. The average percentage of difference
is around fourteen percent. In contrast, a similar comparison for the flash-boiling

condition reveals a different trend, as shown in Fig. 9.29. For the flashing condition,
the momentum of the plumes are found to be higher compared to the respective
nozzle’s exit. This might be indicative of the vapor acceleration of the downstream

plumes due to flashing. In the plume-based Lagrangian simulations, these momentum

discrepancies can be accounted for by introducing a source term in the Eulerian N-S

159
Figure 9.28: Total momentum comparison between the Eulerian HRMFoam and the
Lagrangian plume-based simulations for the non-flashing condition. The Eulerian
momentum is evaluated at the exit of the nozzle. The Lagrangian momentum is
calculated based on the parcel momentum

160
Figure 9.29: Total momentum comparison between the Eulerian HRMFoam and the
Lagrangian plume-based simulations for the flashing condition. The Eulerian momen-
tum is evaluated at the exit of the nozzle. The Lagrangian momentum is calculated
based on the parcel momentum

161
equation for the non-condensible gas. However, the source term needs to be modeled
based on different ambient conditions.

9.4 Summary and conclusions

A novel plume-based coupling approach is developed in the current chapter. This
approach couples the Eulerian spray simulation with the Lagrangian parcel simulation
based on the derived plume surfaces as the transition region. This approach is applied

to study the flow of a practical GDI injector operating under different ambient condi-
tions. Predictions from the Lagrangian simulations are further validated against the
experimental measurements obtained from the publicly available ECN database [1].

Key details of the plume-based coupling approach are summarized below:

• The plume-based coupling approach requires running the Eulerian injector flow
simulation with the primary atomization model, i.e., Σ − Y model, a priori. In
the current work, a formulation of the Σ − Y model based on the critical Weber

number is used. Two different W ecr , i.e., one and six are used to predict the
SMD value in the near nozzle regions.

• Transition criteria based on the inter-droplet spacing are developed. These

transition criteria are used to define the plume boundaries in the Eulerian

framework. Subsequent data tabulations on these plume boundaries are ac-

complished. These tabulated data are used in CONVERGE 2.4 to initialize the
parcels in the Lagrangian framework.

• Lagrangian simulations are performed for the ECN [1] prescribed G1(non-flashing)

and G2(flashing) conditions. Qualitatively, the spray morphologies predicted by

the CFD simulations are compared with DBI measurements obtained from the
ECN data set [1] for both non-flashing and flashing conditions, and well agree-

162
 ments are achieved. In flashing, higher axially and radially penetrating liquid
spray morphologies are captured compared to the non-flashing condition.

• Predictions for liquid penetration and vapor penetration by the plume-based

coupling approach under different ambient conditions are compared with the
DBI and Schlieren experiments performed at the Sandia National Lab [81] and
the University of Melbourne [150]. In the non-flashing condition, the percent-

age differences between the CFD predictions and the SNL measurements for
the maximum liquid and vapor penetrations are one and four percents, re-
spectively. Whereas, the same differences are twenty and six percent for the

flash-boiling condition. Predictions for the flash-boiling predictions can be im-

proved by increasing the size of the embedding region, which extends until 50

mm downstream.

• Lastly, the axial gas velocity predictions from the plume-based approach are

compared with the published PIV measurements [210]. CFD simulations are
able to capture all subphases described by Sphicas et al. [210] for the non-

flashing condition.

• The plume-based coupling approach needs additional source term modeling for

different ambient conditions to account for the momentum conservation.

163
 CHAPTER 10

PILOT STUDY AND CONCLUSIONS

The current dissertation is premised on the computational studies of near-nozzle

and external spray behaviors of a 8-hole GDI injector under different ambient con-

ditions (non-flash-boiling and flash-boiling conditions). The research includes the

development of a needle closure algorithm to study the end of injection phenomena.
Subsequently, the primary atomization model of the in-house developed HRMFoam

solver is improved. Finally, a novel plume-based coupling approach is developed for

external spray simulations.
The next section explores a pilot study to include the surface tension models in

the HRMFoam solver. These added features will enable the exploration of the low
Weber number flows, i.e., the tip-wetting flow. The following section will include the
concluding remarks and based on the major findings of this dissertation.

10.1 Pilot study: Dynamically coupled sharp and diffuse in-

terface approach
As discussed in Ch. 3, gasoline sprays are characterized by high Weber num-
ber and Reynolds number. Although a sharp interface exists between the liquid and

gaseous phase, when projected on to a finite sized computational mesh, the sharp
interface becomes diffused [48], [49]. The in-house HRMFoam solver uses this lim-
itation to its advantage and simulates the hydrodynamics of the spray atomization

process using a diffuse interface-based approach. Consequently, it is computationally

less expensive for engineering applications such as automotive sprays. Furthermore,

164
the inherent nature of the surface tension force attempts to make all liquid struc-
tures spherical in nature due to higher stability. However, in the dense spray region,

nonspherical liquid structures like liquid ligaments and distorted large droplets are
observed. Hence, the HRMFoam solver has been preferred in the dense spray region.

Nonetheless, the flow around the injector tip that drives the wetting phenomenon is
associated with low Weber number and Reynolds numbers. The wetting is further in-
fluenced by the solid-liquid contact angle. Hence, surface tension models are required

to correctly resolve and quantify the fuel film, as indicated in Ch. 7. VOF solvers
of OpenFOAM framework, such as interFoam can be applied to simulate such flows.
Thus, modifications in HRMFoam solver which are inspired by the interFoam solver

will be explored in the subsequent subsections.

10.1.1 Governing equation: interFOAM

A volume of fluid (VOF) method to solve multiphase problems was first postu-

lated by Hirt and Nicholas [78]. The method relies on the definition of an indicator
function to resolve the interface. The VOF approach solves a transport equation for
the indicator function in conjunction with other conservation equations. The indica-

tor function determines whether the computational cell is occupied by one fluid or
the other, or a mixture of both fluids. The VOF solver of OpenFOAM, interFoam,
is designed for incompressible flow without phase-change phenomena. Performances

of the solver have been evaluated by many publications [55], [222], [113]. The gov-
erning equations of the interFoam solver include conservation equations for mass and
momentum and a transport equation for the phase fraction (α), as shown below.

∇·U = 0 (10.1)
∂ρU
+ ∇ · (ρU U ) = −∇p∗ + g · x∇ρ + ∇ · τ + fσ (10.2)
∂t
∂α
+ ∇ · (U α) + ∇ · [Uc α (1 − α)] = 0 (10.3)
∂t

165
 In Eqn. 10.2, U is the velocity vector, p∗ is the modified pressure formulated by re-
moving the hydrostatic-(ρg · x) from the total pressure and τ is the shear stress tensor.

Furthermore, ρ is the mixture density calculated based on the expression in Eqn. 10.4
by accounting for the density of individual fluids and the volume fraction (α).

ρ = ρl (α) + ρg (1 − α) (10.4)

Additionally, In Eqn. 10.2, fσ is the volumetric surface tension force which is

estimated based on the continuum surface force model (CSF) proposed by Brackbill
et al. [36]. At first, the model determines the curvature of the interface based on the
gradient of the indicator function, α for the interFoam solver. Then, the CSF model

calculates fσ as per the expression presented in Eqn. 10.6.

∇α
κ = (10.5)
|∇α|
fσ = σκ∇α (10.6)

The transport equation for the volume fraction, Eqn. 10.6, includes an additional

convection term involving the compression velocity (Uc ), as suggested by Weller et

al. [230]. This artificial compression term attempts to compress the free surface
towards a sharper interface. However, this term is not included in the original VOF

formulation proposed by Hirt and Nicholas [78]. The compression velocity (Uc ) is
further determined based on the correlation presented in Eqn. 10.7

∇α
Uc = Cα |U | (10.7)
|∇α|

In Eqn. 10.7, Cα is a user defined constant that determines the effect of the compres-
∇α
sion term, and |∇α|
ensures the direction of the compression velocity is normal to the

166
interface. Additionally, to ensure the boundedness of the volume fraction, interFoam
uses the multidimensional universal limiter with explicit solver (MULES) [223] [182]

which is based on the flux corrected transport (FCT) technique [181].

10.1.2 Governing equations: HRMFoam with additional surface tension

models

Inspired by the surface tension modeling approach of the interFoam solver, the
liquid volume fraction, Ȳ , is chosen as the indicator function for HRMFoam. In HRM-

Foam, Ȳ is a derived quantity rather than a transported variable, and can be deter-
mined based on other transported mass fraction variables, as shown in Eqn. 10.8. Sub-
sequently, the volumetric surface tension force (fσ ) is modeled based on Ȳ . Eqn. 10.9

shows the modified momentum conservation equation of HRMFoam.

ρ
Ȳ = (1 − x) (1 − y) (10.8)
ρl
∂ρu
+ ∇ · (Φu) = −∇p + ∇τ̄¯ + f¯ + f¯σ (10.9)
∂t
f¯σ = σκ∇Ȳ (10.10)
Ȳ
κ = (10.11)
|∇Ȳ |

Unlike interFoam, the transport equation for the NCG mass fraction, i.e., Eqn. 4.5

of HRMFoam is modified to keep the computational cost minimal. An extra anti-

diffusion term inspired by the interface compression term is added to Eqn. 4.5, as
shown in Eqn. 10.12. Unlike interFoam, Ur is modeled based on the constitutive

relation proposed by Desantes et al. [52], as shown in Eq. 10.13. The correlation
proposed by Desantes et al. [52] depends upon the gradient of liquid mass fraction.

However, in the current formulation, the relative velocity is calculated based on the
NCG mass fraction.

167
  
∂ρy µt
+ ∇ · (Φy) + (1 − Ψ)∇ · (ρUr y) = −Ψ∇ · ∇y (10.12)
∂t Scr
1 hµ
t
i
ρUr = ∇y (10.13)
y(1 − y) Sc

In Eqn. 10.12, Ψ correspond to an indicator function which carries the value of

unity in the diffuse interface region, i.e., regions with high Weber number, whereas

it is zero in the sharp interface region, i.e., regions with low Weber number. For the
current pilot study, a zero value is chosen for the indicator function, Ψ.

10.1.3 Test case: Oscillating drop

For model validation, a 3-D oscillating drop is simulated with the proposed surface-
tension based HRMFoam solver. A perturbation of the velocity field vx = x, vy =

y, and vz = −2z is applied to the droplet. According to Lamb [101], the drop
oscillation period for this mode n = 2 perturbation is given by Eqn. 10.14.

 −0.5  −0.5
σ 8σ
T = 2π n(n − 1)(n + 2) 3 = 2π (10.14)
ρr0 ρr03

Figure 10.1: Cutplanes colored by the NCG mass fraction showing the initial droplet
with different mesh resolutions. (a) coarse (b) Fine

Thus, for a gasoline droplet of 10 mm radius, the period of oscillation is calcu-

lated to be 0.4 seconds. The current model is simulated on two meshes of different

168
resolutions, i.e., a coarse mesh with the mesh size of 3.3 mm and a fine mesh of size
0.67 mm, as shown in Fig. 10.1. Dai et al. [192] indicates a CFL criteria for the

free surface numerical stability. According to them, the maximum time step for an
oscillating drop can be evaluated based on Eqn. 10.15.

r
ρ∆x3
∆t < C (10.15)
2πσ

In Eqn. 10.15, C is a constant of order unity and ∆x is the cell size. The current
simulations use this stability criterion based on the surface tension to decide the
time step. The modified HRMFoam solver is able to accurately capture the pressure

difference across the interface, as shown in Fig. 10.2. The predicted pressure difference
(4.5 pa) satisfies the Laplace equation for pressure (4.4 pa). However, the predicted
period of oscillation does not agree with the exact solution calculated by Eqn. 10.14.

Although a finer mesh resolution improves the period prediction, but it is still far
from the ground truth, as seen in Fig. 10.3. Due to the inherent diffusive formulation
of the HRMFoam solver, the transported mass fractions (x and y) are diffusive in

nature. Consequently, the derived field, Ȳ is also diffused over multiple cells. Hence,
the current formulation for the indicator function, i.e., the use of Ȳ to track the
interface needs to be reformulated.

Hence, it can be concluded that a sharper interface formulation is required to

match the exact solution of the oscillating droplet case. This can be achieved by
formulating a new sharp color function inspired by the coupled level set volume of fluid

(CLSVOF) [35] method. This function can be built with either the NCG mass fraction
or the liquid volume fraction. Furthermore, a separate transport equation for the new
color function needs to be included in the in-house solver. The proposed version of

HRMFoam does not account for the wettability of the tip wall, which is required to
quantify the amount of fuel that sticks to the tip. According to Brackbill [36], this

169
can be accomplished by modifying the surface normal of the cells next to the wall
based on the pre-defined static/dynamic contact angle.

Figure 10.2: Mid-plane clip colored by the pressure difference indicating a higher
pressure inside the droplet, i.e. around 4.5 Pa for the fine mesh. The analytical
pressure difference in 4.4 Pa for the same droplet.

10.2 Summary and Conclusions

In this dissertation, a novel sealing algorithm is developed. The algorithm is
designed to get activated based on a prescribed needle displacement value. Once

activated, the algorithm applies a drag force in the sealing region. This force inhibits
any flow through the seat area without changing the mesh topology. The algorithm
is successfully integrated into the in-house developed HRMFoam solver.

170
Figure 10.3: Graphs indicating the oscillations of the droplet radius in z-direction for
two different meshes as predicted by the surface tension based HRMFoam

171
 Consequently, a multiple injection event is simulated in a single simulation cycle.
During this simulation, cavitation is captured at low needle positions. The simulations

also capture indications of dribbling fuel after completion of the injection event.
Furthermore, the end-of injection (EOI) phases of a single injection cycle under

different ambient conditions are analyzed. During this phase, ambient pressure waves
are observed to propagate into the sac region. As a consequence, a mixed sac con-
dition is witnessed for the non-flashing condition, i.e., the sac is filled with both

non-condensible gas and liquid fuel. In contrast, a boiling sac condition is captured
in the flashing condition, where the residual fuel starts boiling after needle closure.
The expanding sac fuel vapors inhibit any flow of ambient gas into the sac region.

Eventually, the sac compositions for both non-flashing and flashing conditions are
quantified. Post-closure, cavitation is also encountered inside the sac region in addi-
tion to the seat regions. A λ2 - criterion analysis confirmed the interactive sac vortices

to be the reason behind these sac bubbles.

A qualitative analysis of the tip-wetting phenomena reveals distinct wetting mech-
anisms for different ambient conditions. During the needle opening phase, the expan-

sion wetting mechanism dominates the wetting process. As a result, more tip wetting
is observed in flashing compared to the non-flashing condition. However, during the
early needle closing phase, i.e., when the needle starts moving downwards, the indica-

tions of fuel mass on the tip around the injector axis are captured in the non-flashing
condition. This wetting phenomenon is driven by the upward moving ambient gas
carrying traces of fuel due to the the dispersive nature of the non-flashing spray.

During the needle closure phase, more wetted tip areas are observed for both simu-
lated ambient conditions. During this phase, the wetting process is dominated by the
dribble wetting mechanism. Finally, post-closure tip wetting is influenced by the sac

conditions. Hence, more tip wetting is observed for the flash-boiling condition due to
the vapor dribbling from the boiling sac.

172
 A quasi-steady state analysis of a single injection event captures hole-to-hole vari-
ations in the rate of injection and the rate of momentum for both flashing and non-

flashing conditions. Furthermore, the predictions of hydraulic coefficients are com-

pared with the experimental observations. The percentage differences between the

computational predictions and the experimental observations are found to be around

four percent for the simulated ambient conditions. Furthermore, in-nozzle variations
are observed to influence the near-nozzle behaviors for both flash-boiling and non-

flashing conditions.
A novel plume-based coupling approach is also developed in this dissertation. This
approach involves running Eulerian near nozzle simulations with the primary atom-

ization model beforehand. Transition criteria to couple the Eulerian simulations with
the Lagrangian framework are developed based on the inter-droplet spacing.Then,
using the transition criteria the plume surfaces with tessellated area are defined at

a regular time interval. Data such as droplet diameter, temperature, TKE, density,
and velocity in these tessellated cells are tabulated to generate a single input file for
the Lagrangian simulations. Finally, in the Lagrangian simulations, the parcels in-

formed by these plume data are generated for the secondary atomization. Predictions
from the plume-based approach are validated both qualitatively and quantitatively
against the available experimental data, and good agreements are observed for both

non-flashing and flash-boiling conditions.

HRMFoam was originally formulated for high Weber number flows where the in-
terface is not sharp when projected on to a finite sized computational mesh. However,

to resolve the tip-wetting flow which is associated with low Weber number, surface
tension models are required. Hence, a pilot study is performed to include the sur-
face tension-based models in HRMFoam to resolve the interface. The study is able

to capture a proper response of the pressure field to the curvature of the interface.
However, when validated by an oscillating droplet simulation, the predictions for the

173
oscillation period are found to be in disagreement with the exact solution. Hence, a
logical conclusion could not be achieved from the pilot study.

174
 CHAPTER 11

FUTURE WORK

In this dissertation, attempts are made to answer some of the pressing questions

of the GDI systems under different ambient conditions. Several models are developed
and validated with the available experimental data. Although these models answer
some of the pressing concerns of the GDI system, many remain unanswered or un-

explored. Hence, there are still a lot of potentials for future exploration of the GDI
sprays under both flashing and non-flashing conditions. Based on the findings of the
current dissertation, brief descriptions about possible future explorations are included

in this section.
Although the current dissertation explores the EOI behaviors of a single injec-
tion event under different ambient conditions, it does not answer the effect of these

behaviors on subsequent injections. Thus, computational explorations need to be

performed while informing the initial sac conditions by the EOI behaviors of a single
injection event. Furthermore, EOI analysis of multiple injection events needs to be

executed for a better understanding of the post-closure near nozzle behaviors.

The current dissertation also discusses the mechanisms driving the tip-wetting
behaviors under different ambient conditions. However, it does not quantify the

wetting area due to the current limitations of the in-house solver. The pilot study
explores the integration of surface tension-based sharp interface models into the in-
house solver without the desired accuracy. Thus, different sharp interface models

need to be explored to resolve the fuel film on the tip of the injector under different
ambient conditions.

175
 The plume-based coupling approach developed in this dissertation shows promis-
ing results in both flashing and non-flashing conditions. However, predictions from

this plume-based approach are driven by the Eulerian primary atomization model,
i.e., the Σ−Y model. This model was initially developed and tuned for diesel injection

applications. Hence, no or a little validation of the same model exists in the literature
for the GDI applications. GDI systems operate in highly vaporizing conditions, i.e.,
flash-boiling conditions, which will definitely impact the interface production during

the primary atomization process. The current interfacial area density models do not
account for such effects, hence need further revisit of the model formulations. More-
over, the modeling constants of these models are heavily tuned for diesel applications.

These modeling constants have not been validated for gasoline application due to the
lack of experimental data set. Recently, experimentalists at Argonne have performed
ultra-small angle X-ray scattering (USAXS) measurements for GDI sprays under dif-

ferent ambient conditions. These experiments measure the path-specific surface area
of the spray, a measure of the interfacial area of the liquid-gas boundaries along the
X-Ray beam path. These measurements are equivalent to the interfacial surface den-

sity predictions of the primary atomization model. Hence, future research can be
performed by validating the primary atomization model, Σ − Y model, with ultra-
small angle X-ray scattering (USAXS) measurements for GDI sprays. The inherent

nature of the plume-based coupling approach ensures mass conservation, whereas it

does not conserve the momentum while coupling the Eulerian internal flow simula-
tions with the plume-based Lagrangian simulations. Momentum conservation can

further be ensured by introducing a source term in the Eulerian gas phase equation
of the plume-based Lagrangian simulation. However, this source term can not be an
empirical constant , as that would introduce unwanted errors. Future research should

be conducted to model this source term for different ambient conditions.

176
 In multihole injectors, for certain operating conditions, individual plumes merge
to form a single collapsing spray. The collapsing sprays tend to have a higher pen-

etration compared to other sprays. Hence, wall impingement is also observed for
collapsing conditions. Furthermore, when the spray collapses, the emission also de-

teriorates. The plume-based coupling approach does not discriminate between the
Eulerian plumes while creating parcels in the Lagrangian simulations. Thus, the
plume-based coupling approach should be employed to study the physics of collaps-

ing sprays, which are often observed in flare flashing conditions.

177
 BIBLIOGRAPHY

[1] Engine Combustion Network. https://ecn.sandia.gov/.

[2] Foam-Extend-3.2. https://github.com/Unofficial-Extend-Project-Mirror/

foam-extend-foam-extend-3.2.

[3] The Outlook for Energy: A View to 2040. http://cdn.

exxonmobil.com/~/media/global/files/outlook-for-energy/2016/
2016-outlook-for-energy.pdf.

[4] Agarwal, A., and Trujillo, M. F. A closer look at linear stability theory in
modeling spray atomization. International Journal of Multiphase Flow 109
(2018), 1–13.

[5] Agarwal, A., and Trujillo, M. F. The effect of nozzle internal flow on spray
atomization. International Journal of Engine Research 21, 1 (2020), 55–72.

[6] Akoi, I. Water flash evaporation under low pressure conditions. Heat Transfer
Japan Res. 23 (1994), 544–555.

[7] Akoi, I. Analysis of characteristics of water flash evaporation under low-pressure

conditions. Heat Transfer Asian Res. 29 (2000), 22–33.

[8] Aleiferis, P. G., Serras-Pereira, J., Augoye, A., Davies, T. J., Cracknell, R. F.,
and Richardson, D. Effect of fuel temperature on in-nozzle cavitation and spray
formation of liquid hydrocarbons and alcohols from a real-size optical injector
for direct-injection spark-ignition engines. International Journal of Heat and
Mass Transfer 53 (2010), 4588–4606.

[9] Allen, J., Hargrave, G. K., and Khoo, Y. C. In-nozzle and spray diagnostic
techniques for real sized pressure swirl and plain orifice gasoline direct injectors.
Tech. Rep. 2003-01-3151, 2003.

[10] Anez, J, Canu, R, Duret, Benjamin, Reveillon, J, and Demoulin, FX. Turbulent
statistical transition from Euler to Lagrange using droplet velocity PDF. In
ICLASS 2018-14th Triennial International Conference on Liquid Atomization
and Spray Systems (2018).

[11] Aori, G., Hung, D. L. S., Zhang, M., Zhang, G. M., and Li, T. Y. Effect of nozzle
configuration on macroscopic spray characteristics of multi-hole fuel injectors
under superheated conditions. Atomization and Sprays 26 (2016), 439–462.

178
[12] Aquino, C., Plensdorf, W., Lavoie, G., and Curtis, E. The occurrence of flash
boiling in a port fuel injected gasoline engine. Tech. Rep. 982522, 1998.

[13] Araneo, L., and Donde′ , R. Flash boiling in a multihole GDI injector - effects
of the fuel distillation curve. Fuel 191 (2017), 500–510.

[14] Araneo, L., Slima, K. Ben, and Dondé, R. Flash boiling effect on swirled injector
spray angle. ILASS-Europe, 2002.

[15] Arcoumanis, C., Badami, M., Flora, H., and Gavaises, M. Cavitation in real-size
multi-hole diesel injector nozzles. Tech. Rep. 2000-01-1249, 2000.

[16] Avedisian, CT. Bubble growth in superheated liquid droplets. Encyclopedia of

Fluid Mechanics 3 (1985), 130–190.

[17] Badawy, T., Attar, M. A., Hutchins, P., Xu, H., Venus, J. K., and Cracknell,
R. Investigation of injector coking effects on spray characteristic and engine
performance in gasoline direct injection engines. Applied energy 220 (2018),
375–394.

[18] Balachandar, S., Liu, K., and Lakhote, M. Self-induced velocity correction for
improved drag estimation in Euler–Lagrange point-particle simulations. Journal
of Computational Physics 376 (2019), 160–185.

[19] Baldwin, E. T., Jr., R. O. Grover, Parrish, S. E., Duke, D. J., Matusik, K. E.,
Powell, C. F., Kastengren, A. L., and Schmidt, D. P. String flash-boiling in
gasoline direct injection simulations with transient needle motion. Int. J. Mul-
tiphase Flow 87 (2016), 90–101.

[20] Bankoff, SG. Entrapment of gas in the spreading of a liquid over a rough surface.
AIChE journal 4, 1 (1958), 24–26.

[21] Bar-Kohany, T., and Levy, M. State of the art review of flash-boiling atomiza-
tion. Atomization and Sprays 26, 12 (2016).

[22] Battistoni, M., Som, S., and Longman, D. E. Comparison of mixture and mul-
tifield models for in-nozzle cavitation prediction. ASME J. Eng. Gas Turbines
Power 136 (2014).

[23] Battistoni, M., Swantek, D. J. Duke A. B., Tilocco, F. Z., Powell, C. F., and
Som, S. Effects of noncondensable gas on cavitating nozzles. Atomization and
Sprays 25 (2015), 453–483.

[24] Battistoni, M., Xue, Q., Som, S., and Pomraning, E. Effect of off-axis needle
motion on internal nozzle and near exit flow in a multi-hole diesel injector. SAE
Int. J. Fuels Lubr. 7 (2014), 167–182.

[25] Beale, J. C., and Reitz, R. D. Modeling spray atomization with the Kelvin-
Helmholtz/Rayleigh-Taylor hybrid model. Atomization and sprays 9, 6 (1999).

179
[26] Beau, P., Me´ nard, T., Lebas, R.and Berlemont, A., Tanguy, S., and Demoulin,
F. Numerical jet atomization: part II—modeling information and comparison
with DNS results. In Fluids Engineering Division Summer Meeting (2006),
vol. 47500, pp. 555–563.

[27] Bhattarai, N., Kim, G-H, and Park, S O. Numerical study of the wake flow of
a single rectangular low-profile vortex generator by using non-equilibrium eddy
viscosity model. In The 4th Asian Symposium on Computational Heat Transfer
and Fluid Flow Hong Kong (2013).

[28] Bianchi, G. M., Minelli, F., Scardovelli, R., and Zaleski, S. 3D large scale
simulation of the high-speed liquid jet atomization. Tech. rep., 2007.

[29] Bianchi, G.M., Pelloni, P., Toninel, S., Scardovelli, R., Leboissetier, A., and
Zaleski, S. Improving the knowledge of high-speed liquid-jets atomization by
using Quasi-Direct 3D simulation. Tech. Rep. 2005-24-089, SAE, 2005.

[30] Bilicki, Z, and Kestin, J. Physical aspects of the relaxation model in two-phase
flow. In Proceedings of the Royal Society of London A: Mathematical, Physical
and Engineering Sciences (1990), vol. 428, The Royal Society, pp. 379–397.

[31] Bird, GA. Perception of numerical methods in rarefied gasdynamics. Progress

in Astronautics and Aeronautics 117 (1989), 211–226.

[32] Blander, M., and Katz, J. L. Bubble nucleation in liquids. AIChE Journal 21
(1975), 833–848.

[33] Blinkov, VN, Jones, OC, and Nigmatulin, BI. Nucleation and flashing in noz-
zles—2. comparison with experiments using a five-equation model for vapor
void development. Int. J. Multiphase Flow 19 (1993), 965–986.

[34] Blokkeel, Gregory, Barbeau, B, and Borghi, R. A 3D Eulerian model to improve

the primary breakup of atomizing jet. SAE transactions (2003), 45–54.

[35] Bourlioux, A. A coupled level-set volume-of-fluid algorithm for tracking mate-

rial interfaces. In proceedings of the 6th international symposium on computa-
tional fluid dynamics, Lake Tahoe, CA (1995), vol. 15.

[36] Brackbill, J. U., Kothe, D. B., and Zemach, C. A continuum method for mod-
eling surface tension. Journal of computational physics 100, 2 (1992), 335–354.

[37] Brennen, C. E. Fundamentals of Multiphase Flows. Cambridge University Press,

2005.

[38] Brown, R., and York, J. L. Sprays formed by flashing liquid jets. AIChE Journal
8 (1962), 149–153.

180
[39] Brusiani, F., Negro, S., Bianchi, G. M., Moulai, M., Neroorkar, K., and Schmidt,
D. Comparison of the Homogeneous Relaxation Model and a Rayleigh Plesset
cavitation model in predicting the cavitating flow through various injector hole
shapes. Tech. Rep. 2013-01-1613, SAE, 2013.

[40] Butcher, A. J., Aleiferis, P. G., and Richardson, D. Development of a real-size

optical injector nozzle for studies of cavitation, spray formation and flash-boiling
at conditions relevant to direct-injection sprak-ignition engines. Internation J
of Enginer Research 14 (2013), 557–577.

[41] Chen, Y., and Zhang, Z. Mechanism of flash boiling and spray analysis with
gasoline, iso-octane, n-pentane and ethanol from a novel heated tip gdi injector.
Applied Thermal Engineering, 115 (2017), 1109–1117.

[42] Chesnel, J., Reveillon, J., Menard, T., and Demoulin, F. Large eddy simulation
of liquid jet atomization. Atomization and Sprays 21, 9 (2011).

[43] Cleary, V., Bowen, P., and Witlox, H. Flashing liquid jets and two-phase droplet
dispersion I: Experiments for derivation of droplet atomisation correlations.
Journal of Hazardous Materials 142 (2007), 786–796.

[44] Colarossi, M., Trask, N., Schmidt, D. P, and Bergander, M. J. Multidimen-

sional modeling of condensing two-phase ejector flow. International Journal of
Refrigeration 35, 2 (2012), 290–299.

[45] Cole, Robert. Boiling nucleation. Advances in heat transfer 10 (1974), 85–166.

[46] Corradini, M. L. Fundamentals of multiphase flow, 1997.

[47] Crua, C., Shoba, T., Heikal, M., Gold, M., and Higham, C. High-speed mi-
croscopic imaging of the initial stage of diesel spray formation and primary
breakup. In SAE Technical Paper 2010-01-2247 (10 2010).

[48] Dahms, R. N., Manin, J., Pickett, L. M., and Oefelein, J. C. Understanding
high-pressure gas-liquid interface phenomena in diesel engines. Proceedings of
the Combustion Institute 34 (2013), 1667–1675.

[49] Dahms, R. N., and Oefelein, J. C. On the transition between two-phase and
single-phase interface dynamics in multicomponent fluids at supercritical pres-
sures. Physics of Fluids 25 (2013).

[50] De Villiers, E, Gosman, AD, and Weller, HG. Large eddy simulation of primary
diesel spray atomization. SAE transactions (2004), 193–206.

[51] Demoulin, F., Beau, P., Blokkeel, G., Mura, A., and Borghi, R. A new model
for turbulent flows with large density fluctuations: Application to liquid atom-
ization. Atomization and Sprays 17, 4 (2007).

181
[52] Desantes, J. M., Garcia-Oliver, J. M., Pastor, J. M., and Pandal, A. A com-
parsion of diesel sprays CFD modeling approaches: DDM versus Σ-Y Eulerian
atomization model. Atomization and Sprays 26 (2016), 713–737.

[53] Desantes, JM, Garcı́a-Oliver, José Marı́a, Pastor, JM, Pandal, A, Baldwin,
E, and Schmidt, DP. Coupled/decoupled spray simulation comparison of the
ecn spray a condition with the Σ-Y Eulerian atomization model. International
Journal of Multiphase Flow 80 (2016), 89–99.

[54] Desderi, U., and Bidni, G. Study of possible optimisation criteria for geothermal
power plants. Energy Convers. Mgmt. 38 (1997), 1681–1691.

[55] Deshpande, S. S, Anumolu, L., and Trujillo, M. F. Evaluating the performance

of the two-phase flow solver interFoam. Computational science & discovery 5,
1 (2012), 014016.

[56] Dhir, Vijay K. Mechanistic prediction of nucleate boiling heat transfer–

achievable or a hopeless task? Journal of Heat Transfer 128, 1 (2006), 1–12.

[57] DOE, US. A workshop to identify research needs and impacts in predictive sim-
ulation for internal combustion engines (presice). In Workshop, US Department
of Energy (2011).

[58] Downar-Zapolski, P., Bilicki, Z., Bolle, L., and Franco, J. The non-equilibrium
relaxation model for one-dimensional flashing liquid flow. Int. J. Multiphase
Flow 22 (1996), 473–483.

[59] Duke, D J, Kastengren, A L, Matusik, K E, Swantek, A B, Powell, C F, Payri,

R, Vaquerizo, D, Itani, L, Bruneaux, G, Grover Jr, R O, Parrish, S, Markle,
L, Schmidt, D, Manin, J, Skeen, S. A., and Pickett, L. M. Internal and near
nozzle measurements of Engine Combustion Network Spray G gasoline direct
injectors. Experimental Thermal and Fluid Science 88 (2017), 608–621.

[60] Duke, D. J., Swantek, A. B., Tilocco, F. Z, Powell, C. F., Kastengren, A. L.,
Fezzaa, K., Neroorkar, K., Moulai, M., and Schmidt, D. X-ray imaging of
cavitation in diesel injectors. Tech. Rep. 2012-01-1404, SAE, 2014.

[61] Dukowicz, J. K. A particle-fluid numerical model for liquid sprays. Journal of

computational Physics 35, 2 (1980), 229–253.

[62] Durbin, Paul A, and Reif, BA Pettersson. Statistical theory and modeling for
turbulent flows. John Wiley & Sons, 2011.

[63] Duret, B, Reveillon, J, Menard, T, and Demoulin, FX. Improving primary

atomization modeling through DNS of two-phase flows. International Journal
of Multiphase Flow 55 (2013), 130–137.

[64] Elias, E., and Chambre, P. L. A mechanistic non-equilibrium model for two-
phase critical flow. Int. J. Multiphase Flow 10 (1984), 21–40.

182
[65] Elias, E., and Chambre, P. L. Bubble transport in flashing flow. Int. J. Multi-
phase Flow 26 (2000), 191–206.

[66] Faeth, GM, Hsiang, L-P, and Wu, P-K. Structure and breakup properties of
sprays. International Journal of Multiphase Flow 21 (1995), 99–127.

[67] Fischer, A., and Thelliez, M. Methodology and tools to predict GDI injector tip
wetting as predecessor of tip sooting. Tech. rep., SAE Technical Paper, 2018.

[68] Froessling, N. On the evaporation of falling drops. Tech. rep., ARMY BIO-
LOGICAL LABS FREDERICK MD, 1968.

[69] Garcia-Oliver, J. M., Pastor, J. M., Pandal, A., Trask, N., Baldwin, E., and
Schmidt, D. P. Diesel spray CFD simulations based on the Σ-Y Eulerian atom-
ization model. Atomization and Sprays 23, 1 (2013).

[70] Gavaises, M, Murali-Girija, M., Rodriguez, C., Koukouvinis, P., Gold, M., and
Pearson, R. Numerical simulation of fuel dribbling and nozzle wall wetting.
International Journal of Engine Research (2021), 1468087420985189.

[71] Gilles-Birth, I., Bernhardt, S., Spicher, U., and Rechs, M. Experimental in-
vestigation of the in-nozzle flow of valve covered orifice nozzles for gasoline di-
rect injection. Proceedings of the Seventh International Symposium on Internal
Combustion Diagnostics (2006), 59–78.

[72] Gopalakrishnan, S, and Schmidt, D P. A computational study of flashing flow

in fuel injector nozzles. SAE International Journal of Engines 1, 1 (2008),
160–170.

[73] Gorokhovski, M, and Herrmann, M. Modeling primary atomization. Annu.

Rev. Fluid Mech. 40 (2008), 343–366.

[74] Guo, H., Ding, H., Li, Y., Ma, X., Wang, Z., Xu, H., and Wang, J. Comparison
of spray collapses at elevated ambient pressure and flash boiling conditions using
multi-hole gasoline direct injector. Fuel 199 (2017), 125–134.

[75] H, Oh. Improvement of combustion characteristics in a spray-guided direct-

injection sprak-ignition engine. PhD thesis, 2013.

[76] Henry, R. B., and Fauske, H. K. The two-phase critical flow of one-component
mixtures in nozzles, orifices and short tubes. Journal of Heat Transfer 93
(1971), 179–187.

[77] Herrmann, M. Refined Level Set Grid method for tracking interfaces. Annual
Research Briefs (2005), 3–18.

[78] Hirt, C W, and Nichols, B D. Volume of Fluid (VOF) methods for the dynamics
of free boundaries. J. Comput. Phys. 39 (1981), 201–225.

183
[79] Hoyas, S, Gil, A, Margot, X, Khuong-Anh, D, and Ravet, F. Evaluation of the
Eulerian–Lagrangian Spray Atomization (ELSA) model in spray simulations:
2d cases. Mathematical and Computer Modelling 57, 7-8 (2013), 1686–1693.

[80] Huh, KY. A phenomenological model of diesel spray atomization. In Proc. of

The International Conf. on Multiphase Flows’ 91-Tsukuba (1991).

[81] Hwang, J, Weiss, L, Karathanassis, I K, Koukouvinis, P, Pickett, L M, and

Skeen, S A. Spatio-temporal identification of plume dynamics by 3D computed
tomography using Engine Combustion Network spray G injector and various
fuels. Fuel 280 (2020), 118359.

[82] Hwang, J, Yasutomi, K, Arienti, M, and Pickett, L M. Numerical investigation

of near nozzle flash-boiling spray in an axial-hole transparent nozzle. Tech. rep.,
SAE Technical Paper, 2020.

[83] Ishii, E., Yoshimura, K., Yasukawa, Y., and Ehara, H. Late-fuel simulation
near nozzle outlet of fuel injector during closing valve. Journal of Engineering
for Gas Turbines and Power 138, 10 (2016), 102801.

[84] Ishii, M. One-dimensional drift-flux model and constitutive equations for rel-
ative motion between phases in various two-phase flow regimes. Tech. rep.,
Argonne National Lab., IL.(USA), 1977.

[85] Iyer, V., and Abraham, J. An evaluation of a two-fluid Eulerian-liquid Eulerian-

gas model for diesel sprays. J. Fluids Eng. 125, 4 (2003), 660–669.

[86] Iyer, V., and Abraham, J. Two-fluid modeling of spray penetration and disper-
sion under diesel engine conditions. Atomization and Sprays 15, 3 (2005).

[87] Janet, J. P., Liao, Y., and Lucas, D. Heterogeneous nucleation in CFD simu-
lation of flashing flows in convering-diverging nozzles. Int. J. Multiphase Flow
74 (2015), 106–117.

[88] Jeong, J., and Hussain, F. On the identification of a vortex. Journal of Fluid
Mechanics 285 (1995), 69–94.

[89] Kalikmanov, V.I. Nucleation Theory. [electronic resource]. Lecture Notes in

Physics: 860. Dordrecht : Springer Netherlands : Imprint: Springer, 2013.,
2013.

[90] Katz, J L, and Wiedersich, H. Nucleation theory without maxwell demons.

Journal of colloid and interface science 61, 2 (1977), 351–355.

[91] Khan, M M, Hélie, J, Gorokhovski, M, and Sheikh, N A. Experimental and

numerical study of flash boiling in gasoline direct injection sprays. Applied
Thermal Engineering (2017).

184
 [92] Kitamura, Y., Morimitsu, H., and Takahashi, T. Critical superheat for flashing
of superheated liquid jets. Ind. Eng. Chem. Fundam. 25 (1986), 206–211.

[93] Knapp, R. T., Daily, J. W., and Hammitt, F. G. Cavitation. McGraw Hill,
1970.

[94] Kocamustafaogullari, G., and Ishii, M. Interfacial area and nucleation site
density in boiling systems. Int. J. Heat Mass Transfer 26 (1983), 1377–1387.

[95] Koci, C P, Fitzgerald, R P, Ikonomou, V, and Sun, K. The effects of fuel–air

mixing and injector dribble on diesel unburned hydrocarbon emissions. Inter-
national Journal of Engine Research 20, 1 (2019), 105–127.

[96] Krämer, M., Kull, E., and Wensing, M. Flashboiling-induced targeting changes
in gasoline direct injection sprays. Int. J. Engine Research 17 (2016), 97–107.

[97] Lacey, J, Poursadegh, F, Brear, MJ, Gordon, R, Petersen, P, Lakey, C, Butcher,

B, and Ryan, S. Generalizing the behavior of flash-boiling, plume interaction
and spray collapse for multi-hole, direct injection. Fuel 200 (2017), 345–356.

[98] Lacey, J, Poursadegh, F, Brear, MJ, Gordon, R, Petersen, P, Lakey, C, Butcher,

B, and Ryan, S. Generalizing the behavior of flash-boiling, plume interaction
and spray collapse for multi-hole, direct injection. Fuel 200 (2017), 345–356.

[99] Lamanna, G. Cryogenic flashing jets:a review. Tech. Rep. AIAA 2016-4787,
2016.

[100] Lamanna, G., Kamoun, H., Weigand, B., and Steelant, J. Towards a unified
treatment of fully flashing sprays. International Journal of Multiphase Flows
58 (2014), 168–184.

[101] Lamb, Horace. Hydrodynamics. University Press, 1924.

[102] Leach, F, Stone, R, Richardson, D, Lewis, A, Akehurst, S, Turner, J, Remmert,

S, Campbell, S, and Cracknell, R F. Particulate emissions from a highly boosted
gasoline direct injection engine. International Journal of Engine Research 19,
3 (2018), 347–359.

[103] Lebas, R, Blokkeel, G, Beau, P-A, and Demoulin, F-X. Coupling vaporization
model with the Eulerian-Lagrangian Spray Atomization (ELSA) model in diesel
engine conditions. Tech. rep., SAE Technical Paper, 2005.

[104] Lebas, R, Menard, T, Beau, P-A, Berlemont, A, and Demoulin, F-X. Numerical
simulation of primary break-up and atomization: DNS and modelling study.
International Journal of Multiphase Flow 35, 3 (2009), 247–260.

[105] Lee, J. M., Fotache, C. R., Gopalakrishnan, S., and Schmidt, D. P. Flashing
flow of superheated jet fuel. Proceedings of the Combustion Institute 32 (2009),
3215–3222.

185
[106] Leick, P, Bork, B, and Geiler, J N. Experimental characterization of tip wetting
in gasoline di injectors. In Proceedings of the 14th triennial international con-
ference on liquid atomization and spray systems, Chicago, IL (2018), pp. 22–26.

[107] Leinhard, J. H., and Karimi, A. Homogeneous nucleation and the spinodal line.
ASME Journal of Heat Transfer 103 (1981), 61–64.

[108] Lemmon, EW, Huber, ML, and McLinden, MO. NIST standard refer-
ence database23: reference fluid thermodynamic and transport properties-
REFPROP, Version 9.1, NIST. Gaithersburgh. MD (2013).

[109] Liao, Y., and Lucas, D. 3D CFD simulation of flashing flows in a converging-
diverging nozzle. Nuclear Engineering Design 292 (2015), 149–163.

[110] Liao, Y, and Lucas, D. Computational modelling of flash boiling flows: A

literature survey. International Journal of Heat and Mass Transfer 111 (2017),
246–265.

[111] Lin, S-P. Breakup of liquid sheets and jets. Cambridge university press, 2003.

[112] Lindgren, R, Skogsberg, M, Sandquist, H, and Denbratt, I. The influence of

injector deposits on mixture formation in a DISC SI engine. SAE transactions
(2003), 852–861.

[113] Lobosco, RJ, Schulz, HE, and Simões, ALA. Analysis of two phase flows on
stepped spillways. Hydrodynamics-Optimizing Methods and Tools (2011), 1–25.

[114] Lothe, J, and Pound, G M. Reconsiderations of nucleation theory. The Journal

of Chemical Physics 36, 8 (1962), 2080–2085.

[115] Luo, H, Nishida, K, Uchitomi, S, Ogata, Y, Zhang, W, and Fujikawa, T. Mi-

croscopic behavior of spray droplets under flat-wall impinging condition. Fuel
219 (2018), 467–476.

[116] Luret, G, Menard, T, Berlemont, A, Reveillon, J, Demoulin, F-X, and Blokkeel,

G. Modeling collision outcome in moderately dense sprays. Atomization and
Sprays 20, 3 (2010).

[117] Mach, Wolfgang, and Bölter, Jens. 13 weniger reibung, aber mehr verschleiß:
Auswirkungen von dieselruß im schmieröl. versuchsresultate, interpretationen,
ursachenmodelle. Ölkreislauf von Verbrennungsmotoren III 3 (2009), 174.

[118] Magnotti, G M. Modeling the influence of nozzle-generated turbulence on diesel

sprays. PhD thesis, Georgia Institute of Technology, 2017.

[119] Magnotti, G. M., and Genzale, C. L. A Novel Spray Model Validation Method-
ology Using Liquid-Phase Extinction Measurements. Atomization and Sprays
25 (2015), 397–424.

186
[120] Maksic, S., and Mewes, D. CFD-calculation of the flashing flow in pipes and
nozzles. Proceedings of ASME Joint US-European Fluids Engineering Division
Conference (2002), 511–516.

[121] Manin, J, Jung, Y, Skeen, S A, Pickett, L M, Parrish, S E, and Markle, L.

Experimental characterization of DI gasoline injection processes. Tech. rep.,
SAE Technical Paper, 2015.

[122] Matusik, K, Duke, D, Sovis, N, Swantek, A, Powell, C, Payri, R, Vaquerizo,

D, Giraldo Valderrama, J S, and Kastengren, A. A study on the relationship
between internal nozzle geometry and injected mass distribution of eight ECN
Spray G nozzles. In Ilass Europe. 28th european conference on Liquid Atomiza-
tion and Spray Systems (2017), Editorial Universitat Politècnica de València,
pp. 313–320.

[123] Medina, M, Alzahrani, FM, Fatouraie, M, Wooldridge, MS, and Sick, V. Mecha-
nisms of fuel injector tip wetting and tip drying based on experimental measure-
ments of engine-out particulate emissions from gasoline direct-injection engines.
International Journal of Engine Research (2020), 1468087420916052.

[124] Menter, F. R. Improved two-equation k-omega turbulence models for aerody-

namic flows. Nasa Sti/recon Technical Report N 93 (1992), 22809.

[125] Mikic, B. B., Rohsenow, W. M., and Griffith, P. On bubble growth rates. Int.
J. Heat Mass Transfer 13 (1970), 657–666.

[126] Miyatake, O., Tanaka, I., and Lior, N. A simple universal equation for bubble
growth in pure liquids and binary solutions with a non-volatile solute. Int. J.
Heat Mass Transfer 40 (1997), 1577–1584.

[127] Mohapatra, C K., Jacobsohn, G L., and Schmidt, D P. The end of injection
study of a multi-hole gasoline direct injector under flash-boiling conditions.
Atomization and Sprays (2021).

[128] Mohapatra, C K, Schmidt, D P, Sforzo, B A, Matusik, K E, Yue, Z, Powell,

C F, Som, S, Mohan, B, Im, H G, Badra, J, et al. Collaborative investigation of
the internal flow and near-nozzle flow of an eight-hole gasoline injector (Engine
Combustion Network Spray G). International Journal of Engine Research,
1468087420918449.

[129] Mojtabi, M., Wigley, G., and Hélie, J. The effect of flash boiling on the atomiza-
tion performance of gasoline direct injection multistream injectors. Atomization
and Sprays 24 (2014), 467–493.

[130] Montanaro, A, and Allocca, L. Flash boiling evidences of a multi-hole GDI

spray under engine conditions by Mie-Scattering measurements. Tech. rep.,
SAE Technical Paper, 2015.

187
[131] Moon, S, Komada, K, Li, Z, Wang, J, Kimijima, T, Arima, T, and Maeda,
Y. High-speed X-ray imaging of in-nozzle cavitation and emerging jet flow of
multi-hole GDI injector under practical operating conditions. In 13th Triennial
International Conference on Liquid Atomization and Spray Systems, Tainan,
Taiwan (2015), vol. 2015.

[132] Moulai, M, Schmidt, D, Grover, R, and Parrish, S. Internal and near-nozzle

flow in a multi-hole gasoline injector under flashing and non-flashing conditions.
Tech. Rep. 2015-01-0944, SAE, 2015.

[133] Mouvanal, S, Burkhardt, A, Bakshi, S, and Chatterjee, D. Numerical study of

purging of a gasoline direct injection nozzle at the end of injection. International
Journal of Engine Research (2020), 1468087420916658.

[134] Nehmer, D A, and Reitz, R D. Measurement of the effect of injection rate

and split injections on diesel engine soot and nox emissions. SAE transactions
(1994), 1030–1041.

[135] Neroorkar, K., Gopalakrishnan, S., Jr., R. O. Grover, and Schmidt, D. P. Sim-
ulation of flash boiling in pressure swirl injectors. Atomization and Sprays 21
(2011), 179–188.

[136] Neroorkar, K, and Schmidt, D. A computational investigation of flash-boiling

multi-hole injectors with gasoline-ethanol blends. In SAE Technical Paper 2011-
01-0384 (2011).

[137] Neroorkar, K, and Schmidt, D. Modeling of vapor-liquid equilibrium of gasoline-

ethanol blended fuels for flash boiling simulations. Fuel 90, 2 (2011), 665–673.

[138] Neroorkar, K, Shields, B, Grover Jr, R O, Torres, A P, and Schmidt, D. Appli-

cation of the homogeneous relaxation model to simulating cavitating flow of a
diesel fuel. In SAE Technical Paper 2012-01-1269 (2012).

[139] Nishimura, A. A model for primary diesel fuel atomization based on cavitation
bubble collapse energy. In Proc. 8th ICLASS (2000), pp. 1249–1256.

[140] Nocivelli, L, Sforzo, B A, Tekawade, A, Yan, J, Powell, C F, Chang, W, Lee,

C, and Som, S. Analysis of the spray numerical injection modeling for gasoline
applications. Tech. rep., 2020.

[141] Nocivelli, Lorenzo, Yan, Junhao, Saha, Kaushik, Magnotti, Gina M, Lee, Chia-
Fon, and Som, Sibendu. Effect of ambient pressure on the behavior of single-
component fuels in a gasoline multi-hole injector. In Internal Combustion En-
gine Division Fall Technical Conference (2019), vol. 59346, American Society
of Mechanical Engineers, p. V001T02A012.

[142] Noh, and Woods. Simple line interface calculation. vol. 59, 5th Int. Conf. Fluid
Dyn.

188
[143] Olsson, E, and Kreiss, G. A conservative level set method for two phase flow.
Journal of computational physics 210, 1 (2005), 225–246.
[144] O’Rourke, P J. Statistical properties and numerical implementation of a model
for droplet dispersion in a turbulent gas. Journal of Computational Physics 83,
2 (1989), 345–360.
[145] O’Rourke, P J, and Amsden, A A. The TAB method for numerical calculation
of spray droplet breakup. Tech. rep., SAE Technical Paper, 1987.
[146] Oza, R. D. On the mechanism of flasing injection of initially subcooled fuels.
Journal of Fluids Engineering 106 (1984), 105–109.
[147] Oza, R. D., and Sinnamon, J. F. An experimental and analytical study of
flash-boiling fuel injection. Tech. Rep. 830590, 1983.
[148] Pai, M G, and Subramaniam, S. A comprehensive probability density function
formalism for multiphase flows. Journal of Fluid Mechanics 628 (2009), 181–
228.
[149] Pakseresht, Pedram, and Apte, Sourabh V. Volumetric displacement effects in
euler-lagrange les of particle-laden jet flows. International Journal of Multiphase
Flow 113 (2019), 16–32.
[150] Paredi, D, Lucchini, T, D’Errico, G, Onorati, A, Pickett, L, and Lacey, J. Vali-
dation of a comprehensive computational fluid dynamics methodology to predict
the direct injection process of gasoline sprays using Spray G experimental data.
International Journal of Engine Research 21, 1 (2020), 199–216.
[151] Park, B. S., and Lee, S. Y. An experimental investigation of the flash atom-
izaiton mechanism. Atomization and Sprays 4 (1994), 159–179.
[152] Park, C, Jung, J, Oh, H, Lee, J, and Bae, C. Injector tip wetting evaluation with
different nozzle outlet configurations. In International Conference on Liquid
Atomization & Spray Systems (2018), UNIVERSITY OF ILLINOIS.
[153] Parker, and Youngs. Two and three dimensional Eulerian simulation of fluid
flow with material interfaces. Tech. Rep. 01/92, UK Atomic Weapons Estab-
lishment, 1992.
[154] Patterson, M A, and Reitz, R D. Modeling the effects of fuel spray charac-
teristics on diesel engine combustion and emission. SAE transactions (1998),
27–43.
[155] Payri, R, Gimeno, J, Martı́-Aldaravı́, P, and Vaquerizo, D. Internal flow char-
acterization on an ECN GDi injector. Atomization and Sprays 26, 9 (2016).
[156] Payri, R, Molina, S, Salvador, FJ, and Gimeno, J. A study of the relation
between nozzle geometry, internal flow and sprays characteristics in diesel fuel
injection systems. KSME International Journal 18, 7 (2004), 1222–1235.

189
[157] Peter, E. M., Takimoto, A., and Hayashi, Y. Flashing and shattering phenom-
ena of superheated liquid jets. JSME International Journal 37 (1994), 313–321.

[158] Peterson, K, Grover, R, and Mitcham, C. Application of optical diagnostics

and simulation to fuel injector tip wetting and soot production. In in 11th
International Symposium on Combustion Diagnostics, Baden-Baden, Germany
(2014).

[159] Pickett, L M, Manin, J, Payri, R, Bardi, M, and Gimeno, J. Transient rate of

injection effects on spray development. In SAE Technical Paper 2013-24-0001
(9 2013).

[160] Pilch, M, and Erdman, CA. 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 (1987),
741–757.

[161] Pinhasi, G. A., Ullmann, A., and Dayan, A. Modeling of flashing two-phase
flow. Reviews in Chemical Engineering 21 (2005), 133–264.

[162] Plesset, M. S., and Prosperetti, A. Bubble dynamics and cavitation. Annual
Reviews of Fluid Mechanics 9 (1977), 145–185.

[163] Pope, SB. An explanation of the turbulent round-jet/plane-jet anomaly. AIAA

journal 16, 3 (1978), 279–281.

[164] Quan, S, Lou, J, and Stone, H A. Interactions between two deformable droplets
in tandem subjected to impulsive acceleration by surrounding flows. Journal of
fluid mechanics 684 (2011), 384.

[165] Rachakonda, S K. Computational Exploration of Flash-Boiling Internal Flow

and Near-Nozzle Spray. PhD thesis, University of Massachusetts Amherst Col-
lege of Engineering, University of Massachusetts Amherst, 2018.

[166] Rachakonda, S K, Wang, Y, Grover Jr, R O, Moulai, M, Baldwin, E, Zhang,

G, Parrish, S, Diwakar, R, Kuo, T-W, and Schmidt, D P. A computational
approach to predict external spray characteristics for flashing and cavitating
nozzles. International Journal of Multiphase Flow 106 (2018), 21–33.

[167] Rahman, MM, K Mohamme, Mohammed, and A Bakar, Rosli. Effects of air-
fuel ratio and engine speed on performance of hydrogen fueled port injection
engine. Journal of Applied Sciences 9, 6 (2009), 1128–1134.

[168] Rayleigh, L. On the pressure developed in a liquid during the collapse of a

speherical cavity. Philos. Mag. 34 (1917), 94–98.

[169] Reiss, HJLER, Katz, JL, and Cohen, ER. Translation–rotation paradox in the
theory of nucleation. The Journal of Chemical Physics 48, 12 (1968), 5553–5560.

190
[170] Reitz, R D. Modeling atomization processes in high-pressure vaporizing sprays.
Atomisation Spray Technology 3, 4 (1987), 309–337.

[171] Reitz, R. D. A photographic study of flash-boiling atomization. Aerosol Science

and Technology 12 (1990), 561–569.

[172] Reitz, R D, and Diwakar, R. Structure of high-pressure fuel sprays. SAE

transactions (1987), 492–509.

[173] Renardy, and Renardy. PROST:A Parabolic Reconstruction of Surface Tesnsion

for the Volume of Fluid Method. J. Comput. Phys. 183 (2002), 400–421.

[174] Reocreux, M. Contribution à létude des debits critiques en ecoulement

diphasique eau-vapeur. PhD thesis, Université Scientifique et Médicale, Greno-
ble, France, 1974.

[175] Richards, K J, Senecal, P K, and Pomraning, E. CONVERGE 2.4 manual,

2017.

[176] Rider, and Kothe. Reconstructing Volume Tracking. J. Comput. Phys. 141
(1998), 112–152.

[177] Rider, W. A marker particle method for interface tracking. In Proceedings of

6th International Symposium on Comput. Fluid Dynamics, 1995 (1995).

[178] Ritcher, H. J. Separated two-phase flow model:application to critical two-phase

flow. International Journal of Multiphase Flow 9 (1983), 511–530.

[179] Riznic, J. R., and Ishii, M. Bubble Number Density and Vapor Generation in
Flashing Flow. Int. J. Heat Mass Transfer 32 (1989), 1821–1833.

[180] R.O., Grover Jr, and S.E., Parrish. Characterization of a gasoline multi-hole
spray under closely-spaced multiple injection operating conditions. ILASS-
America, 2019.

[181] Rudman, M. Volume-tracking methods for interfacial flow calculations. Inter-

national journal for numerical methods in fluids 24, 7 (1997), 671–691.

[182] Rusche, Henrik. Computational fluid dynamics of dispersed two-phase flows

at high phase fractions. PhD thesis, Imperial College London (University of
London), 2003.

[183] Saha, K., Battistoni, M., Li, Y., Pomraning, E., and Senecal, P. K. Numerical
investigation of two-phase evolution of in- and near-nozzle regions of gasoline
direct injection during needle transients. SAE Int. J. Engines 9 (2016).

[184] Saha, K, Quan, S, Battistoni, M, Som, S, Senecal, PK, and Pomraning, E.

Coupled Eulerian internal nozzle flow and Lagrangian spray simulations for
GDI systems. Tech. rep., SAE Technical Paper, 2017.

191
[185] Saha, K., Som, S., Battistoni, M., Li, Y., Quan, S., and Senecal, P. K. Mod-
eling of internal and near-nozzle flow for GDI fuel injector. ASME Internal
Combustion Engine Division Fall Technical Conference.

[186] Santos, E G, Shi, J, Venkatasubramanian, R, Hoffmann, G, Gavaises, M, and

Bauer, W. Modelling and prediction of cavitation erosion in GDi injectors
operated with E100 fuel. Fuel 289 (2021), 119923.

[187] Scardovelli, R., and Zaleski, S. Direct numerical simulation of free-surface and
interfacial flow. Annual review of fluid mechanics 31, 1 (1999), 567–603.

[188] Schmidt, D. P., and Corradini, M. L. The internal flow of diesel fuel injector
nozzles:a review. Int. J. Engine Research 2 (2001), 1–22.

[189] Schmidt, D. P., Gopalakrishnan, S., and Jasak, H. Multi-dimensional simulation

of thermal non-equilibrium channel flow. Int. J. Multiphase Flow 36 (2010),
284–292.

[190] Schmidt, David P, and Bedford, Frederick. An analysis of the convergence

of stochastic Lagrangian/Eulerian spray simulations. International Journal of
Multiphase Flow 102 (2018), 95–101.

[191] Schmidt, David P., Corradini, M. L., and Rutland, C. J. Two-Dimensional,

Non-Equilibrium Model of Flashing Nozzle Flow. 3rd ASME/JSME Joint Fluids
Engineering Conference.

[192] Schmidt, David P, Dai, Meizhong, Wang, Haoshu, and Perot, J Blair. Direct
interface tracking of droplet deformation. Atomization and Sprays 12, 5&6
(2002).

[193] Schmidt, David P, and Rutland, CJ. A new droplet collision algorithm. Journal
of Computational Physics 164, 1 (2000), 62–80.

[194] Schmidt, DP, Gopalakrishnan, S, and Jasak, Hrvoje. Multi-dimensional simu-

lation of thermal non-equilibrium channel flow. International Journal of Mul-
tiphase Flow 36, 4 (2010), 284–292.

[195] Schmitz, I., Ipp, W., and Leipertz, A. Flash boiling effects on the development
of gasoline direct-injection engine sprays. Tech. Rep. 2002-01-2661, 2002.

[196] Sebastian, P., and Nadeau, J. Experiments and modelling of falling jet flash
evaporators for vintage treatment. Int. J. Thermal Sci. 41 (2002), 269–280.

[197] Senecal, PK, Pomraning, E, Richards, KJ, and Som, S. Grid-convergent spray
models for internal combustion engine computational fluid dynamics simula-
tions. Journal of Energy Resources Technology 136, 1 (2014).

192
[198] Serras-Pereira, J., van Romunde, Z., Aleiferis, P. G., Richardson, D., Wallace,
S., and Cracknell, R. F. Cavitation and primary break-up and flash boiling of
gasoline, iso-octane, n-pentane with a real-size optical direct-injection nozzle.
Fuel 89 (2010), 2592–2607.

[199] Shahangian, N, Sharifian, L, Miyagawa, J, Bergamini, S, Uehara, K, Noguchi,

Y, Marti-Aldaravi, P, Martı́nez, M, and Payri, R. Nozzle flow and spray de-
velopment one-way coupling methodology for a multi-hole GDI injector. Tech.
rep., SAE Technical Paper, 2019.

[200] Shen, S, Che, Z, Wang, T, Jia, M, and Sun, K. Numerical study on flash boiling
spray of multi-hole injector. SAE International Journal of Fuels and Lubricants
10, 2017-01-0841 (2017).

[201] Sher, E., Bar-Kohany, T, and Rashkovan, A. Flash-boiling atomization.

Progress in Energy and Combustion Science 34 (2008), 417–439.

[202] Shi, J., Santos, E. G., Hoffmann, G., and Dober, G. Large eddy simulation as
an effective tool for GDI nozzle development. MTZ worldwide 79, 10 (2018),
58–63.

[203] Shikhmurzaev, Y D. Capillary flows with forming interfaces. CRC Press, 2007.

[204] Shin, T. S., and Jones, O. C. Nucleation and flashing in nozzles-1. Int. J.
Multiphase Flow 19 (1993), 943–964.

[205] Shokomand, H., and Atashkadi, P. Performance improvement of single flashing,

binary, combined cycle for geothermal powerplant. Energy 22 (1997), 637–643.

[206] Skripov, V. P. Metastable liquids. New York : J. Wiley, 1973.

[207] Som, S, and Aggarwal, S K. Effects of primary breakup modeling on spray and
combustion characteristics of compression ignition engines. Combustion and
flame 157, 6 (2010), 1179–1193.

[208] Song, H, Xiao, J, Chen, Y, and Huang, Z. The effects of deposits on spray
behaviors of a gasoline direct injector. Fuel 180 (2016), 506–513.

[209] Sphicas, P, Pickett, L M, Skeen, S, Frank, J, Lucchini, T, Sinoir, D, D’Errico,

G, Saha, K, and Som, S. A comparison of experimental and modeled velocity
in gasoline direct-injection sprays with plume interaction and collapse. SAE
International Journal of Fuels and Lubricants 10, 1 (2017), 184–201.

[210] Sphicas, P, Pickett, L M, Skeen, S A, and Frank, J H. Inter-plume aerodynamics

for gasoline spray collapse. International Journal of Engine Research 19, 10
(2018), 1048–1067.

[211] Spielger, K., and El-Sayed, Y. The energtics of desalination processes. Desali-
nation 134 (2001), 109–128.

193
[212] Steimle, F, Kulzer, A, Richter, H, Schwarzenthal, D, and Romberg, C. System-
atic analysis and particle emission reduction of homogeneous direct injection SI
engines. Tech. rep., SAE Technical Paper, 2013.

[213] Strek, P., Duke, D., Swantek, A., Kastengren, A., Powell, C. F., and Schmidt,
D. P. X-Ray radiography and CFD studies of the Spray G injector. Tech. Rep.
2016-01-0858, SAE, 2016.

[214] Subramaniam, S. Lagrangian–Eulerian methods for multiphase flows. Progress

in Energy and Combustion Science 39, 2 (2013), 215–245.

[215] Subramaniam, S, and O’rourke, PJ. Numerical convergence of the KIVA-3 code
for sprays and its implications for modeling. Los Alamos Laboratory Report UR-
98-5465, Los Alamos, NM (1998).

[216] Sussman, M. A second order coupled Level Set and Volume of Fluid Method
for computing growth and collapse of vapour bubbles. J. Comput. Phys. 187
(2003), 110–136.

[217] Tanner, F X. Liquid jet atomization and droplet breakup modeling of non-
evaporating diesel fuel sprays. SAE transactions (1997), 127–140.

[218] Taskinen, P. Effect of fuel spray characteristics on combustion and emission

formation in a large medium speed diesel engine. Tech. rep., SAE Technical
Paper, 1998.

[219] Taylor, G I. The instability of liquid surfaces when accelerated in a direction

perpendicular to their planes. i. Proceedings of the Royal Society of London.
Series A. Mathematical and Physical Sciences 201, 1065 (1950), 192–196.

[220] Theofanus, T., Biasi, L., Isbin, H. S., and Fauske, H. A theoretical study on
bubble growth in constant and time-dependent pressure fields. Chem. Eng. Sci.
24 (1969), 885–897.

[221] Tow, TC, Pierpont, DA, and Reitz, Rolf D. Reducing particulate and N Ox
emissions by using multiple injections in a heavy duty DI diesel engine. SAE
transactions (1994), 1403–1417.

[222] Trujillo, M F, Hsiao, C-T, Choi, J-K, Paterson, E G, Chahine, G L, and Peltier,
L J. Numerical and experimental study of a horizontal jet below a free surface.
In 9th International conference on numerical ship hydrodynamics, Ann Arbor,
MI (2007).

[223] Ubbink, O. Numerical prediction of two fluid systems with sharp interfaces.

[224] Überall, A, Otte, R, Eilts, P, and Krahl, J. A literature research about particle
emissions from engines with direct gasoline injection and the potential to reduce
these emissions. Fuel 147 (2015), 203–207.

194
[225] Vallet, A., Burluka, A., and Borghi, R. Development of a Eulerian model for
the “atomization” of a liquid jet. Atomization and Sprays 11 (2001), 619–642.
[226] VanDerWege, B A., and Hochgreb, S. The effect of fuel volatility on sprays
from high-pressure swirl injectors. Twenty-Seventh Symposium (International)
on Combustion/The Combustion Institute.
[227] Versteeg, HK, Hargrave, GK, and Kirby, M. Internal flow and near-orifice spray
visualisations of a model pharmaceutical pressurised metered dose inhaler. In
Journal of Physics: Conference Series (2006), vol. 45, Iop Publishing, p. 207.
[228] Wallis, G. B. Critical Two-Phase Flow. Int. J. Multiphase Flow 6 (1980),
97–112.
[229] Wang, C, Xu, H, Herreros, J M, Wang, J, and Cracknell, R. Impact of fuel and
injection system on particle emissions from a GDI engine. Applied Energy 132
(2014), 178–191.
[230] Weller, H G. A new approach to VOF-based interface capturing methods for in-
compressible and compressible flow. OpenCFD Ltd., Report TR/HGW 4 (2008),
35.
[231] Westlye, F R., Penney, K, Ivarsson, A, Pickett, L M., Manin, J, and Skeen,
S A. Diffuse back-illumination setup for high temporally resolved extinction
imaging. Appl. Opt. 56, 17 (Jun 2017), 5028–5038.
[232] Wildgen, A., and Straub, J. The boiling mechanism in superheated free jets.
Int. J. Multiphase Flow 15 (1989), 193–207.
[233] Williams, F. A. Spray combustion and atomization. Phys. Fluids 1 (1958),
541–545.
[234] Woods, A., and Bloom, F. Modelling of flash evaporation I: Formulation of the
mathematical model. Math. Comput. Modelling 32 (2000), 1153–1169.
[235] Wu, S, Xu, M, Hung, D LS, Li, T, and Pan, H. Near-nozzle spray and spray col-
lapse characteristics of spark-ignition direct-injection fuel injectors under sub-
cooled and superheated conditions. Fuel 183 (2016), 322–334.
[236] Xiao, D, Qiu, S, Zhang, X, Zhang, Y, Li, X, Hung, D, and Xu, M. Dynamic
behavior and mechanism analysis of tip wetting process under flash boiling
conditions. Fuel 307 (2022), 121773.
[237] Xu, M., Zhang, Y., Zeng, W., Zhang, G., and Zhang, M. Flash boiling: Easy
and better way to generate ideal sprays than the high injection pressure. SAE
Int. J. Fuels Lubr. 6 (2013), 137–148.
[238] Xue, Q., Battistoni, M., Som, S., Quan, S., Senecal, P. K., Pomraning, E., and
Schmidt, D. P. Eulerian CFD modeling of coupled nozzle flow and spray with
validation against X-Ray radiography data. SAE Int. J. Engines 7 (2014).

195
[239] Zalesak, S. T. Fully multidimensional flux-corrected transport algorithms for
fluids. J. Comput. Phys. 31 (1979), 335–362.

[240] Zeng, W., Xu, M., Zhang, G., and Cleary, D. Atomization and vaporization for
flash-boiling multi-hole sprays with alcohol fuels. Fuel 95 (2012), 287–297.

[241] Zeng, W., Xu, M., Zhang, M., Zhang, Y., and Cleary, D. Macroscopic char-
acteristics for direct-injection multi-hole sprays using dimensionless analysis.
Experimental Thermal and Fluid Science 40 (2012), 81–92.

[242] Zhang, G., Hung, D. L., and Xu, M. Experimental study of flash boiling spray
vaporization through quantitative vapor concentration and liquid temperature
measurements. Exp. Fluids 55 (2014).

[243] Zhang, Y, Jia, M, Duan, H, Wang, P, Wang, J, Liu, H, and Xie, M. Numerical
and experimental study of spray impingement and liquid film separation dur-
ing the spray/wall interaction at expanding corners. International Journal of
Multiphase Flow 107 (2018), 67–81.

[244] Zhao, F, Harrington, D L, and Lai, M-C D. Automotive gasoline direct-injection

engines. Warrendale, PA: Society of Automotive Engineers, 2002. 372 (2002).

[245] Zhao, H., Quan, S., Dai, M., Pomraning, E., Senecal, P. K., Xue, Q., Bat-
tistoni, M., and Som, S. Validation of three dimensional internal nozzle flow
model including automatic mesh generation and cavitation effects. Journal of
Engineering for Gas Turbines and Power 136 (2014).

196