Numerical Study of Head-on Binary
Droplet Collisions: Towards
Predicting the Collision Outcomes
Vom Fachbereich Maschinenbau an der Technischen Universität Darmstadt
zur Erlangung des akademischen Grades Doktor-Ingenieur (Dr.-Ing.)
genehmigte Dissertation von Muyuan Liu M.Sc. aus Hebei, China
Tag der Einreichung: 31.05.2017, Tag der Prüfung: 17.10.2017
Darmstadt 2017 — D 17

1. Gutachten: Prof. Dr. rer. nat. Dieter Bothe
2. Gutachten: Prof. Dr.-Ing. Cameron Tropea



Numerical Study of Head-on Binary Droplet Collisions: Towards Predicting the Collision Outcomes

Genehmigte Dissertation von Muyuan Liu M.Sc. aus Hebei, China

1. Gutachten: Prof. Dr. rer. nat. Dieter Bothe
2. Gutachten: Prof. Dr.-Ing. Cameron Tropea

Tag der Einreichung: 31.05.2017
Tag der Prüfung: 17.10.2017

Darmstadt — D 17

Bitte zitieren Sie dieses Dokument als:
URN: urn:nbn:de:tuda-tuprints-70187
http://tuprints.ulb.tu-darmstadt.de/id/user/7018

Dieses Dokument wird bereitgestellt von tuprints,
E-Publishing-Service der TU Darmstadt
http://tuprints.ulb.tu-darmstadt.de
tuprints@ulb.tu-darmstadt.de

Die Veröffentlichung steht unter folgender Creative Commons Lizenz:
Namensnennung – Keine kommerzielle Nutzung – Keine Bearbeitung 4.0 International
http://creativecommons.org/licenses/by-nc-nd/4.0/



Acknowledgment
First of all, I would like to express my sincere gratitude to my supervisor, Prof. Dr. Dieter Bothe,

for introducing me to the world of science, for providing me with useful guidance throughout
my research and for the encouragement while I was hesitating. His kindness and friendly way
of mentoring have made life during my doctoral research more meaningful than research alone
would have. All the things I have learned from him, which are much more than I could have
ever expected at the start, are incredible assets.

I am grateful to Prof. Dr. Cameron Tropea for being my second supervisor and for correcting
my dissertation. His extensive teachings on fluid dynamics are the foundation of my research.

I owe all the former and present colleagues at Mathematical Modeling and Analysis my grati-
tude for their kind help which eased my boarding. We shared interesting and inspiring discus-
sions, many of which turned out to be important for the success of my research. I appreciate
the computational hours provided by the Hochschulrechenzentrum of the Darmstadt University
of Technology and the financial support provided by the German Research Foundation.

Finally, I want to thank my friends for their support and my parents for their love.

Darmstadt, May 2017 Muyuan Liu

i



Publications during the doctoral research

M. Liu, D. Bothe: Numerical study of head-on droplet collisions at high Weber numbers. Journal
of Fluid Mechanics, 289: 785-805, 2016.

C. Planchette, H. Hinterbichler, M. Liu, D. Bothe and G. Brenn: Colliding drops as coalescing and
fragmenting liquid springs. Journal of Fluid Mechanics, 814: 277-300, 2017.

L. Reitter, M. Liu, J. Breitenbach, K. L. Huang, D. Bothe, G. Brenn, K. L. Pan, I. V. Roisman and
C. Tropea: Experimental and computational investigation of binary drop collisions under elevated
pressure. In 28th Annual Conference on Liquid Atomization and Spray Systems, Valencia, Spain,
2017.

M. Liu, D. Bothe: Numerical study of binary droplet collisions at high Weber numbers. Presenta-
tion in ProcessNet Annual Conference, Bingen, Germany, 2016.

M. Liu, D. Bothe: Numerical study of head-on collisions of water droplets at high Weber numbers.
Poster for International Conference on Multiphase Flows, Florence, Italy, 2016.

ii Acknowledgment



Kurzfassung
Binäre Tropfenkollisionen spielen eine wichtige Rolle in der Natur und in vielen technischen

Sprayanwendungen. Die Modellierung der Kollisionsausgänge, nämlich Bouncing, Coalescence,
Tropfenzerfall nach vorläufiger Verschmelzung und Spatter (auch „Shattering “ und „Splashing“
genannt), bildet die Grundlage für die Untersuchung der Zerstäubungsprozesse auf größeren
Skalen. Das Ziel dieser Arbeit richtet sich auf die Entwicklung von numerischen Methoden,
welche für die Prädiktion der Kollisionsausgänge eingesetzt werden, sowie die numerische
Untersuchung derjenigen Phänomena in binären Tropfenkollisionen, welche die Kollisionsaus-
gänge entscheidend beeinflussen. Für die numerischen Simulationen wird der Inhouse-Code
Free Surface 3D (FS3D), der auf der Volume-of-Fluid (VOF) Methode basiert, eingesetzt. Die
numerischen Untersuchungen beschränken sich auf zentrale Kollisionen.

Spatter tritt bei Tropfenkollisionen mit hoher Kollisionsenergie auf, wobei eine dünne Prall-
lamelle entsteht und in standardmäßigen Simulationen unphysikalisch zerfällt. Um Spat-
ter simulieren zu können, wird ein verbesserter Lamellenstabilisierungsalgorithmus entwick-
elt und ausführlich validiert. Mit einem geeignet eingebrachten Verrauschen des Anfangs-
geschwindigkeitsfeldes wird die Instabilität am Rand des Stoßkomplexes ausgelöst und Spat-
ter in der Simulation erfolgreich reproduziert. Die Simulationsergebnisse stimmen sehr gut mit
den Experimenten überein. Basierend auf den Simulationsergebnissen wird die Entwicklung der
Randinstabilität als eine Verstärkung eines Signals durch ein Signalverstärkungssystem, das in
drei seriell verbundene Subsysteme unterteilt wird, betrachtet. Dabei wird festgestellt, dass die
Entwicklung der Randinstabilität in der linearen Phase der Randinstabilität durch die Rayleigh-
Plateau Instabilitätstheorie vorhergesagt werden kann. Der Einfluss der Tropfenviskosität wird
numerisch untersucht und es wird gezeigt, dass der Kollisionsausgang zu Spatter neigt, wenn
die Tropfenviskosität verkleinert wird. Diese Abhängigkeit nimmt während der Abnahme der
Tropfenviskosität ab. Die Tropfenviskosität beeinflusst die Entwicklung der Randinstabilität vor
allem, indem die Basisgeometrie des Randes verändert wird. Eine erfolgreiche Aufklärung des
Mechanismus der Randinstabilität bildet den Grundstein für die Vorhersage des Eintritts von
Spatter und die Vorhersage des Größenspektrums der sekundären Tropfen, die bei Spatter
entstehen. Die Untersuchung des Mechanismus der Randinstabilität im Kontext von binärer
Tropfenkollisionen hat eine generelle Bedeutung, weil der Ausstoß von Sekundärtropfen von
einem instabilen Rand beim Aufprall eines Tropfens auf eine feste Wand oder auf einen Flüs-
sigkeitsfilm ebenfalls auftritt.

Binäre Tropfenkollisionen führen bei relativ kleinen Weberzahlen entweder zu Bouncing oder
zu Coalescence als Kollisionsausgang. Die Simulationen von Bouncing und Coalescence wer-
den durch die Umschaltung der Randbedingungen an der Aufprallwand erfolgreich durchge-
führt. Die Simulationsergebnisse stimmen mit den Experimenten gut überein. Allerdings sind
diese Simulationen nicht prädiktiv, weil der Kollisionsausgang vorgegeben werden muss. Die

iii



Schwierigkeit der Prädiktion des Kollisionsausgangs Bouncing vs. Coalescence liegt darin, dass
der dünne Gasfilm zwischen den kollidierenden Tropfen nicht aufgelöst werden kann und ein
physikalisch sinnvolles Koaleszenzkriterium in der numerischen Methode fehlt. Um die prädik-
tive Simulation zu ermöglichen, wird ein Multiskalen-Simulationskonzept erarbeitet. Neben
dem Hauptlöser FS3D, der die Strömung auf der makroskopischen Skala löst, besteht das
Multiskalen-Simulationskonzept aus drei Anteilen: (1) Ein Subgridskalen(SGS)-Model wird
integriert in den Hauptlöser FS3D. (2) Coalescence wird numerisch unterdrückt, bevor ein
geeignetes Koaleszenzkriterium möglicherweise erfüllt ist. (3) Ein numerisches Koaleszenzkri-
terium wird angewandt.

Basierend auf der Lubrikationstheorie wird das SGS-Model hergeleitet, wobei der Effekt
verdünnter Strömung mit berücksichtigt wird. Das SGS-Model wird in FS3D implementiert
und ausführlich validiert. Zur Kopplung des SGS-Models wird der Druck im Gasfilm, der
durch das SGS-Model gelöst wird, als Druckrandbedingung auf der Aufprallebene aufgeprägt.
Wird die erste Schneidung der PLIC-Flächen mit der Kollisionsebene als Koaleszenzkriterium
angewendet, können die Simulationsergebnisse sowohl Bouncing als auch Coalescence sein.
Allerdings ist der Kollisionsausgang abhängig von der Gitterauflösung. Wird als Koaleszenzkri-
terium eine Gasfilmdicke von null im Rahmen der im Algorithmus verwendeten Toleranzen
angewendet, führen die Simulationen ausschließlich zu Bouncing. Es wird gezeigt, dass auch
weitere mögliche Korrekturen der Geschwindigkeiten, die für den Transport der Flüssigkeit-
sphase eingesetzt werden, nicht ausreichen, den Übergang zwischen Bouncing und Coalescence
vorherzusagen. Als Ausblick wird angedeutet, wie in zukünftigen Forschungsarbeiten z.B. mit
Hilfe der volumengemittelten Volume-of-Fluid (VA-VOF) Methode, die den Geschwindigkeit-
sunterschied innerhalb einer Rechenzelle berücksichtigt, die Genauigkeit des Transports der
Flüssigkeitsphase erhöht werden kann.

Mit Hilfe der Multiskalen-Simulation wird qualitativ gezeigt, dass der Kollisionsausgang bei
hoher Verdünnung in der Gasphase zu Coalescence neigt.



Abstract
Binary droplet collision plays an important role in nature and in many technical processes

involving sprays. The modeling of the collision outcomes, namely bouncing, coalescence,
separation after temporary coalescence, and spatter (also called ‘shattering’ and ‘splashing’),
establishes the basis for the investigation of the atomization processes on larger length scales.
The aim of this thesis is to develop numerical methods that are employed in the prediction of
the collision outcomes and the numerical investigation of the phenomena in binary droplet col-
lisions which affect the collision outcomes. The in-house code Free Surface 3D (FS3D), which
is based on the Volume of Fluid (VOF) method, is employed for the numerical simulations. The
numerical investigations are restricted to head-on collisions.

