FRIEDRICH-ALEXANDER-UNIVERSITÄT ERLANGEN-NÜRNBERG Lehrstuhl für Angewandte Mathematik III
Simulation of Instabilities occurring in Liquid Jets
Master Thesis
Author: Hurmat ul Ain born on 03.08.1980 in Karachi, Pakistan
Supervisors: Prof. Dr. Eberhard Bänsch Prof. Dr. Ulrich Rüde
November 21, 2012
Erklärung:
Ich versichere, dass ich die Arbeit ohne fremde Hilfe und ohne Benutzung anderer als der angegebenen Quellen angefertigt habe und dass die Arbeit in gleicher oder ähnlicher Form noch keiner anderen Prüfungsbehörde vorgelegen hat und von dieser als Teil einer Prüfungsleistung angenommen wurde. Alle Ausführungen, die wörtlich oder sinngemäß übernommen wurden, sind als solche gekennzeichnet.
Erlangen, den November 21, 2012
........................................................
Abstract This work aims at investigating the reliable predictions of droplet disintegration from the liquid jets in Rayleigh regime computationally and the comparison of the numerical simulations for specific parameters with the experimental data obtained from CSIRO (Australia). Discharging a liquid from a nozzle at sufficiently high velocity leads to a continuous jet that due to capillary forces break up into droplets. To study the problem, a system of 2D Navier-Stokes equations with appropriate initial and boundary conditions are solved at finite Reynolds number using the finite element sharp interface mesh moving method. The initial geometry is generated by Delaunay triangulation which evolves with the changing shape of the drop. The algorithm is able to capture fine features of the overall process such as the pinch off length/limiting length of drop breakup and the volume of the primary drop while it is still attached to the liquid before breakup. The accuracy of the work is verified by comparison of the algorithm in the dripping and Rayleigh regime with the existing works from literature. The results are in good agreement with the theory and therefore opens the door of new research in the area of liquid jets for highly viscous fluids at a very high temperature.
i
Acknowledgements This thesis would not have been possible without the guidance and the help of several individuals who in one way or another contributed and extended their valuable assistance in the preparation and completion of this study. First and foremost, my utmost gratitude to my supervisor Prof. Dr. Eberhard Bänsch who offered his continuous support and encouragement throughout this course of work. I am very thankful to him for his systematic guidance, suggestions and great effort he put into training me to produce an organized scientific report of the work. I am grateful to Dr. Mirco Wegener at CSIRO for regular discussions over the results, providing timely feed back and adding valuable comments. I am thankful to all my colleagues at the Applied Mathematics III, who have accepted me as a part of their family and for their moral and technical support. A very special thanks to Jordi Paul for being a very good friend and providing technical assistance even virtually. I would like to express my gratitude to Dr. Nicolas Neuß for his tips about using Latex, to Dr. Michael Fried for administrative matters, to Martin and Stephan for helping me to understand the software, to Frank for his moral supports from time to time. I feel very glad to find a sweet friend Kathrin Bäumler who was always there to take me out of stressful times and for her helpful advise to improve my thesis. My sincere thanks to Prof. Dr. Ulrich Rüde at Informatik 10 for his appreciation of work as a supervisor and also as COSSE coordinator in Germany, to my program coordinator Dr. Michael Hanke at KTH Sweden, for his support and encouragement through out the entire program. Special thanks to Dr. Harald Köstler at Informatik 10 for being helpful in administrative affairs at Erlangen and making my stay smooth at Germany and also for being a great friend. Last but not the least, my 6 years old son Sameer whose cute and lovely actions made me smile even during tough times.
ii
Contents Abstract
i
Acknowledgements
ii
Notations
vi
1 Introduction 1.1 Motivation . . . . . . . 1.2 History of liquid jets . 1.3 Breakup regimes . . . 1.4 Mechanism of pinch off 1.5 Organization of thesis
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
1 1 3 3 6 6
2 Mathematical model 2.1 Formulation of model . . . . . . . . . . . . . . . . 2.2 Governing equations . . . . . . . . . . . . . . . . 2.2.1 Boundary conditions . . . . . . . . . . . . 2.3 Non-dimensionalization . . . . . . . . . . . . . . 2.3.1 Dimensionless group . . . . . . . . . . . . 2.3.2 Navier Stokes equation . . . . . . . . . . . 2.3.3 Boundary conditions . . . . . . . . . . . . 2.4 Navier Stokes equations in cylindrical coordinates
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
7 7 8 8 9 9 10 11 12
3 Numerical method 3.1 Weak formulation . . . . . . . . . . . . . 3.2 Time discretization . . . . . . . . . . . . 3.3 Space discretization . . . . . . . . . . . . 3.3.1 ALE formulation . . . . . . . . . 3.4 NAVIER implementation aspects . . . . 3.4.1 Program structure . . . . . . . . 3.4.2 Rotational symmetry in NAVIER
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
14 14 15 16 16 17 17 19
4 NAVIER validation for dripping regime 4.1 Characterization for dripping mode . . . . . . . . . . . . . . . . . . . . . . . . 4.2 NAVIER results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
20 20 21
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
iii
. . . . .
. . . . .
. . . . .
. . . . . . .
. . . . .
. . . . . . .
. . . . .
. . . . . . .
. . . . .
. . . . . . .
. . . . . . .
4.3
4.4 4.5
4.2.1 Issues with classical NAVIER . . . . . Remeshing . . . . . . . . . . . . . . . . . . . . 4.3.1 Basic concept . . . . . . . . . . . . . . 4.3.2 Effect of remeshing on initial geometry NAVIER results with remeshing . . . . . . . . Comparison with related work . . . . . . . . . 4.5.1 Effect of Bond number . . . . . . . . . 4.5.2 Effect of capillary number . . . . . . .
. . . . . . . .
5 NAVIER Validation for Rayleigh Regime 5.1 Characterization for Rayleigh mode . . . . . . . 5.2 Aspects of jet disintegration from the literature 5.2.1 Growth of capillary instability . . . . . . 5.2.2 Jet breakup length . . . . . . . . . . . . 5.3 Adaption of NAVIER for Rayleigh regime . . . 5.3.1 Establishment of pinch off criteria . . . 5.3.2 Important decisions regarding remeshing 5.3.3 NAVIER main program flow . . . . . . . 5.4 Validation of numerical results . . . . . . . . . 5.5 Factors effecting simulation results . . . . . . . 6 Numerical results for the experimental data 6.1 Experimental description . . . . . . . . . . . . 6.1.1 Test setup . . . . . . . . . . . . . . . . 6.1.2 Test conditions . . . . . . . . . . . . . 6.1.3 Results . . . . . . . . . . . . . . . . . 6.2 Numerical results . . . . . . . . . . . . . . . . 6.2.1 Discussion . . . . . . . . . . . . . . . .
. . . . . .
. . . . . . . .
. . . . . . . . . .
. . . . . .
. . . . . . . .
. . . . . . . . . .
. . . . . .
. . . . . . . .
. . . . . . . . . .
. . . . . .
. . . . . . . .
. . . . . . . . . .
. . . . . .
. . . . . . . .
. . . . . . . . . .
. . . . . .
. . . . . . . .
. . . . . . . . . .
. . . . . .
. . . . . . . .
. . . . . . . . . .
. . . . . .
. . . . . . . .
. . . . . . . . . .
. . . . . .
. . . . . . . .
. . . . . . . . . .
. . . . . .
. . . . . . . .
. . . . . . . . . .
. . . . . .
. . . . . . . .
. . . . . . . . . .
. . . . . .
. . . . . . . .
. . . . . . . . . .
. . . . . .
. . . . . . . .
. . . . . . . . . .
. . . . . .
. . . . . . . .
. . . . . . . . . .
. . . . . .
. . . . . . . .
. . . . . . . . . .
. . . . . .
. . . . . . . .
. . . . . . . . . .
. . . . . .
. . . . . . . .
22 23 23 24 25 28 28 29
. . . . . . . . . .
31 31 32 32 33 34 35 36 38 40 42
. . . . . .
43 43 43 43 44 45 47
7 Conclusion and future work
49
List of Figures
52
List of Tables
53
Bibliography
56
iv
Notations α
the growth rate of disturbance
δ0
the initial amplitude of the perturbation
u
velocity distribution
γ
coefficient of surface tension
ΓD
Dirichlet boundary
Γf
free capillary boundary
D(u) = 12 (∇u + (∇u)T ), strain tensor g
gravity vector
n
the outward normal vector
Ca
=
Ld
= L/d, limiting length non-dimensionalized by diameter
p
pressure distribution
Q
inflow rate of liquid
WeG
Weber number for the ambient gas
WeL
Weber number for the liquid
We
=
µ
dynamic viscosity
Oh
=
Re
=
ρ
fluid density
σ
= 2µD(u) − pI, stress tensor
µU γ ,
capillary number
ρLU2 γ ,
Weber number
√
We Re ,
ρLU µ ,
Ohnesorge number Reynolds number
v
k
the wave number of disturbance
r0
initial radius of the jet
vi
Chapter 1
Introduction A liquid jet can be formally defined as a stream of a liquid or gas emanating from a nozzle/orifice under some force. The stability and disintegration phenomena of liquid jets has been a topic of great interest for over more than one hundred and fifty years due to its wide ranging applications in different industries. Besides its role in several applications, the dynamics of drop formation itself is a very fascinating topic for scientists on either an experimental or theoretical basis. Some of the common examples include sprays, water dripping from a slightly opened tap and ink-jet printers.
(a) Spray: courtesy of wiki
(b) Ink-jet printer: courtesy of cnx.org
Figure 1.1: Applications of Liquid Jets.
1.1
Motivation
The aim of this research is to compare the numerical simulation results with the experimental results obtained from CSIRO (Australia) for the study of the drop formation process during dry slag granulation. Slag is the main component used in blast furnaces, to cover molten metal and to remove impurities like sulphur. When this high temperature (1500◦ C ) slag is no more usable for the 1
furnace then there could be two possible options: either to dispose it in pits or convert it into granules which can be used as an additive for cement production. The later option seems quite reasonable as there is a huge market for it. Today most of the industries are using wet granulation where molten slag is quenched with huge amounts of water. This process has two drawbacks: heat cannot be recovered as water evaporates and energy is needed to dry wet granules. The CPSE High Temperature Processing Department in CSIRO comes up with an idea of dry slag granulation to overcome the drawbacks being faced in wet granulation. Figure.1.2 shows the process of dry slag granulation where slag is tapped onto a spinning disc and transported to the edge of the spinning disc. There it forms ligaments, which consecutively form droplets which detach from the ligaments. Air cools the flying droplets, they become a solid crust, but are still molten within the core. Hot air which is released during cooling can be reused. Solid granules go into a packed bed heat exchanger, the air cools droplets down to 100◦ C and again the hot air can be reused. Particles can be ground and sold as the cement additive.
Figure 1.2: Courtesy of CSIRO.
It is clear that the droplet formation from the ligaments is the decisive step and is driven by capillary forces. There could be differences in composition, temperature and surface tension which leads to Marangoni convection. This means that surface tension is not in equilibrium. The surface tension gradients are caused by temperature gradients tangential to the interface. 2
Therefore, it becomes important to study the time scale of non-equilibrium for surface tension and the time scale of droplet breakup. The aim is to determine the surface tension of molten slags for small time scales. This leads to the idea of studying laminar vertical jets based on Rayleigh stability theory and to measure the growth rate of the instability. Since each experiment is very time and cost consuming, it would be beneficial to simulate the slag jet behavior and predict through a computational approach.
1.2
History of liquid jets
The classical investigations for the liquid jets dates back from the work of Savart(1833), who proposed two laws which quantitatively describes the jet disintegration behaviour. (1) For constant diameter, the length of the continuous part of the jet is proportional to the jet velocity and (2) for constant velocity, the length of the jet is proportional to its diameter [20]. In 1873, Joseph Plateau carried out experimental investigations based on the work of Savart and showed that instabilities in liquid columns arose if its length exceeds the column diameter. Lord Rayleigh’s first analysis [26] in 1878 considered an inviscid jet in a vacuum. He obtained an expression for the growth rate of a symmetrical disturbance and predicted that the breakup length of a liquid jet is subject only to the forces of surface tension and inertia. Rayleigh assumed that the jet breaks up due to the disturbance having maximum growth rate. He obtained a relationship between growth rate and the wavelength of the disturbance. His theory was in total agreement with Savart. However, it did not take into account the important factors of viscosity of the liquid and the possible effect of the ambient medium on jet stability. Later in 1892, Rayleigh [27] added the effect of liquid viscosity to his theory which results in a complicated expression for the growth rate and this made it doubtful. Based on the experimental work of Haenlein [15], Weber was the first one to come up with the quantitative description of viscous jets taking into account the influence of surface tension. The disintegration time of viscous liquid jet proved to be greater than that of inviscid theory. Weber also showed that aerodynamic forces tend to propagate both symmetric and transverse disturbances [14]. Other notable names in this field are that of Grant and Middleman [14], Meister and Scheele [21], McCarthy and Molloy [20], Smith and Moss [25], and many others.
1.3
Breakup regimes
A liquid jet after emanating from the nozzle has to undergo the ultimate phase of breakup into droplets due to hydrodynamic instabilities. The mechanism of the jet breakup has been investigated in numerous theoretical and experimental approaches. According to [6], jet breakup properties are influenced by the following factors: 1. Hydrodynamic forces (viscous and capillary forces), 2. Nozzle internal flow (flow separation, cavitation), 3
3. Jet velocity profile at the nozzle exit, 4. Turbulence at the nozzle exit, 5. Thermodynamic and physical states of both the liquid and the ambient gas Experiments done by Weber and Haenlein in 1931 actually revealed that a jet may breakup in a variety of regimes defined by the forces that dominate the jet breakup under a given set of operating conditions. There are two ways currently established to categorize jet disintegration: Ohnesorge’s classification and the stability curve. Ohnesorge’s classification shown in Figure 1.3(a) has been a subject of study by numerous authors and it gives an idea of different regimes depending upon Re and Oh numbers (refer to the Notations for description). Another popular way of distinguishing between flow regimes is through the stability curve, see Figure 1.3(b) where jet breakup length (L) is plotted against jet velocity (v).
(a) Breakup regimes depending on Reynolds number and Ohnesorge number [11]
(b) General shape of stability curve by GRANT and MIDDLEMAN [14]
Figure 1.3: Two characterizations for identifying liquid jet regimes.
The main regimes identified are the dripping regime, the Rayleigh regime, the first wind induced regime, the second wind induced regime and the atomization regime depending upon surface tension, liquid inertia and aerodynamic forces acting on jet. Figure 1.4 illustrates the liquid jet shape when it disintegrates into different breakup regimes. Dripping regime: characterized by very low velocities such that the liquid only drips out of the tube and cannot form a definable jet breakup length. It is seen as Section ABC on 4
stability curve. Rayleigh regime: also called the varicose or dilational regime. It is characterized by low jet velocities, drop diameters larger than the jet diameter and jet breakup occurs when the surface disturbance reaches the magnitude of the jet radius. For this regime Weber showed that the breakup length was a linear function of velocity. The region ACD on the stability curve accounts for Weber’s theory for the breakup of laminar jet under the influence of surface tension forces. First wind induced regime: also called sinuous regime. It is characterized by jet velocities higher than those in the Rayleigh regime, drop diameters of the order of jet diameters and the jet gives a wavy appearance about its axis and hence the name sinuous regime. Second wind induced regime: characterized by high jet velocities, drops smaller than the jet diameter and jet disintegration occurs only a few jet diameters downstream from the nozzle exit due to highly significant liquid/atmosphere interactions. Aerodynamic forces become more prominent as the jet velocity increases and it can be seen from the stability curve at point D which starts deviating from linearity. In the transition region, flow changes from laminar to turbulent and the breakup length decreases with increasing velocity. Atomization regime: characterized by very high velocities, the jet disintegrates very close to the nozzle exit with drops much smaller than the jet diameter.
Figure 1.4: Jet behavior in different regimes.
5
1.4
Mechanism of pinch off
Pinch off is an important mechanism in liquid jet dynamics and hence needs to be studied in order to hasten or suppress the pinch off process as it may be needed in some cases. In order to explain the process, we consider the case of a vertically falling jet under the force of gravity. The formation of a drop/pendant in this case might be initiated by gravity. However surface tension is a strong factor behind it. The drop remains attached to the liquid if the volume is small. As soon as the volume exceeds a certain threshold, the drop surface cannot balance the body force anymore and starts accelerating downstream from the nozzle. A neck forms between the main drop and the nozzle, see Figure 1.5, giving opportunity to surface tension to constrict the neck. The jet elongates and the neck becomes smaller until eventually a drop is separated from the liquid. This is called that the jet is pinched off. The ligament front separated from the drop at the bifurcation point travels back towards the nozzle and a new bifurcation point is seen to be formed at the upper part of the ligament.
Figure 1.5: Different stages of drop evolution (From Schulkes 1994).
1.5
Organization of thesis
The rest of the thesis is organized in six more chapters. Chapter 2 describes the mathematical model for the liquid jet, i.e. the governing equations and the appropriate boundary conditions describing the system. Chapter 3 gives an introduction to the numerical method used to solve the mathematical model and the implementation aspects of software. Chapter 4 discusses about the first results of software in the dripping regime and its comparison with relevant literature. Chapter 5 gives a detailed insight into the theory of Rayleigh regime and the analysis of numerical results in comparison with the established results by researchers. Chapter 6 is about a comparison done between the experimental data from CSIRO and the numerical results followed by a discussion on the results. Chapter 7 is the last chapter of the report and it provides the concluding remarks about the results achieved during this course of work and suggestions for future work.
6
Chapter 2
Mathematical model 2.1
Formulation of model
We consider a physical situation where an axisymmetric drop of an incompressible Newtonian fluid with dynamical viscosity µ and mass density ρ is injected into a gaseous environment through a vertical, circular capillary tube of constant radius r0 as shown in Figure 2.1. The liquid flows through the capillary tube with a constant inflow rate Q such that the velocity profile approaches a fully developed parabolic profile at z = zi . The gravity vector g, the tube and symmetry axis all coincide and are taken along z-direction. The ambient environment surrounding the fluid exerts a constant pressure.
Figure 2.1: Representation of the domain.
7
2.2
Governing equations
Let Ω = Ω(t) ⊆ Rd , d = 3 be the region occupied by the liquid. Ω is changing with the free capillary boundary Γf . Owing to the situation, it is convenient to use a cylindrical coordinate system (r, z) with radial component r and axial component z. Once the fluid starts flowing from the tube, a drop in pear shape forms and keeps moving in the z-direction where z = L is the length from the tube exit to the tip of the drop. Eventually the drop pinches off point and the length from tube exit to the point where neck radius is minimum is called limiting length/breakup length/pinch off point Ld . At this instant of time, the drop formed is called the primary drop. The system is therefore modeled by the incompressible Navier-Stokes equations, given by the momentum conservation equation, ∂u ρ + (u · ∇) u = −∇p + µ∆u − ρg (2.1) ∂t and the incompressibility condition, ∇ · u = 0.
(2.2)
where, • u is the velocity distribution in • p is the pressure in mN2 • g denotes gravity vector in sm2
m s
The above equations are solved in the domain Ω(t) at times 0 ≤ t ≤ T . The stress tensor is given by σ = 2µD(u) − pI = µ(∇u + (∇u)T ) − pI,
(2.3)
where D(u) is the rate of strain tensor and is given by 1 ∇u + (∇u)T . (2.4) 2 Remark: Bold faced characters denote vectors and matrices whereas plain characters denote scalars. D(u) =
2.2.1
Boundary conditions
We have the following boundary conditions in our system: Along the boundary surface ΓD which defines the inner wall of the nozzle, the fluid should obey the no-slip condition and the velocity profile is parabolic at z = zi .
8
u = 1 − r2
on ΓD (t)
f or
zi ≤ z ≤ 0.
(2.5)
The free capillary surface which is unknown a priori has no mass transfer across it. Therefore the kinematic condition reads u · n = Uf
on Γf (t) ,
(2.6)
where n is the outward normal vector and Uf is the normal velocity of the free surface Γf . Along the fluid/ambient gas interface, we have a balance of linear momentum, n · σn = γH
on Γf (t) ,
(2.7)
with γ the coefficient of surface tension and mean curvature H = (κ1 + κ2 ) , κi are the principal curvatures of Γf . The free capillary surface also has a zero tangential stress, (i = 1, . . . , d − 1) where τ 1 , . . . , τ
2.3
d−1
τ i · σn = 0 on Γf (t) .
(2.8)
is the set of orthonormal tangential vectors on Γf .
Non-dimensionalization
Non-dimensionalization provides a better insight into a complex problem by removing units from the equations involving physical quantities by suitable substitution of variables such that all physical variables are O(1). Hence it is very important to choose correct scaling of dependent and independent variables.
2.3.1
Dimensionless group
We define non-dimensional parameters as follows: ρLU µ ρLU2 We = γ ρgL2 Bo = γ µU Ca = γ Re =
Reynolds number Weber number Bond number Capillary number
9
2.3.2
Navier Stokes equation
To start with the process, we need to identify first independent and dependent variables and then choose correct scaling factor. We choose the following scaling for corresponding scales:
Characteristic Scale
Dimensionless variables
Length L
x∗ =
Time T
t∗ = u∗ =
Velocity U
p∗ =
Pressure P
x L t T
u U,U
p P ,P
=
L T
= ρU 2
Table 2.1: Scaling factors for Navier Stokes Differential operators Due to above scaling factors, we have following dimensionless differential operators: 1 ∗ ∇ L 1 ∆ = 2 ∆∗ L ∂ 1 ∂ = ∂t T ∂t∗ ∇=
Substitutions Replace all the variables and differential operators in (2.1) by their respective non-dimensional equivalents.
∂u ρ + (u · ∇) u + ∇p − µ∆u = −ρgez ∂t 1 ∂ ∗ 1 ∗ µ 1 ∗ ∗ 11 ∗ ∗ ρgez ∗ u U + u U · ∇ u∗ U − ∆ u U+ ∇ p P =− ∗ 2 T ∂t L ρL ρL ρ ∗ 2 U ∂u U U P ∗ ∗ ρgez + (u∗ · ∇∗ ) u∗ − ν 2 ∆∗ u∗ + ∇ p =− T ∂t∗ L L Lρ ρ ∂u∗ U T ∗ T PT ∗ ∗ ρgT ez + (u · ∇∗ ) u∗ − ν 2 ∆∗ u∗ + ∇ p =− ∂t∗ L L LρU ρU
10
(2.9)
Since we have, UT LT = L TL νT 1 1 µ 1 =ν = 2 L UL ρ UL 2 PT ρU T UT = = LρU LρU L ρgT L ρgL2 σ ρgL2 σ ν 1 1 = ρg 2 = = = Bo ρU ρU σ ρU 2 L σ ρU ν U L Ca Re
=1 =
1 Re
=1 =
Bo We
Therefore (2.9) yields ∂u∗ 1 Bo + (u∗ · ∇∗ ) u∗ , − ∆∗ u∗ + ∇∗ p∗ = − ez in Ω∗ . ∗ ∂t Re We ∇ ∗ · u∗ = 0 in Ω∗
(2.10) (2.11)
Since we non-dimensionalize by radius R, we will use the substitution L = R.
2.3.3
Boundary conditions
On the free surface Γ∗f , we have the kinematic condition which is non-dimensionalized as u∗ · n = Uf∗ .
(2.12)
For the other boundary condition on free surface, we choose scaling H = L1 H ∗ which can be basically derived from the scaling of x and then apply some simple transformations.
n · σn = γH n · [2µD (u) − pI] n = γH n · µ(∇u + (∇u)T ) − pI n = γH 1 ∗ ∗ 1 ∗ ∗ T H∗ ∗ ÷ ρU 2 n· µ ∇ u U + (∇ u U ) − p PI n = γ L L L µ P γ n· ∇∗ u∗ + (∇∗ u∗ )T − p∗ 2 I n = H∗ ρLU ρU ρLU 2 1 1 ∗ ∗ ∗ ∗ ∗ T ∗ n· ∇ u + +(∇ u ) − p I n = H Re We 1 ∇∗ u∗ + +(∇∗ u∗ )T − p∗ I then we have the following nonIf we define now σ ∗ = Re dimensionalized form
11
n · σ∗n =
1 ∗ H We
(2.13)
Remark: All the characters represented by star (∗) represents dimensionless quantities. For the sake of brevity we will drop (∗) and relabel the scaled variables again by u, p etc.
2.4
Navier Stokes equations in cylindrical coordinates
Solving a three dimensional system involves several computational complexities. For the systems having rotational symmetry about the longitudinal axis, it is more convenient to use cylindrical coordinates (r, φ, z). To avoid complexity, I will not provide the complete derivation of the Navier Stokes equations from cartesian to cylindrical coordinates. Let U (x, y, z) : R3 → R3 , U = (Ux , Uy , Uz ) and u(r, φ, z) : R3 → R3 , u = (ur , uφ , uz ). Then we can write the Navier-Stokes in cylindrical coordinates as follows: ∂ur ∂p 1 1 2 ∂uφ 1 − ∆cyl ur − 2 ur − 2 + (u · ∇cyl ) ur − u2φ + ∂t Re r r ∂φ r ∂r ∂uφ 1 1 2 ∂ur 1 ∂p − ∆cyl uφ − 2 uφ + 2 + (u · ∇cyl ) uφ − ur uφ + ∂t Re r r ∂φ r ∂φ 1 ∂p ∂uz − ∆cyl uz + (u · ∇cyl ) uz + ∂t Re ∂z 1 ∂rur 1 ∂uφ ∂uz + + r ∂r r ∂φ ∂z
= fr = fφ = fz =0
Due to the assumption of rotational symmetry for our system, geometry will reduce to a simple geometry which is two dimensional and hence easier to solve.
Figure 2.2: Geometry simplification due to rotational symmetry.
Since we have axisymmetric flow so we can take advantage of symmetry by assuming that there are no azimuthal velocities in our system, i.e. 12
uφ = fφ = 0, ∂uφ ∂uz ∂p ∂ur = = = = 0. ∂φ ∂φ ∂φ ∂φ Therefore the final form of the Navier-Stokes system in rotational symmetric form will be the following:
∂t ur −
1 Re
1 ∆cyl ur − 2 ur + (u · ∇cyl ) ur + ∂r p = fr , r 1 ∂t uz − ∆cyl uz + (u · ∇cyl ) uz + ∂z p = fz , Re 1 ∂r (rur ) + ∂z uz = 0, r
(2.14) (2.15) (2.16)
where ∆cyl u = 1r ∂r (r∂r u) + ∂z2 u as defined in [4] is the Laplacian in cylindrical coordinates. On the free capillary surface Γf , we have 1 1 2∂r ur ∂r uz + ∂z ur H. ncyl · − pI ncyl = ∂ u + ∂ u 2∂ u Re We r z z r z z
(2.17)
The kinematic condition reads ur nr + uz nz = Uf .
13
(2.18)
Chapter 3
Numerical method This chapter introduces the numerical method used to discretize the system of equations described in (2.10) to (2.13). We solve a 2D rotational symmetric problem by using the software NAVIER [3] which is a finite element solver for flow problems with free surfaces.
3.1
Weak formulation
We have non-dimensionalized Navier-Stokes equations from Equations (2.10), (2.11) as follows, ∂t u + (u · ∇) u −
1 Bo ∆u + ∇p = − ez , Re We ∇ · u = 0.
(3.1) (3.2)
To derive a weak formulation, we multiply with a test function v ∈ H 1 (Ω)d and q ∈ L20 (Ω) with v|ΓD = 0. Integrating over the domain Ω, we obtain the following: Z
Z ∂t u · v +
Ω
Z (u · ∇) u · v −
Ω
Ω
Z Z 1 Bo ∆u · v + ∇p · v = − ez · v, Re We Ω ZΩ (∇ · u)q = 0.
(3.3) (3.4)
Ω
After applying integration by parts, we have: Z
Z ∂t u · v +
Ω
Ω
1 (u · ∇) u · v + 2Re
Z
Z D(u) : D(v) −
Ω
Ω
1 p(∇ · v) − We
Bo Hn · v = − We
Z Γf
Z q(∇ · u) = 0. Ω
For the detailed derivation, see for instance the diploma thesis of Vogel [29]. 14
Z ez · v, Ω
Now we define the following bilinear and trilinear forms for u, v ∈ H 1 (Ω)d and p ∈ L20 (Ω), Z 1 D(u) : D(v), a (u, v) = 2Re Ω Z b (u, v, w) = (u · ∇v) · w, Ω
Z c (p, v) = −
p∇ · v, Ω
Z s (u, v) =
u · v. Ω
The final form of the weak formulation of our problem reads: Find u, u(t, .) ∈ H 1 (Ω)d with u|ΓD = 0, p ∈ L20 (Ω) such that for all v ∈ H 1 (Ω)d with v|ΓD = 0, q ∈ L20 (Ω), it holds: Z
1 ∂t s(u, v) + a(u, v) + b(u, u, v) + c(p, v) − We
Ω
Z Hn · v = −
(3.5)
Γf
c(q, u) = 0
3.2
Bo s(ez , v) We
(3.6)
Time discretization
NAVIER uses the so-called fractional step θ scheme in an operator splitting variant, see [8, 2]. The benefit of using this method is that it decouples two major numerical difficulties of the Navier-Stokes equations, the solenoidal condition and the non-linearity. Therefore in each time step one has to solve two linear Stokes sub-problems and a viscous Burger equation accounting for the non-linear part. Another important aspect is the semi-implicit treatment of curvature terms in (3.5) which can be represented by, Z Z Hn · φ = − ∇x · ∇φ, (3.7) Γf
Γf
where ∇ is the tangential gradient on Γf , x is the position vector, H is the curvature.
15
Semi-implicit discretization is enforced by using identity xn+1 = xn + τ un+1 , where xn , xn+1 are the position vectors of Γnf , Γn+1 respectively. Therefore we can write f Z Z Z ∇xn+1 · ∇φ = ∇xn · ∇φ + τ ∇un+1 · ∇φ, (3.8) Γn f
Γn f
Γn f
which links the unknown for geometry Ω and the flow variables u, p [3]. The algorithm used is unconditionally stable (in a context of space-time finite elements) and efficient, see [3].
3.3
Space discretization
To discretize in space, NAVIER uses Taylor-Hood elements, i.e. continuous piecewise quadratic functions for velocity and continuous piecewise linear functions for pressure.
Figure 3.1: Representation of the Taylor-Hood element.
In order to have a good resolution for the approximation of the free surface, isoparametric elements are used for the discretization at the free boundary. Then we have boundary elements represented by curvilinear triangles of the same polynomial degree as that of the velocity distribution, i.e. piecewise quadratic elements.
3.3.1
ALE formulation
In our problem, the domain changes its shape due to the movement of the free capillary boundary Γf . In order to determine the curvature at the free boundary edges with high accuracy, it is important to have an explicit representation of the free boundary. This explicit representation of the free boundary can be obtained with the ALE formulation. The concept of ALE formulation is used in general to handle the problems dominated by free surface flows. The arbitrary Lagrangian-Eulerian (ALE) is a finite element formulation in which the computational domain is not a priori fixed in space (e.g. Eulerian-based finite element formulations) or attached to material (e.g. Lagrangian-based finite element formulations). Instead the computational domain may move according to some (arbitrary) grid velocity, which in the simulation will be chosen in a way, that the free boundary coincides with a subset of 16
element edges at all time-steps. Figure 3.2 illustrates the basic concept in ALE formulation. There exists a one to one mapping between the domains, for details see [16].
Figure 3.2: Pictorial representation of ALE.
The effect of ALE on the Navier-Stokes equations is an additional advective term in the momentum equation, see [5].
3.4
NAVIER implementation aspects
NAVIER is a flow solver written in FORTRAN which solves 2D rotational symmetric NavierStokes equations based on unstructured triangular grids with a sharp interface mesh moving method to track free boundary. This section describes the implementation details of NAVIER.
3.4.1
Program structure
Since NAVIER is a software enriched with several subroutines used to perform dedicated tasks. The modular approach of programming makes it easier to integrate new modules with less overhead. For the readers of this report I would like to explain the structure of NAVIER in a simple way and especially keep focus on the highly important features. A graphical illustration of the main program flow can be seen in Figure 3.3 and there is a zoom to the core component, i.e. the Glowinski fractional θ step. In later chapters, the
17
updated flow charts will reflect the changes made in the structure of the code especially because of adding new routines for geometry handling.
Figure 3.3: Main program structure.
18
3.4.2
Rotational symmetry in NAVIER
Rotational symmetry in NAVIER is achieved by scaling for ur such that ur . (3.9) r This is done in order to resolve the singularities of the differential operator at r = 0. The axisymmetric code then solves for the unknown (ur /r, uz ). This cylindrical 2D-ansatz was developed by Patrick Lailly in his PhD dissertation [18]. us :=
19
Chapter 4
NAVIER validation for dripping regime This chapter describes the very first performance check of NAVIER for the liquid jet project by comparison of the numerical results in dripping regime with the existing work from literature followed by a discussion about adaption in software in order to achieve accurate and reliable results.
4.1
Characterization for dripping mode
The jet breakup in dripping occurs right after the nozzle and at relatively low flow rates and the resulting drops are mono-disperse in size as already discussed in Chapter 1. For We 1, drops of constant mass periodically detach from the nozzle at a constant frequency. This regime is often referred as Periodic Dripping (PD) as the droplets detach from the nozzle at a downstream distance of approximately one diameter [10]. A water faucet as shown in Figure 4.1 is one of the most common examples of dripping mode almost known by everyone. The familiarity and simplicity of the system made it ideal for investigating liquid behavior in this regime.
Figure 4.1: Leaky faucet (Wikipedia).
20
The first attempt towards validation of NAVIER was to check the performance of the code in dripping mode. A number of qualitative results from the literature is available for comparison with the results obtained from NAVIER.
4.2
NAVIER results
This section describes the results obtained from NAVIER by studying the effect of variation in Bond number and capillary number on the breakup length. This idea was taken from Wilkes [32] and the same work is used for the validation purpose. NAVIER takes an unstructured grid for geometry as shown below. Figure 4.2 is the typical initial geometry which is then evolved in each time step through front tracking method in the way so as to form a drop and simulation finally stops when the triangulation is totally distorted near pinch off which means that any of the element in the domain exceeds the limit of the greatest/smallest angle criteria. In order to make a distinction between the initial version of NAVIER and the one adapted for the project, I will call the initial version as Classical NAVIER.
Figure 4.2: Typical domain at t=0 and its triangulation.
Table 4.1 gives the listing of results obtained by NAVIER and its visualization can be seen in Figure 4.3. The limiting length L/d has been defined in Chapter 2. The determination of breakup length is done manually in Paraview (data analysis and visualization application). As the free surface moves, triangulation elements get deformed so as to adjust themselves with the new grid calculations until it forms a drop which is near about to detach from the rest of the fluid.
21
(a) At Re = 3.0, Ca = 0.012
Bond Number Bo 0.02 0.06 0.10 0.15 0.20 0.30 0.40
(b) At Re = 3.0, Bo = 0.17
Limiting Length L/d 10.837 8.469 7.691 7.186 6.861 6.485 6.227
Capillary Number Ca 0.02 0.08 0.12 0.20 0.25 0.30 0.40
Limiting Length L/d 7.250 8.263 8.734 9.683 10.193 10.937 11.711
Table 4.1: Classical NAVIER results for Dripping Regime.
Figure 4.3: Classical NAVIER for Re=3.0, We=0.75, Bo=0.17.
4.2.1
Issues with classical NAVIER
The zoom in view of the necking area in Figure 4.4 shows that in order to reach pinch off point there needs to be more constriction in neck and some of the elements flipped due to higher elongation required. Due to this high deformation of mesh, simulation breaks down before a realistic pinch off point against our expectations. This leads to an idea of improving the mesh strategy somehow to have a good quality mesh and to avoid excessive distortion before pinch off.
22
Figure 4.4: Focussed view of neck area.
4.3
Remeshing
In each time step the domain evolves and the mesh is deformed because of the movement of the free capillary surface Γf . However, there already exists a provision of mesh smoothing to cater for mesh deformations. This mesh moving strategy in Classical NAVIER along with smoothing is not enough to guarantee a good quality of mesh especially when the jet is approaching pinch off point. As an improvement, remeshing is introduced in NAVIER to enhance its capability to adjust from time to time the geometry in case of excessive deformation of domain. Therefore remeshing of the whole domain is performed when necessary.
4.3.1
Basic concept
The implementation of remeshing in NAVIER is based on exporting boundary points of the old triangulation and then using of the 2D mesh generator TRIANGLE [24] which is a C program for two-dimensional mesh generation and construction of Delaunay triangulations, constrained Delaunay triangulations, and Voronoï diagrams. Triangle is fast, memory-efficient, and robust; it computes Delaunay triangulations and constrained Delaunay triangulations exactly. Guaranteed-quality meshes (having no small angles) are generated using Ruppert’s Delaunay refinement algorithm. After that, the new mesh is transported back to NAVIER. After each time step, the current mesh is checked to determine if remeshing is necessary. Remeshing can be performed either after every n time steps or in case if one of the following criteria is met: • If the minimum angle of a triangle is less than 10◦ . • If the maximum angle of a triangle exceeds 140◦ . • If the volume of a triangle is changed largely, i.e. approximately by a factor of 10.
23
A single call to remeshing subroutine consists of the following steps: 1. Generate a list of all boundary nodes for the whole domain and write them as planar straight line graph (PSLG). 2. Pass this information to TRIANGLE as an input which will then generate new mesh. 3. TRIANGLE generates only straight simplices, therefore, correction at boundary edges is required as we represent boundary by curvilinear edges. The newly inserted straight edges are projected to their appropriate curved positions as they were in the old mesh. 4. Interpolation is performed to transfer the old grid date i.e. un , pn onto the new mesh by L2 projection scheme, for details see [5].
4.3.2
Effect of remeshing on initial geometry
In NAVIER, remeshing is performed in the very first time step on the simple initial triangulation of domain. The benefit of remeshing at this stage is to have a control over initial mesh quality, see Figure 4.5.
(a)
(b)
Figure 4.5: Initial geometry with remeshing.
24
After this initial remeshing of domain, the remeshing subroutine is called after every 30 time steps or if any of the criteria based on triangle quality is met as discussed above. There are several remeshing techniques discussed in literature but the benefit of using this approach is its simplicity, portability, efficiency and less overhead with the end result of a very fine quality mesh.
4.4
NAVIER results with remeshing
The simulations were tested again with this new feature for the same data set as was used before. Figure 4.6 gives the pictorial illustration of the results obtained and the statistics of the run are listed in Table 4.2. A closer look at the post processing data reveals that improvement in breakup length is observed which is approximately equal to the breakup lengths achieved by Wilkes for their analysis, see next section for detailed analysis. Furthermore the problem of flipping of triangles near pinch off is also cured. (a) At Re = 3.0, Ca = 0.012
Bond Number Bo 0.02 0.06 0.10 0.15 0.20 0.30 0.40
Limiting Length L/d 11.121 8.645 7.864 7.371 7.084 6.760 6.573
(b) At Re = 3.0, Bo = 0.17
Capillary Number Ca 0.02 0.08 0.12 0.20 0.25 0.30 0.40
Limiting Length L/d 7.583 8.869 9.626 11.008 11.897 12.773 14.598
Table 4.2: NAVIER results with remeshing for dripping regime. As can be seen in Figure 4.6, the mesh quality now is better and the mesh is adapted to maintain quality triangles with the evolution of jet after some time intervals. The last stage is the time when the liquid jet is about to pinch off and the droplet is connected to the rest of the liquid through a very thin thread. A very fine neck constriction is observed and the visualization of the results from numerical simulation proved to reflect the expected behavior in this regime. The zoom in view shows the geometry at the necking area represented by very fine triangles, see Figure 4.7 as compare to the triangulation in Figure 4.4. The insertion of the remeshing subroutine in the main program tree is depicted in Figure 4.8. After startup routines, remeshing is called if necessary. In case of remeshing, the new triangulation with projected grid velocities is now given as input to the solver which then performs Glowinski step and gives grid parameters for the next time step. Then in diagnostic routine, a check is added for setting the remesh flag if necessary. 25
Figure 4.6: NAVIER with remeshing for Re=3.0, We=0.75, Bo=0.17.
Figure 4.7: Zoom into neck area.
26
Figure 4.8: NAVIER program flow with remeshing.
27
4.5
Comparison with related work
For validating the results obtained from NAVIER, it is compared with the computational data by Wilkes [32] in dripping regime. The limiting length of the jet is studied by varying Bo and Ca from 0 to 0.4 while keeping Re fixed at 3.0. The limiting length particularly for this validation is the length of the jet from the nozzle exit until the bottom of the drop. Below are some plots showing a comparison of both works.
4.5.1
Effect of Bond number
It is expected that the limiting length L/d should rise with increasing Bo but this also depends on a critical value of Bo, see [32]. This critical value of Bo is different for different Reynolds numbers for example at Re = 0.03, the critical value is 0.3 after which breakup length starts increasing [32]. In our computations, we have not computed until the critical value of Bo after which L/d starts increasing, therefore, it can be observed in Figure 4.9 that L/d decreases as Bo increases for small values of Bo.
Computed variation of drop length when Re=3.0 and Ca=0.012 12
Classical NAVIER Wilkes NAVIER with remeshing 11
L/d
10
9
8
7
6
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
Bo
Figure 4.9: Computed variation of limiting length L/d with the Bond number Bo.
Figure 4.10 shows the visualization of gravitational Bond number (Bo) effect on limiting length while keeping the capillary number and Reynolds number fixed. An important observation made here was that the volume of the primary drop decreases with increase in gravitational Bond number which in turn decreases the limiting length.
28
(a) Bo = 0.02
(b) Bo = 0.06
(c) Bo = 0.10
(d) Bo = 0.15
(e) Bo = 0.20
(f) Bo = 0.40
Figure 4.10: Effect of Bond number variation on limiting length.
4.5.2
Effect of capillary number
Figure 4.11 shows the variation of pinch off length (L/d) with capillary number (Ca) using a fixed Bond number and Reynolds number. It has been observed that increasing the capillary number has a considerable effect on the limiting length. An important aspect to be considered here is how this range of capillary number effects the droplet size and limiting length. The drops in this regime do not change in volume considerably instead the thread connecting primary drop to the rest of the liquid grows in length with increase in capillary number as shown in Figure 4.12. Therefore based on these results one can conclude that Ca has pronounced effect on L/d and the effect of Ca on primary drop volume is quite small as compared to the effect of Bo on primary drop volume. The successful introduction of remeshing feature in Classical NAVIER for analyzing drop formation in dripping mode seems quite promising when compared with data from Wilkes [32].
29
Computed variation of drop length when Re=3.0 and Bo=0.17 18
Classical NAVIER Wilkes NAVIER with remeshing
16
14
L/d
12
10
8
6
4
2
0
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
Ca
Figure 4.11: Computed variation of limiting length L/d with the capillary number Ca.
(a) Ca = 0.02
(d) Ca = 0.25
(b) Ca = 0.12
(e) Ca = 0.30
(c) Ca = 0.20
(f) Ca = 0.40
Figure 4.12: Effect of capillary number variation on limiting length.
30
Chapter 5
NAVIER Validation for Rayleigh Regime This chapter describes the validation of our software NAVIER in Rayleigh regime by comparing its performance with benchmarks from literature. It also includes the details about how the numerics are adapted to accord for the expected performance in this regime. The chapter concludes with a discussion on the results obtained and their correlation with literature.
5.1
Characterization for Rayleigh mode
In contrast to the dripping mode where drops detach near the nozzle exit, the droplet detachment in Rayleigh regime moves downstream from the nozzle exit typically a distance greater than 10 times of the nozzle diameter and a continuous jet is formed as can be seen in Figure 5.1. It is also called varicose or dilational regime.
Figure 5.1: Rayleigh mechanism (from Wiedemann, Oswald 2005).
31
Rayleigh breakup occurs at low velocities which are, however, larger than in dripping. The breakup length increases linearly with the increasing velocity and the drops formed in this regime are larger than the jet diameter [6]. According to the theoretical predictions made by Rayleigh and Weber, it is generally said that the jet breakup occurs when surface disturbances reaches the order of magnitude of jet radius. Each regime is bounded by certain critical velocity which can be expressed in terms of the Weber number. Clanet and Lasheras [10] have described dripping regime for We 1 whereas Birouk and Lekic [6] have described Rayleigh regime for WeL > 8, WeG < 0.4 by taking into account ambient gas also. If WeG > 0.4, we enter the next regime i.e. first wind induced regime. WeL is the Weber number for the liquid and WeG is the Weber number for the ambient gas.
5.2 5.2.1
Aspects of jet disintegration from the literature Growth of capillary instability
The jet instabilities purely arise from surface tension which is also called the capillary phenomenon. Surface tensions between the liquid and the surrounding gas amplify small perturbations under various flow conditions and these perturbations propagate along the interface as surface waves of specific amplitude and wavelength [11]. The main cause of the jet breakup in Rayleigh regime is the amplification of these capillary instabilities. The axisymmetric perturbation may be expressed in terms of Fourier series as: δ(z, t) = δ0 eαt+ikz ,
(5.1)
where δ0 is the initial amplitude of perturbation along the axial coordinate z, α is its growth rate and k is the wave number. The evolution of the jet radius in terms of time and space can be represented in terms of perturbation as: r(z, t) = r0 + δ(z, t).
(5.2)
Free surface perturbations spread along the interface as surface waves [11]. The growth rate 0 of the capillary instability α for a given wave number k = 2πr λ (r0 is the initial radius of capillary tube) is determined by the dispersion relation obtained by Rayleigh as α2 = −
kr0 2 I1 (kr0 ) 1 − (kr ) , 0 τ2 I0 (kr0 )
(5.3)
where I0 , I1 are the modified Bessel functions of first kind of order 0 and 1 respectively and τ is the characteristic time scale of the motion [13] which is expressed as r σ τ= . (5.4) ρr03
32
Dimensionless growth rate α τ
0.4
0.35
0.3
0.25
0.2
0.15
0.1
0.05
0
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Dimensionless wave number k r 0 Figure 5.2: Dispersion curve describing the relationship between dimensionless growth rate and wave number.
A typical dispersion curve is shown in Figure 5.2. The plot between dimensionless growth rate and wave number illustrates that the maximum amplification of the instability occurs at kmax r0 ≈ 0.7 and the corresponding wavelength is given by: kmax =
5.2.2
2πr0 1 =√ , λmax 2
√ ⇒ λmax = 2 2πr0 .
(5.5)
Jet breakup length
Weber’s theory [30] very well describes the Rayleigh regime and the model suggested by him has been extensively used by experimentalists for comparison purpose. He has formulated a quantitative relationship for the prediction of breakup length in this regime which is, L a = ln d δ0
0.5
(We)
3We + Re
.
(5.6)
where L is the limiting length/breakup length non-dimensionalized by diameter d as defined in Chapter 2. According to Weber’s assumption, the interface of a cylindrical jet is perturbed by axisymmetric disturbances which grow in time. The coefficient ln a/δ0 is the relative scale of the initial disturbance where a = r. Based on the experimental results of Haenlein, Weber determined a constant value for ln a/δ0 as 12.
33
A sharp discrepancy was observed by Grant and Middleman [14] between calculated and measured values who found that the critical velocity in the stability curve V ∗ , see Figure 1.3(b), is either underestimated or overestimated by Weber’s theory according to their experimental conditions. They have reported a different constant value for the factor ln a/δ0 in their initial experiments for a range of nozzle design and liquids: L 3We 0.5 . = 13.5 (We) + d Re
(5.7)
But due to the discrepancy between their data and Weber’s theory, they have suggested modifications in Weber’s equation in order to get a best fit for their experimental data. They have derived an empirical relationship between Oh and ln a/δ0 .
ln
a = −2.66ln(Oh) + 7.68, δ0
(5.8)
where Oh is the Ohnesorge number defined by Oh = We1/2 /Re. They also give a further correlation of breakup length for this region by the least squares fit of the data as given by 3We 0.85 L 0.5 = 19.5 (We) + . d Re
(5.9)
Takahashi and Kitamura [17] correlated the coefficient ln a/δ0 for various Newtonian liquids with Oh by L = 6.5We0.5 (1 + 3Oh)/Oh0.15 d L = 12We0.5 (1 + 3Oh) d
for Oh < 0.017,
(5.10)
for Oh > 0.017.
(5.11)
Smith and Moss [25] reported a value of 13 for ln a/δ0 . K. Sallam et al. [23] in their experimental work correlated the mean liquid column breakup length which has been the best fit for their measurements: L = 5.0We1/2 d
5.3
for We < 400.
(5.12)
Adaption of NAVIER for Rayleigh regime
For the validation in dripping regime, remeshing based on the criteria discussed in Chapter 4 was incorporated in classical NAVIER which was enough to give comparable results. However for validation in the Rayleigh regime, few more features were integrated into NAVIER.
34
5.3.1
Establishment of pinch off criteria
In dripping regime, the jet breakup length was calculated using a manual approach as automatic calculation in simulation is difficult to develop. The precise definition of pinch off point is varying in literature and there is no unified criteria formulated until now. However we have come up with two ideas of finding a pinch off point for our simulations which is implemented successfully into NAVIER. It is possible to select either of them for determining pinch off and stop the simulations before it breaks down completely because of low quality triangles. The most obvious benefit for having these criteria implemented is the automation of process with more accuracy. 1. Minimum radius criteria This is the simplest and straight forward criteria which checks the radial distance from symmetry axis to the free boundary edges in pinch off region. Drop breakup is said to have occurred when the radius at some nodes falls below a specified threshold value which in our simulation is set to 0.002r0 , where r0 is the initial radius of the jet at the nozzle exit and the radius at this point is termed as Rmin .
Figure 5.3: The point of minimum radius.
2. Flow separation criteria Another indication of reaching near pinch point could be the separation of flow as shown in Figure 5.4, where the velocity vectors are uni-directional (a) and as the drop is near to be detached from the rest of the liquid, some of the liquid starts flowing upstream shown by velocity vectors (b). We check the vertical component of velocity, i.e. uy against a threshold value of −0.20 for elements describing the pinch region.
35
(a) velocity vector before flow separation.
(b) velocity vector after flow separation.
Figure 5.4: Velocity at flow separation.
5.3.2
Important decisions regarding remeshing
Although remeshing in NAVIER has been discussed quite in detail in Chapter 4, some important aspects of remeshing came into view especially while working in Rayleigh regime. Therefore I would like to summarize few points in this section regarding important decisions taken to adapt NAVIER for Rayleigh regime. The criteria used for remeshing has been discussed in Chapter 4 and they are still valid. 1. Frequent remeshing vs. not so frequent remeshing In dripping we have been using remeshing after every fixed n time steps and also as and when it is required. After having tested simulations without remeshing at frequent intervals, it was observed that simulations still work quite accurately if we do not remesh the triangulation quite often. The frequency of remeshing in this case is approximately 16 times in 1000 time steps which seems to be better if we remesh frequently after every 30 time steps. Therefore in Rayleigh regime, we only perform remeshing to cure for large deformation in any element when the criteria for remeshing is met. 2. Concave edges and bad angle problem One of the remeshing criteria is to check for the minimum angle in a triangle. This minimum angle is termed as the bad angle. It was observed that the bad angle determination for edges at the neck area is wrong, therefore, remeshing is performed earlier than the triangle minimum angle is reached. The edges defining the neck area are concave which makes the determination of minimum angles difficult, as the real angle (as given by the quadratic edge) may be much smaller than the angle just given by the three vertices of the triangle. In NAVIER, only edges at the free boundary are quadratic, so this problem only existed on the free capillary surface in the pinch region. For precise determination of bad angles in case of concave edges, a new strategy is implemented which takes into account the parametrization of a curved simplex.
36
3. Projection scheme In dripping after performing remeshing of the triangulation, we have been interpolating the values of velocity and pressure for new mesh via L2 -projection. In order to make sure that the velocity in the new subspace is discretely divergence free, weighted Stokes projection has been added to enforce the divergence constraint on velocity field. Since the initial results for the Rayleigh regime validation were not as expected, therefore, the addition of Stokes projection was to ensure the divergence free condition for the velocity distrubution in the new mesh. However, it did not have a considerable effect on the simulation result, therefore, we have not repeated the validation for dripping regime with the new projection scheme. 4. Global remeshing vs. adaptive remeshing The most interesting region of study in formation of a drop is the pinch off area which needs to be predicted with maximum possible accuracy to produce reliable results. In this context it might not be necessary to perform remeshing over whole domain. Indeed remeshing of the neck formation region would suffice. Therefore a comparison study was made between global (uniform refinement of whole domain) and adaptive (refinement only in necking region) strategies. Another development during this study was to determine the threshold values for Rmin and uy needed to automate pinch off criteria. For adaptive remeshing, the elements in the neck area need to be refined only. Therefore, it is necessary to first identify the elements in the neck area for which curvature calculation was incorporated into NAVIER. The idea besides calculating curvature for all edges is that the edges in neck area are concave and therefore their curvature is negative. Simulations were allowed to run until triangles are completely distorted exceeding the minimum (2◦ ) or maximum (179◦ ) limit of any triangle element and simulation is terminated. The outcome of test runs is tabulated for comparison in Tables 5.1, 5.2, 5.3.
Parameters Remeshing Limiting length uy at flow separation Rmin
Global remeshing ∆t = 0.005 ∆t = 0.01 ∆t = 0.04 15 times 17 times 14 times 6.177 6.4047 6.7486 -0.205 -0.221 -0.2152 0.0256 0.0158 0.01801
Adaptive remeshing ∆t = 0.005 ∆t = 0.01 ∆t = 0.05 22 times 19 times 12 times 6.396 6.469 6.685 -0.20 -0.210 -0.214 0.00478 0.00655 0.03796
Table 5.1: Comparison of global and adaptive remeshing for Re = 20, Bo = 0.10, We = 4. To compare, Re and Bo are fixed and We is varied from 4 to 12. In terms of remeshing, global and adaptive refinement are almost comparable until We = 8, i.e. the number of remeshing performed during simulation runs are almost equal but for We > 8 it is observed that number of times remeshing is performed in adaptive refinement is almost 3 times than in global refinement, which is computationally expensive. No considerable 37
Parameters Remeshing Limiting length uy at flow separation Rmin
Global remeshing ∆t = 0.005 ∆t = 0.01 ∆t = 0.04 17 times 19 times 17 times 20.924 20.802 22.19 -0.203 -0.20 -0.2103 0.0346 0.0297 0.0042
Adaptive remeshing ∆t = 0.005 ∆t = 0.01 ∆t = 0.05 27 times 16 times 15 times 19.90 20.29 20.298 -0.201 -0.204 -0.3405 0.0379 0.05825 0.04828
Table 5.2: Comparison of global and adaptive remeshing for Re = 20, Bo = 0.10, We = 8. Parameters Remeshing Limiting length uy at flow separation Rmin
Global remeshing ∆t = 0.005 ∆t = 0.01 ∆t = 0.04 17 times 14 times 16 times 38.82 40.299 40.55 -0.282 -0.20 -0.201 0.0155 0.0244 0.0318
Adaptive remeshing ∆t = 0.005 ∆t = 0.01 ∆t = 0.05 57 times 47 times 28 times 32.04 46.41 32.21 -0.20 -0.20 -0.201 0.0122 0.1261 0.0638
Table 5.3: Comparison of global and adaptive remeshing for Re = 20, Bo = 0.10, We = 12. difference has been observed in terms of limiting length between both strategies. Note that limiting length here is not used for correlation with theory instead only mean value of length is used for comparison purpose. Flow separation used to occur mostly around uy = −0.20 which is the finally taken value for threshold in simulations. Adaptive refinement seems to approach more precise formation of neck for smaller time step size in case of We = 4 only but in other cases the behavior is quite arbitrary therefore depending on the visualization, the threshold for Rmin is chosen as 0.002r0 . To arrive at a logical conclusion, it is inferred on the basis of this study that adaptive refinement does not seem to provide a cost-effective solution for remeshing as We is increased and it did not proved to be helpful for increasing the accuracy of results. Therefore it was decided to stay with global refinement as we still have a fine triangulation of domain, less remeshing and stable liquid jets.
5.3.3
NAVIER main program flow
Figure 5.5 shows the program flow of NAVIER with the adaptions required for the project. The highlighted green parts shows the changes made to software while the grey sections shows the original modules which are unchanged.
38
Figure 5.5: NAVIER results in Rayleigh regime
39
5.4
Validation of numerical results
For the purpose of validation in Rayleigh regime, the performance of NAVIER has been tested against the works mentioned in Section 5.2 for the set of parameters specified in Table 5.4. According to the study from available literature, it was clear that it is not possible to determine ln a/δ0 a priori for our simulations so it seems reasonable that a plot n o of NAVIER results 1 3We 2 with breakup length L/d as a dependent parameter and (We) + Re as an independent parameter should yield a straight line of slope ln a/δ0 . Parameter Re We Bo Oh ∆t
Value 60 6.0 - 26.0 0.0 0.04 - 0.08 0.0015
Table 5.4: Test conditions for validation in Rayleigh regime.
50
45
L/d
40
35
30
25
20
15 3
3.5
4
4.5
5
5.5
6
6.5
7
7.5
8
1
(We) 2 + 3We/Re
Figure 5.6: NAVIER results in Rayleigh regime
It seen in n can be o Figure 5.6 that plotting the dimensionless length L/d against the parameter 1 3We (We) 2 + Re yields a straight line of slope 7 which is the derived value for the coefficient ln a/δ0 in case of NAVIER. To establish a correlation with existing literature, the dimensionless jet breakup length of NAVIER is compared with the breakup relations discussed in Section 5.2.2 as shown in Fig40
80
Grant & Middleman Smith & Moss Weber theory Kitamura NAVIER Sallam et al.
70
60
L/d
50
40
30
20
10
0
0
5
10
15
20
25
30
We
Figure 5.7: Dimensionless correlation of breakup data with existing works at Re = 60.
ure 5.7. The data plot shows that the performance by NAVIER satisfies the behavior of liquid jets in Rayleigh regime as described by the stability curve in Figure 1.3(b). The linear behavior and stability of NAVIER gives us the satisfaction about the correctness of numerical method used even with the offset present. At this point it is difficult to establish a one to one correlation with Weber’s theory, but we observe a good correlation with the work of Sallam et al. [23]. Although there exists an offset between our result and existing work from literature but there might be some possible factors behind it which are discussed in detail in the next section. Another measure for validating results is the determination of the amplitude of the dominant wave λmax , which is responsible for the breakup of liquid jet when the amplitude of perturbation is equal to the jet radius. In our case, we observe wavelength of 4.5, see Figure 5.8 while the expected value is √ λmax = 2 2πr0 ⇒
λmax ≈ 9.0r0 .
41
(r0 = 1)
(5.13)
Figure 5.8: Illustration of the wavelength for maximum amplification at Re = 60, We = 16, Bo = 0.
5.5
Factors effecting simulation results
This section enlightens the possible reasons behind simulation behavior in Rayleigh regime and the factors that might be the cause of a smaller factor for logarithmic term as compared to Weber’s equation. • Time step size Time step size could have a significant effect on simulation results. It was an observation that selection of a smaller time step size is really beneficial for simulation and this factor seems to effect more when liquid jet is near to the pinch off time. Currently we use a time step size ∆t = 10−02 and results are reasonable in acceptable computation time. It will be better if adaptive time stepping can be used, i.e. choosing even more smaller time step size when simulation is near to the pinch off point. • Nozzle design Selection of nozzle geometry could have a significant impact on the jet breakup length according to the survey carried out by McCarthy and Molloy [20]. The total velocity distribution at nozzle exit depends on the interaction between physical properties and the geometry of the nozzle. The important parameters of nozzle such as the aspect ratio (ratio between length and diameter), the contraction angle (can vary between 0◦ e.g. a long pipe to 180◦ e.g. a sharp edge orifice) might need to be considered in model.
42
Chapter 6
Numerical results for the experimental data This chapter describes the experimental setup and its results obtained from M. Wegener [31] for their project followed by a comparison of the numerical simulations carried out for the corresponding parameters. They have carried out the initial experiments considering water jet as validation of results is pretty easy as known results are available from the literature. However, the final goal is to carry out experiments with slag jets for the study of interfacial instabilities in slags and its comparison with numerical simulation which is not a part of present work.
6.1 6.1.1
Experimental description Test setup
The initial investigations have been performed for the water jet in Rayleigh mode to visualize the breakup process. Water is issued from a pipette of capacity 20 ml into air at an ambient temperature of 20◦ C. A constant supply of water is used to maintain the level of liquid inside pipette for a constant flow rate of 0.55ml/s and constant velocity. The jet is falling vertically downwards under the force of gravity which finally disintegrates into droplets at some distance downstream from the tube exit. No external excitation is used and the jet is allowed to adapt the natural perturbations. The test duration was 5 seconds and 8000 images were captured per second.
6.1.2
Test conditions
The experiment was carried out under certain operating conditions. Table 6.1 collects the summary of present test conditions. Based on the physical parameters of liquid and choosing the length scale as L= d, we calculate the dimensionless group as follows:
43
Parameter Liquid Initial jet diameter (d) Jet exit mean velocity (U) Flow rate (Q) Liquid density (ρ) Liquid viscosity (µ) Surface tension (γ)
Value Water 0.63mm 1.76m/s 0.55ml/s 998.2 kg/m3 1.01 ×10−03 Pas 0.0728 N/m
Table 6.1: Summary of test conditions.
ρLU , µ 998.2 × 0.63 × 10−03 × 1.76 = , 1.01 × 10−03 = 1085.
Re =
ρLU2 , γ 998.2 × 0.63 × 10−03 × 1.762 = , 0.0728 = 26.4.
We =
ρgL2 , γ 998.2 × 9.81 × (0.63 × 10−03 )2 , = 0.0728 = 0.0533. √ We Oh = , √Re 26.4 = , 1085 = 4.735 × 10−03 . Bo =
6.1.3
Results
The experimental results proved to be in good agreement with theory. The experimental data reported an average length of 40 mm and according to theory of Weber, breakup length for this specific test setup is determined as,
44
3W e L , = 12 (W e)0.5 + d Re L 3 × 26.4 0.5 , = 12 (26.4) + d 1085 = 62.5, ⇒ L = 39.4mm. The jet surface coordinates captured quite accurately are used to determine the wavelength through the distinct necks and swells that develops during the flow. The average wavelength √ was found to be 2.8mm and the theoretical value calculated as per equation λ = 2πd = 2.78mm. The experiment captures accurately the maximum amplification for k = 0.71 as expected, see the maximum value of the dispersion curve as shown in Figure 5.2, plotted against Equation (5.3).
6.2
Numerical results
This section particularly describes the results obtained from our simulation with NAVIER. The objective is to evaluate the performance of the chosen numerical method by simulating a liquid jet of water choosing the same set of parameters as used in the experimental setup. In our simulations, we have made test runs by first considering the effect of gravity and later on neglecting it. Since the Bond number is small, we did not observe a very pronounced effect on results. Figure 6.1 is a frame from the animation we have received from CSIRO and Figure 6.2 shows the visualization of numerical simulation considering the effects of gravity. The purpose of this illustration is to give an idea to the reader about jet evolution and droplet shape in simulation. Since we do not have the time stamps of the animation, it is not possible to correlate both evolutions in time. In our simulation, we have a limitation currently to predict breakup length until primary drop break up and it did not take into account satellite drops and further droplet formation.
45
Figure 6.1: Image capture of water jet during experiment at CSIRO.
(a) tdimensionless = 3.0
(b) tdimensionless = 7.5
(c) tdimensionless = 12.0
(d) tdimensionless = 18.0
(e) tdimensionless = 23.68
Figure 6.2: Evolution of jet at different time steps in simulation.
46
We have compared the results in terms of the jet breakup length. Figure 6.3 illustrates the numerical as well as experimental breakup length against the dimensionless breakup length from Weber, see Equation (5.6) and Sallam et al., see Equation (5.12). The jet breakup length obtained with simulation is 9.5mm which is quite less than the length obtained from experimental results. The expected breakup length from simulations as shown by the black curve in Figure 6.3 is 23.13mm. Note that the length scale is dimensionless in the plot. This expected result for simulations is based on the validation done in Rayleigh regime for simulations where we have obtained the logarithmic coefficient value as per our simulations. Correlation of experimental and numerical result 90
Weber theory Sallam et al. NAVIER expected NAVIER actual Experimental
80
70
L/d
60
50
40
30
20
10
0
0
5
10
15
20
25
30
35
40
45
50
We
Figure 6.3: Comparison of experimental and numerical result with existing works.
This comparison leaves us in a clueless state. As we have determined in last section that experimental data correlates very well with Weber theory and it is also obvious from the plot of data in Figure 6.3. However numerical simulation result could not reach its expected value i.e. a factor of 7 for the logarithmic term and the jet is pinched off quite earlier than expected which is very contradictory with the performance of simulation in Rayleigh regime as dicussed in Chapter 5.
6.2.1
Discussion
The behaviour of numerical simulation in case of experimental data now needs to be studied to provide possible explanation behind it. In my point of view and through the study of relevant literature on liquid jets, there might be following reasons effecting simulation result.
47
• Logarithmic term in breakup equation The coefficient ln(a/δ0 ) in equation (5.6) as a variable needs more attention. The behavior of this coefficient is found conflicting with the literature. Also the result obtained from the numerical simulation for the experimental data did not match with the expected behavior of simulations in Rayleigh regime as discussed in Chapter 5. • Primary drop formation In simulation, we have currently a limitation of calculating the free fall of jet until the formation of the primary droplet. However, in experiment we observe a fully developed jet which seems to be in quasi-steady state and forming droplets at constant length. In real world and also in experiment, the initial state of the jet might differ and it needs some time to achieve this steady state. It might be then difficult to compare experimental and numerical result if the time frame is different for both. • Nozzle parameters We did not take into account the influence of nozzle design parameters on simulation results. The physical characteristics of the nozzle can play a significant role in the disintegration of the jet [20].
48
Chapter 7
Conclusion and future work This chapter re-evaluates the main results obtained during this work and in the light of these results a few suggestions are made for the future work. This study investigated the theoretical aspects of droplet formation in liquid jets and factors influencing this breakup. The goal is to simulate the behavior of liquid jets for the dripping and Rayleigh regime and to verify the numerical results with a real world application specifically in our case the initial experimental results from CSIRO. A finite element method for solving the mathematical model as prescribed in Chapter 2 is presented which is capable of simulating a liquid jet emanating from a nozzle or orifice with predefined inlet conditions. The algorithm is able to capture the gross features of the phenomena such as the limiting length of a drop at breakup, the volume of the primary drop and its features such as the development of a neck in a viscous drop approaching breakup. However the current model does not predict the formation of satellite droplets and secondary droplets. A sharp interface mesh moving method is used to keep track of the free surface. Occasional remeshing is used to address relatively large mesh deformations which possibly occurs due to a moving grid. An automatic estimate of the jet breakup length is difficult to lead in the simulations as there is no fixed criteria for determination of pinch off point from the literature. However we have come up with two different ideas to determine the breakup length based on the change in geometry and change in a vector physical quantity. Several simulations have been performed to show the relative comparison of our numerics with the foundation theory on liquid jets as well as the popular research in this area. In particular, we have first compared in a satisfactory way the breakup length for a range of Bond numbers and capillary numbers in dripping regime. The agreement with the work of Wilkes [32] is quite good and gave us confidence over the numerics especially in this regime. Then we have made investigations in Rayleigh breakup mode for a range of We by changing the inlet velocity and we compare our results with Weber theory, Grant and Middleman, Smith and Moss, Kitamura and Takahashi and Sallam et al. for breakup in this regime. We have found good correlation between our work and Sallam et al. and more over the per49
formance of simulation is linear as expected according to the stability curve. After having reasonable results in Rayleigh regime, we have verified our simulations with the experimental data received from CSIRO which also have parameters in the range of this regime. Although numerical simulation gives stable result, it did not match the expected results for this regime. Some of the possible reasons for this behavior could be the controversial behavior of the logarithmic term in breakup equation. Also we are comparing the primary drop formation with a fully developed steady jet and its difficult to correlate jet behavior at different time stamp. At this point there exists a discrepancy between numerical simulations and experimental results but it is difficult to provide concrete reasons behind this deviation from expected behavior. Numerous extensions of the present work are possible. An important extension of the present study concerns continuous formation of drops or the so-called continuous jet in steady state. It has been observed that time step size has a considerable effect over simulation results. Adaptive time stepping can be used, that means a moderate time step size to be used initially and then switching to a smaller time step size as soon as the simulation is approaching the breakup area could be helpful in achieving accuracy for the breakup length. A useful area of future research would be the development of a model that also accounts for the effect of different nozzle parameters on the breakup process. According to McCarthy and Molloy [20], the nozzle aspect ratio has a highly significant effect on the initial jet velocity profile and subsequent jet surface shape. Another direction of future work might be to study the effect of different initial mesh resolutions and check the convergence of the breakup length of the jet simulations. Long term goals of this project includes comparing numerical simulations against slag jets experimental data. Surface tension which is a function of temperature and concentration could play an important role in this case. It is possible that presence of concentration gradients in molten slag could give rise to significant gradients in surface tension which means that surface tension is not in equilibrium and changes with time which in turn can cause Marangoni convection. In that case, simulation needs to incorporate Marangoni effects as well. To conclude the discussion, the ALE based finite element method with sharp interface mesh moving strategy works quite well in dripping regime. The results in Rayleigh regime satisfy the linear behavior as expected from the stability curve but with a different value for the logarithmic term as compared to the value determined by Weber. The possible reasons for the deviation from expected behavior are discussed but to elucidate if those reasons are the exact reason behind the observed simulation behavior was not possible in the time frame for the present work.
50
List of Figures 1.1 1.2 1.3 1.4 1.5
Applications of Liquid Jets. . . . . . . . . . . . . . . . . Courtesy of CSIRO. . . . . . . . . . . . . . . . . . . . . Two characterizations for identifying liquid jet regimes. . Jet behavior in different regimes. . . . . . . . . . . . . . Different stages of drop evolution (From Schulkes 1994).
. . . . .
1 2 4 5 6
2.1 2.2
Representation of the domain. . . . . . . . . . . . . . . . . . . . . . . . . . . Geometry simplification due to rotational symmetry. . . . . . . . . . . . . . .
7 12
3.1 3.2 3.3
Representation of the Taylor-Hood element. . . . . . . . . . . . . . . . . . . Pictorial representation of ALE. . . . . . . . . . . . . . . . . . . . . . . . . . Main program structure. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
16 17 18
4.1 4.2 4.3 4.4 4.5 4.6 4.7 4.8 4.9 4.10 4.11 4.12
Leaky faucet (Wikipedia). . . . . . . . . . . . . . . . . . . . . . . . . . . Typical domain at t=0 and its triangulation. . . . . . . . . . . . . . . . . Classical NAVIER for Re=3.0, We=0.75, Bo=0.17. . . . . . . . . . . . . Focussed view of neck area. . . . . . . . . . . . . . . . . . . . . . . . . . . Initial geometry with remeshing. . . . . . . . . . . . . . . . . . . . . . . . NAVIER with remeshing for Re=3.0, We=0.75, Bo=0.17. . . . . . . . . . Zoom into neck area. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . NAVIER program flow with remeshing. . . . . . . . . . . . . . . . . . . . Computed variation of limiting length L/d with the Bond number Bo. . . Effect of Bond number variation on limiting length. . . . . . . . . . . . . Computed variation of limiting length L/d with the capillary number Ca. Effect of capillary number variation on limiting length. . . . . . . . . . .
. . . . . . . . . . . .
20 21 22 23 24 26 26 27 28 29 30 30
5.1 5.2
Rayleigh mechanism (from Wiedemann, Oswald 2005). . . . . . . . . . . . . Dispersion curve describing the relationship between dimensionless growth rate and wave number. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . The point of minimum radius. . . . . . . . . . . . . . . . . . . . . . . . . . . Velocity at flow separation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . NAVIER results in Rayleigh regime . . . . . . . . . . . . . . . . . . . . . . . NAVIER results in Rayleigh regime . . . . . . . . . . . . . . . . . . . . . . . Dimensionless correlation of breakup data with existing works at Re = 60. .
31
5.3 5.4 5.5 5.6 5.7
51
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . . . . . . . . .
33 35 36 39 40 41
5.8
6.1 6.2 6.3
Illustration of the wavelength for maximum amplification at Re = 60, We = 16, Bo = 0. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
42
Image capture of water jet during experiment at CSIRO. . . . . . . . . . . . Evolution of jet at different time steps in simulation. . . . . . . . . . . . . . Comparison of experimental and numerical result with existing works. . . . .
46 46 47
52
List of Tables 2.1
Scaling factors for Navier Stokes . . . . . . . . . . . . . . . . . . . . . . . . .
10
4.1 4.2
Classical NAVIER results for Dripping Regime. . . . . . . . . . . . . . . . . . NAVIER results with remeshing for dripping regime. . . . . . . . . . . . . . .
22 25
5.1 5.2 5.3 5.4
Comparison of global and adaptive remeshing for Re = 20, Bo = 0.10, We = 4. Comparison of global and adaptive remeshing for Re = 20, Bo = 0.10, We = 8. Comparison of global and adaptive remeshing for Re = 20, Bo = 0.10, We = 12. Test conditions for validation in Rayleigh regime. . . . . . . . . . . . . . . . .
37 38 38 40
6.1
Summary of test conditions. . . . . . . . . . . . . . . . . . . . . . . . . . . . .
44
53
Bibliography [1] E. Bänsch. Numerical methods for the instationary Navier-Stokes equations with a free capillary surface. Habilitationsschrift, Universität Freiburg, 1998. [2] E. Bänsch. Simulation of instationary, incompressible flows. Acta mathematica Universitatis Comenianae, 67, no. 1:101–114, 1998. [3] E. Bänsch. Finite element discretization of the Navier-Stokes equations with a free capillary surface. Num.Math, 88, no. 2:203–235, 2001. [4] E. Bänsch, C. P. Berg, and A. Ohlhoff. Uniaxial extensional flows in liquid bridges. Journal of Fluid Mechanics, 000:1–27, 2004. [5] E. Bänsch, J. Paul, and A. Schmidt. An ALE FEM for solid-liquid phase transitions with free melt surface. Zentrum für Technomathematik, page 000, 2012. [6] M. Birouk and N. Lekic. Liquid Jet breakup in quiescent atmosphere. Atomization and Sprays, 19(6):501–528, 2009. [7] G. Brenn, Z. Liu, and F. Durst. Linear analysis of the temporal instability of axisymmetrical non-newtonian liquid jets. International Journal of Multiphase Flow, 26:1621 – 1644, 2000. [8] M. O. Bristeau, R. Glowinski, and J. Periaux. Numerical methods for the Navier-Stokes equations. Applications to the simulation of compressible and incompressible viscous flows. Computer Physics Report, 6:73–187, 1987. [9] A. Cervone, S. Manservisi, and R. Scardovelli. Simulation of axisymmetric jets with a finite element Navier-Stokes solver and a Multilevel VOF approach. Journal of Computational Physics, 229:6853–6873, 2010. [10] C. Clanet and J. C. Lasheras. Transition from Dripping to Jetting. Journal of Fluid Mechanics, 383:307–326, 1999. [11] J. Delteil, S. Vincent, A. Erriguible, and P. Subra-Paternault. Numerical investigations in Rayleigh breakup of round liquid jets with VOF method. Computer and Fluids, 50:10–23, 2011. [12] G. Dziuk. An algorithm for evolutionary surfaces. Num. Math., 58:603–611, 1991.
54
[13] J. Eggers and E. Villermaux. Physics of liquid jets. Report on progress in Physics, 71(3):036601, 2008. [14] R. P. Grant and S. Middleman. Newtonian Jet Stability. A.I.Ch.E.J, 12:669, 1966. [15] A. Haenlein. Uber den Zerfall eines Flussigkeitsstrahles. Forshung Gebiete Ingenierw. 2, 4:139, 1931. [16] B. Höhn. Numerik für die Marangoni-Konvektion beim Floating-Zone Verfahren. Phd thesis, Mathematische Fakultät der Albert-Ludwigs-Universität Freiburg, 1999. [17] Y. Kitamura and T. Takahashi. Stability of a Contracting Liquid Jet. Memoirs School of Engineering, Okayama Univ., 7:61–84, 1972. ´ [18] P. Lailly. R´ esolution num´ erique des Equations de Stokes en Sym´ etrie de R´ evolution par une M´ ethode d’El´ ements finis non Conformes. PhD thesis, Universit´ e Paris XI, 1976. [19] S. P. Lin. Breakup of Liquid Sheets and Jets. Cambridge University Press, 2003. [20] M.J. McCarthy and N.A. Molloy. Review of Stability of Liquid jets and the influence of Nozzle design. The Chemical Engineering Journal, 7:1–20, 1974. [21] B. J. Meister and G. F. Scheele. Drop formation at low velocities in liquid-liquid systems. A.I.Ch.E J., 15:700, 1969. [22] Y. Pan and K. Suga. A numerical study on the breakup process of laminar liquid jets into a gas. Physics of Fluids, 18(5):052101–11, 2006. [23] K.A. Sallam, Z. Dai, and G.M. Faeth. Liquid breakup at the surface of turbulent round liquid jets in still gases. International Journal of Multiphase Flow, 28:427–449, 2002. [24] J. R. Shewchuk. Triangle: Engineering a 2D Quality Mesh Generator and Delaunay Triangulator. In Ming C. Lin and Dinesh Manocha, editors, Applied Computational Geometry: Towards Geometric Engineering, volume 1148 of Lecture Notes in Computer Science, pages 203–222. Springer-Verlag, May 1996. From the First ACM Workshop on Applied Computational Geometry. [25] S.W.J. Smith and H. Moss. Proc. Roy Soc., 93A:331, 1917. [26] J.W. Strutt and Lord Rayleigh. On the instability of jets. Proceedings of London Mathematical Society, 10:4–13, 1878. [27] J.W. Strutt and Lord Rayleigh. On the instability of cylindrical fluid surfaces. Philosophical Magazine, 34:177–180, 1892. [28] W. van Hoeve, S. Gekle, J. H. Snoeijer, M. Versluis, and M. P. Brenner et al. Breakup of diminutive Rayleigh jets. Physics of Fluids, 22:122003, 2010. [29] F. Vogel. Modeling and simulation of the liquid metal pinch, 2008. Diploma thesis. [30] K. Weber. Z. angew. Math. Mech., 11:136, 1931. 55
[31] M. Wegener. CSIRO Process Science and Engineering, Australia. collaboration via emails. [32] E. D. Wilkes, S. D. Philips, and O. A. Basaran. Computational and experimental analysis of dynamics of drop formation. Physics of Fluids, 11(12):3577–3598, 1999.
56
