@prefix vivo: . @prefix edm: . @prefix ns0: . @prefix dcterms: . @prefix skos: . vivo:departmentOrSchool "Applied Science, Faculty of"@en, "Engineering, School of (Okanagan)"@en ; edm:dataProvider "DSpace"@en ; ns0:degreeCampus "UBCO"@en ; dcterms:creator "Jalaal, Maziyar"@en ; dcterms:issued "2012-06-11T22:47:26Z"@en, "2012"@en ; vivo:relatedDegree "Master of Applied Science - MASc"@en ; ns0:degreeGrantor "University of British Columbia"@en ; dcterms:description """The work described in the present thesis is related to a series of projects that I worked on toward the better understanding of fragmentation phenomena. In the past decades, the science of fragmentation has attracted many attentions within the researchers due to its wide range of applications. However, because of the complexity of the subject, even its basic concepts need more investigations. This thesis starts with an introduction to fragmentation of droplets using experimental or numerical approaches. It is discussed that the current mathematical and experimental tools are not able to describe all the details. Thus, high performance numerical simulations are the best alternatives to study the breakup of droplets. The introduction is followed by a discussion on the numerical method and the ranges of the non-dimensional groups. It is described that an adaptive, volume of fluid (VOF) method based on octree meshing is used, providing a notable reduction of computational cost. The rest of the thesis basically discusses the obtained results using direct numerical simulations. Two main geometries are investigated: falling droplets and droplets in a stream. For the case of falling droplets, three simulations with different Eötvös numbers are performed. For the case of droplets in a stream, two-dimensional and three-dimensional simulations are performed for a range of Weber number. The results are compared with the available mathematical theories and it is shown that the analysis presented here precisely demonstrates the mechanism of the bag breakup of falling droplets and instability growth over the droplets in an external high-speed flow. The outcomes can significantly assist the development of the secondary atomization and turbulent two-phase flows modelling."""@en ; edm:aggregatedCHO "https://circle.library.ubc.ca/rest/handle/2429/42476?expand=metadata"@en ; skos:note " Direct Numerical Simulation of Fragmentation of Droplets by Maziyar Jalaal B.Sc., University of Tabriz, 2009 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIRMENTS FOR THE DEEGREE OF Master of Applied Science in The College of Graduate Studies (Mechanical Engineering) THE UNIVERSITY OF BRITISH COLUMBIA (Okanagan) May 2012 © Maziyar Jalaal, 2012 ii Abstract The work described in the present thesis is related to a series of projects that I worked on toward the better understanding of fragmentation phenomena. In the past decades, the science of fragmentation has attracted many attentions within the researchers due to its wide range of applications. However, because of the complexity of the subject, even its basic concepts need more investigations. This thesis starts with an introduction to fragmentation of droplets using experimental or numerical approaches. It is discussed that the current mathematical and experimental tools are not able to describe all the details. Thus, high performance numerical simulations are the best alternatives to study the breakup of droplets. The introduction is followed by a discussion on the numerical method and the ranges of the non-dimensional groups. It is described that an adaptive, volume of fluid (VOF) method based on octree meshing is used, providing a notable reduction of computational cost. The rest of the thesis basically discusses the obtained results using direct numerical simulations. Two main geometries are investigated: falling droplets and droplets in a stream. For the case of falling droplets, three simulations with different Eötvös numbers are performed. For the case of droplets in a stream, two-dimensional and three-dimensional simulations are performed for a range of Weber number. The results are compared with the available mathematical theories and it is shown that the analysis presented here precisely demonstrates the mechanism of the bag breakup of falling droplets and instability growth over the droplets in an external high-speed flow. The outcomes can significantly assist the development of the secondary atomization and turbulent two-phase flows modelling. iii Table of Contents Abstract ............................................................................................................................................. ii Table of Contents ............................................................................................................................. iii List of Tables ..................................................................................................................................... v List of Figures .................................................................................................................................. vi Acknowledgements ......................................................................................................................... xii Dedication ...................................................................................................................................... xiv 1. Introduction ................................................................................................................................... 1 1.1 Fragmentation of Droplets ................................................................................................. 1 1.2 Motivation .......................................................................................................................... 7 1.3 Outline of This Thesis ........................................................................................................ 8 2. Methodology & Non-Dimensional Groups ................................................................................... 9 2.1 Geometries ......................................................................................................................... 9 2.2 Governing Equations ........................................................................................................ 11 2.3 Flow Solver (Gerris) ....................................................................................................... 12 2.3.1 Temporal Discretisation ................................................................................... 13 2.3.2 Spatial Discretisation ....................................................................................... 14 2.4 Non-Dimensional Numbers ............................................................................................. 17 3. Fragmentation of Falling Droplets .............................................................................................. 21 3.1 General Description of The Fragmentation Procedure ..................................................... 21 3.2 Drop Deformation and Bag Formation ............................................................................ 24 3.3 Bag Breakup ..................................................................................................................... 30 3.4 Number of Fragments ...................................................................................................... 46 4. Statistical Characteristics of Fragments ...................................................................................... 51 4.1 An Overview on Droplet Deformation and Breakup ....................................................... 51 4.2 Fragment Number Variation ............................................................................................ 52 4.3 Fragment Size and Velocity Distribution ......................................................................... 55 4.4 Variation of Non-Dimensional Groups ............................................................................ 65 5. Deformation of Droplet in a Stream ............................................................................................. 71 5.1 An Overview on Droplet Dynamics in a Stream .............................................................. 71 iv 5.2 Two-Dimensional Results ................................................................................................ 73 5.2.1 Grid Tests ......................................................................................................... 73 5.2.2 Deformation Mechanism .................................................................................. 76 5.2.3 Effect of Weber Number .................................................................................. 78 5.2.4 Comparison with Theories (Shear Instabilities) ................................................ 82 5.2.4.1 Combination of Potential Flow and Very Thin Vorticity Layer ...... 82 5.2.4.2 Non-Zero Vorticity Layer ................................................................ 85 5.3 Three-Dimensional Results .............................................................................................. 87 5.3.1 Deformation Mechanism .................................................................................. 88 5.3.2 Comparison with Theories (Transverse Azimuthal Modulation) ..................... 92 6. Conclusion ................................................................................................................................... 97 6.1 Future Work ..................................................................................................................... 98 Bibliography .................................................................................................................................. 100 v List of Tables 2.1. Non-dimensional numbers studied for the case of falling droplets .......................................... 18 2.2. Dimensionless groups for the case of a droplet in a stream ..................................................... 19 3.1. Values of Reynolds and Weber numbers ................................................................................. 28 4.1. Parameter in log-normal PDF modelling ................................................................................. 60 5.1. Grid properties at 0.87  ....................................................................................................... 74 5.2. Comparison of the present direct numerical simulations and equation 5.6 .............................. 85 5.3. Comparison of the present direct numerical simulations and equations 5.7 and 5.10 .............. 87 5.4. Comparison of the present direct numerical simulations and equation 5.14 ............................ 95 5.5. Comparison of the present direct numerical simulations and equations 5.15 and 5.16 ............ 95 vi List of Figures 1.1. a) Droplet in an external flow. b) Droplet falling in a constant gravity (acceleration). ................. 1 1.2. The algorithm of droplet disintegration and fragments distribution. ............................................. 2 1.3. Main modes of droplet breakup. .................................................................................................... 3 1.4. Breakup modes in terms of Weber number. ................................................................................... 5 2.1. Schematic sketch of the initial condition of simulations for falling droplets. ............................... 10 2.2. Droplet in a stream. Flow direction is from left to right. ............................................................. 11 2.3. Quad-tree grid adaption from level 0 (dotted box) to level 4 (green boxes). ............................... 14 2.4. Grids generated for a case of falling droplet ( 0.761  for 216Eo  , 0.05dOh  , 0.05cOh  , * 10  ). Left: interface of the drop. Right: cross-section of the drop in 0z  plane; each cell is colored by the refinement level-number of the grid. ......................................................... 15 2.5. Iso-surfaces of grids size for fragments after the fragmentation of a falling droplet. a) Dispersed fragments with different shapes for 216Eo  , 0.05dOh  , 0.05cOh  , * 10  . b) Iso-surface for grid size of, * 0 0.5dldl D   . c) * 0.2dl  . d) * 0.1dl  . e) * 0.05dl  . ............. 16 2.6. Iso-surfaces of grids size for a deformed droplet in a stream. a) * 0.1dl  . b) * 0.05dl  . c) * 0.02dl  . d) * 0.005.dl  ............................................................................................................. 16 3.1. The deformation and fragmentation of a liquid drop for: , , , . (a) , (b) 0.1647  , (c) 0.2642  , (d) 0.3575  , (e) 0.4434  , (f) 0.5183  , (g) 0.5924  , (h) 0.6728  . ............................................................................... 22 3.2. Fragmentation of a liquid water in a continuous air jet flow. and . Left: : the green arrow notes one of the liquid bridges generated in the experiment. Purple and orange arrows demonstrate the penetrated core and unstable ligaments, respectively. middle: : green arrow denotes the fragments created by disintegration of liquid columns. The purple and the orange arrows show the deformation and breakup of the core, right: : green arrow points to a liquid column made by the breakup of the core. The orange arrow indicates the fragments made in the fragmentation. The purple arrow signifies a relatively large fragment generated during the core breakup (Cao et al., 2007). ................................................................. 23 288Eo  0.05dOh  0.05cOh  * 10  0.0  29We  0.0015dOh  24.3t ms 35.6t ms 42t ms vii 3.3. Deformation of an initially spherical liquid drop deformation for: . a) 0.0825  , b) 0.137  , c) 0.175  , d) 0.210  , e) 0.244  , f) 0.278  , g) 0.317  , h) 0.358  . Colours demonstrates the vertical velocity of the interface,   * vV p gD     . Temporal radius and height of the drop are illustrated in picture (h). The gravitational acceleration is in –y direction. ......................................................................... 24 3.4. Temporal variation of drop diameter during the falling procedure. denotes the slope of the linear trendlines interpolated (red for the deformation, blue for the bag growth) for different stages. ........................................................................................................................... 25 3.5. Temporal aspect ratio of falling droplets. .................................................................................... 26 3.6. Left: temporal variations of Reynolds number. Right: temporal variation of Weber number. .... 27 3.7. Drag coefficient vs. time. ............................................................................................................. 29 3.8. Experimental results of Bermond and Villermaux (2005) for bursting liquid films for a) (the arrow shows a ligament formed after attaching of two growing holes), b) , c) , scales in each picture show size of a selected hole (Bermond and Villermaux, 2005). ................................................................................................................ 31 3.9. Top views of hole growth for (a-d, ): a) 0.414  , b) 0.427  (the red lines depict the boundaries of the holes), c) 0.447  (the transparent red circle shows the area of the upper torus), d) 0.463  . (e-h, ): e) 0.267  , f) 0.278  , g) 0.289  (the transparent red circle shows the area of the upper torus), h) 0.300  . ............................... 32 3.10. Holes positions for . a) Three dimensional shape of the penetrated drop, 0.447  . b) A cross-section of the drop in the same time: the area denoted by the red circles is the liquid torus. The area noted by the blue torus (circles) is the area which the majority of the holes are initially generated there. Arrows point to the high-curvature partof the interface. ................................................................................................................... 33 216Eo  m 1.03M  1.07M  1.21M  144Eo  288Eo  144Eo  viii 3.11. a) The holes generated on the drop for at 0.427  . Cross-sections are plotted with different colours. b) Expansion of the selected hole. Corresponding times are, 0.411  , 0.419  and 0.431  , respectively. The rounded edge (rim) of the hole is obvious for the thinner part. The sheet has non-uniform, time varying thickness. (The last time step is magnified in a separated box and demonstrates small amplitude capillary waves in the thinner part. These waves, as well as the round edge are not seen in the thicker part of the hole). c) Growth of two neighbourhood punctures and formation of the liquid bridge. Corresponding times are, 0.419  , 0.431  , 0.444  and 0.463  , respectively. The rounded edges again are evident for the plane which connects to holes. Some very weak capillary waves are seen for this plane while they are not observed for the left and right hand side parts that actually show the upper rim. Holes grow until a liquid column, shown by a circular cross-section, is formed. This column in fact connects the upper rim to the core of the drop. d) The development of the hole in the radial direction and separation of the liquid torus. Corresponding times are, 0.411  , 0.419  and 0.431  , respectively. This picture shows the part of the droplet in which the upper rim is disconnected from the rest of the bag. By development of the punctures the bag vanishes and only the liquid bridges that link the core and the torus would remain. (The necking before the hole creation is demonstrated in a box, 0.406  ). .................................................................................................................. 35 3.12. Retraction speed for a selected hole in simulation of . ............................................... 38 3.13. Bag breakup of the falling droplet . a) 0.377  , b) 0.386  , c) 0.394  , d) 0.403  . .............................................................................................................................. 39 3.14. Bag fragmentation for 1) The web of attached ligaments created by several holes. 2) Growth and deformation of the ligaments. 3) Increase in the number of holes. (green and yellow arrows demonstrate the upper rim and flattened core of the drop, respectively). 4-6) Further enlargement of the punctures and detachment of ligaments. (the red box in picture 5 depicts an attached ligament into the upper rim.) 7-9) Breakup of ligaments and creation of droplets. (the red boxes in pictures 8 and 9 point up a ligament attached to the flattened core and a stable droplet, respectively). .............................. 40 144Eo  144Eo  216Eo  216Eo  ix 3.15. Instability and disintegration of the bag for 1-2) Small amplitude waves at the thin sheet. (blue arrows point to some of the waves). Corresponding times are 0.380  and 0.382  , respectively. 3) A hole created on the interface. ( 0.384  ) 4) Growth of the generated hole in the previous step. ( 0.386  ) (The blue arrow shows a noteworthy deformation close to the core. The green arrow highlights the waves generated during the retraction. The yellow arrow shows another significant deformation for the part which is attached to the upper rim.) 5) Detachment of the bag from the flattened core and formation of other holes. ( 0.3882  ) (The blue arrow shoes the gap generated between the bag and the core. The yellow array demonstrates the cross section of a fragment created by growth of two neighboring holes). 6) Further deformation of the interface and enlargement of the gaps. ( 0.3883  ) (Blue arrow depicts the deformation appeared during the retraction) 7-8) Necking and formation of another hole which forms more fragments. Corresponding times are 0.392  and 0.395  , respectively. 9) Generated fragments which will experience more disintegration. ( 0.396  ) (Yellow and green arrows demonstrate a ligament attached to the upper torus and the drop core). .......................................................... 41 3.16. Cross-section view of a ligament showing drop formation, . Domain is coloured by non-dimensional pressure, . a) The capillary waves grow causing a thicker area in the ligament. ( 0.457  ). b-d) Growth of the waves over the ligament. Corresponding times are 0.460  , 0.461  and 0.464  , respectively. e) Neck generation. ( 0.466  ). f) Pinch-off occurs and the droplet detaches from the ligament. ( 0.468  ). ................................................................................................................................ 42 3.17. Rupture of the flattened droplet core for . Pictures are captured from bottom view and a transparent hollow circle is added to separate the core for clarity. a) Flat core after the bag breakup. ( 0.424  ). b) Initial holes and cracks generated on the interface of the core. ( 0.441  ). c) Detachment of the ligaments in the generated web. ( 0.444  ). d-e) Growth of the holes and development of the crumbling. Corresponding times are 0.447  and 0.450  , respectively. f) Separated ligaments. ( 0.454  ). g) Capillary wave growth over the ligaments and drop formations. ( 0.469  ). h) Stable fragments after the core disintegration. ( 0.483  ). (The large volume of the fluid observable in pictures e to d is the upper rim, which can be seen from below after the breakup of the core). ......................................................................................................................................... 43 216Eo  216Eo  * /p pD  216Eo  x 3.18. Deformation of the upper liquid torus and the ring left after the rupture of the flattened core. Corresponding times are 0.462  , 0.520  , 0.551  , 0.582  and 0.610  , respectively. ............................................................................................................................... 44 3.19. Breakup of the liquid bridge made after the deformation of the lower torus. a) the liquid column caused by the deformation of the lower torus. ( 0.666  ). b) the capillary waves grow on the liquid bridge. ( 0.692  ). c) the necks appear due to development of the waves. ( 0.721  ). d) the liquid column disintegrates into several small, stable droplets as well as a few unstable ligaments. These ligaments will be deformed into stable fragments by further breakup or retractions. ( 0.748  ). ......................................................... 45 3.20. Different styles of fragments made after the breakups. a) Stable droplets. b) A small liquid bridge connecting two droplets. c) Deformed ligaments just before the pinch-off and drop formation. d) Larger, deformed fragments. ( 0.886  ). ........................................................... 46 3.21. Variation of number of fragments during the falling prcedure, . ................................. 47 3.22. Variation of number of fragments during the falling prcedure, . ................................. 47 3.23. Variation of number of fragments during the falling prcedure, . ................................ 48 3.24. Droplet size distribution for a) , 32 0.00234665d  , 0.00521aved  b) , 32 0.00187363d  , 0.00444aved  and c) , 32 0.00168808d  , 0.00412aved  . The lines demonstrate the log-normal distribution. .......................................................................... 49 4.1. Droplet disintegration for 144Eo  . a) Initial spherical droplet. b) Fragments generated after the bag breakup and bottom torus deformation. c) Breakup of the droplet core and deformation of the remaining large fragments. d) Cloud of fragments. ...................................... 52 4.2. Fragment number variations for different cases. a) 144Eo  . b) 216Eo  . c) 288Eo  . ................ 53 4.3. Size and velocity distribution of fragments for different stages and Eötvös numbers. a) Diameter distribution for 144Eo  . b) Velocity distribution for 144Eo  . c) Diameter distribution for 216Eo  . d) Velocity distribution for 216Eo  . e) Diameter distribution for 288Eo  . f) Velocity distribution for 288Eo  . ........................................................................... 56 4.4. Variations of Sauter mean diameter, *32d and 43 32d (Square and circle markers, respectively) for different values of Eo . Solid (red) line depicts the empirical relationship. a) 144Eo  b) 216Eo  c) 288Eo  . ................................................................................................................ 58 4.5. Histograms of a) size distributions and b) velocity distributions for last step of each simulations. .................................................................................................................................. 59 4.6. Volumetric fragment size distributions for the simulations. Left axis: V for 144Eo  . Right axis: V for 216Eo  and 288Eo  . ................................................................................... 61 144Eo  216Eo  288Eo  144Eo  216Eo  288Eo  xi 4.7. Scatter plots of velocity-diameter for different simulations and stages. Lines demonstrate the conditional probability for each stage of simulation. Top: 144Eo  , Middle: 216Eo  , Bottom: 288Eo  . ......................................................................................................................... 63 4.8. Eo versus Oh for the last step in different simulations. ............................................................. 65 4.9. We versus Re for different simulations. Diameters of points are proportional to diameter of fragments. Top: 144Eo  , Middle: 216Eo  , Bottom: 288Eo  . ............................................. 67 4.10. We versus Oh for different simulations. Diameters of points are proportional to diameter of fragments. Breakup modes are shown based on the experimental studies of Hsiang and Faeth (1998, 1999). Top: 144Eo  , Middle: 216Eo  , Bottom: 288Eo  . ......................... 69 5.1. Deformation of the droplet in stream for different tested cartesian grids ( 100We  ) at time= 0.97  . Top: Uniform grids of 512 512 (Level 9). Middle: Uniform grids of 1024 1024 (Level 10). Bottom: Adaptive grid based on quad-tree meshing. The maximum level of the adaption is 12. Vorticity field is normalized using * 0 0D U  . Where rot  u u . ....................................................................................................................... 73 5.2. Velocity profiles on the poles for 100We  . h is the height from the interface at pole. ............. 75 5.3. Wave-lengths close to the pole for different tested grids. ............................................................ 75 5.4. Numerical results for a 2D simulation ( 100We  ). Arrows demonstrates the observed scenarios for the ligaments formed due to growth of shear instablities. a) Initiation of the instabilities around the poles due to higher magnitude of the relative velocity. b-f) Amplification of the instabilities in the stream generating ligaments. g-i) Further development (stretching) of the instabilities (ligaments) under the influence of the secondary flows and main stream. j-l) Fragmentation of the stretched ligaments and droplet (fragment) formation. m) Magnified view of d. n) Magnified view of l. Recirculation zones are highlighted by black circles. ................................................................. 77 5.5. Variation of width of the parent droplets versus time for differentWe . ...................................... 79 5.6. Motion of the droplets in stream-wise direction for differentWe . .............................................. 79 5.7. Interface deformation and growth of the shear instabilities for a) 50We  . a) 100We  . and a) 200We  . Snapshots are all in 0.27  and are just shown for the upper half-circles. ............ 80 5.8. Wave-lengths close to the poles for different We . ...................................................................... 81 5.9. Deformation of the droplets at 1.64  for Left: 50We  , Middle: 100We  , and Right: 200We  . ........................................................................................................................... 81 xii 5.10. a) A schematic view of a droplet with diameter of 0 02D r in a stream. To convert the cartesian framework to spherical polar coordinate, let 2 2 2 2r x y z   , so that cosz r  . The area in the dotted box is enlarged in picture b. b) Gas flow around the droplet at the beginning of the streaming before the deformation. c) Possible scenario for the primary breakup of the droplet where the relative velocity, g dU U U   is large enough. . ............................................................................................................................ 83 5.11. Wave-numbers over the droplet from equation 5.6. Left: Influence of Weber number. Right: Effect of Reynolds number (curves are for 3D cases). .................................................. 84 5.12. Linear velocity profile used in Kelvin-Helmholtz instability analysis. ...................................... 86 5.13. Droplet deformation for 100We  and time 1.9  for the 3D simulation. a) Droplet interface and the generated grids. b) Adaptive meshed and the vorticity field. c) Magnified picture of the dotted box in b. ................................................................................................... 88 5.14. Generation and development of instabilities over the droplets in the 3D simulation. Time is from 0  to 1.92  a) Initial spherical droplet. b) Onset of shear instabilities close to the pole. c) growth of shear instabilities in the stream-wise direction. d) Azimuthal modulation in the plane normal to the flow direction. e-f) Growth of shear and azimuthal waves. g) Formation of ligaments and further development of instabilities. h-i) Further stretching of ligaments. ............................................................................................................. 89 5.15. a) Radial (transverse) growth of instabilities for planes perpendicular to the external stream. Time is 0.58  corresponding to figure 5.14-e. ........................................................ 90 5.16. Short-wave-length RTI generating small fragments. The box in the corner demonstrates the magnified view of the same picture. ................................................................................... 91 5.17. 3D patterns on wake of the droplet. Lines show the streamlines at 1.92  . .......................... 92 5.18. a) The growth of KHI over the droplet. Blue dotted lines show the arbitrary cutting planes. b-d) possible deformation of the interface in the radial direction due to Rayleigh-Taylor Instability. A uniform spatial perturbation growth with 4, 6, and 9 spikes are schematically shown here. ............................................................................................................................... 93 5.19. Wave-number distribution around the droplet. a) The effect of amplitude of the primary (shear) instability. b) The effect of upstream velocity. ............................................................. 94 xiii Acknowledgements This thesis arose in part out of two years of research that has been done since I came to The University of British Columbia, Okanagan Campus. This dissertation would not have been possible without the guidance and the help of several individuals. It is a pleasure to convey my gratitude to them all in my humble acknowledgment. First and foremost I offer my sincerest gratitude to my supervisor, Dr Kian Mehravaran. He provided me persistent encouragement and support in various ways. His truly scientist insight has made him as a constant oasis of ideas and passions in science, which exceptionally inspire and enrich my growth as a student, a researcher and a scientist want to be. I am indebted to him more than he knows. It is an honor for me to have Dr. Mina Hoorfar, Dr. Kenneth Chau, and Dr. Sunny Li as my internal advisory committee members and Dr. Rebecca Tyson for their encouraging words, thoughtful criticism, and time and attention during busy semesters. I would like to extend my sincerest appreciations to the School of Engineering that has provided the support and equipment I have needed to produce and complete my thesis. Beyond science, in my daily life, I have been blessed with a friendly and cheerful group of fellow students in UBC Okanagan. I am ever grateful to all my friends who always support me. My parents deserve special mention for their inseparable support and prayers. My Father, Sohrab, who was my first physics teacher showing me the joy of physics and mathematics ever since I was a child. My Mother, Mahrokh, who was my first biology teacher, is the person who put the fundament my learning character. They sincerely raised me with their caring and gently love. Sahar and Donya thanks for being supportive and caring sisters. I am always proud f you. I owe my deepest gratitude to Salma for her unflagging love that means everything for me. Finally, I wish to thank all who contributed in the completion of my Master Thesis. xiv Dedication To my parents and Salma who supported me each step of the way. 1 Chapter 1 Introduction 1.1 Fragmentation of Droplets The deformation and breakup of liquid droplets is encountered in a wide range of industrial applications as well as in natural situations. Engineering examples are found in diesel engines and other types of combustion applications, electro-sprayed paint and cosmetics, ink-jet printers, turbines, micro-fluidics, cooling systems, etc. In the case of geophysical phenomena, examples abound, ranging from volcanic eruption and tephra formation to rain phenomenon. A thorough description of fragmentation phenomena for a variety of applications can be found in Villermaux (2007). Depending on the circumstances a droplet is dispersed in, the mechanism of the two-phase system can be altered from moving in non-stationary surrounding media and free-fall in a quiescent environment due to the constant gravity body force. Schematic views of these two situations are shown in figure 1.1. Fig. 1.1. a) Droplet in an external flow. b) Droplet falling in a constant gravity (acceleration). No amount of experimentation can ever prove me right; a single experiment can prove me wrong. Albert Einstein (1879 –1955) 2 Shock tube and drop tower tests are examples of circumstances a and b in figure 1.1, respectively. The main difference between the two mentioned situations is in velocity change. In contrast to step velocity change in shock tube tests, droplet accelerates more gradually in drop towers due to gravity. However, for the both cases, the mechanism of fragmentation is quite similar. This mechanism can be described using the schematic diagram shown in figure 1.2. Fig. 1.2. The algorithm of droplet disintegration and fragments distribution. Depending on physical properties of the droplet-surrounding media combination, the initial spherical droplets can be stable or not. If the initial droplet is relatively stable, it remains spherical or slightly deforms into spheroid or cap shapes. In the case of unstable droplets, breakup process occurs. However the type of fragmentations might be different depending on physical properties of the system. To understand the droplet dynamics, several experimental, theoretical, and numerical studies have been performed. The majority of the investigations are performed experimentally using different techniques including shock tubes, continuous jets, drop towers, and hybrid methods. Pilch & Erdman 3 (1987) described the breakup modes separately for a single drop experiencing an external flow under different ranges of the Weber number (see Table 2.1). Based on this experimental study (followed by further investigations), droplet breakup can be categorized mainly in five modes: vibrational mode, bag mode, Sheet thinning mode, catastrophic, and Multi-mode breakup. Figure 1.3 demonstrates the schematic view of four main breakup modes. Fig. 1.3.. Main modes of droplet breakup. O’Brian (1961) and Burger et al. (1983), studied the free-fall of oil and gallium drops, respectively and presented the different modes of breakup for liquid-liquid two-phase systems. According to their results, it could be understood that the classification described by Pilch & Erdman (1987) is also valid for liquid-liquid systems. Besides the work of Pilch & Erdman (1987), several experimental investigations are performed to classify the breakup regimes. Krzeczkowski (1980) made a study of deformation and fragmentation of liquid droplets due to an external air stream. Using image processing, Krzeczkowski (1980) introduced four main breakup modes in terms of Weber number consisting of bag, bag-jet, transition, and shear modes. He also studied the breakup duration and the dependency of this parameter on the 4 viscosity ratio. He found that the viscosity ratio does not affect the mechanism and breakup period notably and the most dominant factor is the Weber number. Chou et al. (1997), Chou and Faeth (1998), and Dai and Faeth (2001) carried out a series of experimental studies on temporal properties of shear, bag, and multimode breakup regimes. Several factors such as deformation, droplet velocities, breakup time, critical Weber number and drag coefficient, etc. are systematically presented using the shadowgraph method. The multimode regime in the work of Dai and Faeth (2001), is divided into two main regimes called bag-plume and plume-shear modes. Cao et al. (2007), using high speed camera and shadowgraphs, identified a new mode close to multimode regime called dual- bag breakup for a liquid droplet in a uniform air jet flow. The mechanism of this mode consists of two bag fragmentations of the main drop and the core drop generated in the first breakup. It was postulated that the dual-breakup is different from the observation of Krzeczkowski (1980) and Dai and Faeth (2001) for the similar Weber range obtained in shock tubes. It was argued that the dual- breakup mode is due to dissimilar kinds of disturbances made in uniform jet flow. Cao et al. (2007) mentioned that bag or bag-stamen modes can be seen for the lower limit of dual-bag regime, as well as shear or explosive regimes for higher limit. Recently Guildenbecher et al. (2009) reviewed the technical literature for experimental methods and morphology of droplet atomization for Newtonian and non-Newtonian fluids. They highlighted that ―the mechanism leading to fragmentation of the bag is not well understood‖ while the study of this mode is particularly essential because it ascertains the criteria for secondary atomization. They also report a map/table as a function of Weber number for different breakup modes and presented a relatively comprehensive description on what we know about the mechanism of each mode. Recognizing various breakup modes for droplet fragmentation is the basic step for analysing the mechanism of breakup. However, what is seen in the experiments for highly unstable regions demonstrates that modes are combined and multi-mode-fragmentation happens. Consequently in most cases, shear, bag, wave crest, and other modes are included in the entire process of fragmentation. To demonstrate the involvedness of the topic, the different breakup regimes of five recent studies are demonstrated in figure 1.4. 5 Fig. 1.4. Breakup modes in terms of Weber number. The overlapping of breakup modes occurs for a wide range of Weber numbers. This issue is especially remarkable for the range of 20 80We  , where different types of breakup are reported. It should be noted that the studies mentioned here used shock tube or wind tunnel experiments. In contrast to drop deformation and breakup in wind tunnel or shock tubes, a smaller number of studies can be found investigating falling droplet fragmentation. Reyssat et al. (2007) experimentally explored an unstable range of water drops falling in air. They demonstrated an identifiable, specific class of bag breakup for drops with diameter much larger than capillary length. The distorted shape of such drops consists of a thin aqueous sheet surrounded by air and bounded at the bottom by a torus. During the falling process, as the instability intensifies, the sheet ruptures and a deformed toroidal rim comes into sight. This high-speed fragmentation is followed by destabilization and breakup of liquid bridges due to Savart-Rayleigh-Plateau instability. In the end, a very large drop falling in air decays in several stable fragments of individual size smaller than the capillary length. Reyssat et al. (2007) extended their work for very large water drops, larger than 3 cm, and showed that multiple bags are formed because of Rayleigh-Taylor-like instabilities. A similar scenario including inflation, disintegration and destabilization of liquid bridges occurs for each of the unstable fragments. A comparable observation of large water drop fragmentation can also be found in the recent paper of Villermaux & Bossa (2009). Using a hybrid system including falling tower and ascending stream of air, the authors showed that the flattening of the drop into a pancake shape is followed by bag 6 creation. Afterwards, the similar fragmentation cascade observed in the work of Reyssat et al. (2007) is perceived. Villermaux & Bossa (2009) stated an imperative argument on uncertainty of the experiments performed. They showed that the formation of a bag is not obligatory for the initial drop to break up. In some cases, the falling drop deforms into a pancake shape, the torus is made, and liquid columns break up into stable daughter droplets. This issue shows the chaotic nature of the problem and the concern of how much the deformation modes are related to the operating condition of the experiments. Villermaux & Bossa (2009) scrutinized the atomization process by calculating the bag growth and distribution of drops. They declared that for moderate Weber numbers ( 6We  ), the drop radius oscillates around a mean. However, for large Weber numbers ( 6We  ), the radius increases exponentially, approaching the bursting moment. Subsequently, the distribution of produced drops is evaluated and compared against Marshall-Plamar exponential distribution of rain drops. As a momentous consequence, it is found that ―single-drop fragmentation determines size distribution of raindrops.‖ In the case of numerical simulations, the studies vary from small deformation of droplets to breakup of dispersed phases. Nevertheless, highly unstable situations still need more investigations. Using the Lattice Boltzmann method, Fakhri & Rahimian (2009, 2011) presented numerical simulations of descending droplets for a variety of Ohnesorge, Archimedes and Eötvös numbers (see Table 2.1). All simulations were performed in 2D or axi-symmetric domains. Different shapes of breakup were observed and the results of stable, non-fragmented drops were compared with theoretical drag formulations for disks and spheres. Some previously established facts like the increase of instabilities by enlarging the droplet Archimedes or Eötvös numbers were also shown in the results. At the most unstable condition of the study, four fragments are observed after the breakup. Although the proposed numerical method is able to capture the shear and bag breakups, no results are presented for highly unstable drops, which are essential for a wide range of applications. A comparable attitude is followed in the direct numerical simulations of Ni et al. (2006), while the influences of wall repulsion are also studied. They also simulated the coalescence of two drops but again no information was provided for the breakup of drops. Recently, Feng (2010) studied the steady-state flow fields inside and outside of a deformable liquid drop falling at its terminal velocity using a Galerkin finite element method. Different viscosity and density ratios as well as Reynolds and Weber numbers are examined, and the drag coefficient is calculated. Although the considered drop was deformable, the breakup was not studied. A similar approach can be seen in the work of Gottesdiener (2004) where the velocity of a stable falling drop is obtained for different situations and the breakup mechanism is out of the scope of the study. Han and Tryggvason (1999), using a finite difference/front tracking method, studied the 7 unsteady motion of falling liquid drops for different fluid-fluid combinations. Results are classified for two density ratios, 1.5 and 10, and different Eötvös and Ohnesorge numbers. Bags are shown in the simulations for the range of 144Eo  , and results are summarized in six Eo Oh shape maps. Similar to most of the other numerical simulations, outcomes of Han and Tryggvason (1999) are limited to deformed shapes and drop velocity variations and no results were presented for breakup. 1.2 Motivation Due to the complexity and chaotic nature of the problem, including very small length and time scales, three dimensional patterns, gas-liquid interaction, etc., theoretical studies normally require simplifying assumptions and approximations. Thus, a purely theoretical explanation of the deformation and breakup of the droplets has not, and is not expected to be verified in the near future. The paper of Harper et al. (1971) may be the most advanced theoretical investigation of this field. Using numerical as well as second-order asymptotic methods, the transient response of an externally excited droplet was explored for a range of Bond numbers. Harper et al. (1971) organized a regimes- map for the dynamic response of drops to an external flow as a conclusion. However, the deformation patterns obtained are far from the reality. Theoretical studies are crucial and fundamentally important for interfacial science. Nevertheless the mathematical techniques are still ineffective in solving the chaotic problems such as fragmentation of liquid drops. In the case of experimental investigations, due to three-dimensional occurrences in very small temporal and spatial scale (ms and μm, respectively), the current measurement techniques (from shadowgraph to high-tech Particle Image Velocimetry) are not accurate enough to show the details of the various phenomena during the fragmentation and capture all the characteristics of droplet fragmentation. Thus, empirical results are normally limited to deformation shape and breakup modes. The story of numerical simulations is different. By using a careful and accurate approach considering the properties of the governing equations, even chaotic problems such as turbulent two-phase flows can be resolved. The limits of numerical simulations are defined by computational costs. In other words, an excessively large number of finite discretized volumes and extensive computational times are required to solve these complex problems. Due to limitations described above, the mechanism of liquid droplet breakup in small scales has not yet been clarified. However, as mentioned in the review paper of Guildenbecher et al. (2009), Direct Numerical Simulation (DNS) of Navier-Stokes equations combined with interface tracking is the best known technique for the study of the process. 8 The present study, to our knowledge, is the first attempt for DNS of liquid fragmentation, aiming to clarify the mechanism of fragmentation and give a detailed description of fragment properties. In this thesis, the deformation and breakup of droplets for two cases of falling droplets in constant acceleration and droplets in an external stream are investigated. Direct simulations are performed based on appropriate numerical algorithms in the context of the volume of fluid (VOF) method, combined with high performance computing techniques. 1.3 Outline of This Thesis This thesis primarily discusses the mechanism of fragmentation using DNS. It starts with the details of the numerical method we use for the simulations in chapter 2. The next chapter (chapter 3) contains the details of the results obtained for the deformation and mechanism behind the breakup of falling droplets in constant acceleration. In chapter 4, a comprehensive analysis for characteristics of fragments, generated in the breakup of falling droplets, is presented. Chapter 5 deals with deformation of droplets in a stream. Our studies for this topic are limited to generation and growth of instabilities over the droplets. Subsequently, a summary of the results and conclusions are given. 9 Chapter 2 Methodology & Non-Dimensional Groups 2.1 Geometries As described in chapter 1, two cases are studies in the current thesis: 1) a free-falling droplet and 2) a droplet in a stream. For the first case, a spherical droplet of initial diameter D with density and viscosity of d and d , respectively is considered to fall under the action of gravitational acceleration g , with zero initial velocity,  * * 0 0 t V   . The quiescent carrier fluid has a viscosity of c and density of c . The computational domain is schematically illustrated in figure 2.1. The whole of science is nothing more than a refinement of everyday thinking. Albert Einstein (1879 – 1955) 10 Fig. 2.1. Schematic sketch of the initial condition of simulations for a falling droplet. The dimensions of the domain and initial position of the droplet are placed so that the boundaries exert negligible influence on the fragmentation process. Thus, the conditions of free fall are satisfied. No external flow is exerted in the domain and the droplet descends due to constant gravity body force. All of the simulations for this case are performed in three dimensions. For the second case (single droplet in a stream), a periodic box with the size of ( 0 0 05 5 5D D D  ) is considered and an initially spherical droplet is placed at the center of the box. The physical problem and the computational domain for the 3D simulations are illustrated in figure 2.2. 11 Fig. 2.2. Droplet in a stream. Flow direction is from left to right. For this case, both 2D and 3D simulations are considered, since the instability initiation and growth of instability in the stream-wise direction can also be captured in 2D simulations. For the 2D simulations the dimensions are the same but only z and x directions are considered. Thus, the droplet has a circular geometry for the 2D case. The gas part of the domain is initialized with a uniform initial velocity, 0U , and the droplet is stationary. It should be mentioned that the passage of the shockwave cannot be presented due to incompressible flow solver used here. However, after a very short time (no deformation), a smooth divergence-free velocity field is obtained. 2.2 Governing Equations The incompressible, variable-density, iso-thermal Navier-Stokes equations with surface tension govern the motion of the dispersed and continuous phases. The set of equations can be written for the computational domain as, 12 . 0, u (2.1)  . 0,t   u (2.2)    . . 2t sp          u u u D n g (2.3) where ( , , )u v wu is the velocity vector, and ( , )t  x and ( , )t  x are respectively the local fluid density and the dynamic viscosity.   2ij i j j iD u u   D is the deformation tensor. The Dirac delta ( s ) ensures that the surface tension term is concentrated on the air/droplet interface. In equation 2.3,  is the surface tension coefficient, and  and n are the curvature and normal vector with respect to the interface, respectively. The density and viscosity are calculated based on the volume fraction of the droplet phase,  ,c tx ,    1 ,d cc c c     (2.4)    1 .d cc c c     (2.5) The advection equation for density can be replaced with an equivalent equation for the volume fraction,  . 0.tc c  u (2.6) 2.3 Flow Solver (Gerris) Gerris is an open source code for the solution of the Navier-Stokes equation. The source code is available free of charge under the Free Software GPL license. Gerris is supported by NIWA (National Institute of Water and Atmospheric research) in New Zealand and Institut Jean le Rond d'Alembert in 13 France. The code is relatively new, however It solves the time—dependant incompressible variable- density Euler and Navier-Stokes equations. The main feature of the code (and the main reason we use this code) is the adaptive mesh refinement (i.e. the resolution is adapted dynamically to the features of the flow). This feature is supported with second-order discretisations in space and time. For interface tracking, a Volume of Fluid (VOF) advection scheme is employed in the code. 2.3.1 Temporal Discretisation A second-order accurate time discretisation is obtained as follows, using a staggered in time discretisation of the parameters:      1 111 1 1 1 1 1 222 2 2 2 2 . . ,n n n n s nnn n n n n p t                                 u u u u D D n g (2.7)   1 1 2 2 0, n n n n c c c t       u (2.8) 0.n u (2.9) Using the Chorin classical time-splitting projection method (Chorin 1969) the system simplifies to,      * 111 1 1 1 * 222 2 2 2 . . ,n n s nnn n n nt                              u u u u D D n g (2.10)   1 1 2 2 0, n n n n c c c t       u (2.11) 1 * 1 1 2 2 ,n n n t p        u u 1 0.n u (2.12) This system of equations requires the solution of Poisson equation for pressure: 14 1 * 1 2 2 . n n t p                 u 1 0.n u (2.13) The velocity advection term in equation 2.10 is estimated using the robust second-order upwind scheme of Bell et al. (1989). This scheme is stable for CFL numbers smaller than one. To solve equation 2.11, a piecewise-linear geometrical Volume Of Fluid (VOF) scheme is used. Interface reconstruction and interface advection are the two main steps in the geometrical VOF scheme. The details of these two steps can be found in the recent paper of Popinet (2009). 2.3.2 Spatial Discretisation A quadtree(2D)/octree(3D) spatial discretization method is used. An example of this technique for a 2D box is illustrated in figure 2.3. Fig. 2.3. Quad-tree grid adaption from level 0 (dotted box) to level 4 (green boxes). In figure 2.3, the root cell is shown using dotted lines. The root cell has a Level of zero and the level increases by one for each successive generation (each parent cell can have zero or four children cells for 2D simulations and zero to eight children in 3D) in the tree. The adaption is continued up to the 15 leaf cells (cells without any children). It should be mentioned that the level of neighboring cells (directly or diagonally) cannot differ by more than one. The refinements can be done using different criteria. For the falling droplet simulations, since an accurate prediction of interface deformation and the influence of the carrier phase on droplet dynamics are essential to the breakup process, the radius of curvature and vorticity magnitude are considered as refinement criteria. The vorticity magnitude (  curl u ) is also important for resolution of the local turbulent eddies generated in continues (gas) phase and close to the interface. For the droplets in a stream, the gradient of the stream-wise velocity is also considered for the adaption, as it significantly assists in the resolution of the boundary layer over the droplet. For the current simulations, adaption is followed reaching level 12. The octree mesh permits a saving up to two order of magnitude in computational cost, making these simulations feasible. The adaptive refinement, based on the octree mesh, is relatively inexpensive, and sensible to execute in every time step. An example of grid refinement for the centre-plane of a 3D simulation of a falling droplet is illustrated in figure 2.4. Fig. 2.4. Grids generated for a case of falling droplet ( 0.761  for 216Eo  , 0.05dOh  , 0.05cOh  , * 10  ). Left: interface of the drop. Right: cross-section of the drop in 0z  plane; each cell is colored by the refinement level-number of the grid. As shown in figure 2.4, the grids are refined at the interface of the drop and the region around and on top of the droplet because of the curvature and vorticity gradients in these areas. The iso-surfaces for the grid sizes in a falling droplet simulation (same as figure 2.4, but after the fragmentation) are also illustrated in figure 2.5 to show the quality of the grid adaption in three dimensions. 16 Fig. 2.5. Iso-surfaces of grids size for fragments after the fragmentation of a falling droplet. a) Dispersed fragments with different shapes for 216Eo  , 0.05dOh  , 0.05cOh  , * 10  . b) Iso- surface for grid size of, * 0 0.5dldl D   . c) * 0.2dl  . d) * 0.1dl  . e) * 0.05dl  . As is shown in figure 2.5, the fragment cloud is covered layer by layer with different levels of refinement. The refinement levels are selected such that very small length scales of fragments (up to 1/1000 of the initial droplet diameter) can be captured after fragmentation. Figure 2.6 demonstrates the grid size Iso-surfaces for the case of single droplet in a stream. Fig. 2.6. Iso-surfaces of grids size for a deformed droplet in a stream. a) * 0.1dl  . b) * 0.05dl  . c) * 0.02dl  . d) * 0.005.dl  . 17 A complete description of the discretisation and numerical schemes used can be found in Popinet (2003) and Popinet (2009) for Euler and Navier-Stokes solutions, respectively. The open-source code used, ―Gerris‖, has been accurately validated for conservation errors and against the linear instability theory by Popinet (2003). For the case of surface tension driven fluid flows, the code is validated for different common problems; for example, generating a sinusoidal perturbation over a flat interface between two fluids initially at rest, Popinet (2009) studied the amplitude of the capillary waves. The results were compared with the exact solution of Prosperetti (1981) and a very good agreement was observed. Moreover, comparing the relative errors, the higher efficiency of the method in comparison with front-tracking methods, etc. is shown. The ability of the code is also checked for two other important problems which are essential for the current topic: liquid sheet retraction, Fuster et al. (2009), and the Rayleigh plateau instability of a liquid column Popinet (2009). For the latter case, obtained results of the interface deformation are compared with the theoretical solutions of Rayleigh (1982) and Weber (1931) for inviscid and viscous flows, respectively, and very good agreement was observed. 2.4 Non-dimensional Numbers A wide range of material properties and physical phenomena affect the droplet fragmentation process. In the case of a falling droplet in a motionless medium, the acceleration due to gravity, momentum, surface tension between the phases, and viscosity are the main parameters. Considering these parameters, different non-dimensional groups are presented by the researchers which govern the behaviour of the droplet. A list of important dimensionless numbers can be found in table 2.1. . 18 Table 2.1. Non-dimensional numbers studied for the case of falling droplets. Dimensionless number Symbol Definition Range Eötvös number Eo   2d cg D    144-216 Ohnesorge number (based on droplet properties) dOh   1 2 d d D    0.1 Ohnesorge number (based on continues phase properties) cOh   1 2 c c D    0.1 Weber number We 2 c Dv  - Reynolds number Re c c Dv  - Viscosity ratio * d c   10 Density ratio * d c   10 Depending on the situations, some of the numbers listed in table 2.1 are more meaningful than others. In the circumstance of descending liquid droplets, using Re and We to categorize the droplet behaviour is not suggested due to the time dependent velocity of the droplet and surrounding medium (initially zero) as well as the presence of breakup. In fact, Re and We are useful for the stable droplets moving in a given velocity or when the stream velocity of the surrounding medium can be controlled. Thus, we used these two numbers for statistical descriptions of the fragmentation when all of the drops have a firm shape. The current study aims to clarify the mechanism of breakup and the way fragments are created in the atomization process. Thus, all the non-dimensional numbers are fixed, except the Eötvös number. Subsequently, Eötvös and Ohnesorge numbers are selected as the basic dimensionless factors. These numbers remain constant, as long as the droplet is still intact. Due to applicability in real-world applications, Ohnesorge numbers are selected to be less than 0.1 (low Oh region). As Krzeczkowski (1980) found experimentally, the viscosity ratio is negligible for this region. In the presence of breakup, an individual number will be used for each fragment to determine 19 the distribution variation of different parameters in time. Furthermore, to study the temporal procedure, time is normalized with respect to the time scale 30t D   , ( 0 t t   ). Based on the current non-dimensional numbers for falling droplets, an intense fragmentation in upward bag regime is predicted, Han and Tryggvason (1999). The maximum Weber number for the current simulations is in the range of 37 59We  , corresponding to bag breakup of droplets, based on experimental studies. Considering the air-water system, the current non-dimensional numbers translate to initial drop size and falling velocities in the range of 03 6cm D cm  and max8.5 / 11.5 /m s v m s  . Comparing with previous studies, the systems we are investigating for falling droplets are more unstable. Thus, an intense fragmentation is expected as it occurs in several industrial and natural phenomena. The fragmentation that occurs for this range of non-dimensional groups is the main subject of our study. In aero-breakup of droplets, aerodynamic force deforms the droplet and causes fragmentation, surface tension resists the deformation, liquid viscosity delays the deformation and dissipates the stream energy exerted on the droplet, and gas viscosity dissipates the energy of the stream. The non- dimensional groups describing the behavior of the droplet in a stream, accounting the parameters described, in addition to the ranges used here are listed in table 2.2. Table 2.2. Dimensionless groups for the case of a droplet in a stream Dimensionless number Symbol Definition Range Ohnesorge number Oh   1 2 d d D   0.002-0.01 Weber number We 2 0cDU  50-200 Reynolds number Re 0c cDU  22300 Viscosity ratio * d c  10 Density ratio * d c  10 In table 2.2, Reynolds number can also be presented as 0.51 0.5 * *Re Oh We    . For this case, the dimensional time is also normalized using 0 0tU D  . Again, it is significant to say several experiments have depicted that the breakup modes are independent of Ohnesorge number for the range of 0.2Oh  . The range of dimensionless groups we have, based on the experimental 20 observations, predict a shear breakup of droplets. This mode of fragmentation is the most common regime and most complex to analyze. The experimental observations for this range cannot describe the mechanisms at the early stages of the breakup due to the small length and time scales of deformation and fragmentation. Thus, this regime is selected for our numerical simulations. We aim to clarify the mechanism of shear breakup of droplets in a stream focusing on the initiation and growth of instability at the early stages of fragmentation. 21 Chapter 3 Fragmentation of Falling Droplets 3.1 General Description of The Fragmentation Process The scope of the present chapter is to simulate and analyze the breakup mechanism of a falling droplet in a lighter medium through direct simulations. Therefore, the selected region is chosen based on the criterion for the occurrence of breakup. According to the study by Han and Tryggvason (1999), the deformation of falling droplets in low circumstances ( ) can be compartmentalized into four main styles in order of increasing : steady deformation, backward- facing bag, transient mode, and frontward-facing bag. Although no results were reported for the breakup of the drop, it could be expected that atomization will occur for high instability regimes which in this case mean high . Hence, we decided to cover a more unstable region in comparison with the previous studies, to analyze the different modes and phenomena of breakup. It should be noted that when is small, the surface tension is much more essential than viscous stresses. Thus, has little effect on the breakup and is the controlling factor. Subsequently, is fixed and three different are selected, as shown in table 2.1. A sequential picture of droplet deformation and fragmentation for is shown in figure 3.1. At the preliminary stages, as the droplet falls, we see the deformation of the droplet from the initial spherical shape to an oblate ellipsoid followed by bag formation (figure 3.1, a-c). The reason for this deformation is the non-uniform distribution of the hydrodynamic pressure around the droplet which is set into motion by gravity Han and Tryggvason (1999). A similar fashion of deformation for falling drops can be found in the studies of Han and Tryggvason (1999), Ni et al. (2006) and Feng (2010) for deformable drops. The deformation is resisted by the surface tension. Nevertheless, this amount is not enough to hold the droplet in a stable shape (This balance can be presented by the Weber number and Oh 0.05d cOh Oh  Eo Eo Oh Oh Eo Oh Eo 288Eo  All science is either physics or stamp collecting. Ernest Rutherford (1871 –1937) 22 will be discussed further). Then the bag grows and breaks up, forming the ligaments, liquid bridges and daughter droplets. A very thin disk-like core will be left during the bag breakup which will break up in a short period of time. Various types of instabilities such as Rayleigh- Plateau, Rayleigh-Taylor, and capillary wave instabilities occur throughout the atomization process. The procedure will continue and fragments in each stage break up further until all the fragments reach a stable state. Fig. 3.1. The deformation and fragmentation of a liquid drop for: , , , . (a) , (b) 0.1647  , (c) 0.2642  , (d) 0.3575  , (e) 0.4434  , (f) 0.5183  , (g) 0.5924  , (h) 0.6728  . Similar behavior is seen in experimental investigations in wind tunnels or shock tubes. Figure 2.2 represents an example for the fragmentation of a liquid drop in an air stream observed in the experimental study of Cao et al., (2007). 288Eo  0.05dOh  0.05cOh  * 10  0.0  23 Fig. 3.2. Fragmentation of a liquid water drop in a continuous air jet flow. and . Left: : the green arrow notes one of the liquid bridges generated in the experiment. Purple and orange arrows demonstrate the penetrated core and unstable ligaments, respectively. middle: : green arrow denotes the fragments created by disintegration of liquid columns. The purple and the orange arrows show the deformation and breakup of the core, right: : green arrow points to a liquid column made by the breakup of the core. The orange arrow indicates the fragments made in the fragmentation. The purple arrow signifies a relatively large fragment generated during the core breakup. (the original picture is from Cao et al. 2007.) It should be mentioned here that the external flow imposed on the droplet in the study of Cao et al. (2007) was a continuous air flow, which is very close to the situation of falling drops where the droplet, as well as the generated fragments, fall in a continuous range of velocity. After the formation and breakup of a thin bag with a rim and a core, Cao et al. (2007) observed that the core will remain connected with some ribs and ligaments, which are generated due to the disintegration of the bag. The formed liquid columns disintegrate briefly because of the combination of Rayleigh-Plateau and aerodynamic instabilities. Afterwards, the core begins to break up. In this stage, another (smaller) bag is seen, which is related to a high enough of the core imposed by the continuous air stream. The disintegration of the core continues until small, stable fragments form. The significant points are the similarities found between the mechanisms observed in the experiment of Cao et al. (2007) and the simulation made in this study. In both cases, the drop initially deformed generating the bag. Then, the bag is broken up, forming ligaments and liquid columns, as well as small droplets. The remaining core and ligaments continue to disintegrate and the atomization continues until a universal stable state is established. Through general analysis of the results, we decided to classify the study of the fragmentation process in three main categories: bag formation and growth, bag breakup, and droplet number variations and distributions. These categories will be discussed in detail in the next sections. 29We  0.0015dOh  24.3t ms 35.6t ms 42t ms We 24 3.2 Drop Deformation and Bag Formation The initial deformation of the drop will be followed by formation of the bag during the fall. This deformation process is the origin of the breakup. Thus, a complete understanding of the deformation characteristics is needed. This section aims to clarify the mechanism of deformation as well as bag formation and the influence of on the behaviour of falling drops before the deformation. Figure 3.3 is a sequential picture of drop deformation for . The interface of the liquid phase is coloured by the non-dimensional spatial falling velocity (in the y-direction). The velocity component is normalized by , where and . Fig. 3.3. Deformation of an initially spherical liquid drop deformation for: . a) 0.0825  , b) 0.137  , c) 0.175  , d) 0.210  , e) 0.244  , f) 0.278  , g) 0.317  , h) 0.358  . Colours demonstrate the vertical velocity of the interface,   * vV p gD     . Radius and height of the drop are illustrated in picture (h). The gravitational acceleration is in the y-direction. The deformation of the droplet shown in figure 3.3, can be divided into three main parts: slight deformation (figure 3.3 a-c), bag formation (figure 3.3 c-f) and bag growth (figure 3.3 g-h). At the beginning, the pressure difference causes a smooth deformation starting from the rear part of the droplet. As can be seen in figure 3.3 (b and c), the vertical velocity of the upper side (top) of the drop is larger than the other parts. So, the rear side of the drop turns into a flat surface while the front side remains rounded. The amount of this deformation increases due to the accelerating motion of the drop in constant gravity. Therefore, the drop deforms into a bowl-like shape (figure 3.3 d and e). This stage is the beginning of bag formation. Up to this point, the diameter of the drop, increases and Eo 216Eo    1 2 p gD    d c     d cp     216Eo  a 25 the height, h reduces. By formation of the bag, the drop, which is now a continuous liquid sheet with non-uniform thickness, is stretched. Therefore, both and increase as well as the volume occupied by the droplet. In addition, the interaction of recirculation zones, bulk falling velocity and the surface tension causes a convex shape at the end of the bag. As it can be observed in figure 3.3 (h), the velocity in the y-direction for this area is less than in the other parts of the drop. This will increase the stretching and thinning of the sheet. The whole progression is followed by bag breakup. To investigate the deformation and the behavior of the bag before breakup, the drop characteristics are discussed quantitatively. Figure 3.4 demonstrates the variation of the parent drop diameter, normalized by initial diameter of the spherical drop, during the falling process for three Eötvös numbers. Fig. 3.4. Temporal variation of drop diameter during the falling process. denotes the slope of the linear trendlines interpolated (red for the deformation, blue for the bag growth) for different stages. Deformation decreases with reducing , and occurs at a slower rate. This fact can be simply described by comparing the values of buoyant force and surface tension. In the case of higher , the resistant force made by the surface tension is smaller. Thus, the drop deformation is more a h m Eo Eo 26 pronounced. A similar explanation can be used to elucidate the influence of on the rate of radius variation. As seen in figure 3.4, the rate of changes for each step (deformation and bag growth) increases by increasing . Another point is the amplification of rate of radius growth before and after the bag formation. As depicted in figure 3.4, for a given , the rate of radius variation is notably magnified. This fact is due to the morphology of the bag, which actually is a continuous liquid sheet. This liquid sheet, which has a larger surface area, can be deformed much more easily than a bulk liquid drop. Similar behaviour is seen in the shock tube experiments by Chou et al. (1998) and wind tunnel tests by Cao et al. (2007). To consider the parameter in the deformation rates, the aspect ratio is defined here calculating the ratio of temporal radius to height of the droplet. Figure 3.5 demonstrates the variation of aspect ratio versus time for different Eötvös numbers. Fig. 3.5. Temporal aspect ratio of falling droplets. As it was seen in figure 3.4, the radius of the drop increases monotonically. Before bag formation, this happens while the height of the drop decreases. Thus aspect ratio grows as is shown in figure 3.5. Just after bag generation, the height of the droplet increases notably and this causes a plateau in aspect ratio. After that, by significant expansion of the droplet in the cross-stream direction, the value Eo Eo Eo h 27 of the aspect ratio increases again, reaching a peak. The aspect ratio decreases once breakup begins. The deformation of the drop during the breakup will be discussed later. The temporal Reynolds and Weber numbers are calculated to demonstrate the variations of ratio of inertia, to viscous and surface tension forces, respectively. For these cases, the length scales, which are shown by in table 2.1, can be fixed or time dependent, . Both of these conditions are considered and the results are depicted in figure 3.6. Fig. 3.6. Left: temporal variations of Reynolds number. Right: temporal variation of Weber number. For the case of fixed length scales, both and We vary with the velocity changes. For a given , the velocity increases with time reaching a peak and then reduces. The peak generated here is due to formation and growth of the bag. Creation of the bag significantly augments the area of the drop in the cross-stream direction. It will magnify the drag force and cause a drop in velocity. An important result can be obtained here, comparing the curves of different . Increasing the Eötvös number means adding to the ratio of driving force (buoyant force) to surface tension resistance force, which means a higher falling velocity as well as larger deformation. In the case of high Eötvös numbers, which are investigated here, the effect of driving force is much greater than the surface tension force, especially at the beginning of the motion. By formation of the bag, however, we expect higher velocities when we increase , the drag force increases and as seen in figure 3.6, the velocities drop. The peaks are not observed for the Rea (time variant length scale) because of the notable influence of D a Re Eo Eo Eo 28 radius changes, which compensate for the effect of the velocity drop after bag formation. For the aWe (time variant length scale), the peaks are still seen. This is due to the second order variation of velocity in Weber number, which magnifies the role of falling velocity. The other vital outcome that we can deduce from figure 3.6 is the operating range of and , as well as the values of other parameters at the moment of breakup. These values are summarized in table 3.1. Table 3.1. Values of Reynolds and Weber numbers 59.7 98.1 155 276 37.6 84.2 123 276 81.2 138 180 354 47.8 122 125 354 106 174 207 396 59 152 154 396 The point which should be mentioned here is that deformation and breakup are continuous phenomena. In other words, for falling drops at the moment of , the drop is completely deformed. Thus, calculating the Weber and Reynolds numbers based on a time varying radius will help to include the effect of deformation in the estimation of breakup regimes. To calculate the ratio of buoyant to inertial forces during the falling process, the following coefficient is determined using, d g V C     , 2c v A (3.1) where v is the falling velocity, V is the volume of the initial droplet and A is the area projection in the direction of fall. Figure 3.7 demonstrates the variation of versus time for the three simulations. Re eW max DWe max aWe max DRe max aRe . .B U DWe . .B U aWe . .B U DRe . .B U aRe 144Eo  216Eo  288Eo  . .B U DWe We dC 29 Fig. 3.7. Drag coefficient vs. time. For a better understanding, equation 3.1 can be rewritten in terms of , and the aspect ratio: (3.2) As discussed above, both the Weber number and aspect ratio increase as the droplet falls. Therefore, a drop in with time is expected for all cases. Comparing the different Eo cases, although increases proportionally with Eo , the product of increasing parameters in the denominator (with Eo ) cause an overall drop in the curves with respect to Eo . After the significant deformation discussed so far, the droplet starts to break up, which is the subject of the next section. DWe Eo 2 4 3 d D Eo C a We D        dC dC dC 30 3.3 Bag Breakup The behavior of the droplets before breakup, including deformation and bag generation, was discussed in the previous section. It was seen that the falling droplet significantly deforms into a highly unstable, continuous thin sheet. This unstable, curved liquid sheet disintegrates in a relatively short time period. Generally speaking, the results demonstrate that the breakup procedure consists of three main parts. 1) Initially the bag bursts, generating droplets and filaments as well as the upper rim and a flattened core (figure 3.1 d). 2) The core deforms into another torus (figure 3.1 e). 3) The rounded or straight liquid columns disintegrate during the falling process generating smaller fragments (figure 3.1 g). These three parts are the subjects of this section. The bursting and retraction of a liquid sheet are classical topics, investigated analytically through instability analysis. After the pioneering studies of Rayleigh (1883) and Taylor (1950), who demonstrated the instability of superimposed variable density fluids, Keller & Kolodner (1954) extended the solution of Taylor (1950) by taking the influence of surface tension into account. Keller & Kolodner (1954), using a first order perturbation, calculated the rate of growth of perturbations over the interface. They show that for the case of constant acceleration, this amplification would be exponential for the most unstable modes, which are considered responsible for the disintegration of the sheet. The dispersion relation obtained by Keller & Kolodner (1954), written in a non- dimensional form in the recent study of Bermond and Villermaux (2005) is (3.3) where is the dispersion of the interface, denotes the wave number of the mode and   1/ 2 2 /c dk gh  is the capillary wave number. Explaining equation 3.3 in terms of the shape function of the interface displacement, Bermond and Villermaux (2005) show that (dependent on the thickness and the magnitude of the acceleration) the film modulates its thickness and is subsequently punctured by numerous holes in different positions. Due to capillary force, the holes generated at the interfaces will grow radially, attaching the other punctures in the vicinity. This amalgamation process will result in a network of attached ligaments, further disintegrating into several droplets. These      1/ 2 1/ 2 43 2coth 1 1 1 / tanh ( )ck k k k k                  k 31 processes have been demonstrated in experiments performed in the shock-tube for three different Mach numbers. Figure 3.8 shows the holes generated for different Mach numbers in the experiments by Bermond and Villermaux (2005). Fig. 3.8. Experimental results by Bermond and Villermaux (2005) for bursting liquid films for a) (the arrow shows a ligament formed after attaching of two growing holes), b) , c) , scales in each picture show size of a selected hole. (the original picture is from Bermond and Villermaux, 2005) The number of the punctures increases in time for all cases. The important point here is related to the variation of the number and perforation rate of holes for different Mach numbers (acceleration). As can also be derived by theory, this parameter increases with acceleration or Mach number. However, the mean size of the holes notably reduces by adding to the magnitude of the imposed acceleration. 1.03M  1.07M  1.21M  32 Despite the apparent differences between a liquid sheet subject to sudden acceleration and a falling droplet deformed into a bag, the underlying mechanism is similar.. For the current study, the variation of acceleration can be shown in Eötvös number instead of Mach number. It is expected that phenomena similar to the ones observed in the experiments of Bermond and Villermaux, (hole generations, ligament formation, and breakup into small fragments) happen. Furthermore, a higher number of holes and an increased perforation rate is anticipated for larger Eötvös numbers, as investigated below. The results obtained in the simulations showed that the behaviors of the droplets in the three simulations are different while the mechanisms are the same. In all cases, several holes are generated on the bag. The holes grow forming the ligaments. The ligaments then disintegrate and generate small fragments, where most of them are stable. Remaining unstable fragments will further crumble after some time. By increasing the Eötvös number, the number of generated holes increases, as well as the growth rate. Figure 3.9 demonstrates the generated holes and their growth for and . Fig. 3.9. Top views of hole growth for (a-d, ): a) 0.414  , b) 0.427  (the red lines depict the boundaries of the holes), c) 0.447  (the transparent red circle shows the area of the upper torus), d) 0.463  . (e-h, ): e) 0.267  , f) 0.278  , g) 0.289  (the transparent red circle shows the area of the upper torus), h) 0.300  . 144Eo  288Eo  144Eo  288Eo  33 The first issue which should be mentioned here is that, from this point (hole generation and growth) deformation and disintegration are three dimensional, and two-dimensional or axi-symmetric simulations are not able to correctly predict the growth behavior. As can be seen in figure 3.9, holes are formed in different positions, times and shapes. This process is followed by formation of ligaments with different characteristics. Our analysis shows that the wave lengths generated over bag breakup as well as core breakup (will be shown later) are in good agreement with the range of most unstable wave lengths obtained by equation 3.3. For instance, for the droplet in , where the average capillary wave number is approximately 0.4ck  , the wave lengths observed in simulation are in the range of 0.17 0.31kh  , close to the   max 0.2kh  obtained by theory. By comparing the pictures of and simulations in figure 3.9, it is obvious that the higher Eo case has more holes. To investigate the mechanism of hole formation and growth in more detail, the simpler case of is further studied. Then, the bag breakup in the more complex case of will be studied. As was shown in figure 3.9, fewer holes are generated in the case of . By enlarging and changing the view in figure 3.9-c and figure 3.10, it is seen that the holes are formed mostly closer to the upper rim. Fig. 3.10. Holes positions for . a) Three dimensional shape of the penetrated drop, 0.447  . b) A cross-section of the drop in the same time: the area denoted by the red circles is the liquid torus. The area noted by the blue torus (circles) is the area which the majority of the holes are initially generated there. Arrows point to the high-curvature part of the interface. 216Eo  144Eo  288Eo  144Eo  216Eo  144Eo  144Eo  34 The thickness of the bag in the high curvature area is less than the other regions. Therefore, the holes are formed close to the upper torus. Subsequently, the holes grow until filaments form between them (figure 3.9 –d). By expansion of the punctures, these filaments morph into liquid bridges, connecting the core of the drop to the upper torus. Eventually, these liquid bridges disintegrate into small, stable fragments, leaving the flattened core of the drop, as well as the liquid torus. At this moment, the bag is completely broken up and a new part of the breakup will begin. Before discussing the torus and core fragmentation, the hole expansions and liquid bridge formations, which are the main parts of the bag breakup will be explained for the case. Since the number of punctures in the case of is small, investigating the hole enlargement and ligament formation is easier. Thus, we have selected a hole to show the quality of the expansion. In addition, a neighbouring hole is also considered to show how the liquid bridge is formed. These developments are shown in figure 3.11. 144Eo  144Eo  35 Fig. 3.11. a) The holes generated on the drop for at 0.427  . Cross-sections are plotted with different colours. b) Expansion of the selected hole. Corresponding times are, 0.411  , 0.419  and 0.431  , respectively. The rounded edge (rim) of the hole is obvious for the thinner part. The sheet has non-uniform, time varying thickness. (The last time step is magnified in a separated box and demonstrates small amplitude capillary waves in the thinner part. These waves, as well as the round edge are not seen in the thicker part of the hole). c) Growth of two neighbourhood punctures and formation of the liquid bridge. Corresponding times are, 0.419  , 0.431  , 0.444  and 0.463  , respectively. The rounded edges again are evident for the plane which connects to holes. Some very weak capillary waves are seen for this plane while they are not observed for the left and right hand side parts that actually show the upper rim. Holes grow until a liquid column, shown by a circular cross-section, is formed. This column in fact connects the upper rim to the core of the drop. d) The development of the hole in the radial direction and separation of the liquid torus. Corresponding times are, 0.411  , 0.419  and 0.431  , respectively. This picture shows the part of the droplet in which the upper rim is disconnected from the rest of the bag. By development of the punctures the bag vanishes and only the liquid bridges that link the core and the torus would remain. (The necking before the hole creation is demonstrated in a box, 0.406  ). As it can be seen in figures 3.9 and 3.11, retraction of the liquid sheet by several holes generated on the droplet interface due to the hydrodynamic instabilities can be considered as the main mechanism 144Eo  36 of the bag breakup. The basic observations and discussions on the fluid retraction were carried out by Dupré (1867) and Rayleigh (1891). Nevertheless, the most basic and foremost studies were performed by Taylor (1959) and Culick (1960) who found the constant speed of retraction by balancing the rate of change of rim momentum and the force due to surface energy for an inviscid film: (3.4) where h is the thickness of the sheet. The constant Taylor-Culick velocity, TCu , is in contrast with the accurate experimental results of Debrégeas et al. (1995 & 1998) which exhibit an exponential growth of the hole radius. The exponential growth of the hole is related to viscous dissipation during the retraction procedure. Brenner & Gueyffier (1999) classified the regimes of retraction in three main categories: low viscosity regime ( ), where the rim is followed by capillary waves, and high viscosity regime ( ), in which the rim at the edge as well as the capillary waves disappear. There is also a transient regime ( ), in which the capillary waves vanish but the rim still remains. The above taxonomy was obtained by numerical solution of one-dimensional lubrication equations: (3.5) (3.6) Where the mean curvature of the film, is defined as: 1 2 2 TC d u h          0.1fOh  10 fOh 0.1 10fOh    0t xh hu     4 d t x x x x dd u u u h u h              ,x t 37 (3.7) In equations 3.5 to 3.7, subscript x denotes the differentiation with respect to space and u is the one-dimensional velocity. Brenner & Gueyffier (1999) argued that the one-dimensional nature of equations 3.5-3.7 limit the generality of the solution. Thus, the initial exponential regime found by Debrégeas et al. (1995) cannot be extracted. However they verified that in all regimes the terminal velocities would be close to Taylor-Culick velocity, equation 3.4. Conversely, Savva and Bush (2009) analytically demonstrated that the lubrication equations are able to explain the exponential expansion of a retracting hole. Using the same equations, equations 3.5-3.7, they presented an analytical equation for the case of high , which is able to display the exponential enlargement of the hole: , (3.8) where the tip velocity ( *u ), time ( *t ) and position ( *x ) are normalized by , and , respectively. As the author has mentioned, equation 3.8 is strictly valid for the high regimes. This equation can be reduced to the dimensional form of for the tip velocity. The retraction velocity of the tip for the selected hole in figure 3.11 is extracted and compared with the analytical results of Savva and Bush (2009), equation 3.8. Results are compared in figure 3.12.        3 2 2 1 2 1 1 4 x x x h h           fOh * *2 * * * * * 1 exp 416 4 t x x u x erfc t t                 TCu 2 d h  f Oh h fOh   1 2 3 TC fu tu Oh h 38 Fig. 3.12. Retraction speed for a selected hole in simulation of . At early times after the hole creation, the retraction speeds from the simulation are comparable to values by equation 3.8. After the early stage of the retraction, also called the exponential regime, the theory underestimates speed. There are many differences between the assumptions made in the theoretical solution of Savva and Bush (2009) and the situations we have in the bag breakup. In the theory, the end of the retraction sheet is fixed in a constant height. However, this quantity is not only varying in our case, but it is non-uniform. Thus, we have to calculate a mean value for this amount. Moreover, equation 3.8 was proposed for early stages of the puncture growth as well as the high conditions. These differences may cause the observed underestimation. For the considered hole, using the average value for the thickness of the film, the Ohnesorge number was obtained to be, . According to the results of Brenner & Gueyffier (1999) and Savva and Bush (2009), for the achieved here, the capillary waves vanish and the rim appears around the tip. This is in agreement with what is seen in figure 3.11; where the rim can be seen and the capillary waves are very weak. 144Eo  fOh 0.1931fOh  0.1931fOh  39 As it was observed above, a few number of holes are generated at the interface of the bag for . The number of holes significantly increases with Eötvös number, as it was also expected from the results of Bermond and Villermaux (2005) shown in figure 3.8. Rise in the number of cracks over the bag will not affect the total mechanism of breakup. However, the number of generated ligaments will notably increase. To clarify the differences between the breakup fashions for higher Eötvös numbers and lower ones, the disintegration of the bag for will be discussed in detail. Figure 3.13 demonstrates the bag breakup for . Fig. 3.13. Bag breakup of the falling droplet . a) 0.377  , b) 0.386  , c) 0.394  ,d) 0.403  . Figure 3.13 shows the remarkable similarity between the bag breakup of a droplet with the liquid sheet breakup observed in figure 3.8. The bag is perforated by several punctures, once generated grow and form the mesh of the ligaments. Similar to , the upper rim and the flattened core of the drop will be left after the bag breakup. The main distinction seen in higher Eötvös numbers is the creation of the ligament network that crumbles quickly and forms a number of droplets. This is shown in detail in figure 3.14. 144Eo  216Eo  216Eo  216Eo  144Eo  40 Fig. 3.14. Bag fragmentation for 1) The web of attached ligaments created by several holes. 2) Growth and deformation of the ligaments. 3) Increase in the number of holes. (green and yellow arrows demonstrate the upper rim and flattened core of the drop, respectively). 4-6) Further enlargement of the punctures and detachment of ligaments. (the red box in picture 5 depicts an attached ligament into the upper rim.) 7-9) Breakup of ligaments and creation of droplets. (the red boxes in pictures 8 and 9 point up a ligament attached to the flattened core and a stable droplet, respectively). As it is seen in figure 3.14, the holes which create the attached ligaments enlarge and finally break up, forming the ligaments attached to the upper rim or the flattened core, as well as some stable or unstable fragments between them. The unstable fragments will experience another disintegration generating smaller droplets. To illuminate the deformation of the interface as well as the creation and breakup of the ligaments shown in figure 3.13, a cross-section view of the droplet after the formation of the bag is illustrated in figure 3.15. 216Eo  41 Fig. 3.15. Instability and disintegration of the bag for 1-2) Small amplitude waves at the thin sheet. (blue arrows point to some of the waves). Corresponding times are 0.380  and 0.382  , respectively. 3) A hole created on the interface. ( 0.384  ) 4) Growth of the generated hole in the previous step. ( 0.386  ) (The blue arrow shows a noteworthy deformation close to the core. The green arrow highlights the waves generated during the retraction. The yellow arrow shows another significant deformation for the part which is attached to the upper rim.) 5) Detachment of the bag from the flattened core and formation of other holes. ( 0.3882  ) (The blue arrow shoes the gap generated between the bag and the core. The yellow arrows demonstrates the cross section of a fragment created by growth of two neighboring holes). 6) Further deformation of the interface and enlargement of the gaps. ( 0.3883  ) (Blue arrow depicts the deformation appeared during the retraction) 7-8) Necking and formation of another hole which forms more fragments. Corresponding times are 0.392  and 0.395  , respectively. 9) Generated fragments which will experience more disintegration. ( 0.396  ) (Yellow and green arrows demonstrate a ligament attached to the upper torus and the drop core). As it can be observed in figure 3.15, the liquid sheet that makes up the bag case and then breaks up. The waves discussed for the bag breakup of the drop in the are observable here again. Another point is the flattening process of the core which can be seen in figure 3.15 (6-9), after the detachment of the bag. A phenomenon seen in various stages of the atomization process is the drop formation from a ligament. This will happen for some parts of the fragments formed due to the bag breakup, as well as 216Eo  144Eo  42 to the ligaments attached to the upper torus or the core. Figure 3.16 illustrates this process for a typical ligament selected from the simulation of . Fig. 3.16. Cross-section view of a ligament showing drop formation, . Domain is coloured by non-dimensional pressure, . a) The capillary waves grow causing a thicker area in the ligament. ( 0.457  ). b-d) Growth of the waves over the ligament. Corresponding times are 0.460  , 0.461  and 0.464  , respectively. e) Neck generation. ( 0.466  ). f) Pinch-off occurs and the droplet detaches from the ligament. ( 0.468  ). The mechanism observed in figure 3.16 is very close to droplet formation detected in the numerical study of Shinjo and Umemura (2010), who investigated the dynamics of ligaments and droplets in a liquid jet, using direct numerical simulations. A detailed explanation is provided in their work on wave modes, mode selection and the mechanism of the pinch-off. It was mainly shown that, in most cases, droplet formation occurs by the short-wave modes and the pressure inside the ligament increases by thinning of the ligaments at the neck. In the current study, the mechanism of drop formation from the ligament is the same. Moreover, the multi drop formation from a ligament, or drop pinch-off from an attached ligament (to upper torus or core) are seen, as well. More detailed investigations can also be found in studies of Ashgriz and Mashayek (1995), Mollot et al. (1993) and Li et al. (2008). After the breakup of the bag, three main parts remain: the upper rim, the flattened core, and the fragments created between these two. From here the disintegration of the core, deformation of the 216Eo  216Eo  * /p pD  43 upper rim, and motions of the fragments will be considered. Figure 3.17 presents a series of pictures showing the bursting of the liquid core remaining after the bag breakup. Fig. 3.17. Rupture of the flattened droplet core for . Pictures are captured from bottom view and a transparent hollow circle is added to separate the core for clarity. a) Flat core after the bag breakup. ( 0.424  ). b) Initial holes and cracks generated on the interface of the core. ( 0.441  ). c) Detachment of the ligaments in the generated web. ( 0.444  ). d-e) Growth of the holes and development of the crumbling. Corresponding times are 0.447  and 0.450  , respectively. f) Separated ligaments. ( 0.454  ). g) Capillary wave growth over the ligaments and drop formations. ( 0.469  ). h) Stable fragments after the core disintegration. ( 0.483  ). (The large volume of the fluid observable in pictures e to d is the upper rim, which can be seen from below after the breakup of the core). Subsequent to the bag breakup, a considerable volume of the original droplet would be left in the bottom, named the droplet core here. This core flattens quickly because of the hydrodynamic pressure on the bottom of the core (stagnation point). However, only the middle of the core would be squeezed, surrounded by a thick torus (figure 3.17-a). Due to instabilities discussed above for the bag breakup, the interface of the core begins to rupture (figure 3.17-b). Again a network of ligaments is 216Eo  44 formed and further disconnected rapidly (figure 3.17-c). The cracks made by the holes expand and finally a number of separated ligaments remain, as well as the thick torus (figure 3.17-f). The separated ligaments are highly unstable due to rapid growth of the capillary waves (figure 3.17-g). The pinch-off modes develop over these ligaments and droplets are formed by the same mechanism discussed in figure 3.16. The number of the droplets generated from a ligament depends on the ratio of the length of the ligament to the acting wave-length. As an example, three main droplets are formed by disintegration of the ligament magnified, in figure 3.17-g. Subsequent to the breakup of the ligaments, stable fragments remain. Once these droplets are produced, they will join the ones generated by the bag breakup. At this stage, the domain of fragmentation includes two liquid tori (upper rim and the torus made after the core bursting) as well as the fragments made by the bag and core fractures. In figure 3.18, a sequence of the deformations of the tori is illustrated. Fig. 3.18. Deformation of the upper liquid torus and the ring left after the rupture of the flattened core. Corresponding times are 0.462  , 0.520  , 0.551  , 0.582  and 0.610  , respectively. The tori will fall faster than the small fragments. This happens due to two main reasons; the first one is the low velocity of the droplets made in bag breakup because of the energy lost due to the disintegration. The second one is the higher weight of the tori. Thus, the relative positions of the fragments with respect to the positions of the upper torus change and they will occupy the upper area of the tori. During the fall, the upper and lower tori will go under significant deformation specific to the tori. This type of deformation is a combination of Rayleigh instability for a torus, Pairam and Ferna´ndez-Nieve (2009), and the influence of the gravity (can also be explained by Rayleigh-Taylor instability). Thus the waves generated over the tori will be magnified and tend to move into the direction of gravity. This fashion of deformation generates the liquid bridges, which connect the fragments. These large liquid bridges can be observed in the last steps of figure 3.18. 45 The liquid bridges rapidly disintegrate into fragments. In figure 3.19, a sequential picture of liquid bridge column breakdown is shown. Fig. 3.19. Breakup of the liquid bridge made after the deformation of the lower torus. a) the liquid column caused by the deformation of the lower torus. ( 0.666  ). b) the capillary waves grow on the liquid bridge. ( 0.692  ). c) the necks appear due to development of the waves. ( 0.721  ). d) the liquid column disintegrates into several small, stable droplets as well as a few unstable ligaments. These ligaments will be deformed into stable fragments by further breakup or retractions. ( 0.748  ). As it can be seen in figure 3.19, the liquid column becomes thinner due to stretching. Meanwhile, the capillary waves grow causing the fragmentation. All the liquid columns will break up rapidly into small droplets and ligaments due to Rayleigh instability. The ligaments will further breakup or retract due to surface tension, forming small, stable droplets. At these stages the fragmentation domain includes two main parts. The first one is the stable small fragments created during different stages of the fragmentation. The second part consists of larger fragments. These portions can be ligaments, deformed drops or small liquid bridges connecting two droplets. Figure 3.20 demonstrates the fragments made after the breakup of tori. 46 Fig. 3.20. Different styles of fragments made after the breakups. a) Stable droplets. b) A small liquid bridge connecting two droplets. c) Deformed ligaments just before the pinch-off and drop formation. d) Larger, deformed fragments. ( 0.886  ). Most of the stable droplets formed during the bag or core bursts are located on the upper part of the fragmentation domain. The larger, heavier pieces are in the lower part of this domain. Depending on the size and the velocity of these fragments, they might or might not experience further breakup, producing smaller droplets. As an example figure 3.20-b,c depict two different forms of the larger fragments, generating smaller droplets. Eventually all ligaments and deformed droplets will be in stable conditions with no further fragmentation. At these stages, the fragmentation domain consists of a group of small fragments with different sizes and characteristics. These characteristics and properties are the topic of the next section. 3.4 Number of Fragments The mechanisms of drop fragmentation were described in previous sections. As demonstrated, in each step of the atomization process droplets broke up into a number of smaller fragments with different sizes. In this section, initially we show the variation of fragmentation numbers during the atomization. Afterwards, the fragmentation size distribution will be discussed. 47 Figures 3.21 to 3.23 depict the variation of the number of fragments during the falling procedure for , , and , respectively. Fig. 3.21. Variation of number of fragments during the falling process, . Fig. 3.22. Variation of number of fragments during the falling process, . 144Eo  216Eo  288Eo  144Eo  216Eo  48 Fig. 3.23. Variation of number of fragments during the falling process, . Similar to what was described in the previous sections, the process is divided into four main parts; deformation, bag breakup, upper rim disintegration and core fragmentation. The number of droplets remains one for the deformation part. After the generation and growth of holes, breakup begins. The number of the droplets suddenly grows due to bag breakup. This rise is followed by another augmentation, due to upper torus breakup. At the end, the core disintegrates into several fragments and all the unstable fragments will further break up. This causes another growth in the diagram. The rate of the growth for the number of fragments is relatively larger for the tori breakups, in comparison with the bag breakup period. In some parts of diagrams associated to , and the number of fragments do not grow. These are the periods of the tori deformations, when no fragmentation occurs. For , this time period is very short because of very fast deformation and small breakup time scales. In general, increasing the Eötvös number, the fragmentation phenomena occurs faster. The time for the bag breakup is found to be 0.282B  , 0.0916B  and 0.0769B  for , , and , respectively. These values for upper torus and the core breakups also decrease by Eötvös number. Generally, the number of fragments at the end of the atomization process will increase for higher Eötvös number. Analysing the curves shown in figures 3.21 to 3.23, it is understood that the number of fragments generated in the bag breakup duration is notably less than those produced by the tori breakups. This 288Eo  144Eo  216Eo  288Eo  144Eo  216Eo  288Eo  49 may be related to the volume occupied by each part. In all the simulations performed here, the volume ratio of the bag, upper torus and the core are about , and , respectively. To compare the fragments size generated in each simulation, the histograms of the non-dimensional droplet diameters are shown in figure 3.24. Fig. 3.24. Droplet size distribution for a) , 32 0.00234665d  , 0.00521aved  b) , 32 0.00187363d  , 0.00444aved  and c) , 32 0.00168808d  , 0.00412aved  . The lines demonstrate the log-normal distribution. 0.25 0.32 0.43 144Eo  216Eo  288Eo  50 A log-normal distribution is fitted for the histograms. As it can be seen from the diagrams shown in figure 3.24, the majority of droplets’ diameters in all three cases is less than 0.1 of the initially spherical droplet. However, the size distribution of the fragments are not the same. The Sauter mean diameter, 32d , and the average diameter of the fragments, aved , are decreased by reducing the value of Eötvös number. As seen, the size of the fragments reduces by Eötvös number. On the other hand, the majority of the fragment diameters are in the range of . As a conclusion, the fragmentation time scale and mean fragment size decrease by Eötvös number. The statistical characteristics of the fragments during and after the breakups need more detailed and different tools of investigation which is provided in the next chapter. 0.01 0.1fD D  51 Chapter 4 Statistical Characteristics of Fragments 4.1 An Overview on Droplet Deformation and Breakup As described in figure 1.2, the unstable droplet disintegrates into several fragments to balance the forces exerting on it (surface tension, drag). Due to the asymmetry of these forces, the droplets do not necessarily fragment into spherical ones, and different fashions of topological changes can occur. Mainly, after a significant deformation, a hollow thin film forms connecting a torus (on top) and the droplet core. This thin film is usually called a bag. Depending on the falling conditions, the bag direction can be downward or upward as mentioned in Han and Tryggvason (1999). It is also found that for the high Eötvös region 70Eo  , the bag direction is always downward Han and Tryggvason (1999). Due to growth of instabilities over the bag, it disintegrates rapidly into several, very small droplets leaving the upper torus, and the core of the parent droplets. The torus will break up due to Rayleigh-Taylor instability, generating more fragments. At the end, the droplet core will flatten and deform into another torus (bottom torus) which fragments further. Subsequently, a cloud of small stable fragments with different size are remained in the domain falling down with different velocities. Figure 4.1 illustrates the breakup of a droplet and creation of the fragments for 144Eo  . Every science begins as philosophy and ends as art. Will Durant (1885 –1981) 52 Fig. 4.1. Droplet disintegration for 144Eo  . a) Initial spherical droplet. b) Fragments generated after the bag breakup and bottom torus deformation. c) Breakup of the droplet core and deformation of the remaining large fragments. d) Cloud of fragments. The knowledge of variation of fragments number and their size and velocity distribution are significant for a wide range of applications. Moreover, an accurate presentation of the fragments’ characteristics assists further numerical modeling. These issues will be discussed in this chapter. 4.2 Fragment Number Variation Naturally, the number of fragments increases by time, while the rate of variation depends on the mechanisms of breakup, as well as the properties of the system. Figure 4.2 demonstrates the increase of fragment number, fN , for different Eötvös numbers. 53 Fig. 4.2. Fragment number variations for different cases. a) 144Eo  . b) 216Eo  . c) 288Eo  . 54 The number of fragments remains one for the initial stage of the falling process while the droplet deforms and a bag is generated. The number of fragments increases considerably once the bag disintegrates. This growth is followed by another increase due to breakup of upper and lower tori. Basically, three different types of fragments can be observed at this stage; small stable droplets, ligaments and unstable fragments. Ligaments retract or disintegrate producing one or more droplets, respectively. Furthermore, the large fragments, usually generated by tori breakup, further break up producing more fragments. This cascade continues until all the generated fragments of the cloud reach a stable condition. Reaching the completely stable cloud of droplets using the numerical simulations requires an extremely large number of grids and large CPU times. Using the current adaptive grids, we are able to follow very small scales of fragments in our simulations. Simulations are stopped when no considerable changes occur for fN and the average diameter of the fragments satisfies * 10 0 1 0.045 fN i f i d d D N    . As demonstrated in figure 4.2, the time of fragmentation decreases with increasing Eo . Moreover, the value of fN at the end of the process increases with Eo . These observations can be explained using Rayleigh-Taylor Instability where increasing the acceleration (Eötvös number), a shorter time is required for dispersion of perturbation. Furthermore, smaller wave-lengths are amplified with increasing the acceleration. Thus, smaller droplets (larger fN ) are generated for larger Eo . Although the number of fragments, fN is an important factor by itself, no information can be obtained about the size, velocity and other characteristics of fragments using figure 4.2. Therefore, three different points are selected for each case (also shown in figure 4.2) to analyze the temporal variations of fragment characteristics. The first point is selected close to occurrence of bag breakup. The second point is chosen after the breakup of the droplet core. Finally, the last stages of the simulations are selected as the third points. Fragment characteristics at these points are extracted and then compared to show the development of cloud of fragments by time. 55 4.3 Fragment Size and Velocity Distribution The distribution of size and velocity of droplets are perhaps the most important factors in fragmentation studies. Using experimental techniques, usually a considerable number of droplets are missed in measurements. However, all the fragments can be followed easily in the current numerical simulations. Thus, an accurate calculation can be performed for distribution of fragment characteristics. Figure 4.3 demonstrates the histograms of non-dimensional diameter and velocity for different stages of simulations. Since the falling velocity is the characteristic velocity for the fragments, only this value is presented here. The most important variation to be considered is the peak value in each step. In the case of size distribution, most of the fragments are in the range of 00.01 0.1fD D  and the peak of the curves increases by time, demonstrating the increase in the number of droplets. Comparing the size distribution of first and second points, the number of fragments after the tori breakup is approximately 7 8 times greater than its value after the bag breakup. The peak of probability after the bag breakup for the all three cases is 0/ 0.03fD D  . This value is reasonably close to the maximum stable size of the droplets for the present two-phase system. This peak separates into two peaks after the tori breakups; one close to the previous diameter with considerably more fragments and another one in the range of 00.04 / 0.06fD D  . The lateral range demonstrates the larger fragments generates in the process of upper tori and droplet core breakups. These relatively large fragments (if unstable) will further disintegrate, producing more droplets in different range of diameters. Subsequently, two peaks are remained in size distribution histograms. The significant one corresponds to the stable small droplets close to the maximum stable droplet size. The second one corresponds to the relatively large fragments remained from the tori breakups. In the case of velocities, the behavior of range of variations for each simulation and stage is more complicated. This is due to chaotic dynamics of droplet disintegration caused by significant morphological changes, different fashions of breakup, and carrier phase flow structures. Basically the range of velocity variations expands with Eötvös number due to higher magnitude of the gravitational force. For 144Eo  , 216Eo  and 288Eo  these ranges are approximately 00 30fv v  , 00 33fv v  , and 00 41fv v  , respectively. The peaks of the curves again increase by time demonstrating the cascade of breakup and generation of more fragments in the mentioned ranges. 56 Fig. 4.3. Size and velocity distribution of fragments for different stages and Eötvös numbers. a) Diameter distribution for 144Eo  . b) Velocity distribution for 144Eo  . c) Diameter distribution for 216Eo  . d) Velocity distribution for 216Eo  . e) Diameter distribution for 288Eo  . f) Velocity distribution for 288Eo  . 57 Since non-dimensional diameter of droplets does not explain the variations of surface areas or volume of the fragments, the other types of characteristic diameter representations are obtained using the characteristics diameter formula, Mugele and Evans (1951):       max max min min 1 * 1 0 0 0 f f f f f f p q D D p q pq f f f f D D d D f D dD D f D dD D                              (4.1) Where p and q are positive integers and  0 ff D is the number of fragments probability density function and the characteristic diameter is normalized employing the initial diameter of the droplet. Normalized Sauter mean diameter (SMD), * 3 232 0 1 1 f fN N i i i i d d D d             , and normalized de Brouckere mean diameter (also called mass mean diameter, MMD), * 4 343 0 1 1 f fN N i i i i d d D d             are the common characteristics diameter in atomization investigations. SMD is the diameter of droplet that corresponds to the volume/surface ratio for the entire liquid. MMD is also the average diameter based on the unit volume (mass) of a particle. Experimental investigations of Simmons (1997-a and –b) and Hsiang and Faeth (1992 and 1993) on secondary atomization of industrial sprays and single droplets show a constant rational correlation between SMD and MMD of fragments. (e.g. * *43 32/ 1.2d d  ) This correlation was found to be applicable for different modes of breakups including bag, multimode, and sheet thinning. The values of SMD and the ratio of 4332 /d SMD MMD are illustrated in figure 4.4 comparing with the mentioned empirical constant relationship. 58 Fig. 4.4. Variations of Sauter mean diameter, *32d and 43 32d (Square and circle markers, respectively) for different values of Eo . Solid (red) line depicts the empirical relationship. a) 144Eo  b) 216Eo  c) 288Eo  . 59 The value of *32d as well as * 43d reduce by time demonstrating the disintegration of larger fragments into smaller ones. Moreover, the calculated values of 4332d obtained by simulations are in agreement with the constant relationship achieved in experiments. Since the size and velocity distributions at the end of the breakup process, when an stable cloud of particles are achieved, are more important, the velocities and diameters of fragments at the third stage of all simulations are reproduced in figure 4.5. Fig. 4.5. Histograms of a) size distributions and b) velocity distributions for last step of each simulations. In the case of size distribution, the peaks are close to 0/ 0.03fD D  . These peaks are followed by other (less pronounced) peaks for the range of 0/ 0.05fD D  . These larger fragments are usually generated in the core and tori breakup stage, and are stable enough to resist further fragmentation. In the case of velocity distributions, the peaks are in the range of 02.5 / 7.5fv v  while the tail of the distribution extends by increasing Eötvös number. This means higher falling velocities are obtained in higher Eo due to lager body force exerted on the fragments. For more details, data obtained at the last step of the simulations are modeled using log-normal probability density functions shown in figure 4.5. The validity of log-normal PDF for droplet size 60 distribution is verified and it is usually used in droplet size modeling, Babinsky and Sojka (2002). The log-normal formulation used here is as follows,       2 0 ln1 1 exp 2 lnln 2 f f DD LNf LN D D f D D                (4.2) where D and 0DLN  represent the logarithmic mean size and width of the distribution, respectively. In the case of velocity distribution, D and DLN are replaced by V and v LN , respectively. (Some negligible negative values in falling velocity distributions are obtained due to complex fashion of breakups. These values are ignored in calculation of log-normal distributions for velocities.) Table 4.1 shows the list of parameters in log-normal PDF for different Eo . Table 4.1. Parameter in log-normal PDF modelling Parameter 144Eo  144Eo  144Eo  D 0.0428 0.0358 0.0348 D LN 2.0279 1.8678 1.8956 V 6.6224 7.2748 8.3221 v LN 1.8634 1.8591 2.1442 Tables 4.1 represent the logarithmic mean size and width of the log-normal distributions. As it is clear, the mean size, D , is reduced by increasing Eo . Moreover, the width of the size distributions, D LN , for 144Eo  , is larger than two other simulations denoting the larger fragments in simulation of 144Eo  . These fragments remained in the cloud of droplets the surface tension resistance in this case is large enough to avoid breakup. No sensible difference was observed between the widths of 216Eo  and 288Eo  . For velocity distributions, the logarithmic mean V is increased by increase of 61 Eo . This shows the growth of fragments’ velocities by adding to Eötvös numbers due to larger amount of gravitational force. Moreover the width of log-normal distributions, vLN , is larger for 288Eo  which demonstrates the larger range of velocities obtained for this simulation. No sensible difference is observed for 144Eo  and 216Eo  , in this case. The data analysis we had so far, included the variation of size and velocity of fragments. However, for a complete understanding of statistical characteristics of fragments we need the contributions of the fragments in volume as well as the correlation of size and velocities. Figure 4.6 demonstrates the volumetric fragment Size distributions, finding the value of   3 0fD DV for quantiles of the normalized diameters. Fig. 4.6. Volumetric fragment size distributions for the simulations. Left axis: V for 144Eo  . Right axis: V for 216Eo  and 288Eo  . The curve in figure 4.6, in fact demonstrates that how much of the volume of the initial drop is converted to what range of the fragment size. The first peak in curves shown in figure 4.6 62 corresponds to the large number of medium size fragments ( 00.06 0.08fD D  ) in the domain. This peak is approximately 0 0.07fD D  for the both simulations. However, the droplet size distribution curves show 0 0.03fD D  as the peak. This is due to particularly small size of the fragments of the range 00 0.04fD D  that do not contribute to the volume. The second peak in figure 4.6 (about 0 0.145fD D  ) depicts the considerable role of larger fragments despite their small probability as shown in the FSD curve. This shows a considerable volume of the initial drop is converted to few large fragments. The important point in figure 4.6 is the independency of peaks on the Eötvös numbers. 63 Fig. 4.7. Scatter plots of velocity-diameter for different simulations and stages. Lines demonstrate the conditional probability for each stage of simulation. Top: 144Eo  , Middle: 216Eo  , Bottom: 288Eo  . 64 To show the correlation of size and velocity, the scatter plots of these parameters as well as the conditional distribution of velocity with respect to diameter:       0 | | ji f f ji ji f f f f D D f D D D D f v V D V D N v       (4.3) are illustrated in figure 4.7. In equation 4.3,  jif fD D denotes the range of diameters subsets used to calculate |V D . In the present study, this range is fixed at   0.05jif fD D  and   0.1jif fD D  for 00 / 0.15fD D  , and 00.15 / 1fD D  , respectively. The highest probability area for size-velocity distribution in all simulations approximately appears in the range of 0 0 0.02 / 0.1 0 / 15 f f D D v v          . However, the higher range of velocity 010 / 30fv v  in the same diameter area has a considerable probability as well. For a given simulation (Eötvös number), the manner of changes of |V D can also describe the disintegration of the droplet. At the first stage, (e.g. A1, B1, or C1), a number of small fragments plus a single large fragment 0/ 0.1fD D are observed. The large fragment is actually the parent droplet remaining after the bag breakup, and the small droplets are those generated due to bag breakup. This fashion of distribution makes a three-part curve of |V D including a plateau for small fragments, a peak for the large fragment and a zero value between these two parts which demonstrates the considerable difference between the size of the fragments generated by bag disintegration and the parent droplet. After the breakup of droplet core, a number of large fragments remaine in the domain, expecting to further breakup. These fragments relatively have larger falling velocities due to their larger weights. The peaks observed for |V D curve at this stage, (beside the plateau of smaller fragments) are the consequence of these fragments. By further disintegration of unstable fragments, the number of small droplets increases and large fragments disintegrate into several smaller ones. Due to these supplementary breakups, the |V D curve shows a smooth manner increasing by fragment size and a peak occurs at the end for the largest fragments. Additional disintegration (is not presented here) of fragments will not change the pattern of the curve but slightly the values of the peak (both diameter and velocity). 65 Distributions of size and velocity, individually, do not show the stability of the fragments. This led us to calculate the non-dimensional numbers for each fragment to analyze the stability of droplets at each stage. 4.4 Variation of Non-Dimensional Groups As it was mentioned in section 2.4, the non-dimensional group governing the dynamics of falling droplets consists of Eo , Oh , We , and Re . The values of these groups are calculated for each fragment at each stage and presented in this section. As described before, Eo and Oh represent the charachteristics of a droplet ar rest. For a moving droplet, these values have the same meaning, ignoring the influences of surrounding media. (e.g. characteristics of the equal-volume spherical droplet of each fragment if they begin to fall from the rest). Figure 4.8 illustrates the variations of Eötvös versus Ohnesorge number for all stages. Fig. 4.8. Eo versus Oh for the last step in different simulations. Han and Tryggvason (1999) presented a series of numerical simulations for falling droplets with different parameters. For the case we are investigating ( 0.05Oh  and 10d c   ) they predicted that breakup will start at 50 60Eo  . However, for larger amount of Oh , we need a higher value of Eo for breakup. Thus, ignoring the influences of continuous medium and assuming that all 66 fragments are spherical and just begin to fall, figure 4.8 can be used to analyze the stability of the fragments. Results show that, based on this criterion, all fragments are almost stable at the last step of the simulations. However the simulation of 288Eo  may still contain a few number of fragments that can be disintegrated to smaller ones. It should be noted that the ratio of Eo and Oh only depends on the size of the droplets (e.g.   5d d fEo Oh g D     ) producing the charactrisitics lines for simulations. These lines are illustrated in a log-log diagram in figure 4.8. The analysis above ( Eo vs. Oh ) does not include the influences of velocities. Thus, to take account of the velocity of the fragments, we also obtained the Re versus We maps for different stages in the simulations. The ratio of Re cWe v  only depends on velocity of the fragments. However, analyzing the scatter plot of (We vs. Re ) helps to demonstrate the range of variations for these non- dimensional parameters as well as their correlation (if any). Figure 4.9 depicts the diagrams of ReWe for different simulations. 67 Fig. 4.9. We versus Re for different simulations. Diameters of points are proportional to diameter of fragments. Top: 144Eo  , Middle: 216Eo  , Bottom: 288Eo  . 68 According to figure 4.9, most of the fragments at the end are in the range of 0 Re 20 0 7We         . This domain is slightly larger for 288Eo  due to stronger body force that generates a lager velocity. Due to the same reason, the maximum values of We and Re at the last steps (A3, B3 and C3) also increase with the initial Eötvös number. Up to now, none of the analysis performed in this section provide a good estimation about the stability of fragments. To include all the influencing factors to predict the stability of moving droplets, we can show the classical diagram of (We vs. Oh ) for the simulations. This analysis considers all the dispersed-phase-continuous-phase combination properties, but the morphological changes of the fragments at that stage. So, assuming all the fragments being spherical at each moment, using equivalent diameter of the same volume sphere, Weber and Ohnesorge numbers are calculated and demonstrated in figure 4.10. 69 Fig. 4.10. We versus Oh for different simulations. Diameters of points are proportional to diameter of fragments. Breakup modes are shown based on the experimental studies by Hsiang and Faeth (1992 and 1993). Top: 144Eo  , Middle: 216Eo  , Bottom: 288Eo  . 70 Diagrams shown in figure 4.10 appropriately explain the fragmentation of droplets in time. For a given Eötvös number, a large fragment (parent droplet) with 0.05 37 Oh We       and a number of small fragments with 0.1 12 Oh We       remain at the first stage (A1, B1 or C1). Based on map of Hsiang and Faeth (1992 and 1993), small fragments will not further disintegrate. However, the large fragment is close to bag and multimode regimes. In reality, the remaining large fragment is considerably deformed and punctured severely. Nevertheless the prediction of the curve on further breakup is correct but its fashion is totally different. (e.g. the multimode or bag can happen for a spherical droplet). At the second stage (A2, B2 or C2), the number of unstable fragments, close to the domain of bag and multimode, as well as the small stable droplets increase. This is due to breakup of the upper torus and the droplet core, which remains some large fragment in the domain. The values of Weber and Ohnesorge numbers for the large fragments in this stage are approximately 0.05 0.1 12 35 Oh We         . Based on the map, these fragments (the number of them is about 4 to 5) will break up again generating more stable droplets. This fact happens and at the third stage (A3, B3 or C3) a cloud of stable droplets are remained. Analyzing the data at the last stage shows that about 99.33%, 99.57% and 99.16% of fragments for 144Eo  , 216Eo  , and 288Eo  are in oscillation or deformation regimes, respectively. This shows that the cloud of droplets created at the end of simulations are quite stable. It should be noted that the lateral analysis can underestimate the values of Weber number, neglecting the deformation of the fragments. However, based on our comparison the prediction is in good agreement with the observations. 71 Chapter 5 Deformation of Droplet in a Stream 5.1 An Overview on Droplet Dynamics in a Stream However, several experimental investigations have tackled the problem using different methodologies including shock tubes (Hsiang & Faeth 1993, 1995), continuous jets (Lee & Reitz 2000, 2001, Zhao et al. 2010), and drop tower (Villermaux & Bossa 2009). In most of these experiments, breakup of the droplets is studied analyzing the pictures from shadowgraphs or high- speed cameras. Empirical results are mostly classified in geometrical categories where the Weber ( 2 0 0gWe U D  ) and Ohnesorge ( 0d dOh D   ) numbers distinguish the breakup regimes. Vibrational, Bag, Multimode, Sheet-thinning, and Catastrophic breakups are the modes observed in the experiments with increasing Weber number. However, as it is also mentioned in the review paper of Guildenbecher et al. (2009), mode transition is essentially a continuous procedure and different studies have reported different geometries as well as values of transition Weber numbers. Nevertheless, a review of the empirical results for breakup of droplets reveals that the details of fragmentation process still remain elusive. This is due to the high-speed three-dimensional occurrences during the deformation and breakup of droplets. Thus, it is extremely problematic to capture all the characteristics of droplet fragmentation in a stream using the current experimental techniques. The alternative approach is numerical simulations where the bottle-neck is the cost. In other words, an excessively large number of finite discretized volumes and extensive computational times are required to solve these kinds of complex problems. Such simulations need solutions to two-phase, unsteady, 3D flow with enough resolution over a wide range of length and time scales. Thus, the majority of numerical simulations that represent a time-dependent deformable interface of a droplet The laws of nature are but the mathematical thoughts of God. Euclid (300 BC) 72 are limited to investigations (topological changes) comparable to the drop size (Quan & Schmidt 2006; Helenbrook & Edwards 2002; Wadhwa et al. 2007; Feng 2010). This corresponds to the low Weber and Reynolds numbers regime. However, a large number of instabilities with different length scales are expected to occur during the droplet fragmentation in high Weber and Reynolds numbers, representing the majority of natural and industrial applications. Another point that has to be considered in numerical simulations is the three-dimensionality of the deformations (as it will be discussed in this article). However, a few studies can be found in the literature that have considered three-dimensional geometry to study the deformation or breakup (Quan & Schmidt 2006; Khosla et al. 2006). Due to limitations of experimental and numerical studies described above, the mechanisms of liquid droplets deformation and breakup in small scales have not been clarified yet and needs more developed techniques for further investigations. Meanwhile, reducing the computational costs, by using special techniques, Direct Numerical Simulation (DNS) of the Navier-Stokes equations combined with interface tracking is the best known technique for the study of the process. In the current study, direct simulations are performed in two and three dimensions based on appropriate numerical algorithms in the context of the volume of fluid (VOF) method combined with high performance adaptive octree meshing methods. High Reynolds number and a range of relatively high Weber numbers are chosen addressing the shear breakup of droplets in a stream. The focus of this work is to identify the instabilities generating and growing over the droplets that lead to breakup, rather than the fragmentation to small droplets. The study of initial deformation, instability generation and growth are enormously important because they are the first stage of all aerodynamically-induced- breakups, considerably affecting the characteristics of the atomization. Therefore, a thorough understanding of these processes is crucial to formulate accurate atomization models. The results are categorized in two main parts. In the first part (sections 5.2), the results of 2D simulations are presented for three Weber numbers. In the next part (section 5.3), a 3D simulation is provided to study the instabilities in radial direction. Results in both categories are compared with theories. 73 5.2 Two-Dimensional Results 5.2.1 Grid Tests First of all, a grid test is provided for three different tests; 512 512 , 1024 1024 , and adaptive grids up to the level 12 (see figure 2.3). The quality of the interface deformations and vorticity domains are depicted in figure 5.1. Fig. 5.1. Deformation of the droplet in stream for different tested cartesian grids ( 100We  ) at 0.97  . Top: Uniform grids of 512 512 (Level 9). Middle: Uniform grids of 1024 1024 (Level 10). Bottom: Adaptive grid based on quad-tree meshing. The maximum level of the adaption is 12. Vorticity field is normalized using * 0 0D U  . Where rot  u u . 74 Three criteria are employed for the adaption: velocity gradient, vorticity magnitude, and interface curvature. Adaption based on these parameters assist to perfectly resolve the boundary layer created over the droplet as well as the deformation of the interface due to the gas stream. Moreover, the influence of the secondary flows (recirculation regions) on the interface and changes in the carrier phase are captured in a better quality, using the adaptive grids. The number of the cells used for these simulations, for a specific time, is given in table 5.1. Table 5.1. Grid properties at 0.87  Test number Type Total cells Leaf cells 1 Uniform 349525 262144 (512 x 512) 2 Uniform 1398101 1048576 (1024 x 1024) 3 Adaptive 487379 365458 Level 9 = 258694 Level 10 = 4195 Level 11= 14627 Level 12 = 88004 As it is clear, the number of the cells is significantly decreased, using the adaptive grids. Nevertheless, more accurate interface tracking is provided due to high resolution of the grids close to the interface and circulations. The grid resolution affects the wave-lengths of the Kelvin-Helmholtz Instability (KHI), ligament formation, and furthermore, it influences the stretching of the ligaments produced by development of the KHI. As shown in figure 5.1, ligaments can be tracked for a longer period of time before the disintegration with increasing the resolution of the meshes. The influence of grid size is not limited to interface deformation. It also plays an important role in boundary layer resolving over the droplet. The velocity profiles on the poles of the droplets and the wave-lengths of shear instabilities around the poles are shown in figures 5.2 and 5.3, respectively. 75 Fig. 5.2. Velocity profiles on the poles for 100We  . h is the height from the interface at pole. Fig. 5.3. Wave-lengths close to the pole for different tested grids. 76 As it is clear in figure 5.2, the boundary layer thickness is significantly reduced, using the adaptive grids. We have 10 finite volume cells in the boundary layer zone for the test with grid adaption, allowing for a reasonable resolution. On the other hand, the wave-length of the shear instability decreases, increasing the mesh resolution, as depicted in figure 5.3. This considerably assists to capture the primary deformation and growth of the instabilities in an acceptable quality. Subsequently, the adaptive grid (up to level 12) is used for the 2D and 3D simulations. 5.2.2 Deformation Mechanism The dimensionless groups for the 2D simulations are listed in table 2.2 Three Weber numbers, (50, 100, 200) are chosen to find the influence of the surface tension on deformation and instability growth. Although changing the surface tension coefficient changes the Ohnesorge number, we always meet the condition of 0.001Oh  , that demonstrates the negligible influence of the viscosity. Figure 5.4 demonstrates the initiation and development of the instabilities for 100We  . 77 Fig. 5.4. Numerical results for a 2D simulation ( 100We  ). Arrows demonstrates the observed scenarios for the ligaments formed due to growth of shear instablities. a) Initiation of the instabilities around the poles due to higher magnitude of the relative velocity. b-f) Amplification of the instabilities in the stream generating ligaments. g-i) Further development (stretching) of the instabilities (ligaments) under the influence of the secondary flows and main stream. j-l) Fragmentation of the stretched ligaments and droplet (fragment) formation. m) Magnified view of d. n) Magnified view of l. Recirculation zones are highlighted by black circles. Kelvin-Helmholtz waves are generated over the poles due to higher relative velocity in that region. The waves grow forming ligaments. Recirculation zones are produced here due to existence of the waves (see figure 5.4-m). The secondary flows at this stage affect the growth of the waves. Basically two different scenarios are seen in the simulations: 1) Development, ligament formation, further stretching, disintegration and fragments formation (arrows 1 and 3 correspond to this possibility). 2) 78 Deformation and merging of the waves. This might reproduce a bigger wave that grows further and makes a ligament that follows the first scenario (arrow 2 shows an example of this scenario). The disintegrated ligaments stretch again repeating the same mechanism. Black arrows, in figure 5.4-l, show two of these ligaments. During the formation and growth of the waves, the droplet flattens and a wake is generated on the rear of the droplet, due to separation of the main stream. This wake affects the rear part of the droplet making more spikes. Circles in figure 5.4-i depict two spikes generated due to wake. The topology of the wake after the disintegration of the ligaments will also be influenced significantly. Different geometries of the waves as well as the different sizes of the deformable moving fragments in that region produce more complex patterns as highlighted in figure 5.4-n with the white circle. The mechanism of the deformation described above is the same for the simulations with different Weber numbers. However, surface tension influences the Kelvin-Helmholtz instability (KHI) and the further disintegration. The influence of the Weber number is discussed in the next subsection. 5.2.3 Effect of Weber Number The influence of Weber number, changing the magnitude of the surface tension, is studied for three cases: 50We  , 100We  , and 200We  . To see the effect of Weber number on geometry and position of the parent droplet, the width (see figure 5.4-g) and position of the center-of-mass are depicted in figures 5.5 and 5.6, respectively. 79 Fig. 5.5. Variation of width of the parent droplets versus time for differentWe . Fig 5.6. Motion of the droplets in stream-wise direction for differentWe . 80 From figures 5.5 and 5.6, no notable change can be described as the effect of change of Weber number. Although, the widths of the droplets are reduced to about half of the initial diameter, no specific changes are observed in the position of the droplets. This is due to significant difference between droplet relaxation time and instability growth time scale. i.e. the required time to change the velocity of the droplet from zero to terminal velocity is much larger than the time needed for initiation and growth of instabilities. However, taking a detailed look on the interfaces, noteworthy changes can be seen due to changes of the Weber number. The deformation of the interface for a given time ( 0.27  ) is shown in figure 5.7. Fig. 5.7. Interface deformation and growth of the shear instabilities for a) 50We  . a) 100We  . and a) 200We  . Snapshots are all in 0.27  and are just shown for the upper half-circles. As it is observed in figure 5.7, the less the Weber number, the less the growth of the shear instabilities. i.e. surface tension damps KHI. The surface tension has a slight effect on the length of the shear instabilities over the droplet. Figure 5.8 demonstrates the values of KHI wave-lengths (  ) for different Weber numbers. 81 Fig. 5.8. Wave-lengths close to the poles for different We . Wave-lengths produced due to shear instabilities are decreased with increasing Weber numbers. This is also seen in Rayleigh limit of shear instabilities (compared with experimental observation of co- axial jets in Marmottant & Villermaux (2004)) for high Weber numbers where the most amplified wave-length is not dependent on the surface tension. As described in section above, growth of the ligaments and their fragmentation make complex patterns on the wake of the droplet. The deformation of the interface plus the vorticity field are illustrated in figure 5.9 to compare the simulations in different Weber numbers. Fig. 5.9. Deformation of the droplets at 1.64  for Left: 50We  , Middle: 100We  , and Right: 200We  . 82 More complex patterns are observed with increasing the Weber number, since smaller fragments and irregular ligaments are produced. However, more study is needed for this part to measure the complexity generated on droplet wakes. Basically, based on the two-dimensional simulations, it is found that the effective shear waves (smaller wave-lengths) are produced about the poles due to higher relative velocity. Moreover, reducing surface tension (increasing Weber number), larger wave-numbers are generated for KHI. In the next subsection, the current results are compared with theories. 5.2.4 Comparison with Theories (Shear Instabilities) 5.2.4.1 Combination of Potential Flow and Very Thin Vorticity Layer Consider a smooth flow past a spherical droplet with the uniform upstream velocity of 0U (see figure 5.10-a). Neglecting the deformation of the droplet at the initial stages as well as the rotation of the carrier phase, the governing equations for the incompressible gas flow reduce to Laplace equation for potential flows: 2 0,   (5.1) where u .The general axisymmetric solution of equation 5.1 obtained by separation of variables in the spherical polar coordinates, considering    ,R r    is as follows,    1 0 cos , nn n n n n A r B r P             (5.2) where nA and nB are arbitrary constants and nP is the Legendre polynomial of degree n . For the current problem, considering a uniform flow at  , we have 0 0lim cos r U r     . Moreover, assuming no normal flow across the droplet interface we have 0r  on 0r r . Assuming the solution in the form of     21 1 cosA r B r    , the boundary conditions described find the constants as 1 0A U and 3 1 0 0 2B r U . Subsequently, the solution for the velocity field can be found as (Summer & Fredsøe 2006), 83   3 3 0 0 0 03 3 1 , , 0 , , 0 cos 1 , sin 1 , 0 2 r r r r u u U U r r r                                   u (5.3) Considering the shear instabilities, the second element of the velocity (tangent to the interface) is the important one for our further analysis. Obviously the maximum amount of the tangent velocity, u , is 01.5U on the poles. Since the droplet velocity is negligible in comparison with the gas velocity, d gU U , the maximum relative velocity is max 01.5g dU U U U    . This value for the 2D domain (circular cylinder) changes to 02U since 2 0 0 2 sin 1 r u U r             . Fig. 5.10. a) A schematic view of a droplet with diameter of 0 02D r in a stream. To convert the cartesian framework to spherical polar coordinate, let 2 2 2 2r x y z   , so that cosz r  . The area in the dotted box is enlarged in picture b. b) Gas flow around the droplet at the beginning of the streaming before the deformation. c) Possible scenario for the primary breakup of the droplet where the relative velocity, g dU U U   is large enough. Now consider a horizontal relative velocity over the surface of the droplet (see figure 5.10-b). It is well-known that two fluids with different velocities are essentially unstable, producing Kelvin- Helmholtz instability (hereafter KHI) (see figure 5.10-c). For small vorticity thickness, negligible gravity, and time-independent-velocities, the dispersion relation of the waves in x-direction  exp ikx i t   is (Chandrasekhar 1961), 84    0.5 2 . d d g g d d d d d d d d u u k k i U k                     (5.4) Where surface tension stabilizes the shear instabilities for: 2 .d d d d U k            (5.5) Subsequently, for negligible droplet velocity d gU U , and large density ratios g d  , the most amplified mode and group velocity are:   22 0 sin2 , 3 d g d m U U k        0 sin . g g g d g g d d d U U U U k                (5.6) The constants , are  3 2,3 2 and  8 3,2 for 3D (spherical coordinate) and 2D (circular cylinder), respectively. Variations of the wave-numbers around the droplet are shown in figure 5.11, changing the surface tension (Weber number) and upstream velocity (Reynolds number). Fig. 5.11. Wave-numbers over the droplet from equation 5.6. Left: Influence of Weber number. Right: Effect of Reynolds number (curves are for 3D cases). 85 Figure 5.11 demonstrates the nonlinear growth of the wave-numbers from zero to maximum at 2  (pole). Moreover, increasing Weber number the peak of the wave-number curve increases. i.e. shorter wave-lengths of KHI are generated. On the other hand, increasing the upstream velocity (Reynolds number), the peak of the wave-number curve is increased due to higher relative velocity over the droplets. These variations are similar to our observations in numerical simulations. The wave-numbers obtained in the simulations (close to the pole) are compared with the prediction made by equation 5.6 as listed in table 5.2. Table 5.2. Comparison of the present direct numerical simulations and equation 5.6. Weber Number 0 mD k (DNS) 0 mD k (Equation 5.6) Ratio 50 90 1335 ~14.9 100 110 2671 ~24 200 146 5343 ~36.59 Equation 5.6, obtained by combining the potential flow and zero-vorticity layer, shows a good qualitative agreement with the numerical simulations. However, assuming zero thickness for the vorticity layer is quite unrealistic (only valid for small Weber numbers) and causes a notable difference between the results found in the simulations and the theory. This difference is more pronounced for larger Weber numbers. Thus, taking the effect of the boundary layer, we compared our results with the theories considering linear velocity profiles. 5.2.4.2 Non-Zero Vorticity Layer Consider linear velocity profiles for the phases as shown in figure 5.12. Based on the instability analysis of Villermaux (1998), followed by Kim et al. (2006), the most amplified mode is a function of densities as well as gas boundary layer thickness: 15 1 , 4 d m g g k F              (5.7) where function F is expressed as, 86     0.5 2 35 13 37 275 1 , 6 6 6 2 d c F                (5.8) Fig. 5.12. Linear velocity profile used in Kelvin-Helmholtz instability analysis. Moreover, for the thick vorticity layer, (Villermaux 1998; Marmottant & Villermaux 2004), expanding the Rayleigh approach, the dispersion relation was found as follows:           1 0.5 1 2 exp 2 1 2 , 1 0.5 1 2 d g d g k k k k                  (5.9) where       2g d g d g dU U k U U U U      . The temporal analysis of equation 5.9 results in a bell-shaped dispersion curve. The dominant mode for g dU U , can be found as, 0.5 1 1.5 g m d g k           (5.10) The group velocity at the wave number represented by equation 5.10, is presented by (Marmottant & Villermaux 2004; Dimotakis 1986): 87 . d d g g g g d g d g U U U             (5.11) A similar instability analysis is provided for co-axial jets, resulting in equations 5.7, 5.10, and 5.11. However, the experiments of Marmottant & Villermaux (2004) show an approximately 70% difference between the results of the stability analysis and the empirical dominant wave-lengths,  . This is due to the linear velocity profile approximation. Thus the wave-number value in the stability analysis is larger than (about three times) the once obtained in the experiment. However, in the DNS of Kim et al. (2006), the value obtained for the dominant KHI wave length is much closer to the stability analysis due to thicker boundary layer ( g is 4 times larger than the experiments). For the current study, we have tried to have a sufficient resolution, using adaptive grids, to capture the gas boundary layer (see figure 5.2). The wave-numbers obtained by the simulation and those achieved using equations 5.7 and 5.10 are listed in table 5.3 for comparison. Table 5.3. Comparison of the present direct numerical simulations and equations 5.7 and 5.10. 0 mD k (DNS) 0 mD k (Equation 5.7) – (Ratio) 0 mD k (Equation 5.10) – (Ratio) ~(90-146) 314 – (2.15-3.48) 315.61 – (2.16-3.5) As it can be obtained from table 5.3, in comparison with zero-vorticity theory, the wave-numbers predicted by equation 5.7 and 5.10 are much closer to those achieved in the simulations. However, the difference between experimental data and linear stability analysis are observed again. The ratio of stability analysis to experimental wave-numbers is reported to be ~3-4 by Marmottant and Villermaux (2004), for co-axial jets. This is very close to what we found in our direct simulations. In addition to the 2D simulations, a quite expensive three-dimension (3D) simulation is also provided to investigate the instabilities behavior having the third dimension (span-wise). The following sections include the results obtained for the 3D direct numerical simulation. 5.3 Three-Dimensional Results The 3D simulation is performed with the same initial and boundary conditions of the 2D tests. The geometry is like figure 2.2 and the spatial discretization (algorithm and criteria) is similar to 2D 88 simulations. An example of the octree meshing for 3D simulation in addition to the vorticity field is shown in figure 5.13. Fig. 5.13. Droplet deformation for 100We  and time 1.9  for the 3D simulation. a) Droplet interface and the generated grids. b) Adaptive meshed and the vorticity field. c) Magnified picture of the dotted box in b. A resolution similar to the 2D simulations is seen for the boundary layer. For the 3D simulation, instabilities in the radial direction are also expected, in addition to KHI. The deformation of the droplet and growth of instabilities are discussed in the following sections. 5.3.1 Deformation Mechanism The generation of instabilities due to shear is similar to the observations in 2D simulations, where smaller wave-lengths with faster growth rates are generated close to the pole, due to higher magnitude of the relative velocity. However, here in 3D, generated waves form circular geometries in the transverse plane, causing the transverse azimuthal modulation due to the acceleration. This can also be called the Rayleigh-Taylor Instability (RTI), since a fluid with higher density is accelerated in a lighter one. Figure 5.14 demonstrates the deformation of the droplet, 89 Fig. 5.14. Generation and development of instabilities over the droplets in the 3D simulation. Time is from 0  to 1.92  a) Initial spherical droplet. b) Onset of shear instabilities close to the pole. c) growth of shear instabilities in the stream-wise direction. d) Azimuthal modulation in the plane normal to the flow direction. e-f) Growth of shear and azimuthal waves. g) Formation of ligaments and further development of instabilities. h-i) Further stretching of ligaments. The number of spikes generated due to RTI depends on many factors such as the angle  and the amplitude of the wave initially produced from shear instability. Moreover, the interactions of the waves including the recirculation zones are more complex than the 2D simulations. To show the variation of the number of the spikes, transverse cross-sections of the droplet at different locations are illustrated in figure 5.15. 90 Fig. 5.15. a) Radial (transverse) growth of instabilities for planes perpendicular to the external stream. Time is 0.58  corresponding to figure 5.14-e. As seen in figure 5.15, the number of spikes increases approaching the polar cross-section. This shows the smaller wave-length (larger wave-numbers) of the RTI for the region with higher relative velocities. The higher shear causes a faster growth of the waves that translates into stronger acceleration, causing amplification of shorter wave-lengths. By growth of the stream, and span –wise instabilities over the droplet, wave-numbers change due to variation of the amplitude of the waves. i.e. smaller wave-lengths are observed with increasing the amplitude of the shear waves. These short wave-lengths instabilities generate small ligaments and subsequently tiny droplets. An example of these waves is shown in figure 5.16. 91 Fig. 5.16. Short-wave-length RTI generating small fragments. The box in the corner demonstrates the magnified view of the same picture. Similar to 2D simulations, deformation of the droplet, generation of the irregular geometries and formation of the ligaments make chaotic patterns in wake. Figure 5.17 depicts two streamlines starting from two close points on upstream. (The initial positions 0 0 0, ,x y z for yellow and red particles are  0 02.5 , 0, 0.05D D and  0 02.5 , 0, 0.05D D  , respectively.) 92 Fig. 5.17. 3D patterns on wake of the droplet. Lines show the streamlines at 1.92  . The patterns formed in the wake might generate more deformation and instabilities over the droplet. However, like the 2D test, this topic needs more detailed investigations. Essentially, based on the 3D simulation, it is found that the development of instabilities in span-wise (radial) direction is dependent on the shear magnitude (i.e. the angle between the flow direction and the point over the droplet) and the amplitude of the waves in stream-wise. In the next subsection we try to employ the available theories to find a mathematical description for these changes. 5.3.2 Comparison with Theories (Transverse Azimuthal Modulation) After the production and growth of the KHI, generating axisymmetric waves, the transverse azimuthal modulations are expected to develop in the radial direction due to Rayleigh-Taylor Instability (Marmottant & Villermaux 2004; Kim et al. 2006, Villermaux & Clanet 2002) as illustrated in figure 5.18. 93 Fig. 5.18. a) The growth of KHI over the droplet. Blue dotted lines show the arbitrary cutting planes. b-d) possible deformation of the interface in the radial direction due to Rayleigh-Taylor Instability. A uniform spatial perturbation growth with 4, 6, and 9 spikes are schematically shown here. The Rayleigh-Taylor Instability (RTI) is the instability of an interface between two fluids of different densities, occurring when the heavier fluid is accelerated in the lighter fluid. (Marmottant & Villermaux 2004; Villermaux & Clanet 2002, Villermaux & Rehab 2000) showed that acceleration normal to the interface at the rim causes RTI with wave-length of  . For a planar interface with uniform constant acceleration ( g ), the maximum amplified wave-number and growth rate are:   0.5 , 3 g D m g k              0.250.5 3 2 , 3 3 D RT g                (5.12) for g D  (Chandrasekhar 1961). The maximum gravity in equation 5.12 can be replaced by (Marmottant & Villermaux 2004; Kim et al. 2006): 2 2 .D U g a             (5.13) Where a is the amplitude of the primary wave generated by shear instability. Finding from the instability analysis provided in section 5.2.4.1, the dominant wave number for RTI is found as, 94       0.5 2 0 sin 3 g D d m D a U k U                   (5.14) for the analysis combining the potential flows and very thin vorticity layer. The influences of a and 0U on the RTI wave-number are demonstrated in figure 5.19. Fig. 5.19. Wave-number distribution around the droplet. a) The effect of amplitude of the primary (shear) instability. b) The effect of upstream velocity. Larger wave numbers for RTI are expected around the poles similar to the wave-numbers for KHI. Increasing the upstream velocity, the higher value of wave number is again expected to be distributed around the droplet due to higher acceleration it causes. Moreover, the influence of the wave amplitude on RTI is interestingly shown. Where higher wave-numbers are expected, increasing the height of the wave. The outcomes from theory show a good qualitative agreement with our observations in the 3D simulations. However, neglecting boundary layer generates a notable difference again as shown in table 5.4. 95 Table 5.4. Comparison of the present direct numerical simulations and equation 5.14. Weber Number  0 mD k  (DNS) 0 mD k  (Equation 5.14) Ratio 100 ~ π/2 ~35-40 1627 ~41-46 ~ π/3 ~29-32 1220 ~38-42 ~ π/4 ~20-24 813 ~33-39 As it can be seen in table 5.4, the value obtained by equation 5.14 is much larger than those achieved in direct simulations. This is again due to neglecting the boundary layer over the droplet. To involve the role of boundary layer (non-zero vorticity layer) we replaced the KHI wave-numbers from the linear stability analysis in section 5.2.4.2. After some algebraic manipulation we have:     0.5 2 2 225 481 , d g D D g m g a U F k                                     and (20)     0.5 23 41 , g g D D d m g a U k                                (21) based on equations 5.7 and 5.10, respectively. A comparison between the theoretical equations and the current numerical simulation is presented using the same wave employed in table 5.4. The results are listed in table 5.5. Table 5.5. Comparison of the present direct numerical simulations and equations 5.15 and 5.16. Weber Number 0 mD k  (DNS) 0 mD k  (Equation 5.15) – (Ratio) 0 mD k  (Equation 5.16) – (Ratio) 100 ~35-40 341(8.5-9.75) 342 – (8.5-9.75) As it can be observed in table 5.5, the predictions provided by equations 5.15 and 5.16 are much closer to the simulation results, in comparison with zero-vorticity layer theory. However, a significant 96 difference between the theories and simulations still exist. The difference between the theory above and experimental data (even for low Weber numbers) was shown in Marmottant & Villermaux (2004). Assuming a linear velocity profile to find the wave-number of primary KH waves as well as the assumption of plane interface for equation 5.12 might be two main reasons causing this difference. Nonetheless, the effect of different physical parameters are well accounted in the provided theories. 97 Chapter 6 Conclusion Deformation, disintegration, and dispersion of fragments of a single falling droplet as well as the growth of instabilities over a droplet in an external fluid flow are studied using direct numerical simulation. For the case of falling droplets (Chapter 3), outcomes are divided into three main parts; at the beginning, deformation and bag generation are investigated and the influences of Eötvös number on variation of parent drop diameter, aspect ratio, and temporal Weber and Reynolds numbers are indicated. Afterwards, the mechanism of bag breakup is studied in detail. The growth of waves over the bag and generation of the holes are shown. It is demonstrated that the number of punctures increase as the Eötvös number is increased. The generated holes retract by different mechanisms depending on the circumstances. The hole generation and retraction speeds are compared with the available theories. Studying the mechanism of bag or flattened core, it is shown that the mechanism generally includes: punctures generation, retraction of holes, formation of ligaments network, detachment of ligaments, disintegration of ligaments, and subsequently stable fragments creation. Breakup of bag and core are followed by deformation and disintegration of upper and lower tori due to Rayleigh instabilities. Subsequently, stable clouds of fragments are obtained for the simulations. The variation of fragments number, distribution of drop size and velocity, correlations between the size and velocity of fragments, and distribution of non-dimensional numbers are also studied in details as the statistical characteristics of fragments (Chapter 4). It is shown that the average size of the fragments is reduced by increasing Eötvös number. Moreover, based on the maps (for non- dimensional groups) it is described that the fragments are in a stable range. This part of results can considerably help to develop the statistical models for secondary atomization. For the case of droplets in a stream (Chapter 5), three two-dimensional and one three-dimensional study is performed for a range of Weber number corresponding to shear breakup mode. It is shown that Kelvin-Helmholtz Instability plays the main role in shear breakup. The wave-lengths of 98 instabilities are obtained for different Weber number and it is shown the wave-length of shear instabilities decreases with increasing the Weber number. In the three-dimensional simulation, the transverse azimuthal modulation is studied where the Rayleigh-Taylor Instability acts in the radial direction. Results for the both types of instabilities are compared with mathematical theories for zero and non-zero vorticity layer. The values obtained from the zero vorticity layer theory correctly demonstrate the trends of distribution of wave-numbers and their variations versus Reynolds and Weber numbers. However, the quantities obtained are far from those achieved in the direct simulations. These values are significantly improved using a linear velocity profile (considering the boundary layer). Nevertheless, a notable difference between the simulation and n-n-zero vorticity layer theory still exists. This difference was seen in the previous experimental studies (with the same ratio we obtained). To the best knowledge of our knowledge, the results presented in this thesis is the most precise simulation studying the deformation and fragmentation of droplets while previous investigations- mostly studied the deformation of the droplets in 2D or axisymmetric domains. The obtained results can be used in a wide range of natural and engineering applications, from rain phenomenon to combustion. 6.1 Future Work  Fragmentation of droplets in a turbulence medium happens in several industrial applications. As a future study, breakup of droplets in an isotropic homogeneous turbulence can be studies, using the same methodology.  For the case of a droplet in a stream, we focused on initiation and growth of instability. The fragmentation of droplets in a high speed flows can be investigated as a future work.  Evaporation is a vital phenomenon in combustion science. As a future study, evaporation of the liquid phase and its influence on hydrodynamics instabilities can be studied for falling droplets or droplets in a stream.  Collision of droplets occurs in several natural and industrial applications. As a future study, interaction of two or more droplets in a stream or in free-fall can be studied using the same methodology.  The obtained results demonstrate a detailed mechanism behind the fragmentation of droplets. The current secondary atomization models can be improved using the present findings. 99  Droplet tracking is crucial for development of atomization models. Finding an algorithm to track the fragments in Eulerian framework and its implementation by developing Gerris significantly assist for the further modelling.  Employing a high performance parallel computation technique, more complex problems can be studied. Thus, improvement of the current parallel computing abilities of Gerris is suggested as a future work. 100 Bibliography Ashgriz, N., Mashayek, F. 1995. Temporal analysis of capillary jet breakup, J. Fluid Mech. 291: 163- 190. Babinsky, E. and Sojka, P. E. 2002. Modeling drop size distributions, Prog. in Energy Combustion Sci. 28: 303-329. Bell, J. B., Colella, P. and Glaz, H. M. 1989. A second-order projection method for the incompressible Navier–Stokes equations, J. Comput. Phys. 85: 257–283. Bermond, N. and Villermaux, E. 2005. Bursting thin liquid films, J. Fluid. Mech. 524: 121-130. Brenner, M. P. and Gueyffier, D. 1999. On the bursting of viscous sheet, Phys. Fluids. 11: 737-739. Burger, M., Schwalbe, W., Kim, D. S., Unger, H., Hohmann, H. and Schins, H. 1983. Two phase description of hydrodynamic fragmentation processes within thermal detonation waves, J. Heat Transfer. 106(4): 728-735. Cao, X. K., Sun, Z. G., Li, W. F., Liu, H. F. and Yu, Z. H. 2007. A new breakup regime of liquid drops identified in a continuous and uniform air jet flow, Phys. Fluids. 19 (5): 057103. Chandrasekhar, S. 1961. Hydrodynamic and Hydromagnetic Stability, Oxford University Press, London. Chou, W. H. and Faeth, G. M. 1998. Temporal properties of secondary drop breakup in the bag breakup regime, Int. J. of Multiphase Flow. 24 (6): 889-912. Chou, W. H., Hsiang, L. P. and Faeth, G. M. 1997. Temporal properties of drop breakup in the shear breakup regime, Int. J. of Multiphase Flow. 23 (4): 651-669. Chorin, A. 1969. On the convergence of discrete approximations to the Navier–Stokes equations, Mathematics of Computation. 23 (106): 341-353. Culick, F. E. C. 1960. Comments on a ruptured soap film, J. Appl. Phys. 31: 1128. 101 Dai, Z. and Faeth, G. M. 2001. Temporal properties of secondary drop breakup in the multimode breakup regime, Int. J. of Multiphase Flow. 27(2): 217-236. Debrégeas, G., Martin, P. and Brochard-Wyart, F. 1995. Viscous bursting of suspended films, Phys. Rev. Lett. 75: 3886-3889. Debrégeas, G., De Gennes, G. P. and Brochard-Wyart, F. 1998. The life and death of viscous ―bare‖ bubbles, Science. 279: 1704-1707. Dimotakis, P. 1986. Two-dimensional shear layer entrainment, AIAA J. 24: 1791-1796. Dupré, M. A. 1867. Sixième memoire sur la theorie méchanique de la chaleur, Ann. Chem. Phys. 4(11): 194-220. Fakhri, A. and Rahimian, M. H. 2009. Simulation of falling droplet by the lattice Boltzmann method, Commun. Nonlinear Sci. Numer. Simulat. 14: 3046-3045. Fakhri, A. and Rahimian, M. H. 2011. Investigation of deformation and breakup of a falling droplet using a multiple-relaxation-time lattice Boltzmann method, J. Comp. Fluids. 40: 156-171. Feng, J. Q. 2010. A deformable liquid drop falling through a quiescent gas at terminal velocity, J. Fluid Mech. 658: 438-462. Fuster, D., Agbaglah, G., Josserand, C., Popinet, S. and Zaleski S. 2009. Numerical simulation of droplets, bubbles and waves: state of the art, Fluid Dyn. Res. 41: 065001. Gottesdiener, L., Gueyffier, D., Abdelouahab, M., Gatignol, R. and Zaleski, S. 2004. Numerical simulations of large falling drops, Int. J. Num. Methods Fluids. 45 (1): 109-123. Guildenbecher, D. R., López-Rivera, C. and Sojka, P. E. 2009. Secondary atomization, Experiments in Fluids. 46 (3): 371-402. Han J. and Tryggvason G., 1999. Secondary breakup of axisymmetric liquid drops. I. Acceleration by a constant body force, Phys. Fluids. 11(12): 3650. Harper, E. Y., Grube, G. W. and Chang I. D. 1972. On the breakup of accelerating liquid drops, J. Fluid Mech. 52 (3): 565-591. Helenbrook, B. T. and Edwards, C. F. 2002. Quasi-steady deformation and drag of contaminated liquid drops. Int. J. Multiphase Flow. 28: 1631–1657. 102 Hsiang, L. P. and Faeth, G. M. 1992. Near-limit drop deformation and secondary breakup, Int. J. of Multiphase Flow. 18 (5): 635-653. Hsiang, L. P. and Faeth, G. M. 1993. Drop properties after secondary breakup, Int. J. of Multiphase Flow. 19 (5) 721-735. Hsiang, L. P. and Faeth, G. M. 1995. Drop deformation and breakup due to shock wave and steady disturbances, Int. J. Multiphase Flow. 21: 545-560. Keller, J. B. and Kolodner, I. 1954. Instability of liquid surfaces and the formation of drops, J. Appl. Phys. 25: 918-921. Kim, D., Desjardins, M., Hermann, M. and Moin, P. 2006. Toward two-phase simulation of the primary breakup of a round liquid jet by a coaxial flow of gas, Center for Turbulence Research – Annual Research Briefs. Khosla, S., Smith, C. E. and Throckmorton, R. P. 2006. Detailed understanding of drop atomization by gas cross flow using the volume of fluid method, ILASS Americas, 19th annual conference on liquid atomization and spray systems, Toronto, Canada. Krzeczkowski, S. A. 1980. Measurement of liquid droplet disintegration mechanism, Int. J. Multiphase Flow. 6: 227-239. Lee, C. H. and Reitz, R. D. 2000. An experimental study of the effect of gas density on the distortion and breakup mechanism of drops in highs speed gas stream, Int. J. Multiphase Flow. 26:229-244. Lee, C. H. and Reitz, R. D. 2001. Effect of liquid properties on the breakup mechanism of high-speed liquid drops, Atomization Spray. 11:1-19. Li, R., Ashgriz, N., Chandra, S. and Andrews, J. R. 2008. Contraction of Free Liquid Ligaments, AICHE J. 54(12): 3084-3091. Marmottant. P. and Villermaux. E. 2004. On spray formation, J. Fluid Mech. 498:73-111. Mugele, R. A. and Evans, D. H. 1951. Droplet size distribution in sprays, Ind. Engr. Chem. 43 (6): 1317-1324. Mollot, D. J., Tsamopoulos, J., Chen, T. –Y. and Ashgriz, N. 1993. Nonlinear dynamics of capillary bridges: experiments, J. Fluid Mech. 255: 411-435. 103 Ni, M. –J., Komorib, S. and Morley, N. B. 2006. Direct simulation of falling droplet in a closed channel, Int. J. Heat Mass Transfer. 49 (1-2): 366-376. O’Brian, V. 1961. Why raindrops break up – vortex instability, J. Met. 18: 549-552. Pairam E. and Ferna´ndez-Nieves A. 2009. Generation and stability of toroidal droplets in a viscous liquid, Phys.Rev. Letter. 102: 234501. Pilch, M. and Erdman, C. A. 1987. Use of breakup time data and velocity history data to predict the maximum size of stable fragments for acceleration-induced breakup of a liquid drop, Int. J. Multiphase Flow. 13(6): 741-757. Popinet, S., 2003. Gerris: a tree-based adaptive solver for the incompressible Euler equations in complex geometries, J. Comput. Physics. 190(2): 572-600. Popinet, S., 2009. An accurate adaptive solver for surface-tension-driven interfacial flows, J. Comput. Physics. 228: 5838-5866. Prosperetti, A., 1981. Motion of two superposed viscous fluids, Phys. Fluids. 24: 1217. Rayleigh, L. 1882. On the instability of a cylinder of viscous liquid under capillary force, Philosophical Magazine, 34(207): 145–154. Rayleigh, L. 1883. Investigation of the character of the equilibrium of an incompressible heavy fluid of variable density, Proc. R. Soc. Lond. 14: 170-177. Summer, B. M. and Fredsøe. J. 2006. Hydrodynamics Around Cylindrical Strucures. Advanced series on ocean engineering, vol 26, World Scientific Publishing. Quan, S. and Schmidt, D. P. 2006. Direct numerical study of a liquid droplet impulsively accelerated by gaseous flow, Phys. Fluids 18: 102103. Rayleigh, L. 1891. Some applications of photography, Nature. 44: 249-254. Weber, C. 1931. Disintegration of liquid jets, Zeitschrift für Angewandte Mathematik und Mechanik, 11(2): 136-159. Reyssat, E., Chevy, F., Biance, A. L., Petitjean, L. and Quere, D. 2007. Shape and instability of free- falling liquid globules, EPL. 80: 34005. Savva, N. and Bush, J. W. M. 2009. Viscous sheet retraction, J. Fluid. Mech. 626: 211-240. 104 Shinjo, J. Adn Umemura, A. 2010. Simulation of liquid jet primary breakup: Dynamics of ligament and droplet formation, Int. J. Multiphase Flow. 36: 513-532. Simmons, H. C. 1977-a. The correlation of drop-size distribution in fuel nozzle sprays Part I: The drop-size/volume-fraction distribution, J. Engr. Power. 99: 309-314. Simmons, H. C. 1977-b. The correlation of drop-size distribution in fuel nozzle sprays Part II: The drop-size/number distribution, J. Engr. Power. 99: 315-319. Taylor, G. I. 1950. The instability of liquid surfaces when accelerated in a direction perpendicular to their plane.i, Proc. R. Soc. Lond. A. 201: 192-196. Wadhwa, A. R., Magi, V. and Abraham, J. 2007. Transient deformation and drag of decelerating drops in axisymmetric flows, Phys. Fluids, 19:113301. Villermaux, E. 1998. On the role of viscosity in shear instabilities, Phys. Fluids. 10:368-373. Villermaux, E. and Clanet, C. 2002. Life of a flapping liquid sheet, J. Fluid Mech. 462: 341-363. Villermaux, E. and Rehab, H. 2000. Mixing in coaxial jets, J. Fluid Mech. 425: 161-185. Villermaux, E. 2007. Fragmentation, Annu. Rev. Fluid Mech. 39: 419-446. Villermaux, E. and Bossa, B. 2009. Single-drop fragmentation determines size distribution of raindrops, Nature Phys. 5: 697-702. Zhao, H., Liu, H. F., Li, W. F. and Xu, J. L. 2010. Morphological classification of low viscosity drop bag breakup in a continuous air jet stream, Phys. Fluids. 22: 114103. "@en ; edm:hasType "Thesis/Dissertation"@en ; vivo:dateIssued "2012-11"@en ; edm:isShownAt "10.14288/1.0072826"@en ; dcterms:language "eng"@en ; ns0:degreeDiscipline "Mechanical Engineering"@en ; edm:provider "Vancouver : University of British Columbia Library"@en ; dcterms:publisher "University of British Columbia"@en ; dcterms:rights "Attribution-NonCommercial-NoDerivatives 4.0 International"@en ; ns0:rightsURI "http://creativecommons.org/licenses/by-nc-nd/4.0/"@en ; ns0:scholarLevel "Graduate"@en ; dcterms:title "Direct numerical simulation of fragmentation of droplets"@en ; dcterms:type "Text"@en ; ns0:identifierURI "http://hdl.handle.net/2429/42476"@en .