Spatter occurs at high energetic collisions, resulting in a thin liquid lamella that ruptures
artificially in standard numerical simulations. In order to simulate spatter, an improved lamella
stabilization algorithm has been developed and extensively validated. By means of properly
chosen white noise disturbances of the initial velocity field, the instability of the rim of the
collision complex is triggered and the spatter is successfully reproduced in the simulations. Very
good agreements between the simulation results and the experiments are achieved. Based on
the simulation results, the development of the rim instability is considered as an amplification of
disturbances via a signal amplification system that is subdivided into three sequential connected
subsystems. It is confirmed that the development of the rim instability in the linear phase of
the instability can be predicted by the Rayleigh-Plateau instability theory. The influence of the
droplet viscosity is studied numerically and it is shown that the collision outcome tends to be
spatter when the droplet viscosity is reduced. This dependency decreases with the decrease
of the droplet viscosity. The droplet viscosity influences the development of the rim instability
mainly through varying the geometrical evolution of the rim. A successful elucidation of the
mechanism of rim instability builds the foundation for the prediction of the occurrence of spatter
and the prediction of the size distribution of the secondary droplets arising in spatter. The
investigation of the mechanism of the rim instability in the context of binary droplet collisions
is of general importance because the ejection of secondary droplets from an unstable rim also
emerges in collisions of a droplet on a solid substrate or on a liquid film.

Binary droplet collisions result in bouncing or coalescence at relatively small Weber numbers.
The simulations of bouncing and coalescence have been successfully conducted by switching the
boundary conditions on the collision plane. The simulation results are in good agreement with
corresponding experiments. However, the simulations are not predictive because the collision
outcome must be specified in advance. The difficulty of the prediction of bouncing versus coa-
lescence lies in the fact that the thin gas film between the colliding droplets cannot be resolved
in feasible simulations and that a physically meaningful coalescence criterion is missing in the
numerical method. In order to facilitate the predictive simulation, a multi-scale simulation con-

v



cept has been developed. In addition to the main solver FS3D, which solves the flow on the
macroscopic scale, the multi-scale simulation concept consists of three parts: (1) A sub-grid-
scale (SGS) model is integrated within the main solver FS3D. (2) Coalescence is numerically
suppressed before a suitable coalescence criterion is contingently satisfied. (3) A numerical
coalescence criterion is applied.

Based on the lubrication theory, the SGS model is derived which accounts for the rarefied flow
effect. The SGS model is implemented in FS3D and extensively validated. For the integration
of the SGS model, the pressure in the gas film, which is solved by the SGS model, applies as
a pressure boundary condition on the collision plane. Employing the first intersection of PLIC-
surfaces with the collision plane as coalescence criterion, the collision outcome in the simulation
can be both bouncing and coalescence. The predicted collision outcome, however, depends on
the grid resolution. Employing zero gas film thickness (in algorithm tolerance) as coalescence
criterion, the simulations result only in bouncing. It is shown that various possible corrections of
the velocity field, which decides the transport of the liquid phase, have not led to a meaningful
prediction of the transition between coalescence and bouncing. Further developments, e.g.
the volume-averaged Volume of Fluid (VA-VOF) method, which takes into account the velocity
difference within a computational cell, shall be implemented in future work to increase the
accuracy of the transport of the fluid phase.

By means of the multi-scale simulation it is qualitatively shown that the collision outcome
tends to be coalescence at higher rarefaction in the gas phase.



Contents

1. Introduction 1
1.1. Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2. Fundamentals of binary droplet collisions . . . . . . . . . . . . . . . . . . . . . . . . . 2

1.2.1. Phenomenological description and literature review . . . . . . . . . . . . . . 5
1.3. Organization of the thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

2. Mathematical Modeling of Two-Phase Fluid Systems 11
2.1. Governing equations in incompressible two-phase fluid systems . . . . . . . . . . . 11

2.1.1. Mass conservation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.1.2. Transport equation of the phase indicator . . . . . . . . . . . . . . . . . . . . . 13
2.1.3. Momentum conservation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.1.4. Differential equations governing the flow in two-phase fluid systems . . . . 14

2.2. Interfacial instabilities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
2.2.1. Plateau-Rayleigh instability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
2.2.2. Rayleigh-Taylor instability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16

2.3. Continuum mechanics and rarefied gas flow regimes . . . . . . . . . . . . . . . . . . 17

3. Numerical Methods 19
3.1. Finite Volume discretization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
3.2. Volume of Fluid method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
3.3. Modeling of the surface tension force . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
3.4. Time discretization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
3.5. Typical setup for head-on collisions of two identical droplets . . . . . . . . . . . . . 27
3.6. Computation of the energy contributions . . . . . . . . . . . . . . . . . . . . . . . . . . 27

4. Lamella Stabilization 31
4.1. Stabilization of a liquid lamella . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31

4.1.1. Validation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
4.1.2. Lamella stabilization in terms of the balanced-CSF model . . . . . . . . . . . 34
4.1.3. Correction of the computation of the surface energy . . . . . . . . . . . . . . 35

4.2. Stabilization of a gas lamella . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
4.2.1. Dependence on the grid resolution . . . . . . . . . . . . . . . . . . . . . . . . . 38
4.2.2. Correction of the computation of the surface energy with a present under-

resolved gas lamella . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41

vii



5. Numerical Simulations of Spatter Phenomenon and the Mechanism of the Rim In-
stability 45
5.1. Domain adjustment . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
5.2. Numerical setups for water droplet collisions . . . . . . . . . . . . . . . . . . . . . . . 46

5.2.1. Role of imposed disturbances . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
5.3. Simulation results for water droplet collisions and quantification of the rim de-

velopment . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
5.3.1. Comparison with experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
5.3.2. Quantification of the rim development . . . . . . . . . . . . . . . . . . . . . . 51

5.4. Division of the collision process . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
5.5. The rim instability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57

5.5.1. Initial unstable phase: RT instability? . . . . . . . . . . . . . . . . . . . . . . . 58
5.5.2. Linear unstable phase: PR instability? . . . . . . . . . . . . . . . . . . . . . . . 59
5.5.3. The nonlinear unstable phase . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60

5.6. Influence of liquid viscosity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
5.6.1. Evolution of the rim geometry depending on viscosity . . . . . . . . . . . . . 61
5.6.2. Input signals and output signals in the linear unstable phase . . . . . . . . . 65
5.6.3. Growing rim as a signal amplification system . . . . . . . . . . . . . . . . . . 66
5.6.4. Conclusion on the viscosity effect . . . . . . . . . . . . . . . . . . . . . . . . . . 69

6. Flow in the Gas Layer between the Colliding Droplets 71

7. Multi-Scale Simulations aiming at Predictions of the Transition between Bouncing
and Coalescence 77
7.1. Derivation of the SGS model starting from the lubrication theory . . . . . . . . . . . 77

7.1.1. Analytical solutions of two idealized cases . . . . . . . . . . . . . . . . . . . . 81
7.1.2. Rarefied lubrication equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
7.1.3. Calibration of the slip coefficients . . . . . . . . . . . . . . . . . . . . . . . . . 86

7.2. Implementation of the SGS model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
7.2.1. Approximation of the quantities related to the interface . . . . . . . . . . . 87
7.2.2. Equation system . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88

7.3. Validation of the SGS model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
7.3.1. Validation through the scenario of two approaching solid spheres . . . . . . 89
7.3.2. Validation through a virtual scenario based on a retracing ellipsoid . . . . . 90

7.4. Integration of the SGS model into the main solver FS3D . . . . . . . . . . . . . . . . 93
7.5. Prediction of the collision outcome in terms of coalescence versus bouncing . . . . 94

7.5.1. Preliminary results of predictions . . . . . . . . . . . . . . . . . . . . . . . . . . 95
7.5.2. Modified f -transport . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96

7.6. Rarefaction effect . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98
7.7. Role of the compressibility . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100

viii Contents



8. Summary and Outlook 101
8.1. Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101
8.2. Outlook . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103

A. Simulations of Spatter Conducted at ITLR of the University of Stuttgart 105

B. Lubrication Equation in Cylindrical Coordinates 107
B.1. Transformation between Cartesian coordinates and polar coordinates . . . . . . . . 107
B.2. Lubrication equation in cylindrical coordinates . . . . . . . . . . . . . . . . . . . . . . 107
B.3. Pressure in the gas layer between two approaching parallel plates in cylindrical

form . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108

C. Setups 109

Bibliography 113

Contents ix





Nomenclature

Latin symbols Unit

n normal vector m

u velocity m/s

T stress tensor kg/(ms2)

g gravitation m/(s2)

t time s

x , y, z Cartesian coordinates m

f volume fraction -

r radius /radius coordinate m

R radius -

p pressure Pa

A area m2

V volume m3

L, l b, h length m

B offset of a collision m

E energy J

m mass kg

X impact parameter -

C Courant number -

m instability mode -

a amplitude -

k wave number -

Greek symbols Unit

λ mean free path m

ρ density kg/m3

κ curvature m−1

η dynamic viscosity Pa·s
σ surface tension N/m

Γ interface -

xi



∆ difference -

Φ extensive quantity -

φ density of an extensive quantity -

Ω domain -

δ Dirac operator -

ϕ density ratio -

ψ viscosity ratio -

χ phase indicator -

ω growth rate -

Indexes

0 initial value

a value on interface

l liquid phase

g gas phase

± phase

b ambient

Abbreviations

Re Reynolds number
We Weber number
Oh Ohnesorge number
Kn Knudsen number
DNS Direct Numerical Simulation
FV Finite Volume
CV Control Volume
CFL Courant-Friedrichs-Lewy
VOF Volume of Fluid
FS3D Free Surface 3D
SGS Sub-Grid-Scale
VdW Van der Waals force

DFG German research foundation
TRR Transregio
ITLR Institut für Thermodynamik der Luft- und Raumfahrt
MMA Mathematical Modeling and Analysis

xii Nomenclature



1 Introduction

1.1 Motivation

Droplet collisions are very common phenomena observed in nature: some people even enjoy
staring at the splash of droplets resulting from heavy rains hitting the ground, lost in thought.
The earliest study of droplet collisions dates back to Lord Rayleigh [48] who noted that small
rain droplets bounce upon collision with a pool of water rather than dive into it directly. The
majority of early works on the collision of two droplets dating from the 1960’s focused on water
droplets essentially due to droplet growth relating to precipitation. Gunn [21] pointed out that
binary collisions of water droplets result in four typical types of outcome: (I) bouncing, (II)
coalescence, (III) drop disruption after coalescence, (IV) drop spatter; see Figure 1.1 for a first
impression. He also noted that approximately seven collisions occur for every kilometer of free
fall in typical heavy rain, which is frequent enough to influence the size of rain droplets.

In technical applications, droplet collision is an elementary process in the liquid atomization
processes that convert bulk fluid into a spray composed of a dispersion of small droplets. An
example of a spray is shown in Figure 1.2. Fluid atomization processes and sprays are of impor-
tance in many industrial processes, for example in fuel injection in internal combustions, spray
painting, food processing, pharmaceuticals, environmental protection, and many others. The
main purpose of the atomization process is to increase the contact area between the gas phase
and the liquid phase. An increase in the gas-liquid contact area in a spray augments momen-
tum, heat, and mass-transfer between the gas phase and the liquid phase as well as the related
chemical processes like the combustion of fuel droplets. Apart from the primary atomization,
in which the supplied liquid is demolished into a cloud of droplets, the collision of the disinte-

Figure 1.1.: Water droplet collisions imaged by Gunn [21]. From left to right: bouncing, coales-
cence, disruption after coalescence and drop spatter.

1



Figure 1.2.: Diesel spray imaged by Wang et al. [67]

grated droplets also affects the droplet size spectrum in a spray through the different collision
outcomes.

It is difficult to investigate the whole spay process by means of Direct Numerical Simulations
(DNS), which resolve all the relevant length scales and time scales of a physical problem, due
to the large span of the length scale: for example in typical diesel combustion engines, the
droplet size is in the region of 1− 10µm while the length of the injected fuel spray is 2− 3 cm
[4]. Generally, computational models for simulating sprays are based on Sub-Grid-Scale(SGS)
models which include the averaged effect of the unsolved scales. The large-scale behavior of the
spray is then described by reduced mathematical models, rather than detailed solution of the
Navier-Stokes equations [44]. An accurate prediction of the characteristics of a spray requires
that the SGS model accounts for the effect of droplet collisions in terms of the collision outcomes
due to the collective effect on the flow property. The derivation of models for determining the
collision outcomes is usually based on experimental observations and measurements, which are
restricted by the spatial and temporal resolutions of the obtained photographs recording the
collision process of two droplets. In principle, Direct Numerical Simulations of the collision pro-
cess can provide information that is not accessible by experimental measurements, for example
the detailed local flow properties within a droplet, which facilitates the modeling of the colli-
sion rules. As long as the numerical simulation is able to reproduce the collision process, more
strictly saying predictively, spectacular phenomena like the spatter of droplets can be studied by
means of numerical ‘experiments’.

1.2 Fundamentals of binary droplet collisions

The collision of two droplets, say, the binary droplet collision, is a very complex phenomenon
in practice. The participating droplets may have different sizes and/or different materials which
are miscible or immiscible. The possible complex rheological behavior of the liquid material
makes the study of the collision process even more difficult. In real systems far away from
thermodynamic equilibrium like within diesel combustion chambers, the Marangoni-effect, i.e.
the variation of the surface tension along the droplet surface, the evaporation and burning of
the droplets would also affect the collision process significantly. Looking beyond all of these
complexities, the study in this work is restricted to the collision of two Newtonian droplets of
the same size and the same material under isothermal conditions without chemical reactions.

2 1. Introduction



In the simplified collision system consisting of two identical droplets and environmental gas,
which is schematically described in Figure 1.3, the dimensional analysis yields that the collision
process is determined by the following five dimensionless parameters [38]:

• Weber number

We=
ρl U

2
r D0

σ
(1.1)

• Reynolds number

Re=
ρl Ur D0

ηl
(1.2)

• Density ratio

ψ=
ρl

ρg
(1.3)

• Viscosity ratio

ϕ =
ηl

ηg
(1.4)

• Impact parameter

X =
B
D0

(1.5)

with the parameters in the collision system listed as follows:

ρl density of the liquid phase, kg ·m3

ρg density of the gas phase, kg ·m3

ηl dynamic viscosity of the liquid phase, Pa · s
ηg dynamic viscosity of the gas phase, Pa · s
U0 relative velocity of the droplets, m/s

D0 initial diameter of the droplets, m

σ surface tension, N ·m
B offset of the collision, m

The Weber number describes the relative importance of the droplets’ inertia compared to
the surface tension and the Reynolds number describes the relative importance of the droplets’
inertia compared to the viscous effect. The impact parameter is a measure of the offset of a
collision. A collision is called a head-on collision if the impact parameter is zero.

Depending on the system parameters, several different collision outcomes can occur. The
typical collision processes are illustrated in Figure 1.4. Following the convention in previous
studies [45] [16], the collision outcomes are summarized in the collision diagram, see Figure
1.5. Note that the collision diagram in Figure 1.5 is a qualitative description and will differ with

1.2. Fundamentals of binary droplet collisions 3



D0

D0

Ur

B

Figure 1.3.: Collision of two identical droplets.

Figure 1.4.: Typical collision processes imaged by Jiang et al. [25] (first to fifth column) and by
Roth et al. [56] (sixth column). From left to right: coalescence (sector I), bouncing
(sector II), coalescence (sector III), near head-on separation (sector IV), off-center
separation (sector V) and spatter (sector VI).

4 1. Introduction



We

X

(I
)

C
o
a
le

sc
en

ce
(II)Bouncing

(III)Coalescence

(IV)Near head-on separation

(V)Off-center separation

(VI) Spatter

Figure 1.5.: Collision diagram.

different Re, ψ and ϕ, since each boundary curve between two sectors of different collision
outcomes, without considering other effects, is described as

X = X (We, Re,ψ,ϕ). (1.6)

The Reynolds number in equation (1.6) can be substituted by the Ohnesorge number defined
as

Oh=
ηl

p

D0 ·ρl ·σ
=
p

We
Re

. (1.7)

The advantage of using the Ohnesorge number instead of Reynolds number is that the collision
outcome is then only dependent on We and X , as long as the materials in the collision system
and the droplet size are fixed. Another dimensionless parameter used for the description of the
temporal evolution of a collision process is the dimensionless time defined as

t∗ = t · Ur/D0. (1.8)

The time instant when the distance between the sphere center of the droplets is D0, i.e. when
the droplets are just touching each other, is defined as t∗ = 0.

1.2.1 Phenomenological description and literature review

Most studies of binary droplet collisions have been driven by the prediction of the collision
outcome due to its important applications in the study of systems on larger length scales. The
literature review is subdivided according to the transitions between the collision outcomes,
along with the description of the significant phenomena related to the collision process in each
sector in the collision diagram. Sectors I to VI refer to the sectors I to VI in the collision diagram
in the present thesis.

1.2. Fundamentals of binary droplet collisions 5



Transition between coalescence and bouncing

In relatively low Weber number regime (sector I, II and III), the collision outcome undergoes a
non-monotonic coalescence-bouncing-coalescence transition. In sector I, the droplets approach
each other with very small inertia and coalesce into a single larger droplet. The droplets move
so slowly that the gas phase hardly resists the approach of the droplets and the coalescence
will occur as long as the distance between the colliding surfaces gets so small that the force
of attraction between the molecules of the two droplets, i.e. the Van der Waals force, becomes
dominant. The critical distance for the onset of the Van der Waals force is in the order of 10 nm
[32]. In sector II, the droplets collide with higher inertia resulting in a high pressure in the
gas phase between the colliding droplets due to the viscous resistance against the drainage of
the gas film. The high pressure in the gas phase prevents a further approach of the droplets.
As a result, the liquid in the droplets flows outward near the collision plane and the colliding
surfaces are enlarged, forming a gas layer in-between. In the next stage, the droplets recede due
to the surface tension force and bounce apart from each other without material exchange. If the
inertia is further elevated, the colliding surfaces will approach each other more closely. When
the local gas layer thickness is small enough for the Van der Waals force to become effective, the
droplets coalesce (in sector III). For a simplification of the discussion, the gas phase between
the colliding droplets at very low Weber number regime (in sector I) is also called a gas layer,
though its resisting effect is almost nil. Based on this phenomenological description, it is the
competition between the resisting effect of the gas layer and the intermolecular forces that
ultimately determines whether the droplets merge into a single droplet or bounce apart. This
phenomenological description is partly based on the work of Zhang and Law [71].

Zhang and Law [71] pointed out that before the possible coalescence of the droplets, the dis-
tance of the colliding surfaces must necessarily pass through a region where it is comparable
with the mean free path of the gas molecules. Aiming at the prediction of the non-monotonic
coalescence-bouncing-coalescence transition for head-on collisions and assuming Poiseuille flow
between two expanding flattened interfaces in the gas film, Zhang et al. [71] established an an-
alytical model that accounts for the motion and deformation of the droplets and the interaction
between the droplets and gas. The resisting effect of the gas layer acts in the form of a pressure
force acting on the droplet’s surface. The pressure in the gas layer is derived based on the gas
drainage, the rarefaction effect and the Van der Waals force. This model is able to predict the
non-monotonic coalescence-bouncing-coalescence transition, though the accuracy of the predic-
tion of the critical Weber numbers especially between the sector II and III needs to be further
improved. Estrade and Biscos [15] conducted an experimental investigation of binary collision
of ethanol droplets and developed a theoretical model for predicting the critical Weber number
between sector II and III which derivation is based on the energy conservation. However, the
influence of the gas phase is not included in their model, although it plays an important role in
the transition between coalescence and bouncing: Qian and Law [45] observed that the sector
II (bouncing) expands with increasing ambient pressure.

6 1. Introduction



Numerical simulation results concerning the transition between coalescence and bouncing
are rare. Most numerical studies focused on the reproduction of the collision process with pre-
scribed collision outcome and/or detailed description of the flow field, see [38, 36, 37, 41, 40].
The difficulty for a predictive simulation in terms of the collision outcome is that the gas layer
becomes so thin that it cannot be resolved by feasible simulations. Trying to tackle this problem,
Mason et al. [33] conducted a multi-scale simulation: the macro-scale simulation is conducted
by means of the Volume of Fluid (VOF) method, while the flow in the gas layer is modeled
by means of the classical lubrication theory. A critical gas layer thickness of 40 nm is used as
a coalescence criterion. The simulation can reproduce the coalescence process well, while the
transition between coalescence and bouncing is not mentioned. Murad and Law [35] conducted
molecular simulations yielding both bouncing and coalescence predictively. The droplet diam-
eter in the simulation is restricted to a few nanometers which is much too small compared to
the typical droplet sizes in technical applications. Resolving the droplet collision process with
molecular simulations in technical applications is not realistic, since within one cubic centimeter
air under standard conditions there are already roughly 1019 molecules, the solution of which
is far beyond the capability of modern high-performance computing technology.

Transition between coalescence and separation

In near head-on collisions in sector III and IV, the coalesced collision complex expands first out-
wards after the onset of the coalescence and then retracts due to the surface tension force tem-
porarily forming a cigar-like shape. The cigar-like shape collision complex breaks into two main
droplets with contingent secondary droplets at relatively high inertia of the colliding droplets
(sector IV). At relatively smaller Weber numbers in sector III, the collision complex stays as a
single droplet after an oscillation, through which a part of the kinetic energy is dissipated. At
relatively high inertia and high impact parameter in sector V, the collision complex falls apart
after the collision complex is elongated due to the inertia and centrifugal force.

Many previous studies have been devoted to modeling the boundary curves between sector
III and IV and between III and V. Ashgriz and Poo [5] developed correlations for the boundary
curves III-IV and III-IV based on the inviscid assumption. Qian and Law [45] showed that the
critical Weber number between the sector III and IV is linearly dependent on the Ohnesorge
number in case of head-on collisions. Jiang et al. [25] proposed a model for the boundary
curve III-V which accounts for the viscous effect. Sommerfeld et al. [62] developed a composite
universal model which is an extension and combination of the model of Ashgriz and Poo [5]
and the model of Jiang et al. [25]. This model is able to predict the boundary curves III-
IV and III-V in a large span of viscosity variation. Planchette et al. [42] deduced the critical
velocity distinguishing coalescence and separation in head-on collisions based on an energy
balance of the expansion phase and the retracting phase of the collision process combined with
a Rayleigh-like fragmentation criterion.

The numerical studies have been able to reproduce the boundary curves between coalescence
and separation. Rieber and Frohn [52] predict the boundary curves by means of a Volume of
Fluid (VOF) code called Free Surface 3D, which is also applied in the numerical studies in this

1.2. Fundamentals of binary droplet collisions 7



work. These curves are in good agreement with corresponding experiments of Brenn and Frohn
[11] and Ashgriz and Poo [5]. Pan and Suga [41] reproduced the collision outcomes in section
III, IV and V by means of the Level-Set method, which locates the interfacial position by means
of transporting a signed distance function, without giving the collision outcome in advance.
Saroka et al. [57], using the VOF method, found that the boundary between sector III and IV
in case of head-on collision depends on Reynolds number and the dependence decreases with
increasing Reynolds number. One issue arising in the simulation is that the liquid lamella of
the collision complex, which emerges in the middle of the collision complex in the expanding
phase after coalescence of the colliding droplets at relatively high Weber numbers, ruptures in
the numerical simulations [41, 37, 17, 18]. The rupture of the liquid lamella deteriorates the
accuracy while predicting the collision outcomes. Focke and Bothe [17] and Focke and Bothe
[18] dealt with this problem by correcting the computation of the surface tension force in the
lamella region. The advantage of using numerical methods in order to study the transition
between coalescence and separation is that it is very easy to separately modify each of the
parameters in the collision system, which facilitates a comprehensive and systematic study of
the dependency of the boundary curves on the system parameters and the validation of the
above mentioned theoretical models predicting the boundary curves III-IV and III-V.

Spatter in sector VI

Spatter is a significantly different phenomenon emerging at collisions with large kinetic en-
ergy. Roth et al. [56] conducted experimental studies on head-on binary collisions of iso-
propanol droplets at Weber numbers ranging from 1030 to 2876. Their results show that with
increasing Weber number, the rim of the collision complex becomes increasingly unstable. If
the Weber number is high enough, secondary droplets are spattered out from the rim of the
collision complex. They also conducted numerical simulations but only qualitatively and at rel-
atively low Weber numbers due to the high requirement on the computational effort resulting
from high collision energies. Despite the lower Weber number in their numerical simulations,
the instability of the rim is too strong compared to the experimental results. Pan et al. [39]
conducted an experimental study of head-on binary droplet collisions at Weber numbers up to
about 5000. They presented clear images for the deformation history of water collision com-
plexes. Their results show that the instability of the rim of the water collision complex emerges
at much smaller Weber numbers compared to iso-propanol droplets. Kuan et al. [27] conducted
a numerical study of head-on collisions of water droplets at high Weber numbers, employing a
parallel, adaptive interface tracking method, and compared their results with the experimental
work of Pan et al. [39]. Their simulations capture the unstable rim of the collision complex
and show good agreement with corresponding experimental results of Pan et al. [39] up to
We = 442. However, their comparisons of the collision complex shapes for We > 442 are ei-
ther incomplete regarding the length of the physical simulation time or do not show a good
agreement.

Understanding the spatter phenomenon, especially the mechanism of the rim instability of the
collision complex, is not only a scientific desire but also a precondition for predicting the bound-

8 1. Introduction



ary curve between the section V and VI. It is even more desirable to predict the size spectrum of
secondary droplets ejected from the unstable rim due to its important role in determining the
characteristics of a spray. Moreover, the study of the mechanism of the rim instability is of more
general importance, since a similar phenomenon of spattering of secondary droplets from an
unstable rim emerges also in collisions of a droplet on a solid substrate or on a liquid film.

1.3 Organization of the thesis

The present thesis originates within the framework of the Transregio-75 (TRR-75) of the Ger-
man Research Foundation (DFG) “Droplet Dynamics Under Extreme Ambient Conditions”. The
driving force of this thesis is the realization of the prediction of the collision outcome in terms
of bouncing versus coalescence and in terms of the onset of spatter as well as the prediction of
the size spectrum of the secondary droplets resulting from the spatter phenomenon.

In order to seek the final goal, namely the prediction of the collision outcomes, the following
sub-goals have been set:

• Explanation of the rim instability that is present on the spatter phenomenon due to its
importance in the onset of spatter and in the size spectrum of the secondary droplets.

• Development of numerical methods that enable the reproduction of the spatter phe-
nomenon.

• Development of numerical methods that enable the reproduction of the bouncing and
coalescence phenomena, at first with prescribed collision outcome.

• Acquisition of the flow properties in the gas flow between the colliding droplets before
possible coalescence.

• Development of a SGS model that solves the flow in the gas layer and that is integrated
into the main solver FS3D.

• Numerical simulation by means of the integrated SGS model and the corresponding as-
sessment in terms of the capability of predicting bouncing versus coalescence.

The developments and the numerical investigations are restricted to the head-on collisions.
The structure of this work is as follows. In Chapter 2, the governing equations of two-phase
incompressible flows are presented. Two most important interfacial instabilities and the con-
tinuum hypothesis as well as its transition to the rarefaction effect are addressed. The basic
numerical methods used in this work are introduced in Chapter 3. The models for the compu-
tation of the surface tension force are discussed in terms of their advantages and limitations.
Chapter 4 describes the problem of the artificial interaction while computing the surface ten-
sion force of a thin lamella, which can be either a liquid lamella in high energetic collisions or
a gas film at lower collision energy, and delivers numerical methods that deal with this prob-
lem. In Chapter 5, the spatter phenomenon is in depth investigated based on the numerical
results. One highlight is the rigorous theoretical/computational confirmation of the mechanism

1.3. Organization of the thesis 9



of the rim instability in the context of binary droplet collisions. In Chapter 6, the results of
a 2D simulation with very high grid resolution focusing on the flow in the gas layer between
the colliding droplets are presented. The derivation, implementation and validation of a SGS
model, which solves the flow in the gas layer, are presented in Chapter 7. Preliminary results
on predictive simulations by means of multi-scale simulations are discussed. At the end of this
thesis, a summary and an outlook are given.

10 1. Introduction



2 Mathematical Modeling of Two-Phase
Fluid Systems

In continuum mechanics, a two-phase fluid system is characterized by an interface separating
two immiscible substances. Due to the jump of material properties over the interface, the differ-
ential equations describing the two-phase fluid system differ from those for a single phase flow.
These equations are given first in this chapter. Then, two kinds of interfacial instabilities are
introduced, which serve as a basis for further investigation of the instabilities occurring in the
context of binary droplet collisions. At the end of this chapter, the continuum hypothesis and its
transition to the rarefaction effect are addressed.

2.1 Governing equations in incompressible two-phase fluid systems

We consider a material control volume V (t) moving with the flow, see Figure 2.1. The control
volume possesses a normal vector n pointing outwards and contains two phases separated by a
free interface Γ (t). The two phases are denoted by Ω+(t) and Ω−(t) and the normal vector of
the interface nΓ points to Ω+(t). Furthermore, the intersection of Γ (t) with the surface of V (t)
is bounded by a curve C(t) with a normal vector N pointing outwards and lying tangential to
Γ (t).

Assuming no-slip condition between the tangential velocity components at the interface, i.e.

u+tan = u−tan = uΓ tan, (2.1)

the conservation equation for an extensive quantity Φ with density φ is given as [9]:

d
d t

∫

V (t)

φdV =−
∫

∂ V (t)

jmol · ndA+

∫

V (t)

PdV

−
∫

C(t)

jmol
Γ ·Nds+

∫

V (t)∩Γ (t)

PΓ dA.

(2.2)

The density φ can be either a scalar or a vector. The first term on the right-hand side of equation
(2.2) describes the change of φ through transport. Since V (t) is a material control volume,
only diffusive (molecular) transport comes into play. The third term describes the change of
φ by surface diffusion through the curve C(t). The second and the fourth terms are volume-
specific and surface-specific source terms, respectively. Given φ, P, PΓ , jmol and jmol

Γ properly,
the conservation equations for mass and momentum can be derived.

11



Figure 2.1.: Description of a two-phase fluid system.

2.1.1 Mass conservation

Mass is an extensive quantity with ρ denoting its density. We assume that the interface does
not carry mass, hence jmol

Γ = 0. With φ = ρ, jmol = 0, P = 0 and PΓ = 0, the mass conservation
equations are given as

∂ ρ

∂ t
+∇ · (ρu) = 0 in Ω+ ∪Ω−, (2.3)

[[ρ(u− uΓ )]] · nΓ = 0 at Γ (t), (2.4)

where uΓ denotes the interface velocity and u the velocity in the bulk. The notation [[·]] denotes
the jump of a quantity across the interface, which is defined as

[[Ψ]](x) = lim
h→0+

(Ψ(x+ hnΓ (x))−Ψ(x− hnΓ (x))). (2.5)

For incompressible fluids of constant density without phase change, equations (2.3) and (2.4)
can be further reduced to

∇ · u= 0 in Ω+ ∪Ω−, (2.6)

u+ · nΓ = uΓ · nΓ = u− · nΓ at Γ (t). (2.7)

The equation (2.6) is the corresponding volume conservation equation for incompressible
flow. The equation (2.7) describes that the velocity component perpendicular to the interface is
continuous. With the assumption (2.1), the jump condition of the velocity field is summarized
as

[[u]] = 0 at Γ (t). (2.8)

12 2. Mathematical Modeling of Two-Phase Fluid Systems



2.1.2 Transport equation of the phase indicator

We define a phase indicator χ with value 1 in Ω+ and 0 in Ω−. Then the local density and
dynamic viscosity are given as

ρ = χρ+ + (1−χ)ρ−, (2.9)

η= χη+ + (1−χ)η−. (2.10)

For incompressible fluids of constant density without phase change, the transport equation of
the phase indicator is given as

∂ χ

∂ t
+∇ · (χu) = 0 in Ω+ ∪Ω−. (2.11)

The transport equation (2.11) is the central equation solved in the VOF method, which details
are described in Chapter 3.

2.1.3 Momentum conservation

Given φ = ρu, jmol = S− pI, jmol
Γ = σ(I− nΓ ⊗ nΓ ), P = ρf, PΓ = 0 for incompressible flows

without phase change [51], the momentum conservation equations are written as

∂

∂ t
(ρu) +∇ · (ρu⊗ u) =∇ · T+ρf in Ω+ ∪Ω−, (2.12)

[[− T]] · nΓ = σκnΓ at Γ (t), (2.13)

whereσ is the surface tension assumed to be constant, κ is the interface curvature, and f denotes
the body force. For Newtonian fluids, the stress tensor T is given as

T= −pI+ S, (2.14)

with S= η(∇u+(∇u)T) in incompressible flow. The momentum equations (2.12) for Newtonian
fluids are called the Navier-Stokes equations, which are valid in both phases.

For a two-phase fluid system at rest, i.e. a static droplet or bubble in a surrounding fluid, the
momentum jump condition can be further reduced to the Young-Laplace equation, giving the
Laplace pressure jump across the interface:

∆p = σ(
1
R1
+

1
R2
), (2.15)

where R1 and R2 are the curvature radii of the interface.

2.1. Governing equations in incompressible two-phase fluid systems 13



Figure 2.2.: Breakup of a water jet due to the capillary effect [66].

2.1.4 Differential equations governing the flow in two-phase fluid systems

Since the surface tension force only acts at the interface, the momentum jump condition
can be added to the momentum conservation equation as a singular interface term by using a
δ-function yielding the one-field formulation of the Navier-Stokes equations. The differential
equations governing the hydrodynamics of incompressible two-phase flows, including the one-
field formulation of the Navier-Stokes equations, are summarized below:

∇ · u= 0, (2.16)
∂ ρ

∂ t
+∇ · (ρu) = 0, (2.17)

∂

∂ t
(ρu) +∇ · (ρu⊗ u) = −∇p+∇ · S+ρf+σκnΓδ. (2.18)

2.2 Interfacial instabilities

In two-phase incompressible flows, the interface between the two phases is unstable in cer-
tain conditions, meaning that a small disturbance on the interface will be magnified by the
fluid system. Two kinds of interfacial instabilities are addressed in this section, namely the
Plateau-Rayleigh (PR) instability and the Rayleigh-Taylor (RT) instability, which are the candi-
date theories for the explanation of the rim instability emerging in spatter phenomenon that
will be discussed in Chapter 5.

2.2.1 Plateau-Rayleigh instability

A straight free liquid jet in a gaseous surrounding is unstable due to the capillary effect and
breaks into a chain of small droplets; see Figure 2.2 for an example. This kind of instability is
referred to as the Plateau-Rayleigh (PR) instability.

This phenomenon has been explained analytically by Lord Rayleigh [46], assuming an inviscid
liquid. The instability of the jet begins with a tiny perturbation that is in form of, e.g. , surface
displacement, pressure or velocity fluctuations in the fluid system. The disturbances in the axial
direction of the jet are magnified with wave numbers less than a cut-off wave number kc, and are
decaying otherwise. The wave number with the maximum growth rate decides the number and
the size of the resultant droplets. In order to go into more detail, we consider a harmonic axial

14 2. Mathematical Modeling of Two-Phase Fluid Systems



a

r0

R2

R1

λ

Figure 2.3.: Schematic representation of an unstable jet.

disturbance on a straight fluid jet with initial radius r0 and an infinite length (see Figure 2.3)
and conduct a dimensional analysis of the Navier-Stokes equations. The disturbance is described
by wavelength λ, frequency ω̄, and amplitude a. According to Ashgriz and Nasser [4], the
characteristic length and the characteristic time of the problem are λ and 1/ω̄, respectively. The
characteristic fluid velocity due to the interface motion is estimated as U ≈ aω̄. Therefore, the
time derivative term, the convective term and the Laplacian term of the Navier-Stokes equations
are estimated as ρ∂ u/∂ t ∼ ρUω̄∼ ρaω̄2,∇·(ρu⊗u)∼ ρU2/λ∼ ρa2ω̄2/λ, η∆u∼ ηU/λ2 ∼
ηaω̄/λ2, accordingly. Comparing the time derivative and the convective term shows that the
non-linear convective term can be neglected if a � λ. The Reynolds number is given as Re =
ρUλ/η≈ ρaλω̄/η. Comparing the time derivative and the viscous term shows that the viscous
term can be neglected when ρλ2ω̄/η= Reλ/a� 1.

Neglecting the viscous effect and linearizing the Navier-Stokes equations written in cylindrical
coordinates, the growth rate of a small-amplitude harmonic disturbance is given as (for detailed
derivation see [4] or the original work [46])

ω2 =
σk
ρr2

0

(1− k2r2
0)

I1(kr0)
I0(kr0)

, (2.19)

where ω = ωr + iωi is the growth rate, k is the wave number, In are modified Bessel functions
of the first kind. If the real part of the growth rate ωr is positive, the disturbance grows expo-
nentially in time. In the parentheses in equation (2.19), k2r2

0 comes from the jet axial curvature
and ‘1’ from the jet sectional curvature which represents the capillary pinching. Since the ra-
tio of the modified Bessel functions I1/I0 is positive for all conditions, the jet is unstable if the
capillary pinching is dominant:

1> k2r2
0 , (2.20)

or equivalently:

λ > 2πr0. (2.21)

The cut-off wave number for the instability is therefore given as kc = 1/r0.The spectrum of
the growth rate for an inviscid jet is plotted in Figure 2.4. The fastest growing perturbation
dominating the jet evolution occurs at kr0 = 0.697, which can be used to estimate the break
length of the fluid jet and, furthermore, the size of the generated droplets.

2.2. Interfacial instabilities 15



kr0

0 0.2 0.4 0.6 0.8 1

ω
/
(σ

/
ρ
r3 0
)0

.5

0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

Inviscid

Oh=0.01

Oh=0.1

Oh=1

Figure 2.4.: Non-dimensional growth rate of instabilities on inviscid and viscous jets with respect
to the wave number.

Weber [68] and Chandrasekhar [12] extended Rayleigh’s inviscid theory concerning viscosity
effect. The corresponding characteristic equation for the perturbation spectrum is given as

ω2 +
2ηk2

ρI0(kr0)

�

I ′1(kr0)−
2kl

k2 + l2

I1(kr0)
I1(l r0)

I ′1(l r0)
�

ω

=
σk
ρr2

0

(1− k2r2
0)

I1(kr0)
I0(kr0)

l2 − k2

l2 + k2
,

(2.22)

where l2 = k2 + ρω/η. For η = 0, this expression reduces to equation (2.19). The influence
of viscosity on the growth rate is shown in Figure 2.4 regarding Ohnesorge numbers defined
as Oh = η/pr0ρσ. One sees that the maximum growth rate and the dominant wavenumber
decrease with increasing viscosity, while the cut-off wave number remains kc = 1/r0.

2.2.2 Rayleigh-Taylor instability

The Rayleigh-Taylor (RT) instability is an instability of interface between two fluids with dif-
ferent densities. It occurs when the lighter fluid supports the heavier fluid under gravity or
equivalently, under an acceleration of the fluid system in the direction towards the denser fluid.
Figure 2.5, obtained from own numerical simulations, illustrates the development of the RT
instability under gravity.

The classical modeling of the Rayleigh-Taylor instability given by Taylor [63] assumes two
inviscid fluids with density ρ1 and ρ2 originally at rest with ρ1 > ρ2. The gravitational field g
points to the lighter fluid. Derived on the basis of potential theory, the growth rate of a small
harmonic perturbation at the interface is given as

ω(k) =
p

kgA , (2.23)

16 2. Mathematical Modeling of Two-Phase Fluid Systems



Figure 2.5.: Development of the Rayleigh-Taylor instability under gravity.

where A= (ρ1 −ρ2)/(ρ1 +ρ2). With A> 0, the instability exists for every positive wavenumber
k. For large wave numbers, it is natural to expect that the surface tension force will stabilize
the instability, which is indeed the case. According to Bellman et al. [6], who extended the RT
instability theory concerning the effect of the surface tension force, the interface perturbation is
dampened for

k <

√

√(ρ1 −ρ2)g
σ

. (2.24)

The fastest growth rate is given as

ωmax =

√

√ 2A

3
p

3σ
g
Æ

(ρ1 −ρ2)g, (2.25)

occurring at wavenumber

k =
1
p

3σ

Æ

(ρ1 −ρ2)g. (2.26)

2.3 Continuum mechanics and rarefied gas flow regimes

It is well known that gas is composed of discrete molecules or atoms on the microscopic scale.
If the length scale of a given volume is large enough, the average property of the gas is not
influenced by the concrete number of molecules in this volume. For example, there exist 3×1010

molecules in 10−9 cm3 air under standard temperature and pressure conditions ensuring the
average sense of the gas property [60]. On the basis of this conception, continuum mechanics
assumes that the substance of the object completely fills the space it occupies and the property
at a point is continuous (even twice continuously differentiable) function of coordinates and
time allowing the description of the motion of the continuum by partial differential equations.

2.3. Continuum mechanics and rarefied gas flow regimes 17



We introduce the Knudsen number for the gas medium defined as

Kn=
λ

L
, (2.27)

where λ is the mean free path of the gas and L is the characteristic length scale of the flow.
The continuum hypothesis is valid for very small Knudsen numbers, i.e. Kn� 1. For a gas flow
with the mean free path of the gas comparable to the characteristic length scale of the flow, i.e.
Kn > 0.01, the flow is regarded as rarefied. According to the rarefaction degree, the rarefied
gas flow can be divided into three regimes [65]:

0.01<Kn< 0.1 Slip flow regime,

0.1<Kn< 10 Transitional regime,

Kn> 10 Free molecular regime.

In the free molecular regime, the collisions between molecules and the collision of the molecules
with the surface of a body have to be taken into account, which makes the problem more
complicated and put it beyond the scope of the present work. In the slip flow regime, the
continuum hypothesis can still be used; but it is necessary to introduce some modifications to
the boundary conditions. For particular cases, it is also possible to extend the method used in
the slip flow regime to the transitional regime. This issue will be further addressed in Chapter
7.

18 2. Mathematical Modeling of Two-Phase Fluid Systems



3 Numerical Methods

3.1 Finite Volume discretization

The numerical simulations in this thesis are performed with the Computational Fluid Dynamic
(CFD) program Free Surface 3D (FS3D), which was originally developed at the ITLR (University
of Stuttgart) [51] and has been extended both at the ITLR and in the Mathematical Modeling
and Analysis (MMA) group of Technische Universität Darmstadt. FS3D is based on the Finite
Volume (FV) method. The basic concept of the Finite Volume method is described by dealing
with the conservation equation for an extensive quantity Φ with density φ in a fixed control
volume:

d
d t

∫

V

φdV =

∫

V

PdV +

∫

S

F · ndS, (3.1)

where P and F summarize all the source terms and the fluxes, respectively. Discretizing the
domain into finite control volumes in cuboidal form, applying the above equation for each
control volume Ω with surfaces l = 1, ..., 6 and integrating over time step δt = tn+1 − tn yield

∫

Ω

φ(tn+1)dV −
∫

Ω

φ(tn)dV =

∫

Ω

tn+1
∫

tn

Pd tdV +
l=6
∑

l=1

∫

Sl

tn+1
∫

tn

Fd t · nl dS. (3.2)

Approximation of the spatial average values by the values of the center of the control volumes
or by the values of the center of the cell faces constructs the foundation of the Finite Volume
method. For unsteady problems, the time averaged values are approximated by the values at
the new or the last time step.

In FS3D, a cuboidal box serving as the computational domain is discretized into cuboidal
cells. Specifically within this work, the edge lengths of the computational cells are restricted
to be equidistant. Two layers of dummy cells are settled around the domain, which are used
to employ the boundary conditions as suggested by Blazek [7]. The variables are arranged in
a staggered way according to the MAC scheme [23], whereby the pressure and other scalar
variables like the volume fraction are stored in cell centers and the velocity components are
stored at cell faces separately. This staggered grid arrangement prevents an oscillating pressure
field [58].

19



x

y

1 1 0.7 0

1 1 0.5 0

1 0.9 0.2 0

0.7 0.3 0 0

n

i− 1
2 i i+ 1

2

δtu
i+1

2
δtu

i− 1
2

Figure 3.1.: Schematic description of the PLIC algorithm. Left: Reconstruction of the Interface.
Right: Computation of the flux through cell faces. The volumes in the dark area are
the advected volumes of fluid in one time step. The velocities in x-direction at cell
faces are positive in this example.

3.2 Volume of Fluid method

Within this thesis, the Volume of Fluid method, originally developed by Hirt and Nichols [23],
is employed for locating and tracking the interface position. In the VOF method, the transport
equation

∂ χ

∂ t
+∇ · (χu) = 0 (3.3)

is solved for the phase indicator χ, which has the value one in the liquid phase and zero other-
wise. Integrating the phase indicator function over the volume of a computational cell and then
dividing the integral by the cell volume gives the volume fraction of the liquid phase in finite
volume discretizations. The volume fraction is denoted by f and has the value







1− ε < f ≤ 1 in the liquid phase,
0≤ f < ε in the gas phase,
ε≤ f ≤ 1− ε in cells containing interface,

(3.4)

where the numerical threshold ε is set to 10−6 in FS3D.
The transport equation (3.3) is solved by means of the Piecewise Linear Interface

(Re-)Construction (PLIC) algorithm [50]. First, the interface is reconstructed by fulfilling two
conditions simultaneously:

1. The interface is locally approximated by a plane possessing a normal vector n given as

n = −
∇ f
‖ ∇ f ‖

. (3.5)

20 3. Numerical Methods



For computing ∇ f , the components of ∇ f are evaluated at cell faces by central differenc-
ing. The obtained values are then averaged to obtain∇ f located at cell vertices. Averaging
∇ f located at all vertices of one cell yields an approximation of ∇ f which is located at the
cell center. In total, a 27 (33) cells stencil is used for computing ∇ f , so as for computing
n.

2. The approximated interface and the wetted cell faces enclose the volume f∆x∆y∆z.

The reconstructed interface is exemplarily shown in Figure 3.1. Instead of solving the volume
transport equation (3.3) directly, this advection equation is split into three one-dimensional
transport equations that are solved successively. Furthermore, the volume transport equation
is complemented with a divergence corrector, which ensures that the cells in the interior of the
liquid phase are not confused with interface cells. The flux of the liquid phase through cell faces
for one time step, e.g. in x-direction, is then the volume of the liquid phase located in the volume
uδt∆y∆z; see Figure 3.1 for a schematic description in 2D. The procedure of the computation
of the volume flux in each direction is the same; the sequence of computing the three split
equations is symmetric between two subsequent time steps, which guarantees a second-order
procedure. After one of the three one-dimensional transports, the f -values could be larger than
one or smaller than zero in badly conditioned cells. The f -values are restricted between zero
and one through a post-correction. For more details, we refer to [51].

Having the updated volume fraction in each cell, the local material properties can be approx-
imated:

ρ = f ρL + (1− f )ρG, (3.6)

η= f ηL + (1− f )ηG, (3.7)

where the subscripts L and G denote the liquid phase and the gas phase, respectively. With given
volume fractions, the so-called height function can be computed, which describes the distance
of the interface to a given reference surface. The height function h can be written based on an
elementary volume between the interface and a projected area serving as reference:

h= lim
dA→0

∫∫

hdA

dA
, (3.8)

where
∫∫

hdA is the volume of fluid between the interface and the projected area dA. It is
possible to find the height function in three space directions. Taking z direction for example,
one can approximate

∫∫

hdA by summing up the volume fractions in a cell column having the
bottom area ∆x ·∆y:

V =
n
∑

k=1

( f (k) ·∆z(k)) ·∆x ·∆y, (3.9)

3.2. Volume of Fluid method 21



where k indicates the cell and runs from the cell lying on the reference surface to the cell
containing only gas phase in z direction. The height function can then be computed:

h=
V

∆x ·∆y
=

n
∑

k=1

f (k) ·∆z(k). (3.10)

The equation (3.10) computes the height function based on the liquid phase and the reference
surface must be within the liquid phase. Reversing the f -values, the height function with respect
to a reference surface that is located in the gas phase can be computed with equation (3.10) as
well.

3.3 Modeling of the surface tension force

In order to employ the one-field formulation of the Navier-Stokes equations, the surface ten-
sion force fΓ has to be approximated. In the context of the Volume of Fluid method, the following
approximations are performed:

δ ≈ ‖ ∇ f̃ ‖, (3.11)

ñ ≈−
∇ f̃

‖ ∇ f̃ ‖
, (3.12)

κ≈∇ · (−ñ), (3.13)

yielding the Continuum-Surface-Force (CSF) model [10] for computing the surface tension
force:

fΓ = σκδn ≈ σκ∇ f̃ . (3.14)

The surface normal ñ and the delta function δ are computed on the basis of a smoothed
volume fraction field f̃ . The smoothing is conducted by means of a second-order B-Spline
function. As a result, the surface tension force is computed not only on the interface but in a
small region around the interface, acting as a volumetric force.

The CSF model suffers from the well-known parasitic currents: the velocity field around the
interface oscillates unphysically. As a result, a spherical droplet with zero initial kinetic energy in
resting environment, which is expected to be silent, will oscillate and move in an unpredictable
way [2]. The parasitic currents can be significantly suppressed by discretizing∇ f instead of∇ f̃
in the same manner as discretizing the pressure gradient∇p, since the surface tension force can
then be balanced by the pressure jump due to the compatible discretization [8]. This balanced
version of the discretization of ∇ f leads to the balanced-CSF model for computing the surface
tension force.

Aside from the discretization of ∇ f , it is necessary to compute the interface curvature accu-
rately in the balanced-CSF model. The well-known variants of the balanced-CSF model are from
Renardy and Renardy [49], Francois et al. [20] and Popinet [43], whereby the main difference

22 3. Numerical Methods



n

Figure 3.2.: Left: Computation of the interface curvature by means of height functions in direc-
tion of the maximum component of the interface normal. Right: A fragment of two
approaching droplets as an example scenario.

lies in the estimation of the interface curvature. Renardy and Renardy [49] estimate the in-
terface curvature through computing the curvature of a fitted paraboloid obtained on the basis
of volume conservation. Francois et al. [20] estimate the interface curvature based on height
functions constructed in a fixed region around the interface cells. In this thesis, the implementa-
tion of Popinet [43] is adopted. The algorithm estimates the interface curvature through three
hierarchical methods. The interface curvature is preferably computed by means of height func-
tions constructed in cell columns that align with the direction of the maximum component of
the interface normal; see a 2D example in Figure 3.2. The curvature in 2D can be determined
by means of the height functions:

κ=
hx x

(1+ h2
x)(3/2)

. (3.15)

If not all necessary height functions for computing the surface curvature can be found, the
construction of the height function is conducted in the other space directions. In case the com-
putation of the surface curvature by means of height functions is not possible, the computation
is conducted with a fallback strategy, namely paraboloid fitting of the interface by means of
height functions found in all the three space directions. If this is still not possible, the fitting is
conducted by means of the barycenters of the PLIC surfaces. The fitting in the fallback strategies
is conducted in a 27 (33, and 32 in 2D) cells stencil. The total algorithm is rather sophisticated
and the details can be found in the original work of Popinet [43].

It should be noted that in the cell columns where the volume fractions are summed for the
computation of the height functions, an empty cell and a full cell must be identified, since the
summation of the volume fractions is conducted between them. This leads to a limitation of the
application of the implemented balanced-CSF method especially in simulating binary droplet

3.3. Modeling of the surface tension force 23



Figure 3.3.: Unphysical interface oscillation near contact area in a binary droplet collision. Three
symmetry planes are employed in the simulation.

collisions, as the construction of a proper height function is not possible when a change in
topology, e.g. the coalescence of droplets, is supposed to occur. Take the contact area schemati-
cally shown in Figure 3.2 (right) for example, a necessary empty cell for constructing the height
function in direction perpendicular to the collision plane cannot be found; in tangential direc-
tion, a full cell can not be found. For this scenario, the fallback strategies will be employed in
the simulation; however, only unpredictable results will be obtained, since no meaningful fitting
is possible for computing the interface curvature in a 32 cells stencil in this 2D scenario. In Fig-
ure 3.3, a simulation result of a binary droplet collision shows an unphysical oscillation of the
interface due to the failure while computing the interface curvature in the contact area. Another
limitation of the balanced-CSF model is that it requires higher grid resolution for computing the
curvature of sub-structures such as the rim of the collision complex, which is very thin at the
early stage of a high energetic collision. The disadvantages of the balanced-CSF model restrict
its use for certain scenario.

Another model for computing the surface tension force is the CSS model [28] which computes
the surface tension force as

fΓ =∇ · T , (3.16)

with

T = σ ‖ ∇ f̃ ‖ (I − ñ ⊗ ñ). (3.17)

The CSS model is a conservative model since it is written as a divergence of a tensor of
second-order. The CSS model suffers from parasitic currents as well; however, it does not have
the problems of the balanced-CSF model described above. In addition, the parasitic currents
have relatively less effect at highly dynamic process, which extends the area of applicability of
the CSS model. The comparison between the CSS model and the balanced-CSF model is sum-
marized in Table 3.1. In this thesis, the decision of employing the CSS model, the balanced-CSF
model or, if needed, a combination of the two models is made based on the advantages and the
limitations of the two models in concrete application scenarios. Without explicit indication in

24 3. Numerical Methods



advantages limitations

CSS model 1. stable for small sub-structure

2. stable with topology change

3. less problematic in highly dy-
namic processes, since the par-
asitic currents are small com-
pared to the velocity field near
the interface

1. parasitic currents

balanced-CSF
model

1. much smaller parasitic cur-
rents

1. high grid resolution required for
small sub-structures

2. unstable with topology change

Table 3.1.: Comparison between the CSS model and the balanced-CSF model.

this thesis, the CSS model is employed. The unbalanced CSF model is not used in this thesis.

3.4 Time discretization

The time discretization is conducted by means of a first-order explicit Euler scheme, where
the key semi-discrete equations are given as:

∇ · un+1 =0, (3.18)

ρn+1un+1 −ρnun

δt
=−∇ · [(ρu⊗ u)]n

+
ρn+1

ρ( f n+1)
[−∇pn+1 +∇ · S(µn+1,un) + f n+1], (3.19)

ρn+1 −ρn

δt
=−∇ · (ρu)n. (3.20)

The density ρ( f n+1), the viscosity µn+1 and the body force f n+1 in equation (3.19) are com-
puted using the volume fraction at tn+1, which is updated by the VOF method. The semi-discrete
momentum equation (3.19) can be rewritten for the velocity at the new time as:

un+1 = ũ−
δt

ρ( f n+1)
∇pn+1, (3.21)

3.4. Time discretization 25



with

ũ=
1
ρn+1

[ρnun −δt∇ · [(ρu⊗ u)]n]+

δt
ρ( f n+1)

[∇ · S(µn+1,un) + f n+1].
(3.22)

The provisional velocity ũ is updated by the convective and the non-convective acceleration
that correspond to the first and second term on the right-hand side of equation (3.22), respec-
tively. The convective acceleration and the mass conservation equation (3.20) are computed
simultaneously, whereby an operator split procedure and a divergence corrector are employed.
Inserting equation (3.21) together with the updated velocity field ũ in equation (3.18) yields
the discrete pressure-Poisson equation:

∇ · [
1

ρ( f n+1)
∇pn+1] =

∇ · ũ
δt

. (3.23)

Complemented by homogeneous Neumann conditions or Dirichlet conditions, the pressure-
Poisson equation is solved by means of a multigrid equation solver. The computation of the
velocity field un+1 is then finished with equation (3.21).

Due to the explicit time discretization, the time step is limited according to the following three
criteria to ensure a stable numerical scheme [52] :

Courant-Friedrichs-Lewy condition:

δt ≤
∆x

max{|u|, |v |, |w|}
(3.24)

Constraint due to surface tension force:

δt ≤

√

√(ρl +ρG)(∆x)3

4πσ
(3.25)

Constraint due to fluid viscosity:

δt ≤
(∆x)2

2max{ηL
ρL

, ηL
ρL
}

(3.26)

The constraints due to surface tension force and due to fluid viscosity are constant during a
simulation. Noticing equation (3.24), the time step in simulations can be dynamically controlled
by limiting the maximum allowable Courant-number:

C =
max{|u|, |v |, |w|} ·δt

∆x
. (3.27)

26 3. Numerical Methods



Figure 3.4.: Typical setup for head-on collisions of two identical droplets by means of three sym-
metry planes.

3.5 Typical setup for head-on collisions of two identical droplets

The typical setup for head-on collisions of two identical droplets, which are the main object
of this study, is described in this section. The simulation is initialized with a quarter of a droplet
lying on two symmetry planes. The droplet collides with an initial velocity towards a symmetry
plane, which is equivalent to the head-on collision of two identical droplets, see Figure 3.4. The
cuboidal domain is decomposed into equidistant cuboidal cells. On the symmetry planes the
slip conditions are prescribed; the homogeneous Neumann boundary conditions for the velocity
and zero pressure are prescribed on all the other boundary planes (outer boundaries).

The setup described in this section serves as a template, from which the setups of almost all
of the simulations conducted in this work are derived.

3.6 Computation of the energy contributions

For an isolated fluid system without gravity, the total amount of the kinetic energy KE plus the
surface energy SE is decreasing with the dissipation rate ΦE. The different contributions can be
calculated based on the local flow field. The local kinetic energy in a computational cell is given
as

ke =

∫

Cell

1
2
ρ||u||2dV ≈

1
2
ρ||u||2∆x∆y∆z, (3.28)

3.5. Typical setup for head-on collisions of two identical droplets 27



where the averaged square of the absolute value of the velocity is computed on the basis of
assuming a linear velocity profile within a cell, which ensures a second order accuracy. The
local kinetic energy is then given as

ke =
1
6
ρ(u2

i− 1
2 , j,k
+ u2

i+ 1
2 , j,k
+ ui− 1

2 , j,k · ui+ 1
2 , j,k+

v 2
i, j− 1

2 ,k
+ v 2

i, j+ 1
2 ,k
+ vi, j− 1

2 ,k · vi, j+ 1
2 ,k+

w2
i, j,k− 1

2
+w2

i, j,k+ 1
2
+wi, j,k− 1

2
·wi, j,k+ 1

2
)∆x∆y∆z,

(3.29)

where i, j and k denote the index of cells and ±1
2 denote the cell faces. A graphic example for

the notations is found in Figure 3.6.
The surface energy within a cell is given as

se = σA= σ

∫

Cell

||∇ f ||dV ≈ σ||∇ f ||∆x∆y∆z, (3.30)

where A is the surface area of the interface. ∇ f is computed by means of a 27 cells stencil, as
it is described in section 3.2. The values of ke and se are summed up over all the computational
cells to obtain KE and SE, respectively.

Alternatively, the interface area can also be computed by adding the areas of all PLIC-surfaces.
Rieber [51] states that computing the surface area by means of summing ‖ ∇ f ‖ (as suggested
in equation (3.30)) is a better choice, because it corresponds closely to the discretization of the
surface tension force. This statement is verified through a simulation, whereby a quarter of a
droplet lying at the intersection of three symmetry planes is initialized with an initial velocity.
After dissipation of the kinetic energy, the droplet is expected to retain its original form, i.e.
original surface area. The comparison between the surface area evolutions computed by means
of the two above mentioned methods does support the statement made by Rieber [51]; see
Figure 3.5. However, the deviation between the two profiles shown in Figure 3.5 is not very
large, especially when the oscillation is still relatively strong.

The dissipation rate for a cell is given as

φe = 2η

∫

Cell

D : D dV ≈ 2ηD : D∆x∆y∆z, (3.31)

where D = 1
2(∇u+∇uT ) is the deformation tensor of the flow field. The scalar product A : B

between tensors is defined as A : B =
∑

i, j ai, j bi, j. For computing the local energy dissipation,
the velocity gradients have to be evaluated. Taking the velocity component u as an example, the
velocity gradients are given using central differencing:

∂ u
∂ x
=

ui+ 1
2 , j,k − ui− 1

2 , j,k

∆x
, (3.32)

∂ u
∂ y
=

ui− 1
2 , j+1,k + ui+ 1

2 , j+1,k − ui− 1
2 , j−1,k − ui+ 1

2 , j−1,k

4∆y
, (3.33)

∂ u
∂ z
=

ui− 1
2 , j,k+1 + ui+ 1

2 , j,k+1 − ui− 1
2 , j,k−1 − ui+ 1

2 , j,k−1

4∆z
. (3.34)

28 3. Numerical Methods



t∗
0 2 4 6 8 10

A
/A

0

0.98

1

1.02

1.04

1.06

1.08

1.1

1.12

1.14

CSF

PLIC

Figure 3.5.: Comparison between the surface area evolutions computed by means of ‖ ∇ f ‖ (de-
noted in the Figure by CSF due to their similar discretization) and PLIC-surfaces. The
balanced-CSF model is employed in the simulation. The surface area is normalized
by the original surface area given analytically. The time is normalized by the time
when the droplet finishes the first oscillation in the sense that the solid line reaches
its minimum. The setup of the simulation is listed in Table C.1 in appendix C.

ui+ 1
2 ,j−1,k

ui+ 1
2 ,j+1,k

ui+ 1
2 ,j,k

(i, j, k)

(i, j − 1, k)

(i, j + 1, k)

ui− 1
2 ,j+1,k

ui− 1
2 ,j−1,k

ui− 1
2 ,j,k

Figure 3.6.: Computation of ∂ u
∂ y by means of a 3 cells stencil. The velocities are interpolated from

cell faces to the cell centers. The differencing is then conducted over 2∆y .

3.6. Computation of the energy contributions 29



It should be noted that the computation of ∂ u
∂ y and ∂ u

∂ z is conducted in a 3 cells stencil, see the
schematic illustration in Figure 3.6. The dissipated energy can be computed by summing the
local dissipation rate over all cells and integrating the obtained total dissipation rate ΦE over
time.

30 3. Numerical Methods



4 Lamella Stabilization
In direct numerical simulations of binary droplet collisions using the VOF method, a fluid

lamella arises and ruptures during a droplet collision process under certain conditions, which
leads to an unphysical topology change and thus strong deviation of the collision complex shape
in the successive deformation process. This Lamella can be a liquid lamella in droplet collisions
at high Weber numbers (We > 100) and a gas lamella (the gas lamella is also referred to as the
gas layer or gas film in this thesis) in droplet collisions at low Weber numbers (We = O(1) ∼
O(10)), the context of which will be described in detail in this chapter.

4.1 Stabilization of a liquid lamella

Head-on collisions of two identical droplets at high Weber numbers lead to the formation
of an extremely thin liquid lamella, which artificially ruptures in numerical simulations even
at relatively fine grid resolutions. The rupture of the thin lamella leads to a deviation of the
collision complex shape as illustrated exemplarily by Figure 4.1. Similar lamella ruptures are
also noted in many other numerical works [41, 37, 17, 16]. In experimental studies of binary
droplet collisions, this kind of lamella rupture is not observed even when the Weber number is
increased to We = 2876 [56], which confirms that the mentioned rupture of a liquid lamella is
artificial and must be prevented in numerical studies for capturing sound physics.

According to the work of Focke and Bothe [17] and Focke [16] in particular, the interaction
while computing the surface tension forces of both sides of the liquid lamella leads to the ar-
tificial lamella rupture. With respect to the CSS model in FS3D, the surface tension force is
computed from the divergence of the surface stress tensor, which is computed from a smoothed
volume fraction field f̃ . Both the smoothing of the volume fraction and the calculation of the
divergence of the surface stress tensor are carried out in stencils of 27 cells (33 cells). Altogether
the surface tension force in each cell containing a surface is affected by volume fractions in a
stencil of 125 cells (53 cells) [17]. When the lamella is resolved by fewer than four cells, the

Figure 4.1.: Lamella rupture (left) vs. stabilized lamella (right). Reproduced from [30] with per-
mission.

31



block of 125 cells contains interface cells from the opposite side of the lamella, which causes an
unphysical interaction that leads to the lamella rupture.

In addition, the surface reconstruction is affected by the artificial interface interaction as well.
The surface reconstruction is unambiguously determined by the f -field and the surface normal
computed by equation (3.5). The computation of the surface normal is conducted in a stencil
of 27 cells. If this contains additional interface cells from the opposite side of the lamella, the
computation of the surface normal is affected and, hence, the surface reconstruction as well.

Both the problem in computing the surface tension force and the problem with the surface
reconstruction can be prevented by identifying the lamella region and then treating the cells of
the opposite side of the lamella as ‘fully wetted’. Focke and Bothe [17] implemented the lamella
stabilization method for head-on collisions of two identical droplets based on the collision of a
droplet on a symmetry plane, which is equivalent to the head-on droplet collision. The boundary
conditions of a symmetry plane are incorporated using two layers of dummy cells adjacent to
the computational domain [7]. Focke and Bothe [17] sort out the lamella region from the rim
region of the collision complex by detecting the direct neighboring cells of an interface cell (we
denote it here as a ‘target’ cell) in the first layer of the computational grid on the symmetry
plane. If not every direct neighboring cell contains liquid, the target cell is defined as ‘rim’.
Target cells whose neighboring cells are all also interface cells, are defined as ‘lamella’. The
volume fraction is then set to f = 1 in the dummy cells in the lamella region, thereby avoiding
the above-mentioned unphysical interaction. It turns out that, in this way, the lamella is well
stabilized. However, the rim region is not optimally identified since only the direct neighboring
cells of a target cell are detected, which may lead to the confusion of the rim region with the
lamella region.

In order to correctly identify the rim region, some modifications are introduced to the work
of Focke and Bothe [17]. The implementation of this modified lamella stabilization method is
based on a collision complex on a symmetry plane as well. The procedure is as follows (see
Figure 4.2 for a schematic diagram):

1. The region is denoted as ‘empty’ if the volume fraction in the first layer of computational
cells of this region equals 0 ( f (i, j, 1) = 0). The region is denoted as ‘full’ if the volume
fraction in both the first and second layer equals 1 ( f (i, j, 1) = 1 and f (i, j, 2) = 1).

2. The algorithm detects all the neighboring cells of the remaining cells in the first layer as
follows: if one of the volume fractions f (i+3, j, 1), f (i+3, j+3, 1), f (i, j+3, 1), f (i−3, j+
3, 1), f (i−3, j, 1), f (i−3, j−3, 1), f (i, j−3, 1), f (i+3, j−3, 1) equals 0, the region is denoted
as ‘rim’; the remaining region is set to ‘lamella’. The detection distance ‘3’ is chosen, since it
is large enough to avoid the confusion of the rim region with the lamella region addressed
above and small enough to avoid reaching the neighboring finger-structures in spatter
phenomenon while detecting.

3. The volume fractions of dummy cells in the lamella region are set to 1 and are mirrored in
the remaining regions.

32 4. Lamella Stabilization



Figure 4.2.: Schematic of lamella stabilization. Reproduced from [30] with permission.

4. The surface tension force, the gradient of f and, hence, the surface normal due to equation
(3.5) are computed from this modified f -field.

5. The modification of the f -field in dummy cells is discarded.

To summarize, the stabilization algorithm for the liquid lamella corrects the computation
of the surface tension force and the reconstruction of the PLIC-surfaces for the stabilization
of the liquid lamella by modifying the f -values in dummy cells temporarily. Meanwhile, the
computation of ∇ f so as the computation of ‖ ∇ f ‖ are also corrected. In addition, the lamella
stabilization of Focke and Bothe [17] keeps the f -value in the lamella from falling below 10−6

during the transport of the f -field. Focke and Bothe [18] extended the work of Focke and Bothe
[17] for simulating off-center droplet collision without a symmetry plane by adding liquid into
the lamella to stabilize the lamella. In contrast, no mass needs to be added in the present work
in order to keep the lamella from rupturing.

4.1.1 Validation

In order to validate the new implementation of the lamella stabilization, the numerical results
are compared with the experimental work of Roth et al. [55]. The fluid is iso-propanol. Figure
4.3 shows that the shapes of the collision complex resulting from the simulation and from the
experiment are in good agreement. The slight corrugation at the rim of the collision complex at
the later stage of the deformation (around 400µs) is also predicted by the simulation.

The developments of the collision complex diameter D for two cases with different Weber
numbers We = 268 and We = 518 are presented in Figure 4.4. In the simulations, D is defined
as the average diameter on the collision plane and is computed based on the barycenter of
the PLIC-surfaces; in the experiment, the measurement of D is indicated in the experimental
part of Figure 4.3 (denoted there as d). The diameter of the collision complex and the time
are presented in a dimensionless form defined by D∗ = D/D0 and t∗ = t · Ur/D0, respectively.
In the case of We = 268 (the same case as presented in Figure 4.3), the development of D
predicted by the simulation is in very good agreement with the experimental results. In the case
of We= 518, the curves for the diameter of the collision complex converge with increasing grid
resolution. However, the limit curve lies above the corresponding experimental results. This

4.1. Stabilization of a liquid lamella 33



Figure 4.3.: Comparison of the predicted shapes of the collision complex (right) with the experi-
mental results (left) of Roth et al. [55]. The two columns of the collision complexes
in the experiment result from double exposure. See Roth et al. [55] for more de-
tails about the experimental setups. The Weber number and Reynolds number are
We = 268, Re = 426. The setup of the simulation is listed in Table C.2 in appendix
C. Post-processing is done by means of the open source ray tracer program POV-Ray
(http://www.povray.org). Reproduced from [30] with permission.

deviation is attributed to both the numerical and experimental inaccuracies. In the experiment,
errors arise and are accumulated measuring the droplets’ material and kinetic parameters. For
example, assuming that the velocity in the experiment might be smaller than measured, an
additional simulation with a smaller Weber number (We = 448, the relative velocity is 7%
smaller than the original) is performed, which achieves a very good accordance. In total, these
comparisons show the ability of the lamella stabilization method in predicting the development
of the collision complex.

4.1.2 Lamella stabilization in terms of the balanced-CSF model

In the lamella region that is not well resolved, as it is exemplarily shown in Figure 4.2, it
is not possible to compute the interface curvature in a meaningful way with the algorithm of
the balanced-CSF model described in Chapter 3, since neither meaningful height functions nor
meaningful fittings can be found. As a result, the surface tension force cannot be computed
correctly in the lamella region and the lamella collapses. In order to compute the interface
curvature in the lamella region accurately though meaningful height functions, it is necessary to
find a full cell and an empty cell in a cell column that aligns with the direction of the maximum
component of the surface normal. By using the lamella stabilization algorithm, the absent
full cells in the lamella region can be artificially guaranteed for computing accurate interface
curvatures and thus the lamella can be stabilized. In spite of its functionality, the balanced-CSF

34 4. Lamella Stabilization



t · Ur/D0

0 5 10 15

D
/D

0

0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

5

5.5

Experiment We=518

Mesh 384x384x128 We=518

Mesh 288x288x96 We=518

Mesh 192x192x64 We=518

Mesh 384x384x128 We=448

Experiment We=268

Mesh 288x288x96 We=268

Figure 4.4.: Comparison of the predicted diameter of the collision complex with experimental
results. The case We= 268 is the same case as presented in Figure 4.3. The Reynolds
number of the case We = 518 is Re = 600. The Reynolds number of the case We =
448 is Re= 558. The setups of the corresponding simulations are listed in Tables C.2,
C.3 and C.4 in appendix C. Reproduced from [30] with permission.

model is not used together with the stabilization algorithm within this thesis, since it cannot
deal with the topology change at binary droplet collisions in the contact area. Although it is
possible to switch the model for computing the surface tension force after a topology change
from the CSS model to the balanced CSF model, the complexity is not compensated because
parasitic currents do not play an important role in droplet collisions at high Weber numbers.

4.1.3 Correction of the computation of the surface energy

The lamella stabilization algorithm corrects the computation of ‖ ∇ f ‖ for an under-resolved
lamella and therefore also corrects the computation of the surface energy in the domain com-
puted by means of summing ‖ ∇ f ‖ (CSF-like approach). Nonetheless, one further step of
correction has to be taken for a valid computation of the integral surface energy by means of
the CSF-like approach, i.e. by equation (3.30). This correction is conducted by extending the
region for summing ‖ ∇ f ‖ from the domain to the layer of dummy cells adjacent to the symme-
try plane based on the fact that the value of ‖ ∇ f ‖ is not zero in dummy cells when the lamella
is under-resolved, as it is exemplarily shown in Figure 4.5. This correction is only employed in
the identified lamella region.

The diagram in Figure 4.6 shows the corrected surface energy evolution of a droplet collision
process compared to the uncorrected one. The collision process can be roughly divided into
an expanding phase, when the surface energy increases due to the expansion of the collision
complex and a receding phase, when the surface energy decreases due to the retraction of the
collision complex. At around t∗ = 4, the deviation between the original result and the corrected

4.1. Stabilization of a liquid lamella 35



f = 0.2 f = 0.2 f = 0.2

domain

symmetry plane

dummy cells

(a)

‖ ∇f ‖= 0.1 0.1 0.1

0.5 0.5 0.5

0.4 0.4 0.4

domain

symmetry plane

dummy cells

(b)

Figure 4.5.: A fragment of an under-resolved lamella which is ideally planar. The volume frac-
tions are constant in tangential direction. The cell width is 1. (a) f -values are set to
1 by the lamella stabilization algorithm in dummy cells. (b) summing ‖ ∇ f ‖ in the
extended region yields the surface area exactly.

one begins to become more obvious, since the lamella is then increasingly under-resolved. At
around t∗ = 14.8 the deviation vanishes since the lamella vanishes due to the retraction of
the collision complex. The curve representing the energy evolution computed by summing the
area of the PLIC-surfaces provides an additional evidence for the necessity of correcting the
computation of the surface energy conducted by means of summing ‖ ∇ f ‖.

4.2 Stabilization of a gas lamella

Similar to the liquid lamella of a collision complex, a gas lamella between two liquid-gas inter-
faces also ruptures in VOF-simulations conducted with the CSS model as long as the interfaces
from both sides of the gas lamella are located in a 53 cells stencil due to the artificial interaction
while computing the surface tension force. In binary droplet collisions, this kind of gas lamella
rupture always leads to coalescence of droplets. We consider two droplets as numerically coa-
lesced if the liquid phase allows to connect the droplet centers without leaving the liquid. In
the low Weber number regime (We = O(1) ∼ O(10)), binary droplet collisions result in either
coalescence or bouncing due to the presence of the gas layer between the colliding droplets. In
standard VOF simulations conducted with feasible resolutions, the thickness of the gas layer is
on the sub-grid scale. Thus, the simulations cannot reproduce bouncing due to the unphysical
gas lamella rupture.

Making use of the symmetry plane in the context of head-on binary droplet collisions, the gas
lamella can be simply stabilized by setting the volume fractions in dummy cells to zero while
computing the surface tension force, as it is schematically shown in Figure 4.7. In this way,
the computation of the surface tension force of one droplet is not affected by the other one

36 4. Lamella Stabilization



t · Ur/D0

0 2 4 6 8 10 12 14 16

S
E
/S

E
0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

5

5.5

CSF (original)

CSF (corrected)

PLIC

Figure 4.6.: Surface energy evolution computed by means of the CSF-like approach with and
without correction and by means of summing the area of the PLIC-surfaces in a
binary droplet collision process run with We = 416, Re = 530. The setup of the
simulation is listed in appendix C.5. The surface energy is normalized with the initial
surface energy. The time is normalized with D0/Ur .

.

Figure 4.7.: Stabilization of a gas lamella through a modification of volume fraction in dummy
cells on the symmetry plane. Left: symmetry condition. Right: stabilization of the
gas lamella.

and the gas layer is therefore prevented from rupture. By means of this lamella stabilization
algorithm through modification of the boundary condition of the volume fractions, the bouncing
phenomenon can be reproduced in simulations. However, the coalescence of droplets is then
not possible, since the two droplets do not affect each other at all. In fact, more physics, i.e.
the molecular forces as well as the rarefied flow effect have to be involved for modeling the
coalescence process, as it is discussed in the introduction. Despite the lack of additional physics
in the numerical modeling on the microscopic scale, a prescribed collision outcome ‘bouncing’
or ‘coalescence’ can be simulated by using standard and/or modified boundary conditions of the
f -field. The simulation methods and corresponding results are described in the following.

4.2. Stabilization of a gas lamella 37



Imposed bouncing

The gas lamella stabilization can be applied to both the CSS model and the balanced-CSF
model, in order to impose bouncing. It has been described in Chapter 3.3 that the height
function cannot be found in direction perpendicular to the collision plane near the impact in
head-on droplet collisions due to the lack of necessary empty cells. This problem is however
not present in imposed bouncing, since the the lacking empty cells near the impact can then be
found in dummy cells while applying the gas lamella stabilization. In all simulations of imposed
bouncing in this work, the balanced CSF model is employed to reduce the parasitic currents
that have a prominent effect at low Weber numbers. The comparison between simulation result
and experiment conducted by Pan et. al [40] shows that the evolution of the collision complex
shapes with the outcome bouncing can be well predicted by the simulation; see Figure 4.8.

Coalescence

In experiments, the coalescence of the colliding droplets (especially in sector III of the colli-
sion diagram) does not occur immediately when the droplets are supposed to get in touch (at
t∗ = 0) due to the presence of the gas layer. Instead, the coalescence is delayed to a later time
after significant deformation of the droplets. It has been suggested by Pan et al. [40] that the
abrupt smoothing of the cuspy contour of the collision complex indicates the occurrence of co-
alescence. Based on this approach, Pan et al. [40] studied the coalescence process by removing
the interface that is modeled by the Front Tracking method at a prescribed time obtained from
the experimental observation. In their experimental results shown in Figure 4.9, the coalescence
instant is identified to be between t = 0.366 ms and t = 0.370 ms. Using the interpolated value
t = 0.368 ms as the time of coalescence, the macroscopic phenomena of the collision process
are well reproduced in their simulations. In order to simulate the coalescence phenomenon, the
idea of specifying coalescence at a prescribed time is followed. The simulation is started with
the balanced-CSF method due to its advantage in reducing the parasitic currents. Meanwhile,
the gas lamella is stabilized. At the time of coalescence, the lamella stabilization is discarded by
switching the boundary condition of the volume fractions in dummy cells to standard boundary
condition at a symmetry plane. Meanwhile, the method for computing the surface tension force
is switched to the CS