<%BANNER%>

Three-Dimensional Marker-Based Multiphase Flow Computation Using Adaptive Cartesian Grid Techniques


PAGE 1

THREE-DIMENSIONAL MARKER-BASED MULTIPHASE FLOW COMPUTATION USING ADAPTIVE CARTESIAN GRID TECHNIQUES By RAJKESHAR SINGH A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLOR IDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2006

PAGE 2

Copyright 2006 by Rajkeshar Singh

PAGE 3

To my family and teachers

PAGE 4

ACKNOWLEDGMENTS I would like to express my sincere gratit ude to Dr. Wei Shyy for providing me the opportunity and flexibility to perform this wo rk. I can not thank him enough for being so patient and understanding for years and pushing us to learn more. I would like to thank Dr. Renwei Mei, Dr. James F. Klausner and Dr. Samim Anghaie for agreeing to serve in my thesis committee. Since the surroundings dictate the quality of life in several ways, I thank my research group members and Dr. Siddharth Thakur for helping with academic aspects while providing a memorable company for the past four years. I thank my wife who has been the best fr iend for more than half a decade and has been extremely patient and understanding duri ng all these years. Finally I thank my parents who have been behind me every st ep of the way providing their unconditional support. iv

PAGE 5

TABLE OF CONTENTS page ACKNOWLEDGMENTS .................................................................................................iv LIST OF TABLES ...........................................................................................................viii LIST OF FIGURES .............................................................................................................x LIST OF SYMBOLS .........................................................................................................xv ABSTRACT ....................................................................................................................xvii 1 INTRODUCTION..........................................................................................................1 1.1 Objective and Contributions .................................................................................3 1.2 Outline of the Dissertation....................................................................................7 2 LITERATURE REVIEW...............................................................................................9 2.1 Interface Tracking.................................................................................................9 2.1.1 Lagrangian..................................................................................................9 2.1.2 Eulerian.....................................................................................................10 2.1.3 Mixed Eulerian-Lagrangian......................................................................11 2.1.4 Challenges and Recent Advances.............................................................11 2.1.5 Present Approach: Markers with or without Connectivity.......................16 2.2 Interfacial Dynamics Modeling..........................................................................17 2.2.1 Sharp Interface..........................................................................................18 2.2.2 Continuous Interface................................................................................19 2.2.3 Challenges and Recent Advances.............................................................19 2.2.4 Present Approach: Sharp or Continuous Interface...................................21 2.3 Adaptive Grid Computation...............................................................................22 2.3.1 Cartesian Grid Data Structures.................................................................24 2.3.2 Present Approach: Unstruct ured Adaptive Cartesian Grid......................26 3 MULTIDIMENSIONAL IMMERS ED BOUNDARY COMPUTATION..................27 3.1 Immersed Boundary Method..............................................................................27 3.1.1 Interfacial Conditions...............................................................................29 3.1.2 Material Property Smoothing...................................................................31 3.1.3 Momentum Source Term Computation....................................................33 v

PAGE 6

3.1.3 Interface Advection..................................................................................36 3.2 Conservative Restructuring for Interface Resolution.........................................37 3.3 Interface Reconstruction for Topology Change.................................................40 3.3.1 Reconstruction Criteria.............................................................................41 3.3.2 Reconstruction Procedure and Difficulties...............................................41 3.3.2.1 Two-dimensional interface.............................................................42 3.3.2.2 Three-dimensional interface...........................................................44 3.4 Summary.............................................................................................................48 4 ADAPTIVE CARTESIAN GRID BASED COMPUTATION....................................50 4.1 Data Structure.....................................................................................................50 4.2 Grid Adaptation..................................................................................................52 4.2.1 Interface Geometry based Adaptation......................................................53 4.2.2 Flow Solution based Adaptation..............................................................54 4.3 Immersed Boundary Solution Procedure............................................................55 4.3.1 Staggered Grid Arrangement and Approximations..................................57 4.3.2 Spatial Discretizati on of Velocity Equation.............................................58 4.3.3 Spatial Discretization of Temperature Equation......................................60 4.3.3.1 Sharp interface treatment...............................................................61 4.3.3.2 Interface temperature treatment.....................................................61 4.4 Summary.............................................................................................................63 5 RESULTS AND DISCUSSION...................................................................................65 5.1 Interface Restructuring: Interface in a Vortex Field...........................................65 5.2 Interface Reconstruction: Effect on Mass Conservation....................................68 5.3 Single Phase Navier-Stokes Solver....................................................................69 5.3.1 Decaying Vortex Field.............................................................................70 5.3.2 Lid-driven Cavity.....................................................................................72 5.3.3 Diffusive Heat Transfer between Concentric Cylinders..........................74 5.3.4 Natural Convection in a Square Cavity....................................................76 5.4 Spurious Velocity Currents.................................................................................77 5.4.1 Effect of Fluid Property Jump and Grid Resolution.................................78 5.4.2 Effect of Interface Reconstruction............................................................80 5.5 Buoyancy Driven Single Rising Bubble.............................................................81 5.5.1 Grid Convergence Test.............................................................................82 5.5.2 Rising Bubbles in Different Terminal Shape Regimes............................84 5.6 Coalescence of Two Rising Bubbles..................................................................87 5.6.1 Head-on Coalescence...............................................................................87 5.6.2 Off-axis Coalescence................................................................................90 5.7 Binary Drop Collision.........................................................................................91 5.7 Phase Change Computation................................................................................97 5.7.1 1-D Phase Change....................................................................................98 5.7.2 Stationary Bubble Growth........................................................................99 6 SUMMARY AND FUTURE WORK........................................................................113 vi

PAGE 7

6.1 Summary...........................................................................................................113 6.2 Future Work......................................................................................................115 LIST OF REFERENCES.................................................................................................118 BIOGRAPHICAL SKETCH...........................................................................................126 vii

PAGE 8

LIST OF TABLES Table page 1-1 Some industrial and natural oc currences of two-phase flows....................................2 1-2 Key marker based 3D tracking computations in literature.........................................5 2-1 Summary of some representati ve methods and related issues.................................14 2-2 Summary of recent advances in interface tracking..................................................15 5-1 Error estimates using Equation (5.2) for interface in a time reversed vortex field test comparing conservative and nonconservative restructuring...................................67 5-2 Volume error for reconstruction of a spherical interface and approximate radial perturbation required for correction.........................................................................69 5-3 Max velocity error and computation time normalized with time taken for 16 x16 uniform grid. The adapted grids are shown in Figure 5-4........................................72 5-4 Lid driven cavity computation using uniform and adaptive grids...........................73 5-5 Natural convection computation..............................................................................76 5-6 Pressure drop for an inviscid fluids with density ratio =1.0.....................................79 5-7 Effect of Laplace number on spurious velocity.......................................................79 5-8 Effect of density ratio on pressu re drop and spurious current on 50x50x50 grid (viscosity ratio = 1.0, La = 250)...............................................................................80 5-9 Effect of viscosity ratio on pressu re drop and spurious current on 50x50x50 grid (density ratio = 1000, La = 250)...............................................................................80 5-10 Grid convergence test for Eo = 1.0, M = 1.0x10-3. The density and viscosity ratios are set to 20.0...........................................................................................................83 5-11 Computation of rising bubbles with different terminal shape regimes....................87 5-12 Computational parameter fo r rising bubble-coalescence tests.................................88 5-13 Properties of tetradecane and nitrogen.....................................................................93 viii

PAGE 9

5-14 Binary drop collision computation with density ratio = 666.081 and viscosity ratio = 179.28....................................................................................................................94 5-15 Computational parameters for one-d imensional phase change computation...........98 ix

PAGE 10

LIST OF FIGURES Figure page 1-1 Illustration of a t ypical multiphase flow....................................................................1 1-2 Mixed Eulerian-Lagrangian tracking.........................................................................5 2-1 Illustration of Lagra ngian interface tracking............................................................10 2-2 Schematic of Eulerian tracking us ing level set and volume of fluid........................11 2-3 Interface in a time reversed vor tex field using level-set method.............................12 2-4 A conceptual interface reconstruction using piecewise linear segments in the volume-of-fluid method...........................................................................................13 2-5 A surgical alteration of the inte rface data to accomplish topology change..............15 2-6 A three-dimensional marker based film boiling computation..................................15 2-7 Illustration of conceptu al differences between mark er-based tracking with and without connectivity information.............................................................................17 2-8 Surface tension, pressure and normal component of shear stress acting on an interface....................................................................................................................18 2-9 Sharp and continuous interface methods..................................................................21 2-10 The effect of fluid property jumps on the computation time for a rising bubble.....23 2-11 A cell-by-cell and block-by-block grid adaptation..................................................25 2-12 A tree-based adaptive Cartesian grid-data...............................................................25 3-1 Schematic of the immersed boundary method.........................................................28 3-2 Probe based technique for temperature gradient computation on the interface.......30 3-3 A triangulated interface and the data-structure........................................................31 3-4 Indicator function and the material flags for cells completely outside or inside the interface....................................................................................................................33 x

PAGE 11

3-5 Momentum source term due to surfa ce tension force on a circular interface..........35 3-6 Computation of the unit normal and ta ngent vectors on interface triangles............35 3-7 Curvature computation for a unit circle using Equation (3.18) and cubic spline with (a) twelve and (b) sixty points..................................................................................36 3-8 Interface restructuring approach...............................................................................38 3-9 Conservative marker deletion for two-dimensional interfaces.................................39 3-10 Conservative marker deletion for three-dimensio nal interfaces...............................39 3-11 Level-contour based interface reconstruction..........................................................40 3-12 Interface reconstruction criterion uses two probes from the interface along the normal vector and bi-linearly inte rpolates the indicator values Ind1 and Ind2 ........41 3-13 Two dimensional reconstructed interface edges ( e 1, e 2 and e 3) and corresponding connectivity data......................................................................................................42 3-14 Illustration of 2D cell-by-cell interface reconstruction............................................43 3-15 A 2D cell with four markers creat ing two locally disconnected edges....................44 3-16 A reconstructed interface segment in a cell.............................................................45 3-17 Some common reconstructed interface segments and integer vertex-flags.............45 3-18 Difficulties due to degeneracy of reconstructed interface data and implemented remedy......................................................................................................................46 3-19 Orientation setting of reconstructed interface..........................................................47 3-20 Examples of 3D interface reconstruction.................................................................47 4-1 Cartesian grid data-structure....................................................................................51 4-2 Adaptive grid refinement.........................................................................................52 4-3 Cartesian cells sharing a face and the corn er cells are not allowed to differ by more than one level of refinement.....................................................................................53 4-4 Geometry based adaptation......................................................................................54 4-5 Staggered grid arrangement.....................................................................................58 4-6 Spatial discretizati on of convection term.................................................................59 xi

PAGE 12

4-7 A one-dimensional interface between i and i+1 ......................................................61 4-8 Linear interface temper ature correction within one cell on each side of the interface....................................................................................................................62 4-9 Presently employed interface temperature correction technique solves a correction equation sharply on each side of the interface.........................................................63 5-1 Interface in a time reversed vortex field...................................................................67 5-2 Error estimate and grid convergence for interface in time reversed vortex field.....67 5-3 Grid convergence test for volume error due to reconstruction.................................69 5-4 Decaying vortex field...............................................................................................71 5-5 Grid convergence test for decaying vort ex field showing maximum velocity error computed on uniform and adaptive grids.................................................................71 5-6 Streamlines for lid-driven cavity simulation for Re = 100.......................................72 5-7 Lid driven cavity......................................................................................................73 5-8 The u velocity profile along a ve rtical line through the center of lid-driven cavity using (a) uniform grids and (b) adaptive grids.........................................................74 5-9 Indicator contour used for modeling the cylindrical computational domain...........75 5-10 Computed temperature and error for h eat transfer in a concentric cylinder............75 5-11 Natural convection with Pr = 0.71 and Ra = 1.0x105..............................................76 5-12 Static bubble computation........................................................................................78 5-13 Grid convergence tests for static bubble..................................................................79 5-14 Effect of reconstruction fre quency on spurious velocities.......................................81 5-15 Computational set-up for a single rising bubble......................................................82 5-16 Aspect ratio computation.........................................................................................83 5-17 Bubble rise with Eo = 1.0, M = 1.0x 10-3, density and viscosity ratio of 20.............84 5-18 Computational domain fo r cases in Table 5-11.......................................................85 5-19 Computed shapes from Table 5-11..........................................................................86 5-20 Computed case for M = 0.971, Eo = 97.1................................................................86 xii

PAGE 13

5-21 Computational set-up for head-on coalescence test.................................................88 5-22 Head-on coalescence................................................................................................89 5-23 Grid and streamlines at t = 0.125 sec for head-on coalescence................................89 5-24 Instance of interface reconstruction during head-on coalescence of two rising bubbles at t = 0.116 s................................................................................................90 5-25 Experimental photographs of bubble coal escence taken at 0.03s time intervals.....90 5-26 Off-axis coalescence................................................................................................91 5-27 Binary drop collision regimes marking th e location of computed parameters with solid circles...............................................................................................................93 5-28 Impact parameter B is defined as B=h/d..................................................................93 5-29 Binary drop collision Case 1 from Table 5-14.......................................................100 5-30 A cross section of the computational grid and the time history of grid data-size for Case 1 from Table 5-14..........................................................................................101 5-31 Binary drop collision Case 2 from Table 5-14.......................................................102 5-32 A cross section of the computational doma in showing the interface and the velocity vectors for Case 2 from Table 5-14........................................................................102 5-33 A cross section of the computational grid and the time history of grid data-size for Case 2 from Table 5-14..........................................................................................103 5-34 Binary drop collision Case 3 from Table 5-14.......................................................104 5-35 A cross section of the computational grid and the time history of grid data-size for Case 3 from Table 5-14..........................................................................................105 5-36 A cross section of the computational doma in showing the interface and the velocity vectors on the left and pressure profile on the right for Case 3 from Table 5-14..106 5-37 Binary drop collision Case 4 from Table 5-14.......................................................107 5-38 A cross section of the computational grid and the time history of grid data-size for Case 4 from Table 5-14..........................................................................................108 5-39 A cross section of the computational doma in showing the interface and the velocity vectors for Case 4 from Table 5-14........................................................................109 5-40 One dimensional phase change setup.....................................................................110 xiii

PAGE 14

5-41 Comparison of computed and theoretical solution for Case 1 in Table 5-15.........110 5-42 One-dimensional phase change comput ation for Case 2 in Table 5-15 showing interface location on the left and temperat ure distribution at t = 0.1 on the right..111 5-43 Stationary bubble growth rate................................................................................112 xiv

PAGE 15

LIST OF SYMBOLS ,, Density, viscosity and surface tension ,, CK Heat capacity, thermal conductivity and latent heat of evaporation m Interfacial mass transfer rate (ve for evaporation) ,, kfn Interface local curvature, surf ace tension force, unit normal vector s F Momentum source term due to surface tension force I Indicator function D Dirac-delta and discre tized Dirac-delta function d Bubble diameter (,,) Uuvw Velocity P Pressure T Temperature g Gravity int A property on interface markers 1 A property of fluid outside the interface 2 A property of fluid inside the interface 1 2 Ratio of a fluid property xv

PAGE 16

s atT Saturation temperature (,,) X xyz Spatial coordinates along x, y and z coordinate axes L Length scale ReUL Reynolds number 2UL We Weber number 2U Fr L g Froude number 1 2CT Ja Jacob number CUL Pe K Peclet number Pr C K Prandtl number 2L La Laplace number 23CgTL Ra K Rayleigh number maxU Ca Capillary number 4 3 g M Morton number 2 gL Eo Etvos number xvi

PAGE 17

Abstract of Dissertation Pres ented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy THREE-DIMENSIONAL MARKER-BASED MULTIPHASE FLOW COMPUTATION USING ADAPTIVE CARTESIAN GRID TECHNIQUES By Rajkeshar Singh August, 2006 Chair: Wei Shyy Major Department: Mechanical and Aerospace Engineering Multiphase flows associated with interfacial dynamics, steep jump in fluid properties and moving boundaries between different phases pose substantial computational challenges in terms of both m odeling as well as computational cost. The present work uses an immersed boundary technique to model the interfacial dynamics, Lagrangian markers to track the moving phase boundaries, and a stati onary Cartesian grid to solve all the flow governing equations. Time dependent triangulated surface meshes are employed to represent the time dependent interfaces shape and location. Based on the solution characteristics, multi-level, three-dimensional adaptive grid techniques are incorporated into the computational framewor k to help meet the resolution requirements. Furthermore, a conservative marker redistribu tion technique is deve loped to maintain a desired marker-spacing, and a connectivity preservi ng level contour-bas ed reconstruction technique is devised to handle topological changes associated with the interfacial dynamics. The flow equations were solved using the projection method with a finite xvii

PAGE 18

volume staggered grid formulation on adaptive grids. For phase change problems, accuracy of the mass transfer computation cr itically affects the overall computational outcome. Efforts have been made to address th is via a sharp interface -based mass transfer mode combined with the immersed boundary method. The capabilities and accuracy of the individual components and overall computat ional system are tested with a range of computations including demonstration of improvements with conservative interface restructuring, reconstruction and its e ffect on immersed boundary solution accuracy, estimation of computational cost saving with adaptive grids, rising bubble coalescence, binary drop collision, and stationary bubble growth in a superheated liquid pool. The method has demonstrated its capability for handling high density ratio, O (1000), multiphase fluid dynamics. xviii

PAGE 19

CHAPTER 1 INTRODUCTION Multiphase flows are characterized by two or more fluids in relative motion. The constituent fluids usually have different physical properties separated by an interface that evolves with the flow. Figure 1-1 depicts some of the defining elements of a typical multiphase flowa moving/defo rming surface that separate s different fluids/phases; discontinuous fluid properties (density, visc osity, conductivity, etc. ) across the interface; interfacial effects such as surf ace tension, phase change, etc. Fluid 1 Fluid 2 Surface tension force Interface Fluid 1 Fluid 2 Surface tension force Interface Figure 1-1. Illustration of a typical multiphase flow. Presence of such flows in every day life can be recognized by simple examples of ocean waves, boiling water, blood flow in the body. Some specific practical applications may be cited as boiling heat transfer in pow er industries, biochemical and metallurgical industries, cavitation and ultr asound technology for medical and industrial applications. The constituents in these flows may be liquidliquid, gas-liquid, liquid or gas with solid. Table 1-1 by Kleinstreuer (2003) gives a summary of some of the industrial or natural occurrences of such flows. Since cons ideration of all the combinations in Table 1-1 is a 1

PAGE 20

2 prohibitively difficult task, the scope of the present work is controlled by focusing only on incompressible liquid-liquid or gasliquid flows for bubble and drop dynamics problems. Table 1-1. Some industrial and natura l occurrences of two-phase flows. Phase combination Phase configuration Occurrence Gas-Solid Dispersion of solid particles in gas Aerosol, particulate matter pollution; filters and particle-collection devices; Pneumatic transport of particles; gasfluidized beds Liquid-Solid Dispersion of particles in liquid Cells and particles in blood; hydraulic transport of particles; liquid-fluidized beds Dispersion of liquid drops in gas Mist, fog, cloud; coalescence of drops and formation of rain; droplet removal devices; fuel droplet injection; spray coating and cooling Gas-Liquid or Liquid Gas Dispersion of gas bubbles in liquid Transport of oil and gas; foam formation; evaporators, boiling flow in pipe Dispersion of liquid drops in liquid Emulsion and creams; break-up of immiscible drops; coalescence of drops and phase separation Liquid-liquid Both bodies of liquids connected Transport of two liquids in horizontal pipe; movements of oil-water interface in porous rock. Note: Adapted from Kleinstreuer (2003) Presence of multiple time and length scales makes numerical simulation of such flows a difficult and expensive affair. Time dependent interfaces with topology changes (bubble coalescence or breakup) and interfacia l mass/heat transfer add further to the computational complexity and cost. Recent ye ars spanning more than a decade have witnessed an ever increasing computational ac tivity with the res earch spanning various aspects ranging from enhancement of computational efficiency to interface tracking and flow modeling. The computational expense an d difficulties in handling interfaces are in part responsible for relatively fewer computati ons in 3D as compared to a vast number of 2D or axisymmetric cases. These limitations are even more visible with phase change

PAGE 21

3 computations. Besides the high computational cost, phase change co mputations present several other difficulties with the mathema tical modeling and com putational algorithms. Some of the modern day phase change com putations are Dhir (2001), Son (2001), Son et al. (2002) and Luo et al (2005) using level-set; We lch and Wilson (2000) using volume-of-fluid; Esmaeeli and Tryggvason (2004 ), Tryggvason et al. (2005), Shin et al. (2005) and their earlier work us ing marker based tracking. Although continuing developments in computational techniques and resources are making the numerical simulation more and more tractable, several algorithmic and modeling challenges remain unresolved. With further discussion, it will be highlighted in the following chapter that some of the co mputational aspects often have conflicting requirements and demand several compromises to be made based on the problem of interest. 1.1 Objective and Contributions The primary goal of the present work is to develop an effective multidimensional interfacial flow computati on-system. The process involve s assembling key ideas from different approaches and techniques into a unified framework along with suitable improvements and modifications wherever necessary. Development of a parallel computing platform is an essential requireme nt for large scale computing. However, the current efforts are part of a larger research group activity and thus do not cover parallel computing aspects in this document. The computational component s are divided into three generic categories: inte rface tracking; flow modeling; flow computation. The general emphasis and objectives in these categories are: Multi-fluid interface tracking o Time dependent triangulated surface grids for interface representation. o Maintain resolution by locally adding /deleting nodes on the surface grid.

PAGE 22

4 o Capability to handle complex topology changes. Interfacial dynamics and flow modeling o Single fluid formulation for the mass and momentum conservation equations. o Surface tension forces modeled as smoothed volumetric source term in the momentum equation. o Accurate computation of the interfacial mass transfer rate for phase change problems. Adaptive grid computation o Grid generation and dynamic adaptation based on interface location and flow features. o Adaptive grid base incompressibl e Navier-Stokes computation. The interface is tracked on a stationary Cartesian grid using time dependent triangulated surfaces ( Figure 1-2 ). This approach is also referred to as a marker based tracking method. Since all th e flow computations are performed on the background stationary grid and the interface is tracke d using Lagrangian markers, the overall approach is of a mixed Eulerian-Lagrangian type. Although this approach has several advantages due to explicit knowledge of th e interface shape/location for computing the geometric properties, the difficulties in data-handling have deterred a wider usage. Table 1-2 presents some of the key usages and techniques demonstrated to handle complex multidimensional interfacial flows with marker based tracking. The present method is similar to that of Tryggvason et al (2001) with primary contributions made in interface handling. The triangulated-interface is restructured by locally adding or deleting markers to avoid undesirably long or short edges. The commonly used restructuring techniques in literature (Tryggvason et al., 2001, 2005; Sousa et al. 2004) introduce some amounts of error in the interface volume for which a conservative algorithm is proposed. Topology changes were handled us ing the level-contour based interface r econstruction of Shin and Juric (2002) with a ppropriate algorithmic consider ations to establish a valid interface connectivity data.

PAGE 23

5 Interface markers Fluid 2 Fluid 1 Interface markers Fluid 2 Fluid 1 (a) (b) Figure 1-2. Mixed Eulerian-Lagrangian track ing. (a) A two dimensional interface using connected markers on a background grid and (b) a triangulated interface. Table 1-2. Key marker based 3D tr acking computations in literature. 3D Explicit Tracking Author Summary of technique A set unconnected triangles Shin and Juric (2002) Easy algorithms for computing geometric properties and dealing with complex topology changes. Reduced flexibility in contro lling interface resolution and dealing with multiple interfaces. A set of unconnected points Torres and Brackbill (2000) Easy to handle complex topology changes but interface geometry computation more involved that Shin and Juric (2002) Reduced flexibility in contro lling interface resolution and dealing with multiple interfaces. Triangulated surface grid. Tryggvason et al. and group (2001, 2005) Flexible control over interface resolution and dealing with multiple interfaces. Uses reconstruction technique of Shin and Juric (2002) to handle topology changes. The interface moves on the background Cart esian grid that is dynamically adapted based on the interface location and flow soluti on to optimize the computational cost. All the flow governing equations are solved on background grid. The data-structure and Cartesian grid generation is based primarily on the work of Aftosmis (1997). The NavierStokes solution is performed using projection method with a staggered grid based finitevolume formulation. The underlyi ng non-uniform nature of the grid makes staggered grid relatively more difficult to implement than the corresponding colloca ted grid (Singh &

PAGE 24

6 Shyy, 2006), but provides a compact pressu re-velocity coupling. As detailed in subsequent chapters, several simplifying appr oximations in the numerical discretization were made for a convenient flow solver im plementation on the sp atially adaptive grid. These approximations are present only on the coarse-fine grid-interfaces which are kept away from the physical location of the in terface by suitable refinement criteria. Most of the interfacial dynamics mode ling and development uses the already available concepts from literature. The unde rlying single fluid formulation for mass and momentum equation with marker based tracking is commonly recognized as the immersed boundary method (Peskin, 1977; Tryg gvason et al., 2001). This method smears the fluid properties across the interface a nd models the surface tension forces as smoothed volumetric source terms in the momentum equation. Although the mass and momentum equations use immersed boundary method, it is possible to avoid such approximations for the temperature equation. Since interfacial mass transfer depends on the jump in temperature gradients across the in terface, an attempt has been made to treat the temperature equations sharply in order to accurately compute the mass transfer rate. The immersed boundary method of Juric and Tryggvason (1998) uses an iterative procedure to maintain the interface temperat ure while the present method uses a simple linear correction si milar to the level-set based computations of Morgan (2005). At this point, the current approach deviates from a conventional immersed boundary method and borrows the philosophy from recent developmen ts in level-set and sharp interface ghost fluid method of Luo et al. (2005) and Gibou et al. (2003). The presented method solves the temperature equation on both sides of the interface and does not require assuming fixed temperature inside the interface as in the work of Son (2001).

PAGE 25

7 The specific contributions of the pr esent work may be listed as: 1. Marker based interface tracking a. Conservative interface restructuring t echnique to maintain resolution without introducing errors in the interface volume. b. Development of interface reconstruc tion algorithm for topology change. Several challenges/difficulties are high lighted along with simple remedies. c. Conservative restructuri ng along with reconstructi on technique offers an enhanced flexibility in handling mu ltiple interfaces as compared to a connectivity-free track ing method of Shin and Juric (2002). 2. Interfacial dynamics modeling a. Implementation of multidimensional immersed boundary method. b. A sharp interface treatment for temperat ure equation for accurate interfacial mass-transfer computation during phase change. Strategies to incorporate level-set based sharp interface techni ques within the present marker based tracking are presented. 3. Adaptive grid computation a. Implementation of an adaptive Cart esian grid generator with dynamic adaptation based on the interface location and flow solution. b. Development of a staggered grid, fi nite volume based Navier-Stokes computation algorithm for velocity and pressure computation. c. Development of a sharp in terface, collocated grid computation procedure for temperature equation in heat-transfer/phase cha nge problems while using marker based tracking. The capabilities of the developed approach and algorithms have been demonstrated via a set of complex multiple rising bubbles an d binary drop collision computations as one of the few known cases in literature usin g marker based tracking. For example, the only known 3D binary drop col lision with marker based tracking are by Nobari and Tryggvason (1996); Shin and Juric (2002) perfor med for moderate density ratios of O(10) while present computations were performed for flows with density ratios of O(1000). 1.2 Outline of the Dissertation The document is organized into six chapters. Chapter 2 presents brief review of alternative approaches in literature. Th e basics of immersed boundary modeling and related key components dealing with interface management are presented in Chapter 3.

PAGE 26

8 Chapter 4 presents the adaptive Cartesian grid generation technique and discusses the Navier-Stokes solution procedure highlig hting various simplifying assumptions in the numerical discretization. The key results are presented in Chapter 5. Tests characterizing and demonstrating interface restructuring, reconstruction, singl e and multiple rising bubbles and binary drop collision are presented along with prelimin ary assessment of the accuracy of phase change computation. Chapter 6 provides a summary of the presen ted work and highlights various issues, approximations following a brief comment on future extensions and refinements.

PAGE 27

CHAPTER 2 LITERATURE REVIEW This chapter provides a brie f review of some of the popular techniques based on the categories outlined in the previous chapter: interface tracking; interfacial dynamics modeling; adaptive grid based computation. Va rious pros and cons of these techniques are visited in an attempt to place the current choices and hi ghlight some of the relevant issues and difficulties to be addressed in following chapters. 2.1 Interface Tracking As reviewed by Shyy et al. (1996), the three main categories are Lagrangian, Eulerian and mixed Eulerian-Lagrangian. Du e to the extensive amount of available literature, the scope of this section is limited to providing a compact background and some of the key issues with these tec hniques along with the state-of-the-art developments. 2.1.1 Lagrangian Lagrangian method discretizes the comput ational domain with the interface as a boundary for creating a body fitted grid ( Figure 2-1 ). This approach requires adjustment or regeneration of the grid in order to accommodate changes in the interface geometry. Some usages of this approach can be found in the works of Ryskin and Leal (1984), Raymond and Rosant (2000), Wasekar and Mang lik (2003) using structured curvilinear grids and work of Perot and Nallapati (2003) using triangulated (t etrahedrons in 3D) grids. 9

PAGE 28

10 Although the body fitted nature of the grid makes it a desirable choice, handling of body fitted grids require sophisticated grid generation capabilities especially when the interfaces undergoes topology change. Observing the frequency of the usage, this is not among the most popular choices and hence it will not be considered further. Fluid 1 Fluid 2 Interface Figure 2-1. Illustration of Lagrangian interface tracking. 2.1.2 Eulerian Eulerian methods use a scalar function to represent the interface implicitly. The level-set (Osher & Fedkiw, 2001, 2003) a nd volume-of-fluid (Youngs, 1982; Rider & Kothe, 1998; Scardovelli & Za leski, 2002) are examples of Eulerian tracking. This method captures the interface using the inform ation contained in a scalar function. This scalar function is a signed distance function in the leve l-set method and volume fraction in the volume-of-fluid method ( Figure 22 ) Unlike the Lagrangian methods, Euleri an methods track the interface on a stationary grid usi ng a transport Equation (2.1) for the scalar func tion evolution that contains the interface information. Due to th e Eulerian nature of the interface, this method does not necessitate any changes in the computational grid as required by Lagrangian methods. 0 V t (2.1)

PAGE 29

11 (a) (b) Figure 2-2. Schematic of Eulerian tracking using level set and volume of fluid. (a) A volume of fluid captures th e interface (solid curve) using volume fraction information; (b) level set recognizes the interface as zero contour of signed distance function 2.1.3 Mixed Eulerian-Lagrangian This class of method derives its name from the fact that it tracks the interface using Lagrangian markers moving on an Eulerian computational grid ( Figure 1-2 ). These markers move on the comput ational grid using a si mple advection Equation (2.2) where the marker velocities are interpolated from the Eulerian grid. Thus, (similar to Lagrangian) any geometric information about interface is explicitly available and (similar to Eulerian methods) the interface evoluti on does not require any modification to the computational grid (Tryggvason et al., 2001; Glimm et al., 1998, 2000; Singh & Shyy, 2006). int intX U t (2.2) 2.1.4 Challenges and Recent Advances The above introduced methods for interface tracking have some associated advantages and disadvantages that make gene ral preference of any particular method difficult. Table 2-1 gives a summary of the key featur es of these methods. The Eulerian methods are among the easiest to implement and some of their most cited advantages lie

PAGE 30

12 in their natural ability to handle interfaces undergoing topology changes. On the flip side, level-set method (LS) suffers from errone ous mass loss that may become unacceptably large and detrimental. Such a mass-loss problem is highlighted by the errors shown in the tests conducted by Enright et al. (2002) where a spherical interface undergoes severe deformation and it is expected to return to original shape when placed in a time-reversed vortex field. Figure 2-3 (a) shows the sever loss in volume of a spherical interface in such a field using level-set method for interface tr acking. This drawback was ameliorated by Enright et al. (2002) using a particle leve l-set method that uses Lagrangian markers around the interface to avoid e rroneous breakup/deletion of th e level-set characteristics ( Figure 2-3 (b)). (a) (b) Figure 2-3. Interface in a time reversed vorte x field using level-set method. (a) Severe mass loss with level set method; (b) im provements produced by particle levelset method (Enright et al., 2002). [Reprinted with permission] Unlike the level-set method, the volumeof-fluid tracks the volume (mass) explicitly and hence does not suffer from mass-loss problems. However, volume-of-fluid methods have been observed to produce s purious flotsam (floating wreckage) and jetsam (jettisoned goods) in the regions of high interface-curvatu re under-resolved by the computational grid (Rider & Kothe, 1998) The most serious difficulty with the VOF is due to errors in computation of in terface curvature required for surface tension

PAGE 31

13 computation. These difficulties arise due to the sharply changing nature of the volumefraction and difficulties in reconstructing a continuous interface using the volumefraction information ( Figure 2-4 ). A more comprehensive detail of various VOF based interface reconstruction and subs equent developments could be found in Rider and Kothe (1998), Scardovelli and Zaleski (2002), Renardy and Renardy (2002), Meie r et al. (2002). Figure 2-4. A conceptual interface reconstruc tion using piecewise linear segments in the volume-of-fluid method. The mixed Eulerian-Lagrangian methods po ssess several desirable features when compared with Lagrangian or Eulerian me thods. Due to explicit definition of interface using markers, the interface does not inaccurately diffuse in time and loose mass as level set methods. It has also been noticed that e xplicit tracking does not require as high grid resolution as Eulerian methods (Glimm et al., 1998). Explicit definitio n of the interface also avoids inaccuracies caused by the e rrors in interface reconstruction for VOF methods. The most serious concern against the Eule rian-Lagrangian methods has been raised due to the complexity of interface data that becomes difficult to handle especially for topology changes. These difficulties arise due to the algorithmic complexity involved in surgical altering of the interf ace data to accommodate topo logy changes. The source of difficulty lies in establishing valid connec tivity information. Although such a procedure

PAGE 32

14 is easily accomplishable for two-dimensional interfaces ( Figure 2-5 ), extension to three dimensions is a much more involved proce dure. The potential complexities in markerbased tracking may be highlighted from the inte rface shapes seen in so me of the 3D film boiling computations of Shin and Juric (2002) ( Figure 2-6 ). Several attempts have been made by re searchers to simplify the marker based algorithm by resorting to elimination of the need to maintain the connectivity information. One of the first such attempts is due to Torres and Brackbill (2000) who used a set of unconnected markers to describe the interface. Lack of the need to maintain explicit connectivity information, this appro ach makes it easier to handle topological changes. The subsequent development of connectivity free track ing using unconnected triangular elements by Shin and Juric (2002) simplified the algorithms further. Table 2-2 presents some key developments by comb ining various indivi dual techniques to exploiting their attractive features. Table 2-1. Summary of some represen tative methods and related issues. Methods Advantages Issues Lagrangian Interface explicitly known and part of the computational grid. Geometric property computed directly on interface. Frequent grid regeneration. Difficult and expensive for complex 3D computations. Eulerian (Level-set) Easy geometry-computation. All computations on Eulerian grid. Easy to handle topology changes. Severe mass loss/gain in underresolved areas of the interface. Eulerian (Volume of Fluid) No mass loss/gain. All computations on Eulerian grid. Easy to handle topology changes. Poor performance in underresolved areas of the interface.. Difficult geometry computation. Mixed EulerianLagrangian Interface explicitly known even at sub-grid levels. Flow computations on Eulerian grid. Complex data-structure and book-keeping for 3D interfaces. Difficult for topological change.

PAGE 33

15 Table 2-2. Summary of recent advances in interface tracking. Base method Improved method Features Level-set (LS) Particle level-set (Enright et al., 2002) Ameliorates the excessive/erroneous mass loss using a cloud of markers around level-set interface. Volume of fluid (VOF) Hybrid marker-VOF (Aulisa et al., 2003, 2004) Employs makrer based interface to avoid interface reconstruction for geometric-property computation. Improved performance in under-resolved areas. LS and VOF Coupled LS -VOF (Sussman, 2003) Excellent mass-conservation due to VOF. Easy geometry-computation using Level-set. Marker tracking The point-set tracking (Torres & Brackbill, 2000) A set of unconnected. Lack of connectivity-infor mation makes it easier to handle topological changes. Geometry-computations require spline-fitting. Connectivity free tracking with levelcontour reconstruction (Shin & Juric, 2002) A set of unconnected triangles. Possible to compute geometry (curvature, normal) without any surface fitting or interpolation. Level-contour reconstruc tion for topology change. merger breakup 1 4 3 2 10 6 9 8 7 5 11 4 3 2 10 6 9 8 7 5 1 merger breakup 1 4 3 2 10 6 9 8 7 5 1 4 3 2 10 6 9 8 7 5 11 4 3 2 10 6 9 8 7 5 1 11 4 3 2 10 6 9 8 7 5 1 Figure 2-5. A surgical alteration of the interface data to accomplish topology change. (a) (b) Figure 2-6. A three-dimensiona l marker based film boiling computation. (a) A snapshot of the computational domain and (b) a triangulated interface before and after pinch-off from Shin and Juric ( 2002). [Reprinted with permission]

PAGE 34

16 2.1.5 Present Approach: Markers with or without Connectivity Since it is evident that the Eulerian type approaches are relatively easier, it makes the choice of the method difficult based on any single considerati on. Since the most serious difficulty with marker based trackin g is due to data-handling, the explicit knowledge of the interface shape/location and r ecent algorithmic advan ces in literature have been the criteria for choosing explicit marker based method. The method of level-contour based rec onstruction with connectivity-free markerbased tracking of Shin and Juri c (2002) offers an attractive choice. However, the lack of connectivity information on the interface makes it difficult to locally alter the markerspacing by marker addition and deletion to maintain a desired marker-spacing on the interface. As seen in Figure 2-7 (a), deletion of an edge requires collapsing of vertices of the neighboring edges (point p2 and p5) but the lack of connectivity makes the identification of edges e1 and e3 difficult i.e., no explicit in formation is available to suggest, for example, that nodes p2 and p3 have the same physical location. For computations such as bubbly flows by Es maeeli and Tryggvason (1998, 1999), several bubbles (~100 or more) move in close proxim ity without merger. In these simulations, the bubbles are not allowed to merge with e ach other but local restructuring of the interface is still requi red to maintain desirable resolu tion. The connectivity-free method resorts to global reconstruction of the interface even to maintain the marker spacing and may inadvertently produce mergers. For such reasons, the current work maintains the connectivity information ( 2 e Figure 2-7 (b)) to allow greater flex ibility in handling multiple interfaces. This decouples the procedure for handling interface resolution and interface reconstruction by using local restructuring when needed and performing reconstruction only to achieve topology changes.

PAGE 35

17 p6 p1 p5 p4 p3 p2 e1 e3 e2 p6 p1 p5 p4 p3 p2 e1 e3 e2 p6 p1 p4 p2 e1 e3 e2 p6 p1 p4 p2 e1 e3 e2 (a) (b) Figure 2-7. Illustration of conceptual differences between marker-based tracking with and without connectivity information. (a) Connectivity free method of Shin and Juric (2002) does not store any connectiv ity information between edges; (b) connectivity based approach maintains such information explicitly. 2.2 Interfacial Dynamics Modeling A moving/deforming interface in the com putational domain presents some obvious difficulties in numerical computation of the flow field. The difficulties arise primarily due the nonlinear coupling of fluids and steep fluid-property jumps across the interface. In the event of phase change, effects such as mass and heat transfer across the interface also need to be considered. For simplicity, the basic introduction is provided here with respect to isothermal flows. Issues and techniques related to phase change modeling will be considered in subsequent sections and chapters. In general, the pressure and viscous stresses show disc ontinuities across the interface that is related to the surface tension force and fluid property jumps. Figure 2-8 shows the pressure (P) and the normal shear stress components ( n n ) across the interface and Equation (2.3) relates the jump in flow pr operties with the surface tension force. Based on the level of detail and accuracy in treatment of th e discontinuities, the existing methods can be classified into tw o categories: sharp interface method (SIM); continuous interface method (C IM). The sharp interface me thod (SIM) is a class of techniques that attempt to satisfy the jump condition (Equation (2.3) ) explicitly while

PAGE 36

18 continuous interface method (CIM) represents techniques that smooth the jump over few cells across the interface. 2121() PPnn (2.3) Figure 2-8. Surface tension, pressure and normal component of shear stress acting on an interface. 2.2.1 Sharp Interface There are three prominent techniques with in this category: sharp interface of Ryskin and Leal (1984), Ye et al. (1999, 2001); immersed interface (LeVeque & Li, 1994; Li & Lai, 2001); ghost fluid (Liu et al., 2000; Kang et al., 2000). The technique of Ryskin and Leal (1984) and Ye et al. (1999, 2001) are the most uncompromising in terms of satisfying the jump conditions across th e interface. The basic approach can be explained with the cut-cell method of Ye et al. (1999, 2001). This method cuts the background Cartesian grid with the interface and creates sub-domains ( Figure 2-9 (a)) that act like body fitted grids with arbitrary shaped cells on the interface. For stability reasons, the cut-cells are merged with suitable nei ghbors to remove arbitrarily small cells. The governing equations are solved separately in each sub-domain w ith the corresponding fluid properties. At every time st ep, an iterative procedure, by perturbing the interface in local normal direction, is followed to find th e shape/location of th e interface such that Equation (2.3) is satisfied along with the mass and momentum conservation equations.

PAGE 37

19 The immersed interface and the ghost-fluid techniques use the jump conditions to construct boundary conditions for flow variab les on the interface. As an example, the ghost-fluid technique solves th e governing equations in each phase separately by creating ghost-cell values using the analytical jump condition. Unlike the sharp interface method of Ye et al. (1999, 2001), these techniques do not solve for the interface shape and flow field simultaneously or create cut-cells. 2.2.2 Continuous Interface In contrast to the sharp interface method, the continuous interface method smoothes the fluid property and models the surface tension forces as a smoothed momentum source term (Peskin, 1977). This facilitates a single fluid formulation for the entire domain with the flow/fluid properties varying smoothly across the interface over a thin zone. This smearing/smoothing region across the interface shown in Figure 2-9 (b) is usually three to four cells in thickness. The continuous inte rface method has been widely used with mixed Eulerian-Lagrangian tracking (i.e., the immersed boundary method) (Peskin, 1977; Tryggvason et al. 2001; Francois & Shyy, 2003 ) as well as Eulerian tracking such as level-set and volume-of-fluid. 2.2.3 Challenges and Recent Advances The iterative procedure of coupling the pha ses across the interface and the complex procedure of cell-cutting makes the extension of the sharp interface technique of Ye et al. (1999, 2001) to three dimensions a tedious a nd challenging task. Although the ghost fluid is among the simplest sharp interface method to use, it has traditionally been used with level-set based interface repres entation. Although successful demonstration of ghost fluid type approach with Eulerian-Lagrangian tracking were made by Hao and Prosperetti

PAGE 38

20 (2004) for some bubble dynamics problems, furthe r research is required to establish the overall stability and performance. Accuracy studies by Ye et al. (2004) and Shyy (2004) comparing the sharp interface and continuous interface methods for a spherical drop in static equilibrium sheds some light on their characteristic differences. With a spherical drop in static equilibrium, the difference in pressure is balanced by surface tens ion force according to the Young-Laplace Equation (2.4) The sharp interface method attempts to capture the pressure jump as a discontinuity whereas the continuous interface method smoothes the pressure jump over few cells ( Figure 5-12 (b)). P (2.4) For a static bubble computation, the theo retical solution produces uniformly zerovelocity in the domain. However, due to certain amount of net numerical stressimbalance across the interface, computations tend to produce a non-zero velocity field termed as parasitic currents or spurious velocities ( Figure 5-12 (a)). These spurious velocities can destabilize the in terface especially if viscous forces are not strong enough (Rider & Kothe, 1995). As an observation, the magnitude of the spurio us velocity can be kept orders of magnitude smaller with the sharp interface treatment. The second order accurate sharp interface technique of Ye et al. (2004) could keep this velocity up to machine level as compared to continuous interface methods produci ng velocities of the order 10 -4 (in terms on Capillary number Ca) (Tryggvason, 2001). The first order accurate ghost fluid technique could also k eep these velocities to the order of 10 -6 or lower (Kang et al., 2000). Several attempts using improved curvature/surface-tension

PAGE 39

21 computation have also been attempted in lite rature to keep such erroneous velocities as low as possible (Renardy & Renardy 2002; Shin et al., 2005). cut by interface Fluid 2 Fluid 1 Fluid 2 Fluid 1 Sub-domain 2 Fluid 2 Sub-domain 1 Fluid 1 (a) Fluid properties smeared Fluid 1 Fluid 2 Fluid properties smeared Fluid 1 Fluid 2 (b) Figure 2-9. Sharp and continuous interface methods. (a) Sharp interface cut-cell method and (b) a continuous interface usi ng the immersed boundary technique. 2.2.4 Present Approach: Sharp or Continuous Interface The continuous interface method provides a numerically stable system due to smoothing of discontinuities but leads to smeari ng of the interfacial effects. This aspect of the continuous interface method poses accuracy concerns for certain types of problems like solidification dynamics at morphological s cales that require hi gher precision in interface treatment (Udaykumar et al., 1999). For the reasons of computational efficiency and the relative difficulties in using the sharp interface method with marker based tracking, the present work employs the continuous interface method via immersed b oundary technique. Despite the accuracy

PAGE 40

22 issues related to smearing of discontinuities, it is still a reasonably accurate method for a large number of drop dynamics problems. Over the years, the immersed boundary technique has been used extensively to study a wide range of physical processes: Peskin (1977) for studying blood flow in heart-valves ; Juric and Tryggvason (1996) for dendritic solidification; Agresal et al (1998) for deformation and a dhesion of circulating cells; Esmaeeli and Tryggvason (1998, 1999) for si mulation of bubbly flows; Francois and Shyy (2002, 2003) for micro-scale drop dynamics. The overall interfacial dynamics treatment in the presence of phase change is of a hybrid nature i.e., the mass and momentum equations are solved using the immersed boundary method, but the temperat ure equation is solved in a sharp fashion. The primary intention, as pointed out in previous chapter, is to accomplish higher accuracy in interfacial mass transfer comput ation. Due to the absence of a surface tension type source term in the temperature equation, sharp tr eatment of the temperature equation can be performed with relative ease as compared to a sharp treatment of the surface tension forces, pressure, etc. Since th e interface has been a ssumed to be at saturation temperature at all times, it simplifies the process by allowing the interface temperature to be used as a boundary condition in the numerical discretiz ation of the thermal diffusion term in temperature equation. This discretization approa ch is similar to the recent developments by Gibou et al. (2003) and Morgan (2005) us ing ghost fluid based sharp interface for phase change computations. 2.3 Adaptive Grid Computation Most of the computational time is devot ed to performing computation/operations on the interface and solution of the flow-govern ing equations. Usually the latter part is responsible for majority of the computational burden. Since the governing equations of

PAGE 41

23 the flow are incompressible, the procedure involves solution of an elliptic equation to compute the pressure field making it one of the primary reasons for the high computational cost. Since the system has multip le fluid components, the jump in the fluid properties, in general, has an adverse bearing on the overall computation-time. Figure 210 from the work of Francois et al. (2004) is a representative case of degrading performance with increasing fluid propert y jump. For brief information, these computations were performed for a buoyanc y-driven rising bubble using the immersed boundary method. The Reynolds number, Weber number and Froude number were set to 100, 4 and 1.0. A W-cycle multigrid technique was used to solve the pressure Poisson equation. The horizontal axis in Figure 2-10 shows the level of multigrid used for computation where level = 1 represents computation without employing any multigrid acceleration. (a) (b) Figure 2-10. The effect of fluid property jumps on the computation time for a rising bubble. (a) Effect of density ratio; (b) e ffect of viscosity ra tio (Francois et al., 2004). [Reprinted with permission]

PAGE 42

24 The performance seems to be affected the most by the density ratio factor: without multigrid, the computation time for density ra tio of 1000 is nearly ten times higher in Figure 2-10 than for a density ratio of 10. Al though these observations come from a specific method and problem, they still repres ent the general trend (F rancois et al., 2004; Lrstad & Fuchs, 2004) and recommend a ppropriate measures to reduce the computational cost as much as possible. The presented work uses adaptive Cartesian grids to resolve the flow features locally and thus attempts to reduce the overall computational data-size. 2.3.1 Cartesian Grid Data Structures The approach for generating adaptive Cartesia n grids is divided into two categories: a cell-by-cell adaptation ( Figure 2-11 (a)); a block-by-block adaptation ( Figure 2-11 (b)). Adaptive mesh refinement known as AMR is an example of a block-by-block adaptation (Berger & Oliger, 1984; Sussman, 1999, 2005) th at interprets the gr id as a union of uniformly refined grid-blocks (rectangles in 2D and hexagons in 3D). Due to the approach of collecting cells in blocks and refining them, AMR tends to produce a larger grid data-size than necessary. As observed by Aftosmis (1997), typically only 70% of the cells in a rectangular patch may actually need refinement. In the sense of reducing the computational data-size, cell-by-cell adap tation is a more flexible approach. A cell-by-cell adaptation uses either a tree-based data (Wang, 1998; Popinet, 2003; Losasso et al., 2004) or an unstructured data format (Aftosmis, 1997; Ham et al., 2002; Singh & Shyy, 2005, 2006) to store the gr id. A tree-based data-structure ( Figure 2-12 ) offers a compact storage system. All the conn ectivity information is implicitly contained in the location of cells in the tree-structure. One of the most common queries in a

PAGE 43

25 numerical simulation is the ce ll-cell connectivity information needed, for example, to compute convection/viscous fluxes. Since no e xplicit information about the identities of the neighbor cells is stored, it needs to be extracted from th e information contained in the tree. As reported by Aftsosmis (1997), a pproximately 10-20% of the CPU time for CFD solvers is dedicated to tree traversals answering grid-connectivity queries. (a) (b) Figure 2-11. A cell-by-cell and block-by-bloc k grid adaptation. (a) A cell-by-cell adaptation in the present work; (b) sche matic of a block-by-block adaptation. 1 1 34 2 1 34 25 76 1 1 34 2 1 34 25 76 1 7652 Level=0 Level=2 Level=1 4321 1 7652 Level=0 Level=2 Level=1 4321 (a) (b) Figure 2-12. A tree-based adaptive Cartesian gr id-data. (a) A cell und ergoing two levels of successive refinement; (b) data stru cture showing the place of cells based on the refinement level.

PAGE 44

26 2.3.2 Present Approach: Unstruct ured Adaptive Cartesian Grid An anisotropically ad aptive grid (Ham et al., 2002; Singh et al., 2005) offers an optimized grid data by using a directionbased refinement strategy. However, the resulting data-structure and solution algor ithms become difficult and tedious to implement. For simplicity and convenience, an isotropic adaptati on is employed where a cell marked for refinement is simultaneously sp lit in all coordinate directions. The flowgoverning equations and numeri cal discretization aspects ha ve been discussed in the related chapters. The grid is stored using unstructured data with cell-by-cell refinement approach. It is believed that the extra memory requireme nts due to explicit storage of connectivity information will be offset in part by compact storage system exploiting the Cartesian nature of the grid and easier accessibility of the data for the flow solver.

PAGE 45

CHAPTER 3 MULTIDIMENSIONAL IMMERSED BOUNDARY COMPUTATION The immersed boundary method originally proposed by Peskin (1997) has been used to model the surface tension force as a smoothed volume source term in the momentum equation. The jump in density a nd viscosity are smoothed within few cells across the interface. The isothermal immers ed boundary computation is similar to the work of Tryggvason et al. (2001). This chap ter presents the key components of the immersed boundary modeling and interface track ing including algorithms for maintaining interface resolution and technique for handli ng topology changes. The numerical solution technique including adaptive grid generation an d equation discretization are presented in the following chapter. 3.1 Immersed Boundary Method The incompressible Navier-Stokes equati ons for mass and momentum and energy conservation represented by Equations (3.1) (3.2) and (3.3) are solved on the stationary background grid. The source terms F s models the surface tension as a smoothed momentum source term using a sm oothed Dirac delta function ( Figure 3-1 (a)).. The fluid properties are smoothed using an indicator function I that varies smoothly from a value of unity outside the interface to zero inside ( Figure 3-1 (b)). The density and viscosity are smoothed using Equation (3.4) (Prosperetti, 2002) where subscripts 1 and 2 refer to the fluid outside and inside the interface, respectively. The source term Q s in the energy equation accounts for the heat source due to phase change. Since the interface is considered at satu ration temperature, the job of this term in 27

PAGE 46

28 the immersed boundary method of Juric and Try ggvason (1998) is essentially to maintain the interface temperature. Despite smeared treatment of the hydrodynamics equations, the temperature equation in the presented work is solved in a sharp fashion to compute the mass-transfer accurately. Since the temperature equation treatment is different than the traditional immersed boundary method, role of the source term Q s and corresponding numerical solution procedur e will be visited in the following chapter. 0U t (3.1) T sU UUPUUFg t (3.2) s T CTUKT t Q (3.3) 212 212I I (3.4) Bound of numerical delta function Region of influence Actual interface I = 0 I = 1 (a) (b) Figure 3-1. Schematic of the immersed boundary method. (a) Interfac ial effects smoothed over a thin zone across interface; (b) indicator function for fluid property smoothing.

PAGE 47

29 3.1.1 Interfacial Conditions The mass conservation equation across th e interface is written as Equation (3.5) The jump in pressure across th e interface is related to the su rface tension force, jump in shear stress and fluid ve locities using Equation (3.6) Using the mass continuity equation, the jump in the normal component of the velocity can be written using Equation (3.7) The tangential component of the velocity is cont inuous across the inte rface and the normal component is also continuous if there is no mass transfer. 11int 22int() ()UUnm UUnm (3.5) 22 2 2 11 1 1 2 int 1 int nn nnUUUPnnUUUPnn n (3.6) 12 1211 () UUnm m (3.7) The mass transfer rate across the in terface is computed using Equation (3.8) for which a probe based technique is used to co mpute the required temperature gradients on the interface. These gradients are computed using Equation (3.9) and Figure 3-2 showing two probes of length 1.2 along the interface norma l direction. The term is the background grid size and T int is interface temperature (set to saturation temperature T sat ). The temperature T 1 and T 2 on the probes are computed usin g bilinear interpolation from the background. 12QKTnKTn (3.8) 1int 1 int2 21.2 1.2 TT KTn TT KTn (3.9)

PAGE 48

30 Using the mass continuity across the interf ace, the interface velocity can be written as: 12 12 int22 UU Unnm (3.10) The interface velocity due to mass-transfer component is computed with the help of Equation (3.8) and unsmoothed fluid properties. The interface velocity component due to fluid velocity is simply interpolated from the background grid using a Dirac delta function described later. The velocities in individual phases are considered inco mpressible and the jump in velocity is smoothed using Equation (3.11) (Juric & Tryggvason, 1998). The effect of mass transfer is modeled using locally non-divergence-free velocities near the interface. The Equation (3.12) relates the divergence of the veloci ty field with the interfacial masstransfer after taking divergence of Equation (3.11) and applying incompressibility condition for velocities in the individual phases where U 1 and U 2 are fluid velocities outside and inside the interface. 212UUUUI (3.11) 12 1211 UUUIm I (3.12) T1T2Vapor phase Liquid phase Figure 3-2. Probe based technique for temper ature gradient computation on the interface.

PAGE 49

31 3.1.2 Material Property Smoothing Figure 3-3 shows a common finite-ele ment type data-structure used to store the triangulated interface. It assi gns a unique identification (i nteger number) to each node (marker) and stores its coordinates along with the triangle-nodes. All the triangles are assigned a unique orientation such that the right hand thumb rule using the ordering of nodes of a triangle gives the local normalvector pointing outside the interface. Node ID x y z 1.. .. .. 2.. .. .. .. .. .. nnode.. .. .. Triangle ID node1 node2 node3 1 ...... 2...... ....... ntria...... Figure 3-3. A triangulated inte rface and the data-structure. The indicator function is computed at Cartesian cell centers and linearly interpolated to cell faces whenever requir ed. It is computed by solving the Poisson Equation (3.13) where ( X X int) is a Dirac-delta function non-zero only at X=X int and D is its discretized form for which Equation (3.14) (Peskin & McQueen, 1996; Griffith & Peskin, 2005) is used. Since the fluid properties are constant in most part of the domain, the indicator function is solved only within few cells (4~5) across the interface. Fo r rest of the cells, the indicator value is set ma nually. The cells completely outside the interface (Fluid 1) are assigned a material flag MatFlag(cell) = 1 and cells completely inside the interface (Fluid 2) are set to MatFlag(cell) = 2 ( Figure 3-4 ). This knowledge is used to manually set the indicator function to zero or one and use it as a Dirichlet boundary condition for Equation (3.13)

PAGE 50

32 f2 (x-x) interfacentriangle ( ) I dA nDrdA (3.13) ()()() 1 32141 : if [0,1) 8 1 52743 : if (1,2] 8 0 : else xxyyzz xyzrrr Dr hhhh r hhh h h (3.14) A ray-tracing algorithm popular in computer-graphics a nd computational-geometry community could be employed for accurate evalua tion of the state of a cell i.e., inside or outside the interface. However, since such accu rate information is required only outside the smearing zone, a simple and efficien t method based on the painters algorithm frequently used in compute r-graphics rendering is used. Unlike the ray-tracing algorithm, the painters algorithm does not require expe nsive computation of three-dimensional linesurface intersection and it is su fficient for the current purposes to deal with simply connected closed topologies. Open surfaces ca n also be dealt with ease following suitable algorithmic considerations. This algorithm stores the state of a cell in an array named MatFlag ( cell ). This approach exploits the fact that a cell outside the interface has all its neighbors either outside or on the interface. The outline of the implemented algorithm is presented below: Set MatFlag ( cell ) = 2 for all the cells contains a marker else set it to zero. Select any cell outside all interfaces, call it SeedCell and put it as the first entry in a temporary array named SeedCellList ( ). Set TotalSeedCell = 1; SeedCellList ( TotalSeedCell ) = SeedCell Set MatFlag ( Seedcell ) = 1 Use a temporary variable NextCheck to march through the array SeedCellList( ). Set NextCheck = 1

PAGE 51

33 Do while ( NextCheck TotalSeedCell ) o SeedCell = SeedCellList ( NextCheck ) o Go through the neighbor cells which share a face with SeedCell If ( MatFlag ( NeighborCell ) > 0)Go to another neighbor cell TotalSeedCell = TotalSeedCell + 1 SeedCellList ( TotalSeedCell ) = NeighborCell MatFlag ( NeighborCell ) = MatFlag ( SeedCell ) ( = 1) o NextCheck = NextCheck + 1 End Do While Loop Set MatFlag ( Cell ) = 2 for all the cells which still have MatFlag ( Cell ) = 0 Fluid 2 MatFlag(cell)=2 I = 0.0 Fluid 1 Fluid 2 MatFlag(cell)=1 I = 1.0 Figure 3-4. Indicator function a nd the material flags for cells completely outside or inside the interface. 3.1.3 Momentum Source Term Computation The surface tension force is computed on th e interface triangles and distributed to the background grid as smoothed source term F s in Equation (3.2) Figure 3-5 shows the smoothed nature of the momentum source fo r a two-dimensional interface. The surface tension force on a discretized in terface element (curves in 2D and triangles in 3D) can be evaluated in several ways: computation with Equation (3.15) where unit normal vector and curvature can be computed using curve fitting for two-dimensional interfaces (Francois & Shyy, 2002, 2003; Ye et al., 2001) and surface fitting for three-dimensional interfaces (Sousa et al., 2004); co mputation using a line integral form shown in Equation (3.16) and fitting curves/surfaces to obtain no rmal and tangent vectors (Al-Rawahi & Tryggvason, 2002; Tryggvason et al., 2001).

PAGE 52

34 A f kndA (3.15) ()AAs f kndAnndAtnds (3.16) There are two important observations to be made here: the net surface tension force on a closed surface should be zero (conser vation); curvature computation using interpolation based methods are numerically sensitive and often requires some form of data smoothing (Francois & Shyy, 2002, 2003; Ye et al., 2001; Sousa et al., 2004). The use of Equation (3.15) does not enforce conservation whereas the line-integral form does not require explicit curvature computation and maintains the conservation. The present work uses the line integral form and comput es the local normal and tangent vectors along the triangle edges using the simple appro ach of Al-Rawahi and Tryggvason (2004) shown in Figure 3-6 To ensure conservation (i.e., the net surface tension force on a closed interface is zero), the cross-pr oduct of normal and tangent on a triangle edge is obtained by averaging its value between the triangles sharing this edge. The surface forces computed on triangle edges are distributed directly to the background grid using Equation (3.19) If required, the curvature can be computed using Equation (3.18) Such a simple technique is seen to produce sufficient accuracy demonstrated by Figure 3-7 comparing curvature of a unit circle using present method and a cubic-spline interpolation. The overall accuracy of the present surface tension force computation and modeling has already been demonstrated by Shin and Juric (2002) for boiling flows and Al-Rawahi and Tryggvason (2004) for dendritic solidification. 1,2,3edge edge f tnds (3.17)

PAGE 53

35 f n dA (3.18) () sedge F fDr edge (3.19) (a) (b) (c) Figure 3-5. Momentum source term due to su rface tension force on a circular interface. (a) The x-component, (b) y-component and (c) the x-component along a horizontal line though th e center of circle. p1 p2 p3 12 12 tt n tt triangle triangle => p1, p2, p3 32 1 32 p p t p p 13 2 13 p p t p p 21 3 21 p p t p p Figure 3-6. Computation of the unit normal and tangent vectors on interface triangles.

PAGE 54

36 5 10 0.995 1 1.005 1.01 1.015 1.02 1.025 1.03 1.035 1.04 1.045 1.05 1.055 1.06Current method SplineCurvatureInterface node 0 10 20 30 40 50 60 1 1.0005 1.001 1.0015 1.002Current method Spline Interface node 0 10 20 30 40 50 60 1 1.0005 1.001 1.0015 1.002Current method Spline Interface node (a) (b) Figure 3-7. Curvature computation for a unit circle using Equation (3.18) and cubic spline with (a) twelve and (b) sixty points. 3.1.3 Interface Advection The interface velocity in Equation (3.10) has a component due to the fluid motion and another due to phase change. The compon ent due to phase change is obtained by computing the mass transfer rate on the in terface markers, as described earlier. The component due to fluid velocity is interpol ated from the backgr ound grid and the final marker velocity expression is: 12 int 12() () 2 UUDrdxdydzmn cell cell (3.20) For simplicity, the marker movement due to background fluid velocity is solved using 2 nd Order Runge-Kutta method but the velocity due to phase change is treated using forward Euler method as shown in Equation (3.21) *() () 12 12* *** ** 1 22nx xn XXtU XXtU n XX n X ntm (3.21)

PAGE 55

37 Error associated with interpolated marker velocity and integration of the advection equation produces some error in the interface volume. These errors are usually small but can accumulate over time and become unacceptably large. To correct such defects, an explicit phase-volume correction step is pe rformed by perturbing the markers in local normal direction and applying bisection met hod to find a location th at keeps the volume error within some prescribed level (~10 -7 %). This correction may be performed either at every time-step or periodically. All the pres ented computations performed this operation at every time step and usually one or two bise ction iterations were found to be sufficient. For flows without phase change, the interf ace volume does not change and hence the amount of volume error and hence the correc tion is always known. For phase change computations, the amount of volume error and thus the correction was estimated based on the mass-transfer rate. 3.2 Conservative Restructuring for Interface Resolution Due to the limited support of the Dirac delta function used for distribution of surface forces and other interpolation operations, excessively long edges (e.g. more than two Cartesian-cell-size) will appear as a hole on the interface. On the other hand, excessively small edges can produce unphysical s ub-grid undulations. To ameliorate such effects, markers are continuously added and deleted based on the following criteria: _ Break edge (marker addition) _ Delete edge (marker deletion) 3 EdgelengthCellsize Cellsize Edgelength (3.22) It is required that the marker addition and deletion causes minimal disturbances on the interface and avoids introducing volum e errors. However, the commonly employed techniques in literature (e.g. Figure 3-8 ) do not preserve the interface volume. More

PAGE 56

38 precisely, it is the marker deletion that is non-conservative. Such non-conservative behaviors are demonstrated by Sousa et al. (2004) who encountered up 5% volume errors due to interface restructuring. (b) (a) p1p2 p3 (b) (a) p1p2 p3 p3 Figure 3-8. Interface restructuring approach. (a ) Long edges are split at the mid-point and (b) a small edge p 1 p 2 is collapsed to its midpoint p 3 The presently employed restructur ing approach is same as in Figure 3-8 with an added correction step to the marker deletion procedure in order to preserve the interface volume. This correction-algorithm is an ex tension of the volume preserving marker relocation of Sousa et al. (2004) who used such an approach for interface smoothing and removal of sub-grid undulations on the inte rface. The underlying procedure is explained using a 2D example shown in Figure 3-9 that shows a marker removal by collapsing the edge formed by point p 2 and p 4 to a point p 3 resulting in a net phase -area loss. To correct this error, first a reference point p ref at p 3 n is selected ( Figure 3-9 (b)) where n is the local unit normal vector at p 3 A reference phase area of the polygon p 1 p 2 p 4 p 5 p ref is computed and point p 3 is relocated after deleting the edge to recover preserve the reference area ( Figure 3-9 (c)). The conservative marker/edge deleti on on 3D surfaces is similar to its 2D counterpart and it is summarized below:

PAGE 57

39 1. Select a point p at the midpoint of the edge to be deleted. A local normal vector n at this point is taken as the average of the normal vector on the triangles connected to this edge. 2. Define a reference point p ref as p-n and compute a reference volume by adding the tetrahedral-volumes made by the reference point and the triangles in 0 Figure 310 (a). 3. Collapse the edge to its mi dpoint and compute a volume by adding the volumes of tetrahedrons made by point p with the triangles in Figure 3-10 (b). 4. Relocate the point p to 0pn for preserving the local phase volume 0 p1p4p5p3p2prefn p3(a) (b) p1p4p5p2n (c) p3 p1p5pref Figure 3-9. Conservative marker deleti on for two-dimensional interfaces. (a) Old interface p 1 p 2 p 4 p 5 and new interface p 1 p 3 p 5 with non-conservative method; (b) define a reference point ( p ref ) and compute a referen ce area; (c) relocate p 3 in local normal direction to pr eserve the reference area. p n pref p n pref (a) (b) Figure 3-10. Conservative marker deletion for three-dimensional interfaces. (a) A reference point p ref is defined at p-n ( p is edge mid-point) and compute a reference, (b) collapse edge to point p and move p along normal vector to preserve the reference volume.

PAGE 58

40 3.3 Interface Reconstruction for Topology Change When two interfaces come within close proximity they are merged together and similarly when sections of the same interf ace are in close proxim ity the interface is broken. The term close proximity at this in stance is intended to be around one grid-cell size and will be addressed further in the doc ument. Since the indicator function reflects the interface shape, it is used to redefine the interface to achieve topology changes. Figure 3-11 shows a 2D interface and indicator contours along with indicator = 0.5 contour exhibiting interface merger. The work of Shin and Juric ( 2002) exploited this idea to reconstruct the in terface by deleting the interface data completely and reconstructing a new interface at indicator = 0.5 contour. Their work tracked the interface using a connectivity free method and needed to perform reconstruction to control interface resolution as well as facilitate topol ogy changes. Such a situation does not arise in the current work as the resolution is co ntrolled via conservative restructuring while reconstruction is performed only to achieve topology change. -0.1 -0.05 0 0.05 0.1 (a) (b) (c) Figure 3-11. Level-contour based interface recons truction. (a) A circle of radius 0.1 and three indicator contours corresponding to le vels 0.2, 0.5 and 0.8, (b) interfaces in close proximity and corresponding (c) indicator = 0.5 level contour.

PAGE 59

41 3.3.1 Reconstruction Criteria The need for reconstruction is checked at prescribed interval s and interfaces are reconstructed based on some indication of a possible topological change. A simple criterion shown in Figure 3-12 sends two probes along the local normal vector from the marker points. The length of each probe is se t to one and half times the background gridcell-size. The indicator functi on value on the probes is inte rpolated from the background grid using bilinear interpolation. If both the probes have indicator function value greater than half or less than half, representing merger or breakup fo r the given grid resolution, reconstruction is carried out. The need for reco nstruction is checked at frequent intervals depending on the problem. When uncertain, it was checked at every few (e.g. 50) time steps in the presented computations. Ind2 Ind1 Figure 3-12. Interface reconstruction criterion uses two probes from the interface along the normal vector and bi-linearly in terpolates the i ndicator values Ind1 and Ind2. 3.3.2 Reconstruction Procedure and Difficulties The approach for reconstruction is to first create a set of globally numbered markers, one on each edge with indicator = 0.5 contour intersecting it. Cells with new markers are visited one at a time and markers on the cell-edges are connected to define interface segments within the Cartesian cells ( Figure 3-13 ).

PAGE 60

42 The following sections present the overall reconstruction procedure using a twodimensional interface and the ideas are extended to three-dimensional interfaces. Some critical difficulties in maintaining the integr ity of the connectivity data are highlighted and simple remedies are provided. e1 => (p1, p2) e2 => (p2, p3) e3 => (p3, p4) p1 p2 p3 p4 e1 e2 e3 Figure 3-13. Two dimensional r econstructed interface edges ( e 1, e 2 and e 3) and corresponding connectivity data. 3.3.2.1 Two-dimensional interface The indicator function values computed at cell centers are linea rly interpolated to Cartesian cell-vertices. The vertices are flagge d and a rule for marker creation on an edge is set as: 1 : indicator0.5 (outside interface) vertex_flag 0 : otherwise (inside/on the interface) Marker is created on an edge only if its two vertices have different flag values. The above rules avoid ambiguities and si mplify the algorithm by eliminating the need to look into markers in neighboring cells during reconstruction in a given cell. As an example, consider Figure 3-14 (a) where an edge e 2 is constructed only while visiting cell 1 and not while in cell 2 even though the corresponding mark ers physically lie on the edges of both these cells but are consider ed only on the vertical edges of cell 1 due to their opposite vertex flag values. This procedure, however, allows creation of edges with zero lengths in order to create a valid connectivity data. Such a situation arises when, for

PAGE 61

43 example, the reconstructed markers p 1 p 2 and p 3 in Figure 3-14 (b) all same physical location i.e., indicator = 0.5 contour is at the Cartesian-cell vertex. No special attention is paid to such cases as a valid connectivity data is still produced while such edges/elements are later deleted by the marker/edge de letion procedure presented earlier. e1e3e2cell2cell1cell3cell4 1 0 0 0 0 0 0 1 1 1 1 1 e1e3e2cell2cell1cell3cell4 1 0 0 0 0 0 0 1 1 1 1 1 e1e2p1p3p2cell2cell1 1 1 0 1 1 e1e2p1p3p2cell2cell1 1 1 0 1 1 (a) (b) Figure 3-14. Illustration of 2D cell-by-cell interface reconstruction. (a) Markers on edges with opposite vertex-flag values (0 and 1) are connected to form edges e 1 e 2 and e 3 in cells cell 3 cell 1 and cell 4 respectively; (b) markers p 1 p 2 and p 3 are connected in cell 1 and cell 2 irrespective of their physical coordinates. A 2D cell can have two, three or four markers. Cells with two markers produces single edge; cell with three markers produce two edges connected within the cell; cell with four markers are used to produce two edge s non-connected inside th e cell. It is last case that causes some difficulty in deciding th e orientation of the two disconnected edges. These four markers are not used to create a cl osed interface within the cell as it is not required since the interface will be closed in some other cell easily deduced by using the vertex flag information. It can be seen that if the interface was to be reconstructed in triangulated cells, such ambiguities with four markers do not arise as a triangular cell can have only two or no markers producing at the most single edge w ithin the triangle. Although, the markers are

PAGE 62

44 created only on Cartesian cell edges, a set of locally triangulated cells are used to decipher the connectivity/orientation of th e edges whenever needed as shown in Figure 315 where the one vertex of the triangles in the cell is at the Cartesian cell-center. The indicator function is known at the cell cente r and hence a vertex flag value can be assigned. 1 1 1 0 0 1 1 1 0 0 1 0 1 0 0 1 0 1 0 0 Figure 3-15. A 2D cell with four markers cr eating two locally disconnected edges. The connectivity of markers is decided based on the vertex flags of the temporarily triangulated cell. 3.3.2.2 Three-dimensional interface The basic approach for three-dimensi onal interfaces is similar to the twodimensional counterpart i.e., create globally numbered markers on cell edges and connect within one cell at a time. Figure 3-16 shows a schematic reconstruction by creating markers on edges that have opposite vertex fl ags and connecting them to create polygons. The reconstructed polygons are broken into triangles by inserti ng a marker at the geometric center of the polygon. Figure 3-17 shows some typical r econstructed polygons with corresponding vertex-fla gs. The last two cases in Figure 3-17 contain six markers each but one has a single six-sided polygon wh ile the other has two triangles. These two cases with six markers are distinguished from each other by inspecting the relative location of edges containing the markers.

PAGE 63

45 p1 0 0 0 0 0 1 1 1 p4p3p2p5 p1p3p2 p4p5 p1 0 0 0 0 0 1 1 1 p4p3p2p5 p1p3p2 p4p5 (a) (b) Figure 3-16. A reconstructed interface segmen t in a cell. (a) Reconstructed markers are connected to produce a polygon (p 1 p 2 p 3 p 4 p 5 ); (b) polygon is broken into triangles with a new marker at the geometric center. 0 0 0 0 0 0 0 0 1 0 1 0 0 1 0 0 1 0 1 0 0 3 markers 6 markers 4 markers 5 markers 0 0 1 1 0 1 1 0 1 0 1 1 1 1 1 0 0 0 1 1 0 1 1 0 1 0 1 1 1 1 1 0 1 1 Figure 3-17. Some common rec onstructed interface segments and integer vertex-flags. Since the indicator function itself is smooth in nature, polygons shown in Figure 317 are observed the most often. Similar to 2D reconstruction, occasionally two or more disconnected polygons within a ce ll are observed. Presence of more than six markers in a 3D cell presents an obvious algorithmic difficu lty in terms of how to define the surface segments using the cloud of markers. As shown in Figure 3-18 (a) eight markers could be connected in at least two ways to define tw o differently oriented and connected polygons.

PAGE 64

46 If the reconstruction were to be conducted in tetrahedrons, only thr ee or four markers and hence just one polygon is present ( Figure 3-18 (b)). Similar to 2D case, whenever required the cells are locally broken into tetrahedrons ( Figure 3-18 (c)) by connecting the cellcenter to cell-vertices in order to decide the connect ivity. Although only the markers created on Cartesian edges are used to cons truct the interface, possible presence of markers on tetrahedron edges (based on vertex flags for tetrahedrons) are used simply as a means to decide how to connect marker on Ca rtesian cell edges in case there are more than six of them in a cell. 1 0 0 0 1 0 0 0 0 1 0 1 1 1 0 0 0 1 0 1 1 1 0 0 (a) (b) (c) Figure 3-18. Difficulties due to degeneracy of reconstructed interface data and implemented remedy. (a) A cell containing eight markers and one of several ways to define reconstructed polygons ; (b) a tetrahedron can have only one polygon; (c) some edges of tetrahedr ons created in a Cartesian cell. The newly created triangles are oriented us ing the gradient of the indicator function ( Figure 3-19 (a)) to point the normal vector in the di rection of increasing indicator value. Triangles with zero or arbitrarily small area can not be reliab ly oriented and may compromise the integrity of th e interface data. The orientati on of all such triangles is recursively corrected by using orient ation of the neighboring triangle ( Figure 3-19 (b)). The triangles on the entire in terface can be correctly or iented starting with proper

PAGE 65

47 orientation of any single triangle. Figure 3-20 shows examples of reconstruction creating a rupture and merger of two interfaces. I np r q tria tria=> p, q, r tria nghtria p q r s nghtria=> r, q, s (a) (b) Figure 3-19. Orientation setting of reconstructed inte rface. (a) Orientati on of triangles is set in the direction of incr easing indicator function va lue; (b) orientation of a triangle (tria) is used to set the or ientation of its neighbor (nghtria). A 3D View cross-sect ion 3D reconstructed (a) (b) Figure 3-20. Examples of 3D interface rec onstruction. (a) A surface undergoing rupture, (b) two spheres merging.

PAGE 66

48 3.4 Summary This chapter presented the immers ed boundary formulation and other computational components including the interface tracking and data-management. The key aspects of the presented me thod are summarized below: Immersed boundary method o Immersed boundary flow equati ons and interfacial conditions. o Fluid property smoothing usi ng an indicator function. o Surface force computation and modeling. Interface tracking o A conservative algorithm for interface restructuring. o A level contour based reconstruction algorithm to handle topology changes highlighting several difficulties and strategies employed to address them. The interface restructuring and reconstruction form the key contributions of the present work while using the well establ ished immersed boundary method to model the interfacial phenomenon except the phase change. The temperature equation is treated sharply and the solution algorithms are presented in the following chapter. The interface recons truction introduces some amount of error in the interfacevolume which is corrected explicitly by pe rturbing the markers in the local normal direction. The magnitudes of volume errors and required correction are estimated in chapter 5 along with the effect of reconstruction on the accu racy of flow computations via measurement of spurious velocities for static bubble computa tions. The effect of conservative restructuring is also demonstrated in chapter 5 via a spherical interface placed in a time-reversed vortex field. The marker tracking by itself is not restricted to any particular choice of interfacial dynamics modeling; it may be used for sh arp interface modeling such as ghost fluid method. Such a sharp interface modeling avoids smearing of fluid pr operties and reduces

PAGE 67

49 the spurious currents. So far, the only known work using marker based tracking with sharp interface ghost fluid method is by Hao a nd Prosperetti (2004). They used indicator function to compute the geometric propertie s of the interface for relatively simple computations not exhibiting topology changes. The indicator functi on was used to assign the fluid properties sharply as done for th e temperature equation described in the following chapter. For such a sharp interf ace method, several stability and accuracy related practical issues cau sed by sub-grid undulations interface reconstruction and errors in curvature computation need to be examined closely.

PAGE 68

CHAPTER 4 ADAPTIVE CARTESIAN GRID BASED COMPUTATION The dynamic grid adaptation is based on th e interface location and some feedback from the flow solution. The grid is stored using an unstructured data format and flow equations are solved using projection method. This chapter presents the grid generation/adaptation method along with the algo rithms for Navier-Stokes computation. 4.1 Data Structure Figure 4-1 shows a face-based data -structure used for storage of an isotropically refined grid. Since detailed information may be located in the work of Aftosmis (1997), Ham et al. (2002), Singh et al (2005), Singh and Shyy (2006), only a brief description of the data-structure is repr esented here. The variable level in Figure 4-1 stores the information about how many times a cell has been split from an initial state. The integers i, j and k together with the level information provide the size as well as the physical location of a cell. Also, all the faces form ing a cell are stored to provide cell-to-face connectivity information. For every face, a single byte integer named orientation gives the orientation of the face: faces normal to x-axis have orientation = 1, faces normal to yaxis have orientation = 2 and faces normal to z-axis have orientation = 3. The variable sideCell1 is identity of a cell sharing this face in the direction of orientation and sideCell2 is the cell on the othe r side of this face. The size (dx dy dz ) and center ( xc yc, zc) of a 3D cell wi th coordinates ( level, i, j, k ) are computed using Equation (4.1) Other relevant informati on such as cell volume and face area can be computed and stored in a ve ry compact way. Since the volume of a cell 50

PAGE 69

51 is depends on its level, it can be pre-computed and stored for up to the maximum allowed level of refinement. Similarly the area of a face is a function of only its orientation and the smallest of level of two side cells and hence it can be pre-computed and stored as well. 11 (,,),, 222 (,,)(0.5),(0.5),(0.5)y x z level level levelL L L dxdydz NxNyNz 1 x cyczcidxjdyjdz (4.1) Integer:: level Cell Data Integer:: i, j, k Integer:: cell-face list Integer:: orientation Face Data Integer:: sideCell1, sideCell2 orientation sideCell2 sideCell1 face orientation sideCell2 sideCell1 face (a) (b) Figure 4-1. Cartesian grid data-structure. (a) Cell and face data; (b) the convention for two side cells for a face with orient ation along x, y or z coordinates. The grid generation process starts by uniformly dividing the entire domain into prescribed number of cells in each coordinate direction. All the cells in this uniform grid are referred to as base cells upon which adaptation is carrie d out. These base cells are the largest cells possible in the com putational domain and have their level coordinate set to zero. To briefly explain the proc ess, consider a 2-D domain in Figure 4-2 (a) uniformly split into four parts in x-di rection and three parts in ydirection. Let the number of uniform cells in x and y directions be Nx and Ny respectively. In two-dimensions, each cell has three integer coordinates ( level, i, j). After the initial uniform splitting of the computational domain, all the cells are marked as ground level cells i.e., level = 0. The adaptation is carried out by split ting cells wherever desired. The final grid is obtained by

PAGE 70

52 recursively splitting the cells to desired level of resolution. When a cell with coordinates ( level, i, j) in Figure 4-2 (b) is isotropically split, coor dinates of the new cells shown in Figure 4-2 (c) are computed as: '1 (',')(21,21) ('','')(2,2) levellevel ijij ijij (4.2) (0,2,2) (0,3,1) Length=Lx, Nx=4 Length=Ly, Ny=3 (0,2,2) (0,3,1) Length=Lx, Nx=4 Length=Ly, Ny=3 (a) ( level, i, j ) ( level, i, j ) ( level, i, j ) ( level, i, j ) ( level, i, j ) ( level, i, j ) ( level, i, j ) ( level, i, j ) ( level, i, j ) ( level, i, j ) (b) (c) Figure 4-2. Adaptive grid refinement. (a) A 4x4 base grid ( level = 0), (b) coordinates of a two-dimensional cell before spli tting and (c) after splitting. 4.2 Grid Adaptation Interface location and shape provides the geometry based adaptation criteria. In addition to geometry, the flow-solution is also used as an indicator for adaptation. The refinement biased algorithm for isotropic adap tation is similar to the work of Ham and et al. (2002) i.e., when a cell flagged for refineme nt is to be broken; any cell blocking the

PAGE 71

53 refinement is also broken. This makes th e adaptation a refinement biased recursive process. The coarsening however is not a recu rsive process; a cell is coarsened only if it is permissible by the constraints due quality considerations and the data-structure. Coarsening is achieved by combining four cells in 2D or eight cells in 3D only if all these candidate cells are marked for coarsening. Fo llowing three restrictions for quality control were imposed on the adaptation process: 1. Up to a prescribed number of cell layers (~10) across the interface are always at maximum resolution to push errors du e to non-uniform cells away from the interface. Uniform grid around the interface also allows simple finite difference solution of the Poisson equation for indicator function in chapter 3. 2. Cells sharing a face are not allowed to diffe r by more than one level of refinement for smooth variation in grid size ( Figure 4-3 ) 3. Corner cells are not allowed to differ by more than one level of refinement ( Figure 4-3 ). This criterion is not a requirement when using collocated grid based computation (Singh & Shyy, 2006). For the implemented staggered grid algorithm, this is an essential requirement to facilitate a convenient programming implementation. Splitting not allowed Figure 4-3. Cartesian cells shar ing a face and the corner cells are not allowed to differ by more than one level of refinement. 4.2.1 Interface Geometry based Adaptation All the cells cut by the interface are flagge d for refinement. To facilitate a smooth transition in the cell size, up to five layers of cells around a flagged cell are also refined. This process of flagging and splitting cells is carried out recursively until the prescribed

PAGE 72

54 grid resolution is obtained. The Figure 4-4 shows the process of geometry based adaptation starting from a base grid and refining up to six levels. For demonstration purposes, only two additional layers of cells around an interfacecell were refined. (a) (b) (c) (d) Figure 4-4. Geometry based adaptation. (a) In itial 5x5 base grid and grid with maximum refinement level of (b) one, (c) three and (d) six. 4.2.2 Flow Solution based Adaptation Cells away from interface (outside ~10 cel ls across the interface) are adapted based on the flow solution. The present implementation uses a curl based adaptation criteria (Wang, 1998) that computes a parameter for each cell as shown in Equation (4.3) The length scale l is estimated as the cubic root of cel l-volume. The decision to refine or

PAGE 73

55 coarsen a cell is made by comparing cell to the standard deviation (Equation (4.3) ) using the criteria in Equation (4.4) 32 2 1,1 'cell i iNcellUl Ncell (4.3) Refine cell 0.1' Coarsen cellcell cell (4.4) During the adaptation procedure, the Cartes ian cell center values such as pressure, temperature and face normal velocities need to be reconstructed for the newly created cells and faces. The cell-centered values are reconstructed linearly while the face-normal velocities are reconstructed using the divergence free reconstruction due to Balsara (2001). Flow variable reconstruction during cell and face coarsening is performed simply by averaging of the corresponding cell -centered or face-centered values. 4.3 Immersed Boundary Solution Procedure The following seven steps present the ove rall structure of the computational algorithm: 1. Interpolate marker velocities from b ackground grid and due to phase change. Compute the mass transfer rate on the markers. m Advect markers using 2 nd Order Runge-Kutta for in terpolated velocity and forward Euler for the velocity due to mass transfer. *() () 12 12* *** ** 1 22nx xn XXtU XXtU n XX n X ntm Estimate interface volume change due to mass transfer and explicitly correct any volume error. 2. Check and restructure interf ace to maintain resolution.

PAGE 74

56 3. Compute indicator function and smoothed fluid propert ies. source term. ervals. e error caused by reconstruction. The he first four steps and the step 7 have b een described in the pr evious chapter. Of the remaining two steps, the geometry and solution based adaptation criteria have already been presented in the earlier sections. Th e following sections present the numerical discretization and Navier-S tokes solution procedure. The incompressible mass and momentum e quations are solved using a projection (Chorin, 1968; Kim & Moin, 1985) with stagge red grid finite volume formulation. The temperature equation is solved using colloca ted grid method. The pressure, temperature and fluid properties are stored at the cell center and face-nor mal velocity is stored on Cartesian cell faces ( Figure 4-5 (a)). The flow computation follows the following sequence of steps: Step 1: Solve for temperature field quation is solv ed using sharp interface method. The 4. Compute surface tension forces and di stribute as momentum 5. Solve Navier-Stokes equation to advan ce flow field to next time level. 6. Check for geometry and solution based grid adaptation at prescribed int 7. Check for interface reconstructi on at prescribed intervals: Perform reconstruction when required. Explicitly correct any global phase vol um magnitude of volume correction in a reconstructed interface is inversely proportional to its volume so that lesser co rrection is applied to smaller entities. If multiple interfaces are present, the reconstruction is selectively performed only for the interfaces that request it ba sed on the probe based criteria presented in the previous section. T The following temperature e spatial discretization of the terms is presented later. The fluid properties and diffusion term are treated sharply. The convection and diffusion terms are discretized using 2 nd order Runge-Kutta and fully implicit method, respectively. 1 11 / 2 nn nnTT Cu T t 1 nK T tep 2: Velocity computa tion using projection method. S ntum for an intermediate velocity field U using Equation (4.5) where all the know values such as surface tension source, gravitation, convection S Predictor-step olve the mome

PAGE 75

57 and old time-step viscous term due to Crank-Nicholson method are lumped into S n The temporal discretization of the convection term uses 2 nd order RungeKutta method. The pressure term is a pproximated using the old time pressure field. Subsequently, remove the effect of pressure term by shifting the velocity field back to obtain another intermediate velocity field U ** using Equation (4.6) ** nnV aUPdAfS vv i s c dAt (4.5) *** 1 n nP UUt Corrector-step Correct the predicte d velocity field (U**) using Equation (4.7) The pressure field for this correction is computed by enfo rcing velocity-divergence condition and (4.6) solving the Poisson Equation (4.8) with conjugate gradient method. The divergence of the new velocity field U n+1 is zero if no phase change occurs else it is computed using Equation (4.9) to account for the interfacial mass transfer described in the previous chapter. 1 1** n nP UUt 1 n (4.7) 1 ** 11ndV (4.8) 1 __ n n Cellface Cellface cellP ndA UndAU t 1211 Um I (4.9) 4.3.1 Staggered Grid Arrangement and Approximations4, 2005) is used to simplify the sp A velocity constraint suggest ed by Losasso et al. (200 atial discretization. With this approach, the velocities a nd pressure gradients on fine faces, face1 and face2 of a coarse cell shown in Figure 4-5 (b) are constrained as U 1 = U 2 and P x1 = P x2 i.e., face1 and face2 have the same control volume. The face normal gradients on these fine faces are approximated using Equation (4.10) so that same amount of velocity correction is applied when using Equation (4.7) It was also shown by Losasso et al. (2004, 2005) that applying the same pr essure gradients on finer faces can produce

PAGE 76

58 second order spatial accuracy fo r the computed pressure fi eld. Although assigning the same velocity to fine faces reduces the sp atial accuracy near non-uniform cells, it is accepted as a compromise for algorithmic si mplicity especially when non-uniform cells can be pushed away from the critical regions (e.g. interface) via suitable refinement criteria. 12 122 2xxPPP PP d 3 (4.10) U2 U1P1P2P3Px1Px2 face1 face2 d P,, U (a) (b) Figure ed arr flow/fluid properties; (b) single control volume used for fine faces of a coarse cell. 4.3.2 Spa 4-5. Staggergrida ngement. (a) Spatial location of tial Discretization of Velocity Equation The convection term of type U is represen ted in finite volume form integrated over tUdA (4.11) The normal velocities are computed appropriately to maintain the divergence free c tro he control-volume faces as: UdVUndA cf CV cfU ondition within each conl volume i.e., 0cfUdA Using Figure 4-6 (a) and the velocity constraints, the control volume f ace-noties are computed using facermal veloci area weighted averaging that produces:

PAGE 77

59 112 2130.5() 0.5()cf cfUUU UUU 314 56 40.5() (12) 12cf cfUUU UAreaUArea U A reaArea 4.12) The momentum flux term ( on control volume faces are computed using linear interpolation of facenormal velocities. To conserve the momentum flux across nonunifo rm grid cells, Figure 4-6 (b) makes another approximation that computes the flux on top control-volume face for velocity U 1 and equally distributes to the corresponding control volumes of faces for U 2 and U 3 at a T-node junction. The flux term in Equation (4.11) is computed using 2 nd or 3 rd order essentially non-oscillatory scheme (ENO) where the grid is locally of uniform size. Since fe w layers of cells around interface are always at uniform resolution, 3 rd order ENO scheme is employed there. For non-uniform cells (away from the interface) where uniform grid cells were not present to construct even a 2 nd order ENO scheme, a simple central differencing was used for simplicity. U1 U1Ucf1 Area1 U2 Ucf2 U3 U4 U6U5 U6 Area2 Ucf3 Ucf4 U1 U1Ucf1 Area1 U2 Ucf2 U3 U4 U 6U5 U Ucf4 f f/2 f/2 6 Area2 Ucf3 U2U3 U1 (a) (b) Figure 4-6. Spatial discretiza tion of convection term. (a) C ontrol volume (CV) and facenormal velocities (Ucf); (b) flux at a T-node junction face (top of CV for U is split eq3. f 1 ) ually (f/2) and passed to the control volumes of U 2 and U

PAGE 78

60 The momentum equation, with the known source terms lumped into (S u ,S v ,S w ), can b using Equation (4.13) The last term in Equation (4.13) com e rewritten es from the strategy outlined by Figure 4-6 while the terms due to are added only around the interface where the visc osity is spatially easily discretized Tu for incompressible flows. The discretization of u requires computation of control-volume face-normal velocity gradients a nd the terms due to require computation of mixed derivatives. Discretization of follows Tthere since the grid around the interface is uniform. Owing to divergence free nature of the velocity and constant viscos ity in most of the domain, th e contribution due to terms of Tu are neglected in cells away from th e interface where the grid may be nonuniform. Tu u u varying and it is ,, uvw P Dt ,, ,,u v wDu Su x x uvw DvP Sv Dty y uvw DwP Sw Dtz z (4.13) 4.3.3 Spatial Discretization of Temperature Equation Since the temperature is located at the cel l-center, a collocated grid discretization ed (Ham et al., 2002; Singh & Shyy, with auxiliary variable creation for non-uniform cells is us 2006). Convection terms are discretized using 3 rd order ENO (essentially nonoscillatory) scheme around the interface. A 2 nd order ENO scheme is used for cells away from the interface if the requi red uniform grid stencil can be easily be found, else a

PAGE 79

61 simple central differencing was used and was not found to produce any instabilities in the computed results. 4.3.3.1 Sharp int erface treatment The fluid properties (C K ) a re used sharply based on the interface location. The interf ace location, for this purpose, is appr oximated by indicator = 0.5 contour as was done for interface reconstruction described in the previous chapter. The term KT is computed using central differencing for ce lls away from the inte rface. For interfacial cells, the discretization is base d on the sharp interface treatme nt of Gibou et al. (2002) and Luo et al. (2005) using level-set method. The discretization in the interfacial cells can be explained using Figure 4-7 where the diffusion term for cell i is computed using the interface temperature as a boundary condition (Equation (4.14) ). 1 sati iiTTTT KK 1/2 1/2 iiT dx K xx dx (4.14) delta dx-delta i i+1 east west interface Figure 4-7. A one-dimen sional interface between i and i+1. tion presented in Chapter 3, the role of the sourc 4.3.3.2 Interface temperature treatment Coming back to the temperat ure equa e term Q s in Equation (3.3) is essentially to maintain the interface at saturation temperature. The approach of Juric and Tr yggvason (1998) treated this source term implicitly using a Newton iteration technique. The work of Son et al. (1999) and Son

PAGE 80

62 (2001) eliminated this source term by fixing the temperature inside the interface at saturation temperature and solv ing only outside the interfac e. Later on, Morgan (2005) suggested another approach by linearly correc ting the temperature wi thin one cell on each side of interface. This linear corre ction approach is summarized as: Step 1. Solve temperature equation: uTKT tC 1 (4.15) Step 2. T orrect temperature in cells immediately across the interface. Consider the erface lying between cell i and i+1 in the following figure. After solving the C int temperature equation in step 1, the temperature at cell centers i and i+1 are modified by linear interpolat ion of interface temperature (saturation temperature Tsat) and neighboring cell-center temperatures as shown in Figure 4-8 Corrected Ti x x x x Interface location Interface temperature = Tsat xTsat Ti T i-1 T i+1 T i+2 xCorrected Ti+1 x Figure 4-8. Linear interface temperature corr ection within one cell on each side of the interface. The proach defines the interf ace at indicator = 0.5 and defines the correc present ap tion as a solution of 0 KT This approach does not require explicit knowledge of the minimum distance of a cell from the interface as required in the method of Morgan (2005) using level-set tracking. Th is correction technique is convenient with the immersed boundary method and can easil y correct the temperature within any prescribed thickness across the interface. This equation is referred to as the interface

PAGE 81

63 temperature correction equation and it is solved sharply within a thickness of 1.2 gridcells on each side of the interface. Using Figure 4-9 this method will correct the temperature of T i and T i+1 while leaving T i-1 and T i+2 untouched. The interface saturation temperature (T sat ) and temperatures of cells outside the correction zone are taken as Dirichlet boundary condition for the correction e quation. Since K is constant within each domain, this equation reduces to 0 T and it is enforced while solving the temperature equation. x x x x Ti+2as BC BC Ti-1as BCKT 0 KT TiTi+1 Interface T as0 t o be corrected Figure 4-9. Presently employed interface te re correction technique solves a correction equation sh each side of the interface. This chapter presented the adaptive Cartesian grid computation technique emplo on are possible. For example, a solution based adaptation requires some measure of local solution error and a more general mperatu arply on 4.4 Summary yed in the current work. The discretizat ion of various terms in the flow equation and corresponding constraints on the grid were presented. The hydrodynamics was solved using immersed boundary method while the temperature equation was treated sharply for accurate interfacial mass-tr ansfer computation. Several improvements in the grid adaptati

PAGE 82

64 criter ion suitable for both low and high Re ynolds number flows based on the truncation or discretization error (Ham et al., 2002; Ferm & Lotstedt 2003; Muzaferija & Gosman, 1997) can be considered. Flow solution proce dure and numerical discretization also need further research for in order to improve th e spatial accuracy by elim inating restrictions placed on the velocity and pressure field.

PAGE 83

CHAPTER 5 RESULTS AND DISCUSSION The tests are conducted in thre e categories: basic tests to establish the working of adaptive grid flow solver and characteriz ation of the interface handling techniques; bubble/drop dynamics tests; phase change computations. The conservative interface restructuring technique and reconstruction are tested for various aspects such as accuracy, mass conservation. A spherical surface in a time reversed vortex field is used to test the effectiveness of conservative restructuring. Interface reconstruction is tested for its accuracy in terms of mass-conservation and so lution disturbances it causes. The spurious velocity currents and effect of reconstruc tion along with a set of multidimensional single and multiple rising bubbles and binary drop collision are computed to establish the accuracy of computational procedure. Computa tions of diffusion heat transfer between concentric cylinders, natural convection, one dimensional phase chan ge and a stationary bubble growth are conducted to gauge the perf ormance of the presently developed phase change computation procedure. 5.1 Interface Restructuring: Interface in a Vortex Field A spherical interface pl aced in a time-reversed vortex field ( Figure 5-1 ) deforms and should return back to the original spherical shape (Aulisa et al., 2004). The vortex field with time period T = 4 was imposed on the Cartesian grid using Equation (5.1) and 65

PAGE 84

66 the markers were advected by interpolating this velocity field from the Cartesian grid. 2 2 2(,,)cos()sin()sin(2)sin(2) (,,)cos()sin()sin(2)sin(2) (,,)cos()sin()sin(2)sin(2) uxyztTxzy vxyztTyxz wxyztTzyx (5.1) During the simulation period, markers are continuously added and deleted from the interface to maintain the reso lution. The interface at the en d of one time period should return back to the spherical shape with some errors due to restructuring and marker advection using the inte rpolated velocity field from the Ca rtesian grid. The errors at the end of one time period were estimated for the final surface area, volume and radius of markers (distance from the expected sphere center) as shown in Equation (5.2) Table 5-1 presents the computed errors showing a better performance by the conservative restructuring as compared to the non-conservative restructuring presented in Chapter 3. Figure 5-1 shows the time history of the volume-error in the interface al ong with the grid convergence studies. The conser vative method performs cons iderably well showing much smaller volume errors with better than qua dratic convergence. No te that the volume errors with conservative restructuring are cau sed solely by the integration of the marker advection equation using the in terpolated velocity field. 2 ,, 1,RMS_Radius_Error Area_Error1001 Volume_Error1001ifinaliinitial iN final initial final initalrr A A V V N (5.2)

PAGE 85

67 Table 5-1.Error estimates using Equation (5.2) for interface in a time reversed vortex field test comparing conservative and non-conservative restructuring. Grid RMS_Radius_error (conservative, nonconservative) % Area_error (conservative, nonconservative) % Volume_error (conservative, nonconservative) 30 x 30 x 30 (0.2435, 1.0254) x 10 -2 (4.2898, -7.9391) (-0.5591, -14.0950) 60 x 60 x 60 (0.0091, 0.4212) x 10 -2 (1.5306, -2.4060) (-0.1178, -4.9882) 120 x 120 x 120 (0.0079, 0.1155) x 10 -2 (0.4169, -0.5978) (-0.0216, -1.0455) t = 0 t = T/4 t = T/2 t = 3T/4 t = T t = 0 t = T/4 t = T/2 t = 3T/4 t = T (a) (b) Figure 5-1. Interface in a time reversed vorte x field. (a) Initial position of the interface with diameter 0.15 placed at (0.5, 0.75, 0.5) in a unit cube with its center at (0.5, 0.5, 0.5); computed shape history for vortex field in Equation (5.1) with period T = 4. time phasevolume 0 1 2 3 4 0.0125 0.0130 0.0135 0.0140Conservative Non-conservative Gridresolution %Volumeerror 40 60 80 100 120 10-1 100 101Conservative Non-conservative (a) (b) Figure 5-2. Error estimate and grid converg ence for interface in time reversed vortex field. (a) Interface volume error on 30 x 30 x 30 grid; (b) final volume error on 30 x 30 x 30, 60 x 60 x 60 and 120 x 120 x 120 grids with dotted quadratic reference rate.

PAGE 86

68 5.2 Interface Reconstruction : Effect on Mass Conservation The interface marker spacing and accuracy of the indicator-function computation and hence the level-contour based reconstruction depe nds on the background grid resolution. The aim is to measure interface volume-error caused by reconstruction i.e., the difference between the original and recons tructed interface volumes. The tests were performed using a spherical interface of radius r = 0.2 at the center of a domain with dimensions (0.4, 0.4, 0.4). The interface wa s reconstructed with different background grid resolution to conf irm a quadratic convergence of the volume error ( Figure 5-3 ) observed by Shin and Juric (2002). The volume errors produced by reconstruction are corrected explicitly by perturbing the markers in local normal direction. It is important to assess the approximate amount of perturbation a reconstruction may produce. Denoting the initial surface area and volume of the interface as A and V an approximate radial perturbation to correct a volume defect V may be estimated using Equation (5.3) This shows a second order convergence behavior for the radi al perturbation. To get furt her insight, this perturbation should be measured in terms of the grid resolution. If the diameter of the interface is resolved with grid cells, the radial pe rturbation can be rewritten in terms of the grid cell size as shown in Equation N (5.4) The last column of Table 5-2 shows the behavior of this quantity and reveals that this perturbation is approximately one tenth of the grid cell even for the coarsest grid used. This e rror too exhibits a quadratic convergence with grid refinement. Considering the fact that fl uid properties are smeared over nearly four to five cells across interface making the imme rsed boundary method a first order accurate method, the amount of radial disturbance cau sed by even moderately resolved grids ( > 20) is reasonably small. N

PAGE 87

69 243 VVV r r ArV (5.3) 36VV r r VV N (5.4) Numberofcellsperdiameter %Volumeerror 20 40 60 80 10 0 10-1 100 101 1 2 Figure 5-3. Grid convergence test for volume error due to reconstruction. Table 5-2. Volume error for reconstructi on of a spherical interface and approximate radial perturbation required for correction. Grid cell per-diameter (d/ ) % Volume error Radial perturbation (Equation (5.4) ) 10 -7.5995 0.1266 20 -1.9887 0.0663 40 -0.4881 0.0325 80 -0.1233 0.0164 5.3 Single Phase Navier-Stokes Solver A decaying vortex field and li d driven cavity tests are presented to establish the accuracy of the hydrodynamic solver. The temperature equation is tested for the discretization of the diffusion term using the diffusive heat transfer between concentric cylinders maintained at constant temperatures. The specified temperature on the cylinder geometries mimics the specified interface temperature and measures the error in the diffusion term discretization.

PAGE 88

70 5.3.1 Decaying Vortex Field A vortex field (Kim & Choi, 2000) shown in Equation (5.5) was introduced in a square domain of 1x1 with free-slip condition on all the walls. The center of the computational domain is located at the origin of the coordinate system. The computed decay of the velocity field at Reynolds number Re = 10 with (reference velocity and length scales taken as 1.0) was compared with the theoretical value at time t = 0.3. Figure 5-4 (a) shows the streamlines. 2 22/Re 2/Re(,,)cos()sin() (,,)sin()cos()t tuxytxye vxytxye (5.5) The simulations were carried out on a se ries of adaptive Cartesian grids and compared with uniform-grid computati on with resolutions corresponding to the maximum allowed resolution on the adaptive grid The adaptive grids were generated by recursively refining a 16x16 base grid and the computations were started with prescribed maximum allowed refinement in the entire domain. The solution based adaptation was performed at every 20 time step with the co mputational time step set at 0.001. An adapted grid with two refinement levels at time t = 0.3 is shown in Figure 5-4 (b). The maximum error in numerically comput ed u-component of the velocity field and the corresponding computation-time is presented in Table 5-3 with the time normalized by the time taken for 16x16 uniform-grid computation. As seen in Figure 5-5 uniform grid computations show quadrati c convergence but adaptive grid shows an average convergence rate of ~1.5. Since it was observed that the maximum error occurs along in the boundary cells which have non-uniform neighbor cells ( Figure 5-4 ), the computed error is a reasonable reflection of the accuracy of spatial discretization over

PAGE 89

71 non-uniform Cartesian grids. Although the accura cy from uniform grid has been observed to be superior, it is the computational cost that gives adaptive grid computations an edge. Since most of the cells shown in adapted gr ids are of uniform sizes, these computations may not offer the most attractive picture of co mputation-time savings that are expected to be more pronounced for 3D computations. (a) (b) Figure 5-4. Decaying vortex fi eld. (a) Streamlines; (b) two refinement levels on a 16x16 base grid. S q ua r e r ootofnumbe r ofcells MaxUerror 100 200 300 10-4 10-3 10-2Uniform Adaptive 2ndorderreference Figure 5-5. Grid convergence test for decayi ng vortex field showing maximum velocity error computed on uniform and adaptive grids.

PAGE 90

72 Table 5-3. Max velocity error and computation time norma lized with time taken for 16x16 uniform grid. The adapted grids are shown in Figure 5-4 Uniform grid computation Adaptive grid base computation Max. grid resolution Max. U error Normalized time cells Max. U error Normalized time 16 x 16 2.966 x 10 -3 1.000 256 2.966x 10 -3 1.000 32 x 32 7.184 x 10 -4 4.440 964 1.105x 10 -3 1.946 64 x 64 1.752 x 10 -4 21.602 3616 4.111x 10 -4 10.725 128 x 128 4.344 x 10 -5 104.810 14548 1.611x 10 -4 90.64 5.3.2 Lid-driven Cavity A lid driven cavity flow with Re = 100 was computed where the effect of grid adaptation is more visible than the decay ing vortex field in previous section. Figure 5-6 shows the streamlines of the flow computed on a 120 x 120 grid. Since the adaptation criterion is based on curl of the velocity fi eld, computed vorticity contours are presented in Figure 5-7 (a) and they are in accordance with the co mputations of Ghia et al. (1982). Figure 5-6. Streamlines for lid -driven cavity simulation for Re = 100 The computations were carried out on th ree grids with maximum resolution of 30 x 30, 60 x 60 and 120 x 120 obtained by one, two and three re finement levels as shown in Figure 5-7 Figure 5-8 shows the u-velocity profile along a vertical line at the center of the domain for the uniform and adaptive gr id computations showing no significant difference among the computed results on adapti ve and uniform grids. The exact values of minimum u-velocity in Figure 5-8 are presented in Table 5-4 This table also provides

PAGE 91

73 a comparison of the total numbe r of grid cells and computational time for uniform and adaptive grids. The second-last column of th is table shows ratio of the total number of cells from adaptive and uniform grid comput ations. The last column is the ratio of computational time between adaptive a nd uniform grids indicating up to ~70% computational time savings on large grids. Table 5-4. Lid driven cavity computation using uniform and adaptive grids. Comparison of the number of cells, computation time and minimum u-velocity on the vertical line through center. Uniform grid Adaptive grid Max. grid resolution #Cells Time Min u #Cells Time Min u Cell ratio Time ratio 30 x 30 900 1.000 -0.203 447 0.587 -0.204 0.497 0.587 60 x 60 3600 2.933 -0.211 1263 1.157 -0.210 0.351 0.394 120 x 120 14400 10.910 -0.212 4407 3.264 -0.212 0.306 0.299 (a) (b) (c) (d) Figure 5-7. Lid driven cav ity. (a) Vorticity contour and final adaptive grids with maximum refinement level of (b) one, (c) two and (d) three.

PAGE 92

74 U Y -0.2 0 0.2 0.4 0.6 0.8 1 -0.4 -0.2 0 0.2 0.430x30 60x60 120x120 U Y -0.2 0 0.2 0.4 0.6 0.8 1 -0.4 -0.2 0 0.2 0.4Uniform120x120 Adaptive30x30 Adaptive60x60 Adaptive120x120 (a) (b) Figure 5-8. The u velocity profile along a vertical li ne through the center of lid-driven cavity using (a) uniform grids and (b) adaptive grids. 5.3.3 Diffusive Heat Transfer between Concentric Cylinders The diffusion Equation (5.6) is solved for steady state solution with Dirichlet boundary condition of constant temperature on the cylinder walls. The density, heat capacity and conductivity are set to 1.0 The inner and outer surfaces of the cylinder are maintained at temperature T i = 2.0 and T o = 1.0, respectively. With the inner and outer radii denoted by R i and R o the theoretical steady stat e temperature distribution T ( R ) in circular coordinate syst em is shown in Equation (5.7) The computations were performed on a 2D Cartesian grid. The geometry of concentric cylinders (concentric circles in 2D) was modeled using an indicator function as shown in Figure 5-9 The computations are performed where indicator >= 0.5 (space between the concentric circles). Figure 5-10 shows the computed normalized temperature distribution and grid converg ence studies showing convergence rate between linear and quadratic as expected base d on the discretization of the diffusion terms near the boundary. T CK t T (5.6)

PAGE 93

75 log logi oioiTRTRR TTRR i (5.7) -0.4 -0.2 0 0.2 0.4 0 0.2 0.4 0.6 0.8 1IndicatorRadius RiRo (a) (b) Figure 5-9. Indicator co ntour used for modeling the cylindr ical computational domain. (a) Computations are performed only within the domain with Indicator >= 0.5; (b) cross-section of (a) show ing outer and inner radii as R o = 0.3 and R i = 0.1 respectively. (a) (b) Figure 5-10. Computed temperatur e and error for heat transfer in a concentric cylinder. (a) temperature distribution (T-T i )/(T o -T i ); (b) grid convergence study with resolution 30 (30 x 30 grid), 60 (60 x 60 grid) and 120 (120 x 120 grid) showing L 2 and L infinity error norms with referen ce quadratic and linear lines.

PAGE 94

76 5.3.4 Natural Convection in a Square Cavity The temperature equation solver was test ed for a natural convection case with Prandtl number (Pr ) = 0.71 and Rayleigh number ( Ra ) = 1.0 x 10 5 using Boussinesq approximation. The computations were performed in a 2D unit square ( L = 1) with top left and right walls at hot and cold temperatures (difference = T), respectively. The top and bottom walls were assigned adiabatic boundary conditions. Results in Table 5-5 show good agreement with the computations by Davis and Jones (1983) for the minimum and maximum Nusselt numbers on the left wall ( 0 xNuTx with temperature and lengths non-dimensionalized with T and L ) and v max as the maximum vertical velocity on a horizontal line through th e center of the domain. Figure 5-11 shows the computed streamlines and temperature isotherms. (a) (b) Figure 5-11. Natural convection with Pr = 0.71 and Ra = 1.0 x 10 5 (a) Streamlines and (b) the temperature isotherms. Table 5-5. Natural convection computation. Grid resolution Measured parameter 60 x 60 120 x 120 Davis and Jones (1983) u max 68.29 68.52 68.59 Nu x=0 min. 0.721 0.726 0.729 Nu x=0 max. 7.952 7.781 7.717

PAGE 95

77 5.4 Spurious Velocity Currents Spurious velocities or parasitic currents ( Figure 5-12 (a)) are unphysical velocity fields created due to numerical errors causi ng imbalance of interfacial stresses. These currents are best seen for static bubble simulations where the theoretical velocity is identically zero everywhere and the pressure drop across interface balances the surface tension forces dictated by the Young-Laplace Equation (5.8) ( Figure 5-12 (b)). There are several reasons ranging from the accuracy of surface tension computation to errors in computation of pressure and gr id anisotropy effects. A further detailed insight into the origins and nature of these currents can be found in the work of Torres and Brackbill (2000) and Harvie et al. (2005). Since the magnitude of these currents for continuous interface methods are found to be directly proportional to the surface tension and inversely proportional to the viscosity (Laf aurie, 1994), they have the potential to destabilize a computation and restrict the accu racy and the range of computational flowparameters (fluid properties, bubble sizes) that can be successfully simulated. The expected magnitude of non-dimensional spurious velocities (Capillary number) based on the observations of Shin and Juric (200 2) and Sousa et al. (2004) is of the order of 10 -4 and it has been hypothesized by Lafaurie (1994) that this magnitude remains largely unaffected by the Laplace number. To offer a perspective, the volume of fluid computations of Lafaurie(1994) observed spurious current of the order of 10 -2 that limited their simulations to Laplace num ber (La) < 10 while the curre nt code was tested up to Laplace number 12000 and found to be stable. The properties of the fluid outside the interface and bubble diameter were taken as relevant reference scales. A spherical bubble of diameter d = 0.2 was placed in a cubic domain of size (3 d 3 d 3 d ) with slip-wall boundary condition. The surface tension was set to 0.1 producing

PAGE 96

78 a theoretical pressure jump of 2.0. The computations were performed on 25 x 25 x 25, 50x50x50 and a 100 x 100 x 100 grid. All the time axes in the plots for the spurious velocity computations are non-dimensionalized with the capillary time scale of 1d p (5.8) -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0 0.5 1 1.5 2 Pressure (a) (b) Figure 5-12. Static bubble computation. (a ) Spurious currents of magnitude Ca~10 -4 (b) computed pressure for th eoretical jump of 2.0. 5.4.1 Effect of Fluid Property Jump and Grid Resolution A Laplace number = 250 case was performed on 25 x 25 x 25 and 50 x 50 x 50 grids. The magnitude of spurious current reduced by an order of magnitude on finer grid as shown in Figure 5-13 (a) but was found to exhibit no concrete improvements on finer grids. Since the pressure drop across the interface is captured over a few cells, an error in computed pressure drop was measured with an error norm in Equation (5.10) used by Brackibill et al. (1992) where N in is the number of cells inside the interface. The pressure drop p i in Equation (5.10) is the difference in the pressu re computed at the cell center and average pressure outside the interface ( ) as shown in Equation outp (5.9) for computing the average pressure drop nump The observed convergence rate of ~0.3 by Brackbill et al. (1992) is suppor ted by the computed data in Table 5-6 and Figure 5-13 (b).

PAGE 97

79 The simulated Laplace numbers in the range of 250 to 12000 were found to have no appreciable effect on the Capi llary number as shown in Table 5-7 High density and viscosity ratios (1000 and 100 re spectively) were also tested and found not to cause very serious degradation as shown Table 5-8 and Table 5-9 1111in outNN numinout i j ij in out p pppp NN (5.9) 2 11 Pressure drop errorinN iexact i in exactpp Np (5.10) Time Capillarynumber 0 10 20 30 40 0 0.0005 0.001 0.001525x25x25 50x50x50 ln21 0.38 ln4 EE Slope (25, E1)0.12 0.160.2040 60 80100 (100, E2 ) Grid resolutionError in pressure drop ln21 0.38 ln4 EE Slope (25, E1)0.12 0.160.2040 60 80100 (100, E2 ) Grid resolutionError in pressure drop (a) (b) Figure 5-13. Grid convergence tests for static bubble. (a) Capillary number evolution, (b) error in pressure drop (Equation (5.10) ) using data in Table 5-6 Table 5-6. Pressure drop for an inviscid fluids with density ratio =1.0. Data was taken after three non-dimensional time step of 5.0x10 -3 Grid resolution Pressure drop error (Equation (5.10) ) 25x 25 x 25 0.197 50x 50 x 50 0.151 100x 100x 100 0.117 Table 5-7. Effect of Laplace number on spurious velocity. Laplace number ( 2Lad ) Capillary number ( maxCaU ) 250 5.32x 10 -4 12000 5.35x 10 -4

PAGE 98

80 Table 5-8. Effect of density ratio on pressu re drop and spurious current on 50x50x50 grid (viscosity ratio = 1.0, La = 250) Density ratio = 12 Pressure drop= num exactp p Capillary number maxU 1 0.90 5.32x 10 -4 10 0.89 5.27x 10 -4 100 0.87 5.15x 10 -4 1000 0.87 6.35x 10 -4 Table 5-9. Effect of viscosity ratio on pressure drop and spurious current on 50x50x50 grid (density ratio = 1000, La = 250) Viscosity ratio = 1 2 Pressure drop= num exactp p Capillary number maxU 10 0.87 6.52x 10 -4 100 0.86 6.61x 10 -4 5.4.2 Effect of Interface Reconstruction The current level-contour based interface r econstruction is expected to have the same characteristics as shown in the tests of Shin and Juric (2002). The effect of reconstruction on mass conservation showing quadr atic convergence with grid refinement has already been established. Here the impact of reconstruction on sp urious velocities is assessed. The tests were performed for Laplace number 2500 on a 50 x 50 x 50 grid. Temporary spikes in spurious velocities were observed with every reconstruction ( Figure 5-14 (a)) as also seen by Shin and Juric (2002). Th ese spikes tend die down in time unless the reconstruction is performed too frequently depending on the time scale of the problem. Figure 5-14 (a) shows that if reconstruction is pe rformed every 10 Capillary time scale, the temporary increase in spurious velocity te nds to die down to the profile exhibited if no reconstruction were performed. However, if the reconstruction is performed at every Capillary time scale, the spurious currents ar e an order of magnitude larger and continue to grow. Further observation can be made from Figure 5-14 (b) which shows the interface

PAGE 99

81 normal velocity at a selected marker point. It shows that the oscill atory behavior caused by the surface tension effects quickly die dow n if no reconstruction is performed during the simulation. With reconstruction every 10 Capillary time scale, the amplitude of velocity does not decay but remains stable with mean value of nearly zero. This is to be expected based on the Laplace number of 2500 representing viscous time scales 2500 times slower than the Capillary time scales. But the reconstruction every Capillary time scale clearly shows a completely different mean value that is severa l orders of magnitude higher. These observations furt her strengthen the belief th at reconstruction should be performed as sparingly as possible. timeCapillary numberNo reconstruction Every Capillary time scale Every 10 Capillary time scale No reconstruction Every Capillary time scaletimeV velocity at top pointEvery 10 Capillary time scale (a) (b) Figure 5-14. Effect of reconstruction frequency on spurious velocities. (a) Capillary number evolution, (b) normal velo city at a selected marker point. 5.5 Buoyancy Driven Single Rising Bubble A set of rising bubbles are presented next to establish the accuracy of the present multidimensional immersed boundary proced ure and implementation for handling complex interfacial flows.

PAGE 100

82 5.5.1 Grid Convergence Test A single bubble of diameter d is placed in a computational domain of size and orientation indicated in Figure 5-15 (a). The bubble rises in the positive Y-coordinate direction under the influence of gravity. Th e non-dimensional parame ters used in the simulations are Reynolds number ( Re ), Morton number ( M ), number ( Eo ), density and viscosity ratios. The propertie s of fluid outside the interface and bubble diameter are taken as the required scales. Eotvos 9d3d 3d d z y x 2 0 2 4 6 8 1 2 4 (a) (b) Figure 5-15. Computational set-up for a sing le rising bubble. (a) Computational domain, (b) initial adaptive grid with 3 levels of refinement over 7 x 21 x 7 base grid. Simulations of Eo = 1.0 and M = 1.0x10 -3 representing a high surface tension and low Reynolds number case were conducted. To test the grid conve rgence, simulations were conducted on two grids with thr ee levels of refinement over a 7 x 21 x 7 base grid and a 10 x 30 x 10 base grid. The 7 x 21 x 7 base grid produces maximu m resolution equivalent to

PAGE 101

83 nearly 19 cells per bubble-diameter and the 10 x 30 x 10 base grid produces maximum resolution of nearly 27 cells per diameter. The density and viscosity ratios are both set to 20.0. The aspect ratio is computed as shown in Figure 5-16 and the average bubble rise velocity is computed using Equation (5.11) (Esmaeeli & Tryggvason, 1998). With these parameters, the bubble is expected to remain almost spherical with an estimated rise Reynolds number of 1.6 using the shape di agram of Clift et al. (1978) (Esmaeeli & Tryggvason, 1998). The computed results shown in Figure 5-17 and Table 5-10 confirm a spherical terminal shape for the rising bubble. The aspect ration of the bubble is large than 0.99 for the entire duration of the comput ation. Similar computations by Sousa et al. (2004) using uniform Cart esian grids of size 16x32x16 and 32x64x32 for a domain of size 2dx4dx2d produced an average rise Reynolds number of 1.4. interface1 vc v x nds Vol (5.11) Table 5-10. Grid convergence test for Eo = 1.0, M = 1.0 x 10 -3 The density and viscosity ratios are set to 20.0. Maximum adaptive grid resolution Rise Reynolds number Aspect ratio 10 x 30 x 10 1.53 0.992 56 x 168 x 56 1.51 0.996 b a Aspect ratio = b/a Figure 5-16. Aspect ratio computation

PAGE 102

84 (a) (b) (c) Figure 5-17. Bubble rise with Eo = 1.0, M = 1.0 x 10 -3 density and viscosity ratio of 20. (a) streamlines, (b) aspect ratio, (c ) rise Reynolds number with time ( tgd ). 5.5.2 Rising Bubbles in Different Terminal Shape Regimes A large body of experimental data on shapes and rise velociti es was condensed by Clift et al. (1978) in the fo rm of a single diagram using number, Morton number and Reynolds number as the non-di mensional parameters. This diagram has been used for validation by several researchers (Haario et al., 2004; Annaland et al., 2005) and serves as a general guideline for assessing the shape and rise velocities. Eotvos

PAGE 103

85 A series of Morton number (M) and Eo number (Eo) cases from tvos Table 5-11 were computed and compared with the expect ed terminal bubble-shape and rise velocity based on the shape diagram of Clift et al. (1978) along with comparison with the Volume of fluid computations of Annaland et al. ( 2005). The side walls and bottom surface of the domain ( Figure 5-18 (a)) were assigned free-slip bounda ry condition and outlet condition (specified pressure and zero face normal veloc ity gradients) was used on the top surface. The initial computational grid was created in a domain of size = (5d, 10d, 5d) with three levels of refinement over 12x24x12 base grid ( Figure 5-18 (b)). The computed shapes shown in Figure 5-19 are in good agreement with the expected shape from Table 5-11 and the computed rise Reynolds numbers agree well with computations of Annaland et al. (2005). The skirted bubble surface and streamlines showing several vortices are presented in Figure 5-20 (b) and compared with the comput ation of Haario et al. (2004) for M =1.3 and Eo =125. 10d5d 5d d z y x (a) (b) Figure 5-18. Computational domain for cases in Table 5-11 (a) Computational domain, (b) A 12 x 24 x 12 base grid with refinement levels.

PAGE 104

86 (a) (c) (d) (b) wobbling spherical ellipsoid spherical-cap skirted Dimpled ellipsoidal-cap Eotvos number R eyno ld s num b er LOG M0(a) (b) (c) (d) (e) Figure 5-19. Computed shapes from Table 5-11 (a) Spherical, (b) ellipsoidal, (c) dimpled ellipsoidal and (d) skirted regimes; (e ) shape diagram of Clift et al. (1978) with the corresponding locations of computed shapes indicated. (a) (b) (c) Figure 5-20. Computed case for M = 0.971, Eo = 97.1. (a) Interface and its cross-section, (b) streamlines, (c) streamlines for M=1.3, Eo=125 by Haario et al. (2005).

PAGE 105

87 Table 5-11. Computation of rising bubbles with different terminal shape regimes. Computational parameter Rise Reynolds number M Eo Terminal Shape (Clift et al., 1978) Clift et al. (1978) Annaland et al. (2005) Computed 1.26x 10 -3 0.971 CASE (a): Spherical 1.7 1.6 1.6 1.0x 10 -1 9.71 CASE (b): Ellipsoidal 4.6 4.3 4.6 1.0x 10 -3 97.1 CASE (c): Dimpled ellipsoidal 1.5 1.7 1.7 9.71x 10 -1 97.1 CASE (d): Skirted 20 18 17.8 5.6 Coalescence of Two Rising Bubbles Drop collision simulations involving inte rface reconstruction to merge two bubbles have been conducted and compared with th e volume of fluid base d computations of Annaland et al. (2005). The computational parameters of M = 2 x 10 -4 Eo = 16, density and viscosity ratio = 100 are used for both the simulations. Bubble diameter and property of fluid outside the interface are used as relevant scale to define non-dimensional parameters. 5.6.1 Head-on Coalescence The computational domain and initial grid shown in Figure 5-21 has dimensions (4d, 8d, 4d) discretized by three levels of refinement over 10x20x10 base grid resulting in a maximum resolution of 80x160x80 equivalent to twenty grid cells per bubble diameter. The initial centers of the top and bottom bubbles are at (2d, 2.5d, 2d) and (2d, d, 2d) with the origin located at the bottom left corner of the domain. For comparison of the time history of the bubble shape with literature, the dimensional para meters used in the current simulation and that of Annaland et al. (2005), are presented in the following table. The trailing bubble at the bottom rises fast er in the wake of the leading bubble resulting in an eventual colla pse. The time history of the development is presented in Figure 5-22 It shows the leading bu bble acquiring a spherical cap shape (at t ~ 0.75s) which is a consistent with respect to the shape diagram for the chosen computational

PAGE 106

88 parameters. The trailing bubbl e has a significant influence of the leading bubble and acquires a different shape. The computed results are in good agreement with the simulations of Annaland et al. (2005) and follow the developments indicated by the experiments of Brereton a nd Korotney (1991) shown in Figure 5-25 (a). Figure 5-23 show the computational grid and streamlines at 0.125s containing ~125K grid cells. The interface reconstructi on was performed at time t = 0.116s ( Figure 5-24 ) producing nearly 1.4% interface-volume error that was explicitly corre cted. The computational time was 142 hours on an Intel Xeon 1.8 GHz processor with time step ( tgd ) = 1.0x10 -4 Table 5-12. Computational paramete r for rising bubble-coalescence tests. Bubble diameter 0.01m Computational domain (0.04, 0.08, 0.04)m Grid resolution 80 x 160 x 80 Bubble fluid property Density ( 2 ) = 10 kg/m 3 Viscosity( 2 ) = 0.001 kg/ms Outer fluid property Density ( 1 ) = 1000 kg/m 3 Viscosity( 1 ) = 0.1 kg/ms Surface tension ( ) 0.1N/m 8d4d 4d d 2 0 2 4 6 8 1 2 (a) (b) Figure 5-21. Computational set-up for headon coalescence test. (a ) Domain and (b) a cross-section of the initial grid generated by three levels of refinement over a 10 x 20x10 base grid.

PAGE 107

89 0.000 sec 0.025 sec 0.050 sec 0.075 sec 0.100 sec 0.125 sec (a) (b) Figure 5-22. Head-on coalescen ce. (a) Computed time histor y of the two rising bubbles and (b) computations of Annaland et al (2005). [Reprinted with permission] (a) (b) Figure 5-23. Grid and streamlines at t = 0. 125 sec for head-on coal escence. (a) Interface shape and a cross-section view of the gr id, (b) streamlines w ith respect to the topmost point on the bubble.

PAGE 108

90 0 1 5 0 1 0 0 5 0 0 0 5 0 1 0 1 5 0.6 0.65 0.7 0.75 0.8 0.85 -0.15 -0.1 -0.05 0 0.05 0.1 0.1 5 0.6 0.65 0.7 0.75 0.8 0.85 (a) (b) Figure 5-24. Instance of interface reconstr uction during head-on coalescence of two rising bubbles at t = 0.116 s. (a) Cr oss-section of two interfaces before reconstruction and (b) after reconstruction. (a) (b) Figure 5-25. Experimental photographs of bubble coalescence taken at 0.03s time intervals. (a) A head-on coalescence and (b) an off-axis coalescence. (Taken from Brereton and Korotney, 1991) 5.6.2 Off-axis Coalescence All the computational parameters used fo r off-axis collision are the same as for head-on collision except the initia l bubble centers that are located at coordinates (2.8d, d, 2d) and (2d, 2.5d, 2d). The computed time histories are presented in Figure 5-26 (a). The computational results are in good agreement with the volume of fluid based computations of Annaland et al. (2005) ( Figure 5-26 (a)) and exhibit the behavior similar to the experimental photographs by Brereton and Korotney (1991) shown in Figure 5-25 (b).

PAGE 109

91 0.000s 0.025s 0.050s 0.000s 0.025s 0.050s 0.075s 0.100s 0.125s 0.075s 0.100s 0.125s 0.150s 0.175s 0.200s 0.150s 0.175s 0.200s (a) (b) Figure 5-26. Off-axis coalescence. (a) Comput ed, (b) computations of Annaland et al. (2005). [Reprinted with permission] 5.7 Binary Drop Collision Bubble coalescence tests for two rising bubbl es demonstrated the potency of the developed algorithms for handling comple x interface shapes and topology changes. These cases however were relatively be nign in terms of the dynamics of collision/coalescence along with relatively lower density ratios of 100. Next we consider collision of equal sized tetradecane (C 14 H 30 ) drops in nitrogen and evaluate the

PAGE 110

92 performance by comparing with the experi ments of Qian and Law (1997) for the dynamics of collision and volumes of the sat ellite droplets from the experimental observations of Estr ade et al. (1999). A typical binary drop collision process r esults in either bouncing or coalescence that may be followed by separation that may result in on or more satellite droplets. Figure 5-27 shows various collision regimes depend ing on the Weber number and impact parameter that distinguishes between off-cen tre and head-on collisions. Some of the experimental works on drop co llisions can be found in Ashg riz and Poo (1990), Jiang et al. (1992), Qian and Law (1997), Estrade et al. (1999), Brenn et al. (2001) and Willis and Orme (2003). In comparison, full three-dimens ional numerical simulations are relatively sparse and some of the know n computations are by: Rieber and Frohn (1997) using volume of fluid; Pan and Suga (2005) and Tanguy and Berlemont ( 2005) using level set; Nobari and Tryggvason (1996) and Shin and Juric (2002) using immersed boundary method; Mashayek et al. (2003) using finite element method. Besides computations with modest density ratios (20~30) by Nobari and Tryggvason (1996) and Shin and Juric (2002) no other 3D immersed boundary com putation could be found in literature. The non-dimensional parameters used for de fining the collision of two equal size drops (diameter = d) are impact parameter B (= h/d from Figure 5-28 ), Weber number (We) and Reynolds number (Re). The density and viscosity ra tios are used to define the properties of the surrounding medium. The parameter h in Figure 5-28 is the distance between the velocity vectors of the two drops and U impact is the impact velocity. The impact parameter defines the eccentricity of collision with B = 0 representing a head-on

PAGE 111

93 collision. All the non-dimensional parameters are computed with the fluid properties, drop diameter and the impact ve locity as the relevant scales. Due to excessive computational time a nd grid resolution required for attempting bouncing regimes at low Weber numbers, only regimes 3, 4 and 5 in Figure 5-27 are attempted. The fluid properties shown in Table 5-13 are taken from Pan and Suga (2005). The chosen fluid properties produce density and viscosity ratios of 666.081 and 179.28, respectively. The parameters for the computed cases are shown in Table 5-14 (1) Coalescence (2) Bouncing (3) Coalescence (4) Near head-on separation (5) Off-center separationWe B 60 20 40 0.4 0.8 Figure 5-27. Binary drop coll ision regimes marking the loca tion of computed parameters with solid circles. [Adapted from Qian and Law (1997)] h Uimpact/2 Uimpact/2 d Figure 5-28. Impact parameter B is defined as B=h/d. Table 5-13. Properties of tetradecane and nitrogen. Medium Density (Kg/m 3 ) Viscosity (Kg/ms) Surface tension (Kg/s 2 ) Tetradecane (C 14 H 30 ) 758.0 2.128 x 10 -3 Nitrogen 1.138 1.187 x 10 -5 0.026 (in nitrogen) Note: Adapted from Pan and Suga (2005)

PAGE 112

94 Table 5-14. Binary drop collis ion computation with density ratio = 666.081 and viscosity ratio = 179.28. Test Case We Re B Expected outcome from Qian and Law (1997) 1 32.8, 210.8, 0.08 Coalescence 2 37.2, 228.0, 0.01 Near head-on separation without satellite 3 61.4, 296.5, 0.06 Near head-on separation with satellite 4 60.1, 302.8, 0.55 Off-axis collision with satellites. The cases 1, 2 and 3 in Table 5-14 were performed on a 7.5d x 7.5d x 7.5d domain with four levels of refinement over 11 x 11 x 11 base grid resulting in approximately 24 cells per diameter. Case 4 was performed on a 9d x 4.5d x 4.5d domain with four levels of refinement over 16 x 8 x 8 base grids resulting in approximately 29 cells per diameter. The collisions take place in x-y plane with different impact parameter and impact velocity producing the desired non-dimensional numbe rs. The initial horizontal separation between the drop-centers in x-y plane is appr oximated from the figures of Qian and Law (1997) and it is set to one diameter (surf aces touching) for Case 1 and 1.25 diameters for the other cases. All the figures showing liquid drop shapes, cross-section of the computational grid and velocity vectors are in x-y plane through the center of the initial drops. The initial velocity conditions are obt ained by imposing the impact velocity field in the cells which have indicator function le ss than 1.0. The imposed velocity field is projected in a divergence-free space and used as the initial velocity condition. All the computational domains are assigned outle t condition imposing ze ro normal velocity gradient and constant atmospheric pressure on the boundaries. Outcome of the collision is strongly dependent on the interaction of the kinetic energy of collision competing with the su rface energy due to tension forces. The eccentricity of collision is an important factor in the overall collision outcome due to introduction of rotation energy in the coalesced drop resulting in production of

PAGE 113

95 centrifugal forces responsible for a stretching type separation. At Weber numbers below a critical value of 34 estimated in the work of Qian a nd Law (1997) the collision is expected to result in a permanent coalescence as the surface te nsion forces have a dominating influence. Figure 5-29 has Weber number = 32.8 (Case 1) demonstrating permanent coalescence due to high surface tension forces pulling the long cylindrical surface (t = ~ 1.22 ms) back towards a spherica l shape consistent with the experimental observations. Figure 5-30 shows the computational grid at various time instances along with time history of the grid data size e xhibiting the relationship between the grid data size and instantaneous drop shapes due to solution and geometry based adaptation. Case 2 was computed with the Weber number of 37.2 with a near head-on collision (B = 0.01). The surface forces are still st rong enough to create a dimpled disc shape converting the entire kinetic energy into surf ace energy. High curvature at the rim of the disc creates a high pressure zone as compared to the almost flat disc-center (t = 0.33ms in Figure 5-31 and Figure 5-32 ) causing contraction that eventually moves towards producing a cylindrical shape. Howeve r, the axial velocities shown in Figure 5-32 at t = 1.26ms and 1.71ms are strong enough to overcome the surf ace energy and continue the outward motion of the blobs at cylinder ends resulting in eventual breakup into two equal size droplets. Figure 5-33 shows a snapshot of the grid at several time instances and the corresponding time evolution of the grid data size. Further increase in the Weber number (61. 4 in Case 3) produces stronger kinetic energy as compared to the surface forces and the resulting dynamics is shown in Figure 5-34 The disc creation following the initial coalescence is similar to the earlier case with We = 37.2. However, an impact parameter of 0.06 was seen to pr oduces a large rotation

PAGE 114

96 of nearly 45 degrees seen in Figure 5-34 and Figure 5-35 Most of the rotation takes place around the instance of disc formation that subsequently begins elongation into a cylindrical shape taking place approxima tely between t = 0.23ms to 0.87ms. No significant rotation is observed during the separation process leading to formation of a satellite. Figure 5-36 shows the computed velocity profile in x-y plane across the bubbles along with the pressure distribut ion to highlight the collisio n mechanism. The separation and satellite formation occurs via elongation of the cylinder creating bulbous shape owing to surface tension forces. The separation creates two prim ary droplets of equal size and a satellite with an equivalent diameter th at is 25.8% of the in itial diameter of the individual drops. The comput ed satellite drop size is in good agreement with the experimental observation of Es trade et al. (1999) reporting a va lue of 26% for collision of ethanol droplets that have very simila r fluid properties as tetradecane. Similar computations by Pan and Suga (2005) produced a visually larger satellite drop although no quantitative information was made available. The magnitude and direction of the ro tation computed for Case 3 seem to be visually different than the experimental observations of Qian and Law (1997). Also experimental figures seem to exhibit most of the rotation after formation of the cylindrical shape leading to separation. The re asons for this apparent difference is still unknown and requires further investigation re garding the orientation of the collision plane and photographs in the experimental figures. This is besides the apparent differences in the symmetry between co mputation and experiment seen in Figure 5-34 around t = 0.18ms and 0.23ms for the impact parameter of 0.06.

PAGE 115

97 Case 4 is an off-center collision with hi gh impact parameter of 0.55 shown in Figure 5-37 As shown in the computation as well as the experiment, most of the rotation takes place with the coalesced mass before breakup process due to stretching begins. The experiment shows creation of two satellite droplets at t 1.92ms that quickly merge into a single satellite. However, the computed r esults show a single satellite droplet. A possible cause may be the relatively coarse computational grid but grid refinement studies could not be conducted due to exceedin gly large computational time on a serial computing platform. Further research evaluating grid refinement and the sensitivity of the computational parameters on the collision dyna mics need to be conducted. However, the equivalent diameter of the computed satellite droplet is found to be 15.7% of the initial individual bubble diameter a nd it is in good agreement w ith Estrade et al. (1999) observing a value of 17%. Figure 5-38 shows the computational gr id and grid data size evolution in time and Figure 5-39 provides the velocity fiel d in x-y plane across the interface. 5.7 Phase Change Computation A one dimensional phase change test pr oblem used by Son and Dhir (1998), Welch and Wilson (2000) is performed to evaluate the phase change computation technique. The temperature equation was solved sharply with out considering any smearing of the fluid properties in the temperature equation. The cap ability of the technique to capture Stefan condition accurately is demonstrated. Another test case in 3D for a stationary bubble growth in a superheated fluid is considered to examine the as ymptotic growth rate in time as compared to the theoretically predicted radial growth that is proportional to the square root of time (Prosperetti & Plesset, 1978).

PAGE 116

98 5.7.1 1-D Phase Change The left wall shown in Figure 5-40 was maintained at a constant temperature T wall with the interface location denoted by () t separating the vapor and liquid phases on the left and right sides respectivel y. The right side of the domai n is kept open to allow the liquid to flow out. The liquid is kept at the saturation te mperature. The vapor velocity is zero and as evaporation takes place, liqui d is pushed out from the right side. The theoretical interface location and temperature distribut ion are given by Equation (5.12) where and C are diffusivity and heat capacity of the vapor phase and is the latent heat of evaporation. Parameter is obtained by solving the transcendental Equation (5.13) The computations are started with th e initial interface at x = 0.2 and two cases with conductivity rati os of 1 and 10 from Table 5-15 are computed. The parameter =1.05968701 is used based on the solution of transcendental equation for the chosen fluid properties. Computational domain is fr om x = 0 to x = 1 w ith 500 grid cells. The computations are performed from the time when ()20.2tt with the theoretical solution imposed as the initial condition. Table 5-15. Computational parameters for one-dimensional phase change computation. Parameter Case 1 Case 2 T wall 1.0 1.0 T sat 0.0 0.0 Liquid and vapor density ,lv (1.0, 0.2) (1.0, 0.2) Liquid and vapor heat capacity ,lvCC (1.0, 5.0) (1.0, 5.0) Liquid and vapor thermal conductivity ,lvKK (1.0, 1.0) (10.0, 1.0) ()2 (,) 2wallsat walltt TTx TxtT erf erf t (5.12)

PAGE 117

99 2expwallsatTT erfC (5.13) The computed and theoretical interface location and temperature distribution for Case 1 assigning the same conductivity for both the fluids are shown in Figure 5-41 Also, the density ratio for th is case was increased to 1000 a nd tests were conducted with smeared treatment of the density and heat capacity in the temperature equation. In this case both sharp and smeared treatments produced results within 0.1% of each other. The second case assigning different conductivity to the two fluids is used to assess the differences between the sharp and smeared treatment for the thermal conductivity in the temperature equation. For this case, the smoothed thermal conductivity in the temperature was computed using Equation (5.14) The computed results with smoothed and sharp treatment are presented in Figure 5-42 with the sharp treatment following theoretical profile much better than the smoothed treatment. 1111vlv I KKKK (5.14) 5.7.2 Stationary Bubble Growth : Growth rate of a statio nary vapor bubble in one degree superheated liquid was computed and compared with the theoretical predicted growth rate A high density ratio of 1000 and viscosity ratio of 22, resembling a water-vapor combination, is used. The thermal conductivity and heat capacity ratios were set to 27 and 2, respectively. The Weber number (We) is 2.5e-5, Reynolds number (Re) is 0.58 and Jacob number (Ja) is set to 3.0. Since the diameter of the bubbl e changes in time the length scale is defined arbitrarily as R 0 = 0.5. The reference velocity was defined as 0R where is the

PAGE 118

100 thermal diffusivity of the liquid. The computa tion grid was resolved with ~ 25 grid cells per diameter for the initial bubble radius of R 0. Figure 5-43 shows the computed growth rate ex hibiting the asymptotic growth rate proportional to t Since the temperature equation is solved sharply, the interface temperature is well maintained as shown in Figure 5-43 where a thermal boundary layer develops in the liquid phase while the vapor phase remains at saturation temperature. (a) (b) Figure 5-29. Binary drop collision Case 1 from Table 5-14 (a) Computed time history, (b) experimental photographs of Qian and Law (1997). [Reprinted with permission]

PAGE 119

101 t = 0.00 ms t = 1.22 ms t = 2.26 ms Figure 5-30. A cross section of the computational gr id and the time hist ory of grid datasize for Case 1 from Table 5-14

PAGE 120

102 (a) (b) Figure 5-31. Binary drop collision Case 2 from Table 5-14 (a) Computed time history, (b) experimental photographs of Qian and Law (1997). [Reprinted with permission] 0.25 ms 0.33 ms 0.77 ms 1.26 ms 1.70 ms Figure 5-32. A cross section of the computa tional domain showing the interface and the velocity vectors for Case 2 from Table 5-14

PAGE 121

103 0.0 ms 0.33 ms 1.26 ms 1.60 ms 1.81 ms Figure 5-33. A cross section of the computational gr id and the time hist ory of grid datasize for Case 2 from Table 5-14

PAGE 122

104 (a) (b) Figure 5-34. Binary drop collision Case 3 from Table 5-14 (a) Computed time history, (b) experimental photographs of Qian and Law (1997). [Reprinted with permission]

PAGE 123

105 0.0 ms 0.37 ms 1.95 ms 2.21 ms 2.68 ms Figure 5-35. A cross section of the computational gr id and the time hist ory of grid datasize for Case 3 from Table 5-14

PAGE 124

106 t = 0.37 ms t = 1.42 ms t = 1.65 ms t = 1.95 ms Figure 5-36. A cross section of the computa tional domain showing the interface and the velocity vectors on the left and pressur e profile on the right for Case 3 from Table 5-14

PAGE 125

107 (a) (b) Figure 5-37. Binary drop collision Case 4 from Table 5-14 (a) Computed time history and (b) the experimental photographs of Qian and Law (1997). [Reprinted with permission]

PAGE 126

108 0.0 ms 0.52 ms 0.88 ms 1.00 ms 1.45 ms 2.24 ms Figure 5-38. A cross section of the computational gr id and the time hist ory of grid datasize for Case 4 from Table 5-14

PAGE 127

109 0.0 ms 0.52 ms 0.88 ms 1.00 ms 1.45 ms 2.24 ms Figure 5-39. A cross section of the computa tional domain showing the interface and the velocity vectors for Case 4 from Table 5-14

PAGE 128

110 Figure 5-40. One dimensi onal phase change setup. Time Inter f acelocation 0.02 0.04 0.06 0.08 0.1 0.3 0.4 0.5 0.6Theory Computed x T emperature 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1Theory Computed (a) (b) Figure 5-41. Comparison of computed a nd theoretical solution for Case 1 in Table 5-15 (a) Interface location and (b) temp erature profile at time = 0.1.

PAGE 129

111 Time Inter f acelocation 0.02 0.04 0.06 0.08 0.1 0.2 0.3 0.4 0.5 0.6Theory SmoothedK x Temperature 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1Theory SmoothedK (a) Time Inter f acelocation 0.02 0.04 0.06 0.08 0.1 0.2 0.3 0.4 0.5 0.6Theory SharpK x Temperature 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1Theory SharpK (b) Figure 5-42. One-dimensional phase change computation for Case 2 in Table 5-15 showing interface location on the left a nd temperature distribution at t = 0.1 on the right. (a) Computation with sm oothed and (b) sharp treatment of thermal conductivity in th e temperature equation.

PAGE 130

112 tR/R010-510-410-310-210-11 2 3 X -4 -2 0 2 4 0 0.5 1 1.5T Fluidflag (a) (b) Figure 5-43. Stationary bubble growth rate. (a) Non-dimensional radius and time with dashed line showing the asympto tic growth rate proportion to t ; (b) nondimensional temperature distribution along a line through the center of bubble at t = 0.18 with the inside and outside of the interface marked by fluid flags of one and zero, respectively.

PAGE 131

CHAPTER 6 SUMMARY AND FUTURE WORK 6.1 Summary The primary objective of the present work was to develop an effective computational framework that can handle co mplex interfacial flow computations. The main components of the present work can conc eptually be placed into the following four sections: Triangulated surface grid tracking to represent the multiphase interface. Cartesian grid generation and dynamic ad aptation based on interface location and solution features. Development of a staggered grid incompressible Navier-Stokes pressure based solver for spatially va riable density, viscosity flows. Immersed boundary modeling of interfaci al dynamics to account for fluid/flow property changes across interface and surface tension effects. Modeling of mass transfer due to phase change across interface. Brief comments are provided below to summarize the individual aspects and present status of the above four components. A conservative interface restructuring was used to add and delete markers (interface restructuring) w ithout introducing any error in the interface volume. The volume errors while using conservative rest ructuring with immersed boundary method occur solely due the errors in the marker a dvection using the veloc ity field interpolated from the background Cartesian grids. Since any interface volume error is corrected artificially at the end of ever y time step, the purpose of usin g conservative restructuring is to reduce the magnitudes of such artificial corrections. The magnitude of correction 113

PAGE 132

114 required with conservative restructuring is much smaller and hence reduces solution disturbance buildup due to artif icial relocation of markers to correct the volume errors. The interface reconstructi on is performed as sparingly as possible to minimize errors caused by reconstruction. Since the interface reconstruc tion causes artificial smoothing of the interface, it t oo comes with some error in the interface volume. These errors are corrected explicitly and the re quired corrections were shown to exhibit quadratic grid convergence. Using connectivity free tracking of Shin an d Juric (2002) could have simplified the interface tracking but would require frequent global reconstruction of the interface to maintain marker spacing as well as to accommoda te topology changes. It was shown that a frequent reconstruction may produce undesi rable effects for fast moving time scales of the flow. In this regards, it is believed that the conserva tive restructuring along with level-contour based reconstruction provides an enhanced flexibility by maintaining the interface resolution at all times and performing recons truction only if topology changes are expected. The single phase flow solver was tested for a number of cases to examine the accuracy of the hydrodynamics and temperature equation solvers with the staggered grid method. Since the velocity fields on fine faces of a large cell are constrained to be equal for algorithmic simplicity, a second order accurate solution is obtainable only in the uniform grid zones. In practice, the grid is uniformly adapted in the zones of interest (around interface), such a practi ce cause little degradation in actual performance of the computation.

PAGE 133

115 The spurious velocity currents in immer sed boundary method were tests and shown to be of the order of 10 -4 The overall method was found to be exhibit a spatial accuracy of the order of ~1.3 typical to a conti nuous interface method where the fluid properties are smeared over few cells across the interface. The interface tracking and its ability to handle complex topology changes and large fl uid property jumps were demonstrated via a set of rising bubble coalescence and bina ry drop collision computations. Although the hydrodynamics was solved us ing immersed boundary method with smoothed treatment of the fluid properties, the temperatur e equation itself was solved utilizing a sharp interface treatment without property smearing. The idea was borrowed from the well known sharp interface leve l-set method called ghost fluid method and combined with the immersed boundary trackin g to solve the mass transfer phenomena as accurately as possible. The idea was primarily to let the diffusion term see the same temperature gradients as seen by the numeric ally computed mass-transfer term. As a consequence, a locally non-conservative discr etization of the diffusion term was used near the interface to incorpor ate the interface temperature as a boundary condition. Such a treatment along with a local in terface correction made it easier to maintain the interface at saturation temperature and thus avoided any need to assume constant temperature inside the interface as done in the work of Son (2001). 6.2 Future Work The work conducted during this period is pa rt of a collaborative work aimed at developing an efficient comput ational capability exploiting the advantages of parallel computing. In view of the work needed for three dimensional multiphase fluid flow computations, parallel computing is highl y desirable. However, with dynamically evolving interface geometry and speed, an d possible topological changes, it is

PAGE 134

116 challenging to attain desirabl e parallel efficiency. Several methods and tools developed during this work form a base for fu rther efforts in this direction. The current capabilities are confined to computations performed in a box-shaped computational domain. Immersed boundary method to model moving/stationary solid boundaries has been used by sev eral researchers for a wide range of computations ranging from turbulent flows (Yang & Elias, 2006), fluid-structur e interaction (Zhu & Peskin, 2002; Gilmanov et al., 2003) to multiph ase flow (Liu et al., 2005; Udaykumar et al., 2001; Ye et al., 1999). Among these, the me thod of Ye et al. (1999) and Udaykumar et al. (2001) uses cut-cell approach and is us ually more involved than the others. In this regard, the finite difference b ased method of Liu et al. (2005) is one of the simplest to use. It uses finite difference method near inte rface and computes the gradients of the flow using the interface conditions (e.g., no slip wall) as a boundary condition. The technique used to discretize the temperature equation in the presented work is similar in nature as it also modeled the interf ace as a boundary condition and used finite difference for interface cells and finite volume elsewhere. The computational procedure developed fo r phase change involv es several aspects that need further research for assessing the numerical accuracy and seeking necessary enhancement. An example of this is the artific ial correction of temperature in cells across the interface. Its effect with respect to grid resolution, thermal boundary layer thickness and convection dominated flows needs a cl oser scrutiny although such a correction technique was employed and shown to perf orm well for a number of nucleate boiling computations by Morgan (2005). The overall phase change computation with explicit marker tracking was shown to perform well for cases considered, further work on more

PAGE 135

117 complex cases such as rising bubbles need to be conducted to evaluate the numerical accuracy of this approach. In addition, for morphological level of phase changes, such as crystal growth problems, when the geometry is highly irre gular and the curvatures and surface tension/energy effects more difficult to accurately capture, significant efforts will be required before truly threedimensional physics of practical interest can be simulated.

PAGE 136

LIST OF REFERENCES Al-Rawahi N. and Tryggvason G., Numerical simulation of dendritic solidification with convection: Three-dimensi onal flow, Journal of Com putational Physics, 2004, 194:677-696. Al-Rawahi N. and Tryggvason G., Numerical simulation of dendritic solidification with convection: Two-dimensional geometry, J ournal of Computati onal Physics, 2002, 180:471-496. Aftosmis M. J., Solution adaptive Cartesian grid method for aerodynamic flows with complex geometries, Computational Fluid Dynamics VKI Lectures Series, Belgium, March 3-7, 1997. Agresar G., Linderman J. J., Tryggvason G. a nd Powell K. G., An adaptive, Cartesian, front-tracking method for the motion, de formation and adhesion of circulating cells, Journal of Computati onal Physics, 1998, 143:346-380. Annaland M. V. S., Deen N. G. and Kuiper s J. A. M., Numerical simulation of gas bubbles behavior using a thre e-dimensional volume of fluid method, Chemical Eng. Science, 2005, 60:2999-3011. Ashgriz N. and Poo J. Y., Coalescence and separa tion in binary collisions of liquid drops, Journal of Fluid Mechanics, 1990, 221:183-204. Aulisa E., Manservisi S., Scar dovelli R. and Zaleski S., A geometrical area-preserving volume of fluid method, Journal of Co mputational Physics, 2003, 192:355. Aulisa E., Manservisi S. and Scardovelli R., A surface marker algorithm coupled to an area-preserving marker redistribution method for three-dimensional interface tracking, Journal of Computati onal Physics, 2004, 197:455. Balsara D. S., Divergence-free adaptive mesh refinement for magneto-hydrodynamics, Journal of Computational Physics, 2001, 174:614-648. Berger M. J. and Oliger J., Adaptive mesh re finement for hyperbolic partial differential equations, Journal of Computati onal Physics, 1984, 53:482-512. Brackbill J.U., Kothe D.B. and Zemach C ., A continuum method for modeling surface tension, Journal of Computat ional Physics, 1992, 100:335-354. 118

PAGE 137

119 Brenn G., Valkovska D. and Danov K. D., The formation of satellite droplets by unstable binary drop collisions, Physics of Fluids, 2001, 13:2463-2477. Brereton G. and Korotney D., Coaxial and obl ique coalescence of two rising bubbles, in Dynamics of Bubbles and Vortices near a Free Surface, edited by I. Sahin and G. Tryggvason, Vol. 119, ASME, New York, 1991. Chorin A. J., Numerical solution of the Navier-Stokes equations, Mathematics of Computation, 1968, 22:745-762. Clift R., Grace J. R. and Weber M, Bubbles, Drops and Particles, Academic Press, New York, 1978. De Vahl Davis G. and Jones I. P., Natura l convection in a squa re cavity: a comparison exercise, International Journal for Nume rical Methods in Fl uids, 1983, 3:227-248. Dhir V. K., Numerical simulation of poolboiling heat transfer, AIChE Journal 2001, 47:813-293. Enright D., Fedkiw R. P., Ferzig er J. and Mitchell I., A hybr id particle level set method for improved interface capturing, Journal of Computational Physics, 2002, 183:83116. Esmaeeli A. and Tryggvason G., Direct numeric al simulations of bubbly flows, Part I. Low Reynolds number arrays, Journal of Fluid Mechanics, 1998, 377:313-345. Esmaeeli A. and Tryggvason G., Direct numeric al simulations of bubbly flows, Part II. Moderate Reynolds number arrays, Jour nal of Fluid Mechanics, 1999, 385:325358. Esmaeeli A. and Tryggvason G., Co mputations of film boiling. Part I: numerical method, International Journal of Heat and Mass Transfer, 2004, 47:5451-5461. Estrade J.-P., Carentz H., Lavergne G. a nd Biscos Y., Experiment al investigation of dynamic binary collision of ethanol dr oplets-a model for droplet coalescence and bouncing International, Journal of Heat and Fluid Fl ow, 1999, 20:486-491. Ferm L. and Lotstedt P., Anisotropic grid adaptation for Navier-Stokes equations, Journal of Computational Physics, 2003, 190:22-41. Francois M. and Shyy W., Mi cro-scale drop dynamics for heat transfer enhancement, Progress in Aerospace Sciences, 2002, 38:275-304. Francois M. and Shyy W., Computations of drop dynamics with the immersed boundary method; Part 1Numerical algorithm a nd buoyancy induced effect, Numerical Heat Transfer, Part B, 2003, 44:101-118.

PAGE 138

120 Francois M. and Shyy W., Computations of drop dynamics with the immersed boundary method; Part 2Drop impact and heat tran sfer, Numerical Heat Transfer, Part B, 2003, 44:119-143. Francois M., Uzgoren E., Jackson J. and Shyy W., Multigrid computations with the immersed boundary technique for multiphase flows, International Journal of Numerical Methods for Heat a nd Fluid Flow, 2004, 14:98-115. Ghia U., Ghia K. N. and Shin C. T., HighRe solution for incompr essible flow using the Navier-Stokes equations and a multigri d method, Journal of Computational Physics, 1982, 48:387-411. Gibou F., Fedkiw R., Caisch R. and Osher S, A level set approach for the numerical simulation of dendritic growth, Journal of Scientific Computing, 2003, 19:183-199. Gilmanov A. and Sotiropoulos F., A hybrid Cartesian/immersed boundary method for simulating flows with 3D, geometrically complex, moving bodies, Journal of Computational Phys ics, 2005, 207:457-492. Glimm J., Graham M. J., Grove J., Li X. L., Smith T. M., Tam D., Tangerman F. and Zhang Q., Front tracking in two and thr ee dimensions, Computers and Mathematics with Applications, 1998, 35:1:11. Glimm J., Grove J. W., Li X.L. and Tan D.C., Robust computational algorithms for dynamic interface tracking in three dimens ions, SIAM Journal of Scientific Computing, 2000, 21:2240-2256. Griffith B. E. and Peskin C. S., On the order of accuracy of the immersed boundary method: Higher order convergence rates for sufficiently smooth problems, Journal of Computational Physics, 2005, 208:75-105. Ham F. E., Lien F. S. and Strong A. B ., A Cartesian grid method with transient anisotropic adaptation, Journal of Co mputational Physics, 2002, 179:469-494. Hao Y. and Prosperetti A., A numerical me thod for three-dimensional gas-liquid flow computations, Journal of Comput ational Physics, 2004, 196:126-144. Haario H., Korotkaya Z., Luukka P. And Sm olianski A., Comput ational modeling of complex phenomena in bubble dynamics: vortex shedding and bubble swarms, Presented at European C ongress on Computational Met hods in Applied Sciences and Engineering ECCOMAS, 24-28, July 2004. Harvie D. J. E., Davidson M. R. And Rudm an M., An analysis of parasitic current generation in volume of fluid simulations, Journal of Australian and New Zealand Industrial and Applied Ma thematics, 2005, 46:133-149. Jiang Y. J., Umemura A. and Law C. K., An experimental investigation on the collision behavior of hydrocarbon dr oplets, Journal of Fluid Mechanics, 1992, 234:171-190.

PAGE 139

121 Juric D. and Tryggvason G., Co mputations of boiling flows, International Journal of Multiphase Flows, 1998, 24:387-410. Juric D. and Tryggvason G., A fr ont-tracking method for dendrit ic solidification, Journal of Computational Physics, 1996, 123:127-148. Kang M., Fedkiw R. and Liu X.-D., A boundary condition capturing method for multiphase incompressible flow, Journal of Scientific Computing, 2000, 15:323360. Kim D. and Choi H., A second-order time-accu rate finite volume method for unsteady incompressible flow on hybrid unstructured grids, Journal of Computational Physics, 2000, 162:411-428. Kim J. and Moin P., Application of a frac tional step method to incompressible NavierStokes equations, Journal of Com putational Physics, 1985, 59:308-323. Kleinstreuer C., Two-phase flow: theory and applications, Taylor a nd Francis, New York, 2003. Lafaurie B., Nardone C., Scardovelli R. Zal eski S. and Zanetti G., Modeling merging and fragmentation in multiphase flows with SURFER, Journal of Computational Physics, 1994, 113:134-147. LeVeque R. J. and Li Z., The immersed interface method for ellip tic equations with discontinuous coefficients and singular sources, SIAM Journal of Numerical Analysis, 1994, 31:1019-1044. Li Z. and Lai M-C, The immersed interface method for the Navier-Stokes equations with singular forces, Journal of Com putational Physics, 2001, 171:822-842. Liu X.-D., Fedkiw R. and Kang M., A boundary condition capturing method for Poissons equation on irregul ar domains, Journal of Co mputational Physics, 2000, 160:151-178. Liu H., Kishnan S., Marella S. and Udaykum ar H. S., Sharp interface Cartesian grid method II: A technique for simu lating droplet interactions with surfaces of arbitrary shape, Journal of Computational Physics, 2005, 210:32-54. Luo X.-Y., Ni M.-J., Ying A. and Abdou M. A., Numerical modeling for multiphase incompressible flow with phase change, Numerical Heat Transfer-Part B, 2005, 48:425-444. Losasso F., Gibou F., Fedkiw R., Simulating water and smoke with an octree data structure, ACM Transactions on Graphics, 2004, 23:457-462. Losasso F., Gibou F., Fedkiw R., Spatially adap tive techniques for level set methods and incompressible flow, Computers and Fluids, 2005, in press.

PAGE 140

122 Mashayek F., Ashgriz N., Minkowycz W. J. and Shotorban B., Coalescence collision of liquid drops, International Journal of Heat and Mass Transfer, 2003, 46:77-89. Meier M., Yadigaroglu G. and Smith B. L ., A novel technique for including surface tension in PLIC-VOF methods, European Journal of Mechanics B/Fluids, 2002, 21:61-73 Morgan N. R., A new liquid-vapor phase tran sition technique for the level set method, Ph.D. thesis, Georgia Institute of Technology, 2005. Muzaferija S. and Gosman D., Finite-volum e CFD procedure and adaptive error control strategy for grids of arbitr ary topology, Journal of Co mputational Physics, 1997, 138:766-787. Nobari M. R. H. and Tryggvason G., Numerical simulations of three-dimensional drop collisions, AIAA Journal, 1996, 34:750-755. Osher S. and Fedkiw R. P., Level set met hods: an overview and some recent results, Journal of Computational Physics, 2001, 169:463-502. Osher S. and Fedkiw R. P., Level set met hods and dynamics implic it surfaces, Applied Mathematical Science, Sp ringer-Verlag, NewYork, 2003. Pan Y. and Suga K., Numerical simulation of binary liquid droplet co llision, Physics of Fluids, 2005, 17:1-14. Perot B., Nallapati R., A movi ng unstructured staggered mesh method for the simulation of incompressible free-surface flows, J ournal of Computational Physics, 2003, 184:192-214. Peskin C. S., Numerical analysis of blood flow in the hear t, Journal of Computational Physics, 1977, 25:220-252. Peskin C. S. and McQueen D. M., Fluid dynamics of the heart and its valves, in Case studies in Mathematics Modeling: Ecology, Physiology, and Cell Biology, edited by Othmer H. G., Adler F. R., Lewis M. A., Dallon J. C., Prentice-Hall, Englewood Cliffs, NJ, 1996. Popinet S., Gerris: a tree-based adaptive solver for the incompressibl e Euler equations in complex geometries, Journal of Co mputational Physics, 2003, 572-600. Prosperetti A. and Plesset M.S., Vaporbubble growth in a superheated liquid, Journal of Fluid Mechanics, 1978, 85:349-368. Prosperetti A., Numerical algorithms for free-surface flow computations: An Overview, in Drop-Surface Interactions edited by Rein M., Springer, 2002.

PAGE 141

123 Qian J. and Law C. K., Regi mes of coalescence and separat ion in droplet collision, Journal of Fluid Mechanics, 1997, 331:59-80. Raymond F. and Rosant J.-M., A numerical a nd experimental study of the terminal velocity and shape of bubbles in viscous liquids, Chemical Engineering Science, 2000, 55:943-955. Renardy Y. and Renardy M., PROST: a parabolic reconstruction of surface tension for the volume-of-fluid method, Journal of Computational Physics, 2002, 183:400-421. Rieber M. and Frohn A., Navier-Stokes simu lation of droplet collision dynamics, in Proc. 7 th Int. Symp. On Comp. Fluid Dynamics, edited by Zhuang F., Beijing, China, 1997. Rider W. J. and Kothe D. B., Stretching and tearing interface tracking methods, AIAA Paper 1717, presented at 12 th AIAA CFD Conference, San Diego CA, 1995. Rider W. J. and Kothe D. B., Reconstructi ng volume tracking, Journal of Computational Physics, 1998, 141:112-152. Ryskin G. and Leal L. G., Numerical solution of free-boundary problems in fluid mechanics, Journal of Fluid Mechanics, 1984, 148:1-43. Scardovelli R. and Zaleski S ., Interface reconstruction with least-square fit and split Lagrangian-Eulerian advecti on. International Journal for Numerical Methods in Fluids, 2002, 41:251. Shin S. and Juric D., Modeling Three-dimensional multiphase flow using a level contour reconstruction method for front tracki ng without connectivity, Journal of Computational Phys ics, 2002, 180:427-470. Shin S., Abdel-Khalic S. I., Daru V. and Juric D., Accurate re presentation of surface tension using the level contour reconstr uction method, Journal of Computational Physics, 2005, 203:493-516. Shin S., Abdel-Khalik S. I. and Juric D., Direct three-di mensional numerical simulation of nucleate boiling using th e level contour reconstruc tion method, International Journal of Multiphase Flow, 2005, 31:1231-1242. Shyy W., Udaykumar H. S., Rao M. M. and Smith R. W., Comput ational fluid dynamics with moving boundaries, Taylor & Francis, Washington, DC, 1996. Shyy, W., Multiphase computations using sh arp and continuous interface techniques for micro-gravity applications, Compt es Rendus Mecanique, 2004, 332:375-386.

PAGE 142

124 Singh R. K., N'Dri N., Uzgoren E., Shyy W. a nd Garbey M., Three-dimensional adaptive, Cartesian grid method for multiphase flow computations, AIAA Paper-1390, presented at 43rd AIAA Aerospace Sciences Meeting and Exhibit, Reno, NV, 2005. Singh R.K. and Shyy W., "Three-Dimensional Adaptive Grid Computation with Conservative, Marker-Based Tracking fo r Interfacial Fluid Dynamics", AIAA Paper-1523, presented at 44 th Aerospace Science Meeti ng and Exhibit, Reno, NV, 2006. Son G., A numerical method for bubble moti on with phase change, Numerical Heat Transfer, Part B, 2001, 39:509-523. Son G., Ramanujapu N., Dhir V. K., Numerical simulation of bubble merger process on a single nucleation site during pool nucleate boiling, ASME Journal of Heat Transfer, 2002, 124:51-62. Sousa F. S., Mangiavacchi N., Nonato L. G., Castelo A., Tome M. F., Ferreira V. G., Cuminato J. A. and McKee S., A front-t racking/front-capturi ng method for the simulation of 3D multi-fluid flows with free surfaces, Journal of Computational Physics, 2004, 198:469-499. Sussman M., Almgren A. S., Bell J. B., Colella P., Howell L. and Welcome M., An adaptive level set approach for inco mpressible two-phase flow, Journal of Computational Phys ics, 1999, 148:81-124. Sussman M., A second order coupled leve l set and volume-of-fluid method for computing growth and collapse of va por bubbles, Journal of Computational Physics, 2003, 187:110-136. Sussman M., A parallelized, adaptive algorithm for multiphase flows in general geometries, Journal of Computer s and Structures, 2005, 83: 435-444. Tanguy S. and Berlemont A., Application of a level set method for simulation of droplet collisions, International Journal of Multiphase Flow, 2005, 31:1015-1035. Torres D. J. and Brackbill J. U., The point-set method: front-tracking without connectivity, Journal of Computational Physics, 2000, 165:620-644. Tryggvason G., Bunner B., Esm aeeli A.., Al-Rawahi N., Taube r W., Han J., Nas S. and Jan Y-J., A front tracking method for the computations of multiphase flow, Journal of Computational Physics, 2001, 169:708-759. Tryggvason G., Esmaeeli A. D. and Al-Rawahi N ., Direct numerical simulations of flows with phase change, Computers and Structures, 2005, 83:445.

PAGE 143

125 Udaykumar H. S., Mittal R. and Shyy W., Comp utation of solid-liquid phase fronts in the sharp interface limit on fixe d grids, Journal of Com putational Physics, 1999, 153:535-574. Wang Z.J., A quadtree-based adaptive Cartesian/ Quad grid flow solver for Navier-Stokes equations, Computers & Fluids, 1998, 27:529-549. Wasekar V. M. and Manglik R. M., Short-time transient surfactant dynamics and Marangoni convection around boiling nuclei, Journal of Heat Transfer, 2003, 5: 858-866. Welch S. W. J. and Wilson J., A volume of fluid based method for fluid flows with phase change, Journal of Computationa l Physics, 2002, 160: 662-682. Willis K. and Orme M., Viscous oil droplet collisions in a vacuum, Experiments in Fluids, 2000, 29:347-358. Yang J. and Balaras E., An embedded-boundary formulation for larg e-eddy simulation of turbulent flows interacting with moving boundaries, Journal of Computational Physics, 2006, 215:12-40. Ye T., Mittal R., Udaykumar H. S. and Shyy W., An accurate Cartesian grid method for viscous incompressible flows with co mplex immersed boundaries, Journal of Computational Phys ics, 1999, 156:209-240. Ye T., Shyy W. and Chung J. N., A fixe d-grid, sharp-interf ace method for bubble dynamics and phase change, Journal of Co mputational Physics, 2001, 174:781-815. Ye T., Shyy W., Tai C-F, and Chung J. N., Assessment of sharpand continuousinterface methods for drop in static e quilibrium, Computers & Fluids, 2004, 33:917-926. Youngs D. L., Time-dependent multi-material flows with large fluid distortion, in Numerical Methods for Fluid Dynamics, edited by Morton K. W., Baines M. J., Academic Press, New York, 1982. Zhu L. and Peskin C., Simulation of a flexible flapping filament in a flowing soap film by the immersed boundary method, Journal of Computational Physics, 2002, 179:452468.

PAGE 144

BIOGRAPHICAL SKETCH Rajkeshar Singh was born in 1977 in Uttar Pradesh, Indi a. He was brought up in Mumbai and received his Bachelor of T echnology degree in aerospa ce engineering from the Indian Institute of Technology, in 2000. Thereafter he joined Mississippi State University and graduated in 2002 with Mast er of Science degree in computational engineering. Since then he has been pur suing his Ph.D. degree in mechanical and aerospace engineering at the University of Fl orida. His current resea rch interests lie in computational modeling of two phase flows and adaptive grid based computing. 126


Permanent Link: http://ufdc.ufl.edu/UFE0015703/00001

Material Information

Title: Three-Dimensional Marker-Based Multiphase Flow Computation Using Adaptive Cartesian Grid Techniques
Physical Description: Mixed Material
Copyright Date: 2008

Record Information

Source Institution: University of Florida
Holding Location: University of Florida
Rights Management: All rights reserved by the source institution and holding location.
System ID: UFE0015703:00001

Permanent Link: http://ufdc.ufl.edu/UFE0015703/00001

Material Information

Title: Three-Dimensional Marker-Based Multiphase Flow Computation Using Adaptive Cartesian Grid Techniques
Physical Description: Mixed Material
Copyright Date: 2008

Record Information

Source Institution: University of Florida
Holding Location: University of Florida
Rights Management: All rights reserved by the source institution and holding location.
System ID: UFE0015703:00001


This item has the following downloads:


Full Text












THREE-DIMENSIONAL MARKER-BASED MULTIPHASE FLOW COMPUTATION
USING ADAPTIVE CARTESIAN GRID TECHNIQUES

















By

RAJKESHAR SINGH


A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL
OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT
OF THE REQUIREMENTS FOR THE DEGREE OF
DOCTOR OF PHILOSOPHY

UNIVERSITY OF FLORIDA


2006

































Copyright 2006

by

Rajkeshar Singh

































To my family and teachers















ACKNOWLEDGMENTS

I would like to express my sincere gratitude to Dr. Wei Shyy for providing me the

opportunity and flexibility to perform this work. I can not thank him enough for being so

patient and understanding for years and pushing us to learn more.

I would like to thank Dr. Renwei Mei, Dr. James F. Klausner and Dr. Samim

Anghaie for agreeing to serve in my thesis committee. Since the surroundings dictate the

quality of life in several ways, I thank my research group members and Dr. Siddharth

Thakur for helping with academic aspects while providing a memorable company for the

past four years.

I thank my wife who has been the best friend for more than half a decade and has

been extremely patient and understanding during all these years. Finally I thank my

parents who have been behind me every step of the way providing their unconditional

support.
















TABLE OF CONTENTS



A C K N O W L E D G M E N T S ................................................................................................. iv

LIST OF TA BLES ............. .................... ............ ......................... .. viii

LIST OF FIGURES ............................... ... ...... ... ................. .x

LIST OF SYM BOLS ............. ............... ........... ............. .. ..... ... xv

ABSTRACT ............. .................................................................. xvii

1 IN TR OD U CTION ............................................... .. ....... .... .............. .

1.1 O objective and C contributions ..................................................................... .....3
1.2 O utline of the D issertation .......................................................................... ... 7

2 LITER A TU R E R EV IEW ................................................................ ...................... 9

2.1 Interface T racking ........... .......................................................... ................ .9
2 .1 .1 L ag ran g ian .................................................................................. 9
2 .1 .2 E u le ria n ............................................................................................... 1 0
2.1.3 M ixed Eulerian-Lagrangian ............................ ............. .. ............. 11
2.1.4 Challenges and Recent Advances.......................................................11
2.1.5 Present Approach: Markers with or without Connectivity......................16
2.2 Interfacial Dynam ics M odeling .................................................................. 17
2.2.1 Sharp Interface............ ................................................ .. ............. 18
2.2.2 C ontinuous Interface ........................................ .......................... 19
2.2.3 Challenges and Recent Advances............................. ................19
2.2.4 Present Approach: Sharp or Continuous Interface.............. ...............21
2.3 A daptive G rid Com putation ........................................ .......... ............... 22
2.3.1 Cartesian Grid D ata Structures............................................................ 24
2.3.2 Present Approach: Unstructured Adaptive Cartesian Grid ....................26

3 MULTIDIMENSIONAL IMMERSED BOUNDARY COMPUTATION ................27

3.1 Im m ersed B oundary M ethod .................................................. ..... ............. 27
3.1.1 Interfacial C conditions ........................................ .......................... 29
3.1.2 Material Property Smoothing .......................... .......................... 31
3.1.3 Momentum Source Term Computation....................... ............... 33


v









3.1.3 Interface A dvection ................. ....... .. ............. ........ ................. 36
3.2 Conservative Restructuring for Interface Resolution ......................................37
3.3 Interface Reconstruction for Topology Change ...........................................40
3.3.1 Reconstruction Criteria............................................... 41
3.3.2 Reconstruction Procedure and Difficulties................. ...............41
3.3.2.1 Tw o-dim ensional interface.................................. ............... 42
3.3.2.2 Three-dim ensional interface............................... ............... 44
3.4 Sum m ary ..................................... ................................ ..........48

4 ADAPTIVE CARTESIAN GRID BASED COMPUTATION............................... 50

4 .1 D ata S tru ctu re ............................................................................................... 5 0
4 .2 G rid A adaptation .................... ............... ............ .. ............ ... ..............52
4.2.1 Interface Geometry based Adaptation................... ..................................53
4.2.2 Flow Solution based Adaptation ................................... ............... ..54
4.3 Immersed Boundary Solution Procedure........................................................55
4.3.1 Staggered Grid Arrangement and Approximations..............................57
4.3.2 Spatial Discretization of Velocity Equation...........................................58
4.3.3 Spatial Discretization of Temperature Equation ......................................60
4.3.3.1 Sharp interface treatm ent .................................... ............... 61
4.3.3.2 Interface temperature treatment ............................................. 61
4.4 Summary ............. ...... .. .......... ........................... ......... 63

5 RESULTS AND DISCUSSION ...................................................... ....................65

5.1 Interface Restructuring: Interface in a Vortex Field ....................................65
5.2 Interface Reconstruction: Effect on Mass Conservation ....................................68
5.3 Single Phase N avier-Stokes Solver ........................................ .....................69
5.3.1 D ecaying V ortex Field ........................................ ......... ............... 70
5.3.2 Lid-driven Cavity ...................................... ..........................................72
5.3.3 Diffusive Heat Transfer between Concentric Cylinders ........................74
5.3.4 Natural Convection in a Square Cavity ............................................76
5.4 Spurious V elocity Currents........................... .............. ............... ... 77
5.4.1 Effect of Fluid Property Jump and Grid Resolution ..............................78
5.4.2 Effect of Interface R construction ....................................... ................80
5.5 Buoyancy Driven Single Rising Bubble.................................. ............... 81
5.5.1 G rid C onvergence Test............... ..... .. ................. ............... ... 82
5.5.2 Rising Bubbles in Different Terminal Shape Regimes ............................84
5.6 Coalescence of Two Rising Bubbles ...................................... ............... 87
5.6.1 H ead-on C oalescence ........................................ .......................... 87
5.6.2 O ff-axis C oalescence................................................... ............... ............90
5.7 B inary D rop C ollision................................................ ............................. 91
5.7 Phase Change Com putation.......................................... .......................... 97
5.7.1 1-D Phase Change ............................................................................. 98
5.7.2 Stationary Bubble Growth............................................... .................. 99

6 SUMMARY AND FUTURE WORK ................................................................113









6 .1 S u m m a ry ........................................................................................................... 1 1 3
6 .2 F u tu re W o rk ............................................................................... 1 15

L IST O F R E FE R E N C E S ......................................................................... .......... ......... 118

BIOGRAPHICAL SKETCH .............................................................. ...............126
















LIST OF TABLES


Table page

1-1 Some industrial and natural occurrences of two-phase flows. ............. ................2

1-2 Key marker based 3D tracking computations in literature............... ...................5

2-1 Summary of some representative methods and related issues. .............................14

2-2 Summary of recent advances in interface tracking. .......................................15

5-1 Error estimates using Equation(5.2) for interface in a time reversed vortex field test
comparing conservative and non-conservative restructuring.................................67

5-2 Volume error for reconstruction of a spherical interface and approximate radial
perturbation required for correction. ............................................. ............... 69

5-3 Max velocity error and computation time normalized with time taken for 16x16
uniform grid. The adapted grids are shown in Figure 5-4....................................72

5-4 Lid driven cavity computation using uniform and adaptive grids .........................73

5-5 N natural convection com putation. ........................................ ......................... 76

5-6 Pressure drop for an inviscid fluids with density ratio =1.0..................................79

5-7 Effect of Laplace number on spurious velocity. .....................................................79

5-8 Effect of density ratio on pressure drop and spurious current on 50x50x50 grid
(viscosity ratio = 1.0, La = 250) .................................................................. 80

5-9 Effect of viscosity ratio on pressure drop and spurious current on 50x50x50 grid
(density ratio = 1000, La = 250)........................................................ ............ 80

5-10 Grid convergence test for Eo = 1.0, M= 1.0x10-3. The density and viscosity ratios
are set to 20.0. ........................................................................ 83

5-11 Computation of rising bubbles with different terminal shape regimes ...................87

5-12 Computational parameter for rising bubble-coalescence tests .............................88

5-13 Properties of tetradecane and nitrogen ...... ..............................93









5-14 Binary drop collision computation with density ratio = 666.081 and viscosity ratio
= 1 7 9 .2 8 ....................................................................... 9 4

5-15 Computational parameters for one-dimensional phase change computation .........98
















LIST OF FIGURES


Figure p

1-1 Illustration of a typical multiphase flow ....................................... ............... 1

1-2 M ixed Eulerian-Lagrangian tracking. ............................................. ............... 5

2-1 Illustration of Lagrangian interface tracking........... ..... .................10

2-2 Schematic of Eulerian tracking using level set and volume of fluid...................... 11

2-3 Interface in a time reversed vortex field using level-set method. ..........................12

2-4 A conceptual interface reconstruction using piecewise linear segments in the
volum e-of-fluid m ethod. ........................................... ......................................... 13

2-5 A surgical alteration of the interface data to accomplish topology change .............15

2-6 A three-dimensional marker based film boiling computation ..............................15

2-7 Illustration of conceptual differences between marker-based tracking with and
without connectivity inform ation. ........................................ ........................ 17

2-8 Surface tension, pressure and normal component of shear stress acting on an
in terfa c e ............................................................................................... 1 8

2-9 Sharp and continuous interface methods............................................ .............21

2-10 The effect of fluid property jumps on the computation time for a rising bubble.....23

2-11 A cell-by-cell and block-by-block grid adaptation. ............................................25

2-12 A tree-based adaptive Cartesian grid-data. ................................... ............... 25

3-1 Schematic of the immersed boundary method. ....................................................28

3-2 Probe based technique for temperature gradient computation on the interface.......30

3-3 A triangulated interface and the data-structure. ........................................... ........... 31

3-4 Indicator function and the material flags for cells completely outside or inside the
in te rfa c e .......................................................................... 3 3









3-5 Momentum source term due to surface tension force on a circular interface. .........35

3-6 Computation of the unit normal and tangent vectors on interface triangles. ...........35

3-7 Curvature computation for a unit circle using Equation(3.18) and cubic spline with
(a) tw elve and (b) sixty points ....................................................... ...................36

3-8 Interface restructuring approach .........................................................................38

3-9 Conservative marker deletion for two-dimensional interfaces..............................39

3-10 Conservative marker deletion for three-dimensional interfaces ............................39

3-11 Level-contour based interface reconstruction. ................................. ...............40

3-12 Interface reconstruction criterion uses two probes from the interface along the
normal vector and bi-linearly interpolates the indicator values Indl and Ind2........41

3-13 Two dimensional reconstructed interface edges (el, e2 and e3) and corresponding
connectivity data. ................................................... ................. 42

3-14 Illustration of 2D cell-by-cell interface reconstruction. .........................................43

3-15 A 2D cell with four markers creating two locally disconnected edges ..................44

3-16 A reconstructed interface segment in a cell. ................................. .................45

3-17 Some common reconstructed interface segments and integer vertex-flags. ...........45

3-18 Difficulties due to degeneracy of reconstructed interface data and implemented
re m e d y ........................................................................... 4 6

3-19 Orientation setting of reconstructed interface. ............. ................. ..................... 47

3-20 Examples of 3D interface reconstruction. ...................................... ..................47

4-1 Cartesian grid data-structure. .............................................................................. 51

4-2 A daptive grid refinem ent. ............................................................. .....................52

4-3 Cartesian cells sharing a face and the corner cells are not allowed to differ by more
than one level of refinement.........................................................................53

4-4 Geom etry based adaptation. ............................................. ............................ 54

4-5 Staggered grid arrangem ent. ............................................ ............................ 58

4-6 Spatial discretization of convection term ..................................... ...............59









4-7 A one-dimensional interface between i and i+. .............................................. 61

4-8 Linear interface temperature correction within one cell on each side of the
in te rfa c e .......................................................................... 6 2

4-9 Presently employed interface temperature correction technique solves a correction
equation sharply on each side of the interface. .............................. ......... ...... .63

5-1 Interface in a time reversed vortex field ................. ............. ....... ........... 67

5-2 Error estimate and grid convergence for interface in time reversed vortex field.....67

5-3 Grid convergence test for volume error due to reconstruction..............................69

5-4 D ecaying vortex fi eld. ...................................................................... ...................71

5-5 Grid convergence test for decaying vortex field showing maximum velocity error
computed on uniform and adaptive grids ...................................... ............... 71

5-6 Streamlines for lid-driven cavity simulation for Re = 100....................................72

5-7 L id driven cavity. ............................................. ................... .... ...... 73

5-8 The u velocity profile along a vertical line through the center of lid-driven cavity
using (a) uniform grids and (b) adaptive grids..................................................... 74

5-9 Indicator contour used for modeling the cylindrical computational domain. .........75

5-10 Computed temperature and error for heat transfer in a concentric cylinder. ...........75

5-11 Natural convection with Pr = 0.71 and Ra = 1.0x105. ...........................................76

5-12 Static bubble com putation ............................................... ............................. 78

5-13 Grid convergence tests for static bubble. ..................................... ............... 79

5-14 Effect of reconstruction frequency on spurious velocities ............. ...............81

5-15 Computational set-up for a single rising bubble. .............................................. 82

5-16 A spect ratio com putation .............................................................. .....................83

5-17 Bubble rise with Eo = 1.0, M= 1.0x10-3, density and viscosity ratio of 20.............84

5-18 Computational domain for cases in Table 5-11 .......................................... ........... 85

5-19 Computed shapes from Table 5-11 .........................................................................86

5-20 Computed case forM = 0.971, Eo = 97.1. .................................... .................86









5-21 Computational set-up for head-on coalescence test. ..........................................88

5-22 H ead-on coalescence. ...... ........................... ........................................... 89

5-23 Grid and streamlines at t = 0.125 sec for head-on coalescence..............................89

5-24 Instance of interface reconstruction during head-on coalescence of two rising
bubbles at t = 0.116 s .................................... ................ ................... 90

5-25 Experimental photographs of bubble coalescence taken at 0.03s time intervals. ....90

5-26 O ff-axis coalescence. ...................................................................... ................... 9 1

5-27 Binary drop collision regimes marking the location of computed parameters with
so lid c irc le s ................................................................................................ ..... 9 3

5-28 Impact parameter B is defined as B=hd. ...................................... ............... 93

5-29 Binary drop collision Case 1 from Table 5-14..................................................100

5-30 A cross section of the computational grid and the time history of grid data-size for
Case 1 from Table 5-14........................................ .......... .... ................. 101

5-31 Binary drop collision Case 2 from Table 5-14..................................................102

5-32 A cross section of the computational domain showing the interface and the velocity
vectors for Case 2 from Table 5-14.................................... ........................ 102

5-33 A cross section of the computational grid and the time history of grid data-size for
Case 2 from Table 5-14........................................ .......... .... ................. 103

5-34 Binary drop collision Case 3 from Table 5-14..................................................104

5-35 A cross section of the computational grid and the time history of grid data-size for
Case 3 from Table 5-14........................................ .......... .... ................. 105

5-36 A cross section of the computational domain showing the interface and the velocity
vectors on the left and pressure profile on the right for Case 3 from Table 5-14..106

5-37 Binary drop collision Case 4 from Table 5-14..................................................107

5-38 A cross section of the computational grid and the time history of grid data-size for
Case 4 from Table 5-14........................................ .......... .... ................. 108

5-39 A cross section of the computational domain showing the interface and the velocity
vectors for Case 4 from Table 5-14.................................... ........................ 109

5-40 One dimensional phase change setup ...................................... ........................... 110









5-41 Comparison of computed and theoretical solution for Case 1 in Table 5-15.........110

5-42 One-dimensional phase change computation for Case 2 in Table 5-15 showing
interface location on the left and temperature distribution at t = 0.1 on the right.. 111

5-43 Stationary bubble growth rate. ............ ............................ .......... ......... 112















LIST OF SYMBOLS


p, ,o- Density, viscosity and surface tension

C, K, A Heat capacity, thermal conductivity and latent heat of evaporation


m Interfacial mass transfer rate (-ve for evaporation)

k, f,n Interface local curvature, surface tension force, unit normal vector

F Momentum source term due to surface tension force

I Indicator function

5,D Dirac-delta and discretized Dirac-delta function

d Bubble diameter

U =(u, v, w) Velocity

P Pressure

T Temperature

g Gravity

()mt A property 4 on interface markers

(0), A property 4 of fluid outside the interface

(0)2 A property 4 of fluid inside the interface


)i 2 Ratio of a fluid property 4












(x,y,z)


Re= pUL



/a
We = pU2L


Fr = U


Ja = p CA
/pj


Pe =pCUL/K


Pr = C/K


La = pYL/
/gt

Ra = P2CgfATL3 K


Ca = U ma


M = gt3


Eo = pgL


Saturation temperature

Spatial coordinates along x, y and z coordinate axes

Length scale

Reynolds number


Weber number


Froude number


Jacob number


Peclet number


Prandtl number


Laplace number


Rayleigh number


Capillary number


Morton number


Eotvos number















Abstract of Dissertation Presented to the Graduate School
of the University of Florida in Partial Fulfillment of the
Requirements for the Degree of Doctor of Philosophy

THREE-DIMENSIONAL MARKER-BASED MULTIPHASE FLOW COMPUTATION
USING ADAPTIVE CARTESIAN GRID TECHNIQUES

By

Rajkeshar Singh

August, 2006

Chair: Wei Shyy
Major Department: Mechanical and Aerospace Engineering

Multiphase flows associated with interfacial dynamics, steep jump in fluid

properties and moving boundaries between different phases pose substantial

computational challenges in terms of both modeling as well as computational cost. The

present work uses an immersed boundary technique to model the interfacial dynamics,

Lagrangian markers to track the moving phase boundaries, and a stationary Cartesian grid

to solve all the flow governing equations. Time dependent triangulated surface meshes

are employed to represent the time dependent interfaces shape and location. Based on the

solution characteristics, multi-level, three-dimensional adaptive grid techniques are

incorporated into the computational framework to help meet the resolution requirements.

Furthermore, a conservative marker redistribution technique is developed to maintain a

desired marker-spacing, and a connectivity preserving level contour-based reconstruction

technique is devised to handle topological changes associated with the interfacial

dynamics. The flow equations were solved using the projection method with a finite









volume staggered grid formulation on adaptive grids. For phase change problems,

accuracy of the mass transfer computation critically affects the overall computational

outcome. Efforts have been made to address this via a sharp interface-based mass transfer

mode combined with the immersed boundary method. The capabilities and accuracy of

the individual components and overall computational system are tested with a range of

computations including demonstration of improvements with conservative interface

restructuring, reconstruction and its effect on immersed boundary solution accuracy,

estimation of computational cost saving with adaptive grids, rising bubble coalescence,

binary drop collision, and stationary bubble growth in a superheated liquid pool. The

method has demonstrated its capability for handling high density ratio, 0(1000),

multiphase fluid dynamics.


xviii















CHAPTER 1
INTRODUCTION

Multiphase flows are characterized by two or more fluids in relative motion. The

constituent fluids usually have different physical properties separated by an interface that

evolves with the flow. Figure 1-1 depicts some of the defining elements of a typical

multiphase flow- a moving/deforming surface that separates different fluids/phases;

discontinuous fluid properties (density, viscosity, conductivity, etc.) across the interface;

interfacial effects such as surface tension, phase change, etc.


Fluid 1
Interface

Surface tension force
Fluid 2






Figure 1-1. Illustration of a typical multiphase flow.

Presence of such flows in every day life can be recognized by simple examples of

ocean waves, boiling water, blood flow in the body. Some specific practical applications

may be cited as boiling heat transfer in power industries, biochemical and metallurgical

industries, cavitation and ultrasound technology for medical and industrial applications.

The constituents in these flows may be liquid-liquid, gas-liquid, liquid or gas with solid.

Table 1-1 by Kleinstreuer (2003) gives a summary of some of the industrial or natural

occurrences of such flows. Since consideration of all the combinations in Table 1-1 is a









prohibitively difficult task, the scope of the present work is controlled by focusing only

on incompressible liquid-liquid or gas-liquid flows for bubble and drop dynamics

problems.

Table 1-1. Some industrial and natural occurrences of two-phase flows.
Phase Phase configuration Occurrence
combination
Gas-Solid Dispersion of solid Aerosol, particulate matter pollution; filters and
particles in gas particle-collection devices; Pneumatic transport of
particles; gas-fluidized beds

Liquid-Solid Dispersion of Cells and particles in blood; hydraulic transport of
particles in liquid particles; liquid-fluidized beds

Gas-Liquid or Dispersion of liquid Mist, fog, cloud; coalescence of drops and formation of
Liquid Gas drops in gas rain; droplet removal devices; fuel droplet injection;
spray coating and cooling

Dispersion of gas Transport of oil and gas; foam formation;
bubbles in liquid evaporators, boiling flow in pipe

Liquid-liquid Dispersion of liquid Emulsion and creams; break-up of immiscible drops;
drops in liquid coalescence of drops and phase separation

Both bodies of Transport of two liquids in horizontal pipe; movements
liquids connected of oil-water interface in porous rock.
Note: Adapted from Kleinstreuer (2003)

Presence of multiple time and length scales makes numerical simulation of such

flows a difficult and expensive affair. Time dependent interfaces with topology changes

(bubble coalescence or breakup) and interfacial mass/heat transfer add further to the

computational complexity and cost. Recent years spanning more than a decade have

witnessed an ever increasing computational activity with the research spanning various

aspects ranging from enhancement of computational efficiency to interface tracking and

flow modeling. The computational expense and difficulties in handling interfaces are in

part responsible for relatively fewer computations in 3D as compared to a vast number of

2D or axisymmetric cases. These limitations are even more visible with phase change









computations. Besides the high computational cost, phase change computations present

several other difficulties with the mathematical modeling and computational algorithms.

Some of the modern day phase change computations are Dhir (2001), Son (2001), Son

et al. (2002) and Luo et al. (2005) using level-set; Welch and Wilson (2000) using

volume-of-fluid; Esmaeeli and Tryggvason (2004), Tryggvason et al. (2005), Shin et al.

(2005) and their earlier work using marker based tracking.

Although continuing developments in computational techniques and resources are

making the numerical simulation more and more tractable, several algorithmic and

modeling challenges remain unresolved. With further discussion, it will be highlighted in

the following chapter that some of the computational aspects often have conflicting

requirements and demand several compromises to be made based on the problem of

interest.

1.1 Objective and Contributions

The primary goal of the present work is to develop an effective multidimensional

interfacial flow computation-system. The process involves assembling key ideas from

different approaches and techniques into a unified framework along with suitable

improvements and modifications wherever necessary. Development of a parallel

computing platform is an essential requirement for large scale computing. However, the

current efforts are part of a larger research group activity and thus do not cover parallel

computing aspects in this document. The computational components are divided into

three generic categories: interface tracking; flow modeling; flow computation. The

general emphasis and objectives in these categories are:

* Multi-fluid interface tracking
o Time dependent triangulated surface grids for interface representation.
o Maintain resolution by locally adding/deleting nodes on the surface grid.









o Capability to handle complex topology changes.
* Interfacial dynamics and flow modeling
o Single fluid formulation for the mass and momentum conservation equations.
o Surface tension forces modeled as smoothed volumetric source term in the
momentum equation.
o Accurate computation of the interfacial mass transfer rate for phase change
problems.
* Adaptive grid computation
o Grid generation and dynamic adaptation based on interface location and flow
features.
o Adaptive grid base incompressible Navier-Stokes computation.

The interface is tracked on a stationary Cartesian grid using time dependent

triangulated surfaces (Figure 1-2). This approach is also referred to as a marker based

tracking method. Since all the flow computations are performed on the background

stationary grid and the interface is tracked using Lagrangian markers, the overall

approach is of a mixed Eulerian-Lagrangian type. Although this approach has several

advantages due to explicit knowledge of the interface shape/location for computing the

geometric properties, the difficulties in data-handling have deterred a wider usage.

Table 1-2 presents some of the key usages and techniques demonstrated to handle

complex multidimensional interfacial flows with marker based tracking. The present

method is similar to that of Tryggvason et al. (2001) with primary contributions made in

interface handling. The triangulated-interface is restructured by locally adding or deleting

markers to avoid undesirably long or short edges. The commonly used restructuring

techniques in literature (Tryggvason et al., 2001, 2005; Sousa et al. 2004) introduce some

amounts of error in the interface volume for which a conservative algorithm is proposed.

Topology changes were handled using the level-contour based interface reconstruction of

Shin and Juric (2002) with appropriate algorithmic considerations to establish a valid

interface connectivity data.



















Interface markers
(a)


(b)
(b)


Figure 1-2. Mixed Eulerian-Lagrangian tracking. (a) A two dimensional interface using
connected markers on a background grid and (b) a triangulated interface.

Table 1-2. Key marker based 3D tracking computations in literature.
3D Explicit Author Summary of technique
Tracking

A set Shin and Juric Easy algorithms for computing geometric properties and
unconnected (2002) dealing with complex topology changes.
triangles Reduced flexibility in controlling interface resolution and
dealing with multiple interfaces.

A set of Torres and Easy to handle complex topology changes but interface
unconnected Brackbill geometry computation more involved that Shin and Juric
points (2000) (2002)
Reduced flexibility in controlling interface resolution and
dealing with multiple interfaces.

Triangulated Tryggvason et Flexible control over interface resolution and dealing with
surface grid. al. and group multiple interfaces.
(2001, 2005) Uses reconstruction technique of Shin and Juric (2002) to
handle topology changes.

The interface moves on the background Cartesian grid that is dynamically adapted

based on the interface location and flow solution to optimize the computational cost. All

the flow governing equations are solved on background grid. The data-structure and

Cartesian grid generation is based primarily on the work of Aftosmis (1997). The Navier-

Stokes solution is performed using projection method with a staggered grid based finite-

volume formulation. The underlying non-uniform nature of the grid makes staggered grid

relatively more difficult to implement than the corresponding collocated grid (Singh &









Shyy, 2006), but provides a compact pressure-velocity coupling. As detailed in

subsequent chapters, several simplifying approximations in the numerical discretization

were made for a convenient flow solver implementation on the spatially adaptive grid.

These approximations are present only on the coarse-fine grid-interfaces which are kept

away from the physical location of the interface by suitable refinement criteria.

Most of the interfacial dynamics modeling and development uses the already

available concepts from literature. The underlying single fluid formulation for mass and

momentum equation with marker based tracking is commonly recognized as the

immersed boundary method (Peskin, 1977; Tryggvason et al., 2001). This method smears

the fluid properties across the interface and models the surface tension forces as

smoothed volumetric source terms in the momentum equation. Although the mass and

momentum equations use immersed boundary method, it is possible to avoid such

approximations for the temperature equation. Since interfacial mass transfer depends on

the jump in temperature gradients across the interface, an attempt has been made to treat

the temperature equations sharply in order to accurately compute the mass transfer rate.

The immersed boundary method of Juric and Tryggvason (1998) uses an iterative

procedure to maintain the interface temperature while the present method uses a simple

linear correction similar to the level-set based computations of Morgan (2005). At this

point, the current approach deviates from a conventional immersed boundary method and

borrows the philosophy from recent developments in level-set and sharp interface ghost

fluid method of Luo et al. (2005) and Gibou et al. (2003). The presented method solves

the temperature equation on both sides of the interface and does not require assuming

fixed temperature inside the interface as in the work of Son (2001).









The specific contributions of the present work may be listed as:

1. Marker based interface tracking
a. Conservative interface restructuring technique to maintain resolution without
introducing errors in the interface volume.
b. Development of interface reconstruction algorithm for topology change.
Several challenges/difficulties are highlighted along with simple remedies.
c. Conservative restructuring along with reconstruction technique offers an
enhanced flexibility in handling multiple interfaces as compared to a
connectivity-free tracking method of Shin and Juric (2002).
2. Interfacial dynamics modeling
a. Implementation of multidimensional immersed boundary method.
b. A sharp interface treatment for temperature equation for accurate interfacial
mass-transfer computation during phase change. Strategies to incorporate
level-set based sharp interface techniques within the present marker based
tracking are presented.
3. Adaptive grid computation
a. Implementation of an adaptive Cartesian grid generator with dynamic
adaptation based on the interface location and flow solution.
b. Development of a staggered grid, finite volume based Navier-Stokes
computation algorithm for velocity and pressure computation.
c. Development of a sharp interface, collocated grid computation procedure for
temperature equation in heat-transfer/phase change problems while using
marker based tracking.

The capabilities of the developed approach and algorithms have been demonstrated

via a set of complex multiple rising bubbles and binary drop collision computations as

one of the few known cases in literature using marker based tracking. For example, the

only known 3D binary drop collision with marker based tracking are by Nobari and

Tryggvason (1996); Shin and Juric (2002) performed for moderate density ratios of 0(10)

while present computations were performed for flows with density ratios of 0(1000).

1.2 Outline of the Dissertation

The document is organized into six chapters. Chapter 2 presents brief review of

alternative approaches in literature. The basics of immersed boundary modeling and

related key components dealing with interface management are presented in Chapter 3.









Chapter 4 presents the adaptive Cartesian grid generation technique and discusses

the Navier-Stokes solution procedure highlighting various simplifying assumptions in the

numerical discretization.

The key results are presented in Chapter 5. Tests characterizing and demonstrating

interface restructuring, reconstruction, single and multiple rising bubbles and binary drop

collision are presented along with preliminary assessment of the accuracy of phase

change computation.

Chapter 6 provides a summary of the presented work and highlights various issues,

approximations following a brief comment on future extensions and refinements.














CHAPTER 2
LITERATURE REVIEW

This chapter provides a brief review of some of the popular techniques based on the

categories outlined in the previous chapter: interface tracking; interfacial dynamics

modeling; adaptive grid based computation. Various pros and cons of these techniques

are visited in an attempt to place the current choices and highlight some of the relevant

issues and difficulties to be addressed in following chapters.

2.1 Interface Tracking

As reviewed by Shyy et al. (1996), the three main categories are Lagrangian,

Eulerian and mixed Eulerian-Lagrangian. Due to the extensive amount of available

literature, the scope of this section is limited to providing a compact background and

some of the key issues with these techniques along with the state-of-the-art

developments.

2.1.1 Lagrangian

Lagrangian method discretizes the computational domain with the interface as a

boundary for creating a body fitted grid (Figure 2-1). This approach requires adjustment

or regeneration of the grid in order to accommodate changes in the interface geometry.

Some usages of this approach can be found in the works of Ryskin and Leal (1984),

Raymond and Rosant (2000), Wasekar and Manglik (2003) using structured curvilinear

grids and work of Perot and Nallapati (2003) using triangulated (tetrahedrons in 3D)

grids.









Although the body fitted nature of the grid makes it a desirable choice, handling of

body fitted grids require sophisticated grid generation capabilities especially when the

interfaces undergoes topology change. Observing the frequency of the usage, this is not

among the most popular choices and hence it will not be considered further.


FI id 1
SInterface




Flu d


Figure 2-1. Illustration of Lagrangian interface tracking.

2.1.2 Eulerian

Eulerian methods use a scalar function to represent the interface implicitly. The

level-set (Osher & Fedkiw, 2001, 2003) and volume-of-fluid (Youngs, 1982; Rider &

Kothe, 1998; Scardovelli & Zaleski, 2002) are examples of Eulerian tracking. This

method 'captures' the interface using the information contained in a scalar function. This

scalar function is a signed distance function in the level-set method and volume fraction

in the volume-of-fluid method (Figure 2-2)

Unlike the Lagrangian methods, Eulerian methods track the interface on a

stationary grid using a transport Equation(2.1) for the scalar function evolution that

contains the interface information. Due to the Eulerian nature of the interface, this

method does not necessitate any changes in the computational grid as required by

Lagrangian methods.

S+V.VO-= 0 (2.1)
at










03 0
0.11
o .3 ....- .... --





1 1 0.8

(a) (b)

Figure 2-2. Schematic of Eulerian tracking using level set and volume of fluid. (a) A
volume of fluid captures the interface (solid curve) using volume fraction
information; (b) level set recognizes the interface as zero contour of signed
distance function D.

2.1.3 Mixed Eulerian-Lagrangian

This class of method derives its name from the fact that it tracks the interface using

Lagrangian markers moving on an Eulerian computational grid (Figure 1-2). These

markers move on the computational grid using a simple advection Equation (2.2) where

the marker velocities are interpolated from the Eulerian grid. Thus, (similar to

Lagrangian) any geometric information about interface is explicitly available and (similar

to Eulerian methods) the interface evolution does not require any modification to the

computational grid (Tryggvason et al., 2001; Glimm et al., 1998, 2000; Singh & Shyy,

2006).

aT Unt (2.2)
at

2.1.4 Challenges and Recent Advances

The above introduced methods for interface tracking have some associated

advantages and disadvantages that make general preference of any particular method

difficult. Table 2-1 gives a summary of the key features of these methods. The Eulerian

methods are among the easiest to implement and some of their most cited advantages lie









in their natural ability to handle interfaces undergoing topology changes. On the flip side,

level-set method (LS) suffers from erroneous mass loss that may become unacceptably

large and detrimental. Such a mass-loss problem is highlighted by the errors shown in the

tests conducted by Enright et al. (2002) where a spherical interface undergoes severe

deformation and it is expected to return to original shape when placed in a time-reversed

vortex field. Figure 2-3(a) shows the sever loss in volume of a spherical interface in such

a field using level-set method for interface tracking. This drawback was ameliorated by

Enright et al. (2002) using a particle level-set method that uses Lagrangian markers

around the interface to avoid erroneous breakup/deletion of the level-set characteristics

(Figure 2-3(b)).


"p"i(^ W% c2A


...... -
-fe--i--8


1' 31% 641k


-.me-.M -.


(a) (b)

Figure 2-3. Interface in a time reversed vortex field using level-set method. (a) Severe
mass loss with level set method; (b) improvements produced by particle level-
set method (Enright et al., 2002). [Reprinted with permission]

Unlike the level-set method, the volume-of-fluid tracks the volume (mass)

explicitly and hence does not suffer from mass-loss problems. However, volume-of-fluid

methods have been observed to produce spurious flotsam ('floating wreckage') and

jetsam ('jettisoned goods') in the regions of high interface-curvature under-resolved by

the computational grid (Rider & Kothe, 1998). The most serious difficulty with the VOF

is due to errors in computation of interface curvature required for surface tension









computation. These difficulties arise due to the sharply changing nature of the volume-

fraction and difficulties in reconstructing a continuous interface using the volume-

fraction information (Figure 2-4). A more comprehensive detail of various VOF based

interface reconstruction and subsequent developments could be found in Rider and Kothe

(1998), Scardovelli and Zaleski (2002), Renardy and Renardy (2002), Meier et al. (2002).

0.3 0


1 0.9 \ .4


1 1 0.8


Figure 2-4. A conceptual interface reconstruction using piecewise linear segments in the
volume-of-fluid method.

The mixed Eulerian-Lagrangian methods possess several desirable features when

compared with Lagrangian or Eulerian methods. Due to explicit definition of interface

using markers, the interface does not inaccurately diffuse in time and loose mass as level

set methods. It has also been noticed that explicit tracking does not require as high grid

resolution as Eulerian methods (Glimm et al., 1998). Explicit definition of the interface

also avoids inaccuracies caused by the errors in interface reconstruction for VOF

methods.

The most serious concern against the Eulerian-Lagrangian methods has been raised

due to the complexity of interface data that becomes difficult to handle especially for

topology changes. These difficulties arise due to the algorithmic complexity involved in

'surgical' altering of the interface data to accommodate topology changes. The source of

difficulty lies in establishing valid connectivity information. Although such a procedure









is easily accomplishable for two-dimensional interfaces (Figure 2-5), extension to three

dimensions is a much more involved procedure. The potential complexities in marker-

based tracking may be highlighted from the interface shapes seen in some of the 3D film

boiling computations of Shin and Juric (2002) (Figure 2-6).

Several attempts have been made by researchers to simplify the marker based

algorithm by resorting to elimination of the need to maintain the connectivity

information. One of the first such attempts is due to Torres and Brackbill (2000) who

used a set of unconnected markers to describe the interface. Lack of the need to maintain

explicit connectivity information, this approach makes it easier to handle topological

changes. The subsequent development of connectivity free tracking using unconnected

triangular elements by Shin and Juric (2002) simplified the algorithms further. Table 2-2

presents some key developments by combining various individual techniques to

exploiting their attractive features.

Table 2-1. Summary of some representative methods and related issues.
Methods Advantages Issues

Lagrangian Interface explicitly known and part of Frequent grid regeneration.
the computational grid. Difficult and expensive for
Geometric property computed directly complex 3D computations.
on interface.

Eulerian Easy geometry-computation. Severe mass loss/gain in under-
(Level-set) All computations on Eulerian grid. resolved areas of the interface.
Easy to handle topology changes.

Eulerian No mass loss/gain. Poor performance in under-
(Volume of All computations on Eulerian grid. resolved areas of the interface..
Fluid) Easy to handle topology changes. Difficult geometry computation.

Mixed Interface explicitly known even at Complex data-structure and
Eulerian- sub-grid levels, book-keeping for 3D interfaces.
Lagrangian Flow computations on Eulerian grid. Difficult for topological change.










Table 2-2. Summary of recent advances in interface tracking.


Base method

Level-set
(LS)

Volume of
fluid (VOF)


LS and VOF


Marker
tracking


Improved method

Particle level-set
(Enright et al., 2002)

Hybrid marker-VOF
(Aulisa et al., 2003,
2004)

Coupled LS -VOF
(Sussman, 2003)

The point-set tracking
(Torres & Brackbill,
2000)


Connectivity free
tracking with level-
contour reconstruction
(Shin & Juric, 2002)


Features


Ameliorates the excessive/erroneous mass loss using a
cloud of markers around level-set interface.

* Employs makrer based interface to avoid interface
reconstruction for geometric-property computation.
* Improved performance in under-resolved areas.

* Excellent mass-conservation due to VOF.
* Easy geometry-computation using Level-set.

* A set of unconnected.
* Lack of connectivity-information makes it easier to
handle topological changes.
* Geometry-computations require spline-fitting.

* A set of unconnected triangles.
* Possible to compute geometry (curvature, normal)
without any surface fitting or interpolation.
* Level-contour reconstruction for topology change.


1 6


Figure 2-5. A surgical alteration of the interface data to accomplish topology change.


t~Ep~ a


Figure 2-6. A three-dimensional marker based film boiling computation. (a) A snapshot
of the computational domain and (b) a triangulated interface before and after
pinch-off from Shin and Juric (2002). [Reprinted with permission]


1~









2.1.5 Present Approach: Markers with or without Connectivity

Since it is evident that the Eulerian type approaches are relatively easier, it makes

the choice of the method difficult based on any single consideration. Since the most

serious difficulty with marker based tracking is due to data-handling, the explicit

knowledge of the interface shape/location and recent algorithmic advances in literature

have been the criteria for choosing explicit marker based method.

The method of level-contour based reconstruction with connectivity-free marker-

based tracking of Shin and Juric (2002) offers an attractive choice. However, the lack of

connectivity information on the interface makes it difficult to locally alter the marker-

spacing by marker addition and deletion to maintain a desired marker-spacing on the

interface. As seen in Figure 2-7(a), deletion of an edge e2 requires collapsing of vertices

of the neighboring edges (point p2 and p5) but the lack of connectivity makes the

identification of edges el and e3 difficult i.e., no explicit information is available to

suggest, for example, that nodes p2 and p3 have the same physical location. For

computations such as bubbly flows by Esmaeeli and Tryggvason (1998, 1999), several

bubbles (-100 or more) move in close proximity without merger. In these simulations,

the bubbles are not allowed to merge with each other but local restructuring of the

interface is still required to maintain desirable resolution. The connectivity-free method

resorts to global reconstruction of the interface even to maintain the marker spacing and

may inadvertently produce mergers. For such reasons, the current work maintains the

connectivity information (Figure 2-7(b)) to allow greater flexibility in handling multiple

interfaces. This decouples the procedure for handling interface resolution and interface

reconstruction by using local restructuring when needed and performing reconstruction

only to achieve topology changes.









p2 3 p45 p2 p4
e2 p5

el e3 el e1 e
p6 p6

(a) (b)

Figure 2-7. Illustration of conceptual differences between marker-based tracking with and
without connectivity information. (a) Connectivity free method of Shin and
Juric (2002) does not store any connectivity information between edges; (b)
connectivity based approach maintains such information explicitly.

2.2 Interfacial Dynamics Modeling

A moving/deforming interface in the computational domain presents some obvious

difficulties in numerical computation of the flow field. The difficulties arise primarily

due the nonlinear coupling of fluids and steep fluid-property jumps across the interface.

In the event of phase change, effects such as mass and heat transfer across the interface

also need to be considered. For simplicity, the basic introduction is provided here with

respect to isothermal flows. Issues and techniques related to phase change modeling will

be considered in subsequent sections and chapters.

In general, the pressure and viscous stresses show discontinuities across the

interface that is related to the surface tension force and fluid property jumps. Figure 2-8

shows the pressure (P) and the normal shear stress components (n.r-n) across the

interface and Equation (2.3) relates the jump in flow properties with the surface tension

force. Based on the level of detail and accuracy in treatment of the discontinuities, the

existing methods can be classified into two categories: sharp interface method (SIM);

continuous interface method (CIM). The sharp interface method (SIM) is a class of

techniques that attempt to satisfy the jump condition (Equation(2.3)) explicitly while









continuous interface method (CIM) represents techniques that smooth the jump over few

cells across the interface.

(P, -)- (T T) -n = oK (2.3)


Fluid 1

n-T -n




Fluid 2

Figure 2-8. Surface tension, pressure and normal component of shear stress acting on an
interface.

2.2.1 Sharp Interface

There are three prominent techniques within this category: sharp interface of

Ryskin and Leal (1984), Ye et al. (1999, 2001); immersed interface (LeVeque & Li,

1994; Li & Lai, 2001); ghost fluid (Liu et al., 2000; Kang et al., 2000). The technique of

Ryskin and Leal (1984) and Ye et al. (1999, 2001) are the most uncompromising in terms

of satisfying the jump conditions across the interface. The basic approach can be

explained with the cut-cell method of Ye et al. (1999, 2001). This method cuts the

background Cartesian grid with the interface and creates sub-domains (Figure 2-9(a)) that

act like body fitted grids with arbitrary shaped cells on the interface. For stability reasons,

the cut-cells are merged with suitable neighbors to remove arbitrarily small cells. The

governing equations are solved separately in each sub-domain with the corresponding

fluid properties. At every time step, an iterative procedure, by perturbing the interface in

local normal direction, is followed to find the shape/location of the interface such that

Equation (2.3) is satisfied along with the mass and momentum conservation equations.









The immersed interface and the ghost-fluid techniques use the jump conditions to

construct boundary conditions for flow variables on the interface. As an example, the

ghost-fluid technique solves the governing equations in each phase separately by creating

ghost-cell values using the analytical jump condition. Unlike the sharp interface method

of Ye et al. (1999, 2001), these techniques do not solve for the interface shape and flow

field simultaneously or create cut-cells.

2.2.2 Continuous Interface

In contrast to the sharp interface method, the continuous interface method smoothes

the fluid property and models the surface tension forces as a smoothed momentum source

term (Peskin, 1977). This facilitates a single fluid formulation for the entire domain with

the flow/fluid properties varying smoothly across the interface over a thin zone. This

smearing/smoothing region across the interface shown in Figure 2-9(b) is usually three to

four cells in thickness. The continuous interface method has been widely used with mixed

Eulerian-Lagrangian tracking (i.e., the immersed boundary method) (Peskin, 1977;

Tryggvason et al. 2001; Francois & Shyy, 2003) as well as Eulerian tracking such as

level-set and volume-of-fluid.

2.2.3 Challenges and Recent Advances

The iterative procedure of coupling the phases across the interface and the complex

procedure of cell-cutting makes the extension of the sharp interface technique of Ye et al.

(1999, 2001) to three dimensions a tedious and challenging task. Although the ghost fluid

is among the simplest sharp interface method to use, it has traditionally been used with

level-set based interface representation. Although successful demonstration of ghost fluid

type approach with Eulerian-Lagrangian tracking were made by Hao and Prosperetti









(2004) for some bubble dynamics problems, further research is required to establish the

overall stability and performance.

Accuracy studies by Ye et al. (2004) and Shyy (2004) comparing the sharp

interface and continuous interface methods for a spherical drop in static equilibrium

sheds some light on their characteristic differences. With a spherical drop in static

equilibrium, the difference in pressure is balanced by surface tension force according to

the Young-Laplace Equation(2.4). The sharp interface method attempts to capture the

pressure jump as a discontinuity whereas the continuous interface method smoothes the

pressure jump over few cells (Figure 5-12(b)).

AP = oK (2.4)

For a static bubble computation, the theoretical solution produces uniformly zero-

velocity in the domain. However, due to certain amount of net numerical stress-

imbalance across the interface, computations tend to produce a non-zero velocity field

termed as parasitic currents or spurious velocities (Figure 5-12(a)). These spurious

velocities can destabilize the interface especially if viscous forces are not strong enough

(Rider & Kothe, 1995). As an observation, the magnitude of the spurious velocity can be

kept orders of magnitude smaller with the sharp interface treatment. The second order

accurate sharp interface technique of Ye et al. (2004) could keep this velocity up to

machine level as compared to continuous interface methods producing velocities of the

order 10-4 (in terms on Capillary number Ca) (Tryggvason, 2001). The first order

accurate ghost fluid technique could also keep these velocities to the order of 10-6 or

lower (Kang et al., 2000). Several attempts using improved curvature/surface-tension









computation have also been attempted in literature to keep such erroneous velocities as

low as possible (Renardy & Renardy 2002; Shin et al., 2005).


cut by interface






(a)


Fluid properties smeared
(b)

Figure 2-9. Sharp and continuous interface methods. (a) Sharp interface cut-cell method
and (b) a continuous interface using the immersed boundary technique.

2.2.4 Present Approach: Sharp or Continuous Interface

The continuous interface method provides a numerically stable system due to

smoothing of discontinuities but leads to smearing of the interfacial effects. This aspect

of the continuous interface method poses accuracy concerns for certain types of problems

like solidification dynamics at morphological scales that require higher precision in

interface treatment (Udaykumar et al., 1999).

For the reasons of computational efficiency and the relative difficulties in using the

sharp interface method with marker based tracking, the present work employs the

continuous interface method via immersed boundary technique. Despite the accuracy


Fluid 1


Fluid 2


Fluid 1



Fluid 2









issues related to smearing of discontinuities, it is still a reasonably accurate method for a

large number of drop dynamics problems. Over the years, the immersed boundary

technique has been used extensively to study a wide range of physical processes: Peskin

(1977) for studying blood flow in heart-valves; Juric and Tryggvason (1996) for dendritic

solidification; Agresal et al. (1998) for deformation and adhesion of circulating cells;

Esmaeeli and Tryggvason (1998, 1999) for simulation of bubbly flows; Francois and

Shyy (2002, 2003) for micro-scale drop dynamics.

The overall interfacial dynamics treatment in the presence of phase change is of a

hybrid nature i.e., the mass and momentum equations are solved using the immersed

boundary method, but the temperature equation is solved in a sharp fashion. The primary

intention, as pointed out in previous chapter, is to accomplish higher accuracy in

interfacial mass transfer computation. Due to the absence of a surface tension type source

term in the temperature equation, sharp treatment of the temperature equation can be

performed with relative ease as compared to a sharp treatment of the surface tension

forces, pressure, etc. Since the interface has been assumed to be at saturation temperature

at all times, it simplifies the process by allowing the interface temperature to be used as a

boundary condition in the numerical discretization of the thermal diffusion term in

temperature equation. This discretization approach is similar to the recent developments

by Gibou et al. (2003) and Morgan (2005) using ghost fluid based sharp interface for

phase change computations.

2.3 Adaptive Grid Computation

Most of the computational time is devoted to performing computation/operations

on the interface and solution of the flow-governing equations. Usually the latter part is

responsible for majority of the computational burden. Since the governing equations of










the flow are incompressible, the procedure involves solution of an elliptic equation to

compute the pressure field making it one of the primary reasons for the high

computational cost. Since the system has multiple fluid components, the jump in the fluid

properties, in general, has an adverse bearing on the overall computation-time. Figure 2-

10 from the work of Francois et al. (2004) is a representative case of degrading


performance with increasing fluid property jump. For brief information, these

computations were performed for a buoyancy-driven rising bubble using the immersed

boundary method. The Reynolds number, Weber number and Froude number were set to

100, 4 and 1.0. A W-cycle multigrid technique was used to solve the pressure Poisson

equation. The horizontal axis in Figure 2-10 shows the level of multigrid used for

computation where level = 1 represents computation without employing any multigrid

acceleration.

Performance (10/42 = 10) Performance (p1/p2 = 100)
1 B-D----------------------I 0.4 E---!-----------------!---------------------
1 0.4
X0 0.35 -
1 "-""'-- Il'O2|0.7 0.3 'l/2-=1 --
0.5 0.

.go- 0.3p/--- a -

(a0.2 (b







S Number of Levels Number of Levels
(a) (b)




bubble. (a) Effect of density ratio; (b) effect of viscosity ratio (Francois et al.,
2004). [Reprinted with permission]
1 2 3 1 2 3






2004). [Reprinted with permission]









The performance seems to be affected the most by the density ratio factor: without

multigrid, the computation time for density ratio of 1000 is nearly ten times higher in

Figure 2-10 than for a density ratio of 10. Although these observations come from a

specific method and problem, they still represent the general trend (Francois et al., 2004;

Lorstad & Fuchs, 2004) and recommend appropriate measures to reduce the

computational cost as much as possible. The presented work uses adaptive Cartesian

grids to resolve the flow features locally and thus attempts to reduce the overall

computational data-size.

2.3.1 Cartesian Grid Data Structures

The approach for generating adaptive Cartesian grids is divided into two categories:

a cell-by-cell adaptation (Figure 2-11(a)); a block-by-block adaptation (Figure 2-11(b)).

Adaptive mesh refinement known as AMR is an example of a block-by-block adaptation

(Berger & Oliger, 1984; Sussman, 1999, 2005) that interprets the grid as a union of

uniformly refined grid-blocks (rectangles in 2D and hexagons in 3D). Due to the

approach of collecting cells in blocks and refining them, AMR tends to produce a larger

grid data-size than necessary. As observed by Aftosmis (1997), typically only 70% of the

cells in a rectangular patch may actually need refinement. In the sense of reducing the

computational data-size, cell-by-cell adaptation is a more flexible approach.

A cell-by-cell adaptation uses either a tree-based data (Wang, 1998; Popinet, 2003;

Losasso et al., 2004) or an unstructured data format (Aftosmis, 1997; Ham et al., 2002;

Singh & Shyy, 2005, 2006) to store the grid. A tree-based data-structure (Figure 2-12)

offers a compact storage system. All the connectivity information is implicitly contained

in the location of cells in the tree-structure. One of the most common queries in a









numerical simulation is the cell-cell connectivity information needed, for example, to

compute convection/viscous fluxes. Since no explicit information about the identities of

the neighbor cells is stored, it needs to be extracted from the information contained in the

tree. As reported by Aftsosmis (1997), approximately 10-20% of the CPU time for CFD

solvers is dedicated to tree traversals answering grid-connectivity queries.


(a) (b)

Figure 2-11. A cell-by-cell and block-by-block grid adaptation. (a) A cell-by-cell
adaptation in the present work; (b) schematic of a block-by-block adaptation.

1 2 ---------------------- Level=0
1
3 4

1 -- 2 -- 4------- Level=1
S 25
76
3 4 2 5 7Level=2

(a) (b)

Figure 2-12. A tree-based adaptive Cartesian grid-data. (a) A cell undergoing two levels
of successive refinement; (b) data structure showing the place of cells based
on the refinement level.









2.3.2 Present Approach: Unstructured Adaptive Cartesian Grid

An anisotropically adaptive grid (Ham et al., 2002; Singh et al., 2005) offers an

optimized grid data by using a direction-based refinement strategy. However, the

resulting data-structure and solution algorithms become difficult and tedious to

implement. For simplicity and convenience, an isotropic adaptation is employed where a

cell marked for refinement is simultaneously split in all coordinate directions. The flow-

governing equations and numerical discretization aspects have been discussed in the

related chapters.

The grid is stored using unstructured data with cell-by-cell refinement approach. It

is believed that the extra memory requirements due to explicit storage of connectivity

information will be offset in part by compact storage system exploiting the Cartesian

nature of the grid and easier accessibility of the data for the flow solver.














CHAPTER 3
MULTIDIMENSIONAL IMMERSED BOUNDARY COMPUTATION

The immersed boundary method originally proposed by Peskin (1997) has been

used to model the surface tension force as a smoothed volume source term in the

momentum equation. The jump in density and viscosity are smoothed within few cells

across the interface. The isothermal immersed boundary computation is similar to the

work of Tryggvason et al. (2001). This chapter presents the key components of the

immersed boundary modeling and interface tracking including algorithms for maintaining

interface resolution and technique for handling topology changes. The numerical solution

technique including adaptive grid generation and equation discretization are presented in

the following chapter.

3.1 Immersed Boundary Method

The incompressible Navier-Stokes equations for mass and momentum and energy

conservation represented by Equations(3.1), (3.2) and (3.3) are solved on the stationary

background grid. The source terms F, models the surface tension as a smoothed

momentum source term using a smoothed Dirac delta function (Figure 3-1(a)).. The fluid

properties are smoothed using an indicator function I that varies smoothly from a value of

unity outside the interface to zero inside (Figure 3-1(b)). The density and viscosity are

smoothed using Equation (3.4) (Prosperetti, 2002) where subscripts 1 and 2 refer to the

fluid outside and inside the interface, respectively.

The source term Q, in the energy equation accounts for the heat source due to phase

change. Since the interface is considered at saturation temperature, the job of this term in








the immersed boundary method of Juric and Tryggvason (1998) is essentially to maintain

the interface temperature. Despite smeared treatment of the hydrodynamics equations, the

temperature equation in the presented work is solved in a sharp fashion to compute the

mass-transfer accurately. Since the temperature equation treatment is different than the

traditional immersed boundary method, role of the source term Q, and corresponding

numerical solution procedure will be visited in the following chapter.

Op
S-+V. pU =0 (3.1)
at
dU
p +V*UU V=-VP+V. lV(U+V'u)+F,+pg (3.2)


pC C +v.(TU) = V.(KVT)+Q (3.3)


P = P2 +(1 -2)

p pJ +(p p (3.4)


Bound of numerical delta function

Region of influence

Actual interface = 1


1=0




(a) (b)

Figure 3-1. Schematic of the immersed boundary method. (a) Interfacial effects smoothed
over a thin zone across interface; (b) indicator function for fluid property
smoothing.









3.1.1 Interfacial Conditions

The mass conservation equation across the interface is written as Equation(3.5).

The jump in pressure across the interface is related to the surface tension force, jump in

shear stress and fluid velocities using Equation(3.6). Using the mass continuity equation,

the jump in the normal component of the velocity can be written using Equation(3.7). The

tangential component of the velocity is continuous across the interface and the normal

component is also continuous if there is no mass transfer.

P (U, Umlt)- n = (.
(3.5)
P2(U2 Ut ) n =

P2 ((Ul)2 -(U)nt)-U2+P, --2 n= p, ((U )1 -(Um)-U1 +P,-n- n+o nc

(3.6)


(U, U,2)n= 1 --- (3.7)


The mass transfer rate across the interface is computed using Equation(3.8) for

which a probe based technique is used to compute the required temperature gradients on

the interface. These gradients are computed using Equation (3.9) and Figure 3-2 showing

two probes of length 1.2A along the interface normal direction. The term A is the

background grid size and Tint is interface temperature (set to saturation temperature Tsat).

The temperature T1 and T2 on the probes are computed using bilinear interpolation from

the background.

S= (-KVT n) (-KVT n)2 Ail (3.8)

T Tn
(KVT-n = 1 mt
1.2A (3.9)
(KV T n)2 =nt T
1.2A









Using the mass continuity across the interface, the interface velocity can be written

as:

glnt.H z Il- + U2 i -1 21
U n =( U n ) (3.10)
2 2 (3.10)

The interface velocity due to mass-transfer component is computed with the help of

Equation (3.8) and unsmoothed fluid properties. The interface velocity component due to

fluid velocity is simply interpolated from the background grid using a Dirac delta

function described later.

The velocities in individual phases are considered incompressible and the jump in

velocity is smoothed using Equation(3.11) (Juric & Tryggvason, 1998). The effect of

mass transfer is modeled using locally non-divergence-free velocities near the interface.

The Equation (3.12) relates the divergence of the velocity field with the interfacial mass-

transfer after taking divergence of Equation(3.11) and applying incompressibility

condition for velocities in the individual phases where U1 and U2 are fluid velocities

outside and inside the interface.

U= U +(U U2)I (3.11)

1 (3.1
V.U=(U -U2).VI= -- -V*I (3.12)
SPi P2)


T, 0 Liquid phase


Figure 3-2. Probe based technique for temperature gradient computation on the interface.










3.1.2 Material Property Smoothing

Figure 3-3 shows a common finite-element type data-structure used to store the

triangulated interface. It assigns a unique identification (integer number) to each node

(marker) and stores its coordinates along with the triangle-nodes. All the triangles are

assigned a unique orientation such that the right hand thumb rule using the ordering of

nodes of a triangle gives the local normal-vector pointing outside the interface.

Node ID x y z
1
2

nnode

Triangle ID nodel node2 node3
1
2

ntria

Figure 3-3. A triangulated interface and the data-structure.

The indicator function is computed at Cartesian cell centers and linearly

interpolated to cell faces whenever required. It is computed by solving the Poisson

Equation(3.13) where 6(x-xint) is a Dirac-delta function non-zero only at X=Xint and D is its

discretized form for which Equation (3.14) (Peskin & McQueen, 1996; Griffith & Peskin,

2005) is used.

Since the fluid properties are constant in most part of the domain, the indicator

function is solved only within few cells (4-5) across the interface. For rest of the cells,

the indicator value is set manually. The cells completely outside the interface (Fluid 1)

are assigned a material flag MatFlag(cell) = 1 and cells completely inside the interface

(Fluid 2) are set to MatFlag(cell) = 2 (Figure 3-4). This knowledge is used to manually

set the indicator function to zero or one and use it as a Dirichlet boundary condition for

Equation(3.13).










V V. n-,)dA& VK nD(r)dA (3.13)
interface ) tangle ,


D(r) =
A A A

1(33-2h+ 1+4h(1-h) :ifh e[0,l)
8 (3.14)

h= 1=<(5-2h- -7+4h(3-h)) :ifhe(1,2]

0 :else


A ray-tracing algorithm popular in computer-graphics and computational-geometry

community could be employed for accurate evaluation of the state of a cell i.e., inside or

outside the interface. However, since such accurate information is required only outside

the smearing zone, a simple and efficient method based on the painter's algorithm

frequently used in computer-graphics rendering is used. Unlike the ray-tracing algorithm,

the painter's algorithm does not require expensive computation of three-dimensional line-

surface intersection and it is sufficient for the current purposes to deal with simply

connected closed topologies. Open surfaces can also be dealt with ease following suitable

algorithmic considerations. This algorithm stores the state of a cell in an array named

MatFlag(cell). This approach exploits the fact that a cell outside the interface has all its

neighbors either outside or on the interface. The outline of the implemented algorithm is

presented below:

* Set MatFlag(cell) = 2 for all the cells contains a marker else set it to zero.
* Select any cell outside all interfaces, call it SeedCell, and put it as the first entry in
a temporary array named SeedCellList().
* Set TotalSeedCell = 1; SeedCellList(TotalSeedCell) = SeedCell.
* Set MatFlag(Seedcell) = 1
* Use a temporary variable NextCheck to march through the array SeedCellList().
Set NextCheck = 1









* Do while ( NextCheck < TotalSeedCell)
o SeedCell = SeedCellList(NextCheck)
o Go through the neighbor cells which share a face with SeedCell
If (MatFlag(NeighborCell) > 0)Go to another neighbor cell
TotalSeedCell = TotalSeedCell + 1
SeedCellList(TotalSeedCell) = NeighborCell
MatFlag(NeighborCell) = MatFlag(SeedCell) ( = 1)
o NextCheck = NextCheck + 1
* End Do While Loop
* Set MatFlag(Cell) = 2 for all the cells which still have MatFlag(Cell) = 0

MatFlag(cell)=1
Fluid 1 /= 1.0


Fluid 2

MatFlag(cell)=2
I = 0.0




Figure 3-4. Indicator function and the material flags for cells completely outside or inside
the interface.

3.1.3 Momentum Source Term Computation

The surface tension force is computed on the interface triangles and distributed to

the background grid as smoothed source term Fs in Equation(3.2). Figure 3-5 shows the

smoothed nature of the momentum source for a two-dimensional interface. The surface

tension force on a discretized interface element (curves in 2D and triangles in 3D) can be

evaluated in several ways: computation with Equation (3.15) where unit normal vector

and curvature can be computed using curve fitting for two-dimensional interfaces

(Francois & Shyy, 2002, 2003; Ye et al., 2001) and surface fitting for three-dimensional

interfaces (Sousa et al., 2004); computation using a line integral form shown in Equation

(3.16) and fitting curves/surfaces to obtain normal and tangent vectors (Al-Rawahi &

Tryggvason, 2002; Tryggvason et al., 2001).









f = ckndA (3.15)
5A

f = af kndA= (nx V)x ndA= at x nds (3.16)
3A4 3A

There are two important observations to be made here: the net surface tension force

on a closed surface should be zero (conservation); curvature computation using

interpolation based methods are numerically sensitive and often requires some form of

data smoothing (Francois & Shyy, 2002, 2003; Ye et al., 2001; Sousa et al., 2004). The

use of Equation (3.15) does not enforce conservation whereas the line-integral form does

not require explicit curvature computation and maintains the conservation. The present

work uses the line integral form and computes the local normal and tangent vectors along

the triangle edges using the simple approach of Al-Rawahi and Tryggvason (2004) shown

in Figure 3-6. To ensure conservation (i.e., the net surface tension force on a closed

interface is zero), the cross-product of normal and tangent on a triangle edge is obtained

by averaging its value between the triangles sharing this edge. The surface forces

computed on triangle edges are distributed directly to the background grid using

Equation(3.19). If required, the curvature can be computed using Equation(3.18). Such a

simple technique is seen to produce sufficient accuracy demonstrated by Figure 3-7

comparing curvature of a unit circle using present method and a cubic-spline

interpolation. The overall accuracy of the present surface tension force computation and

modeling has already been demonstrated by Shin and Juric (2002) for boiling flows and

Al-Rawahi and Tryggvason (2004) for dendritic solidification.

sf= 7 -(tfn)edgeds (3.17)
edge=1,2,3










3f~n
K=
cdA

F5 =fedgeD(r)
edge


(3.18)


(3.19)


(b)


Figure 3-5. Momentum source term due to surface tension force on a circular interface.
(a) The x-component, (b) y-component and (c) the x-component along a
horizontal line though the center of circle.

tl1t2
Stl t21


t2 Ip1l-
pl- p3 \p3


Sp3-p2
triangle tl p3 p2I


t3 p2- p
Iv2- 11 p2


triangle => pl, p2, p3


Figure 3-6. Computation of the unit normal and tangent vectors on interface triangles.











1.055
1.05 1.002 -
1.045
1.04 1.0015 -
2 1.035
1.03
S1.025 1.001
1.02
1.015-
1.01 1.0005
Current method
1.005Current method Current method
1 1 -vvvvvvvvvvvvvvvwvvw
0.995 I I I I I I
0 10 20 30 40 50 60
Interface node Interface node

(a) (b)
Figure 3-7. Curvature computation for a unit circle using Equation(3.18) and cubic spline
with (a) twelve and (b) sixty points.

3.1.3 Interface Advection

The interface velocity in Equation (3.10) has a component due to the fluid motion

and another due to phase change. The component due to phase change is obtained by

computing the mass transfer rate on the interface markers, as described earlier. The

component due to fluid velocity is interpolated from the background grid and the final

marker velocity expression is:


Umnt Y llU(cl)D(r)dxdydz-m + P2 n (3.20)
l (c ell) 2pI 2 )

For simplicity, the marker movement due to background fluid velocity is solved

using 2nd Order Runge-Kutta method but the velocity due to phase change is treated using

forward Euler method as shown in Equation(3.21).


X =Xn +AtU,
** *
X =X +AtU(x, (3.21)
n ** 1 __
n+1 Xn+X nat p +p2
2 2X2nAt
2 2p, p2









Error associated with interpolated marker velocity and integration of the advection

equation produces some error in the interface volume. These errors are usually small but

can accumulate over time and become unacceptably large. To correct such defects, an

explicit phase-volume correction step is performed by perturbing the markers in local

normal direction and applying bisection method to find a location that keeps the volume

error within some prescribed level (-10-7%). This correction may be performed either at

every time-step or periodically. All the presented computations performed this operation

at every time step and usually one or two bisection iterations were found to be sufficient.

For flows without phase change, the interface volume does not change and hence the

amount of volume error and hence the correction is always known. For phase change

computations, the amount of volume error and thus the correction was estimated based on

the mass-transfer rate.

3.2 Conservative Restructuring for Interface Resolution

Due to the limited support of the Dirac delta function used for distribution of

surface forces and other interpolation operations, excessively long edges (e.g. more than

two Cartesian-cell-size) will appear as a 'hole' on the interface. On the other hand,

excessively small edges can produce unphysical sub-grid undulations. To ameliorate such

effects, markers are continuously added and deleted based on the following criteria:

Edge_ length > Cell size > Break edge (marker addition)
Cell size (3.22)
Edge _length < >- Delete edge (marker deletion)
3

It is required that the marker addition and deletion causes minimal disturbances on

the interface and avoids introducing volume errors. However, the commonly employed

techniques in literature (e.g. Figure 3-8) do not preserve the interface volume. More









precisely, it is the marker deletion that is non-conservative. Such non-conservative

behaviors are demonstrated by Sousa et al. (2004) who encountered up 5% volume errors

due to interface restructuring.






(a)





Pi (b) -----------


(b)

Figure 3-8. Interface restructuring approach. (a) Long edges are split at the mid-point and
(b) a small edge pp2 is collapsed to its midpointp3.

The presently employed restructuring approach is same as in Figure 3-8 with an

added correction step to the marker deletion procedure in order to preserve the interface

volume. This correction-algorithm is an extension of the volume preserving marker

relocation of Sousa et al. (2004) who used such an approach for interface smoothing and

removal of sub-grid undulations on the interface. The underlying procedure is explained

using a 2D example shown in Figure 3-9 that shows a marker removal by collapsing the

edge formed by point p2 and p4 to a point p3 resulting in a net phase-area loss. To correct

this error, first a reference pointpref atp3-n is selected (Figure 3-9(b)) where n is the local

unit normal vector at p3. A reference phase area of the polygon pipzp4P5Pref is computed

and point p3 is relocated after deleting the edge to recover preserve the reference area

(Figure 3-9(c)). The conservative marker/edge deletion on 3D surfaces is similar to its 2D

counterpart and it is summarized below:









1. Select a point p at the midpoint of the edge to be deleted. A local normal vector n
at this point is taken as the average of the normal vector on the triangles connected
to this edge.
2. Define a reference point pref as p-n and compute a reference volume by adding
the tetrahedral-volumes made by the reference point and the triangles in Figure 3-
10(a).
3. Collapse the edge to its midpoint and compute a volume Q by adding the volumes
of tetrahedrons made by point with the triangles in Figure 3-10(b).
4. Relocate the point p to p + n (Q0/Q) for preserving the local phase volume Q.


Figure 3-9. Conservative marker deletion for two-dimensional interfaces. (a) Old
interface pip2P4P5 and new interface pip3p5 with non-conservative method; (b)
define a reference point (pref) and compute a reference area; (c) relocate p3 in
local normal direction to preserve the reference area.


(a) (b)
Figure 3-10. Conservative marker deletion for three-dimensional interfaces. (a) A
reference pointpref is defined atp-n (p is edge mid-point) and compute a
reference, (b) collapse edge to point and move p along normal vector to
preserve the reference volume.









3.3 Interface Reconstruction for Topology Change

When two interfaces come within 'close proximity' they are merged together and

similarly when sections of the same interface are in 'close proximity' the interface is

broken. The term 'close proximity' at this instance is intended to be around one grid-cell

size and will be addressed further in the document. Since the indicator function reflects

the interface shape, it is used to redefine the interface to achieve topology changes.

Figure 3-11 shows a 2D interface and indicator contours along with indicator = 0.5

contour exhibiting interface merger. The work of Shin and Juric (2002) exploited this

idea to reconstruct the interface by deleting the interface data completely and

reconstructing a new interface at indicator = 0.5 contour. Their work tracked the interface

using a connectivity free method and needed to perform reconstruction to control

interface resolution as well as facilitate topology changes. Such a situation does not arise

in the current work as the resolution is controlled via conservative restructuring while

reconstruction is performed only to achieve topology change.



















in close proximity and corresponding (c) indicator = 0.5 level contour.









3.3.1 Reconstruction Criteria

The need for reconstruction is checked at prescribed intervals and interfaces are

reconstructed based on some indication of a possible topological change. A simple

criterion shown in Figure 3-12 sends two probes along the local normal vector from the

marker points. The length of each probe is set to one and half times the background grid-

cell-size. The indicator function value on the probes is interpolated from the background

grid using bilinear interpolation. If both the probes have indicator function value greater

than half or less than half, representing merger or breakup for the given grid resolution,

reconstruction is carried out. The need for reconstruction is checked at frequent intervals

depending on the problem. When uncertain, it was checked at every few (e.g. 50) time

steps in the presented computations.



Ind1




Ind2




Figure 3-12. Interface reconstruction criterion uses two probes from the interface along
the normal vector and bi-linearly interpolates the indicator values Indl and
Ind2.

3.3.2 Reconstruction Procedure and Difficulties

The approach for reconstruction is to first create a set of globally numbered

markers, one on each edge with indicator = 0.5 contour intersecting it. Cells with new

markers are visited one at a time and markers on the cell-edges are connected to define

interface segments within the Cartesian cells (Figure 3-13).









The following sections present the overall reconstruction procedure using a two-

dimensional interface and the ideas are extended to three-dimensional interfaces. Some

critical difficulties in maintaining the integrity of the connectivity data are highlighted

and simple remedies are provided.


p3
2 e2 e3
e3
/l el
PZ 1 ___ _


el => (pl,p2)
e2 => (p2, p3)
e3 => (p3, p4)


Figure 3-13. Two dimensional reconstructed interface edges (el, e2 and e3) and
corresponding connectivity data.

3.3.2.1 Two-dimensional interface

The indicator function values computed at cell centers are linearly interpolated to

Cartesian cell-vertices. The vertices are flagged and a rule for marker creation on an edge

is set as:

flag : indicator > 0.5 (outside interface)
* vertex{flag =
0 : otherwise (inside/on the interface)
* Marker is created on an edge only if its two vertices have different flag values.

The above rules avoid ambiguities and simplify the algorithm by eliminating the

need to look into markers in neighboring cells during reconstruction in a given cell. As an

example, consider Figure 3-14(a) where an edge e2 is constructed only while visiting cell1

and not while in cell2 even though the corresponding markers physically lie on the edges

of both these cells but are considered only on the vertical edges of cell due to their

opposite vertex flag values. This procedure, however, allows creation of edges with zero

lengths in order to create a valid connectivity data. Such a situation arises when, for









example, the reconstructed markers pl, p2 and p3 in Figure 3-14(b) all same physical

location i.e., indicator = 0.5 contour is at the Cartesian-cell vertex. No special attention is

paid to such cases as a valid connectivity data is still produced while such edges/elements

are later deleted by the marker/edge deletion procedure presented earlier.

-1 1 1 11


cell,

Q P3 P2
1 1~ 0 1
e6 e> e2 e,
el
Pi
cell3 cell2 cell4 cell1 cell2

0- 0 0 0 1
(a) (b)

Figure 3-14. Illustration of 2D cell-by-cell interface reconstruction. (a) Markers on edges
with opposite vertex-flag values (0 and 1) are connected to form edges el, e2,
and e3 in cells cell3, cell1 and cell4, respectively; (b) markers p, p2 and p3 are
connected in cell, and cell2 irrespective of their physical coordinates.

A 2D cell can have two, three or four markers. Cells with two markers produces

single edge; cell with three markers produce two edges connected within the cell; cell

with four markers are used to produce two edges non-connected inside the cell. It is last

case that causes some difficulty in deciding the orientation of the two disconnected edges.

These four markers are not used to create a closed interface within the cell as it is not

required since the interface will be closed in some other cell easily deduced by using the

vertex flag information.

It can be seen that if the interface was to be reconstructed in triangulated cells, such

ambiguities with four markers do not arise as a triangular cell can have only two or no

markers producing at the most single edge within the triangle. Although, the markers are









created only on Cartesian cell edges, a set of locally triangulated cells are used to

decipher the connectivity/orientation of the edges whenever needed as shown in Figure 3-

15 where the one vertex of the triangles in the cell is at the Cartesian cell-center. The

indicator function is known at the cell center and hence a vertex flag value can be

assigned.



0 1 0
-- ~-~ -^ ---- 0 --1 --





1 0 0



Figure 3-15. A 2D cell with four markers creating two locally disconnected edges. The
connectivity of markers is decided based on the vertex flags of the temporarily
triangulated cell.

3.3.2.2 Three-dimensional interface

The basic approach for three-dimensional interfaces is similar to the two-

dimensional counterpart i.e., create globally numbered markers on cell edges and connect

within one cell at a time. Figure 3-16 shows a schematic reconstruction by creating

markers on edges that have opposite vertex flags and connecting them to create polygons.

The reconstructed polygons are broken into triangles by inserting a marker at the

geometric center of the polygon. Figure 3-17 shows some typical reconstructed polygons

with corresponding vertex-flags. The last two cases in Figure 3-17 contain six markers

each but one has a single six-sided polygon while the other has two triangles. These two

cases with six markers are distinguished from each other by inspecting the relative

location of edges containing the markers.










0 1 ,4
P3 P c5 P3 P5
1 1
P2 P2




0 0
(a) (b)

Figure 3-16. A reconstructed interface segment in a cell. (a) Reconstructed markers are
connected to produce a polygon (pl p2 p3 p4 p5); (b) polygon is broken into
triangles with a new marker at the geometric center.

0 0 1 1 0
0 0


0 0 0 0 ---------- -
S00
0 0 0 0 0 0/
3 markers 4 markers 5 markers




0--- -.----


o 1 1
6 markers
Figure 3-17. Some common reconstructed interface segments and integer vertex-flags.

Since the indicator function itself is smooth in nature, polygons shown in Figure 3-

17 are observed the most often. Similar to 2D reconstruction, occasionally two or more

disconnected polygons within a cell are observed. Presence of more than six markers in a

3D cell presents an obvious algorithmic difficulty in terms of how to define the surface

segments using the cloud of markers. As shown in Figure 3-18(a) eight markers could be

connected in at least two ways to define two differently oriented and connected polygons.









If the reconstruction were to be conducted in tetrahedrons, only three or four markers and

hence just one polygon is present (Figure 3-18(b)). Similar to 2D case, whenever required

the cells are locally broken into tetrahedrons (Figure 3-18(c)) by connecting the cell-

center to cell-vertices in order to decide the connectivity. Although only the markers

created on Cartesian edges are used to construct the interface, possible presence of

markers on tetrahedron edges (based on vertex flags for tetrahedrons) are used simply as

a means to decide how to connect marker on Cartesian cell edges in case there are more

than six of them in a cell.







l~ ------ - - -
1 00

0 1 0
(a) (b) (c)

Figure 3-18. Difficulties due to degeneracy of reconstructed interface data and
implemented remedy. (a) A cell containing eight markers and one of several
ways to define reconstructed polygons; (b) a tetrahedron can have only one
polygon; (c) some edges of tetrahedrons created in a Cartesian cell.

The newly created triangles are oriented using the gradient of the indicator function

(Figure 3-19(a)) to point the normal vector in the direction of increasing indicator value.

Triangles with zero or arbitrarily small area can not be reliably oriented and may

compromise the integrity of the interface data. The orientation of all such triangles is

recursively corrected by using orientation of the neighboring triangle (Figure 3-19(b)).

The triangles on the entire interface can be correctly oriented starting with proper







47


orientation of any single triangle. Figure 3-20 shows examples of reconstruction creating

a rupture and merger of two interfaces.

VI
nI
r



P
L--r-:: ---q-- -- ----
q
q

tria => p, q, r nghtria => r, q, s
(a) (b)

Figure 3-19. Orientation setting of reconstructed interface. (a) Orientation of triangles is
set in the direction of increasing indicator function value; (b) orientation of a
triangle (tria) is used to set the orientation of its neighbor (nghtria).


Before reconstruction
An


A 3D View cross-section
(a)


. . .


3D reconstructed


Figure 3-20. Examples of 3D interface reconstruction. (a) A surface undergoing rupture,
(b) two spheres merging.









3.4 Summary

This chapter presented the immersed boundary formulation and other

computational components including the interface tracking and data-management. The

key aspects of the presented method are summarized below:

* Immersed boundary method
o Immersed boundary flow equations and interfacial conditions.
o Fluid property smoothing using an indicator function.
o Surface force computation and modeling.

* Interface tracking
o A conservative algorithm for interface restructuring.
o A level contour based reconstruction algorithm to handle topology changes
highlighting several difficulties and strategies employed to address them.

The interface restructuring and reconstruction form the key contributions of the

present work while using the well established immersed boundary method to model the

interfacial phenomenon except the phase change. The temperature equation is treated

sharply and the solution algorithms are presented in the following chapter.

The interface reconstruction introduces some amount of error in the interface-

volume which is corrected explicitly by perturbing the markers in the local normal

direction. The magnitudes of volume errors and required correction are estimated in

chapter 5 along with the effect of reconstruction on the accuracy of flow computations

via measurement of spurious velocities for static bubble computations. The effect of

conservative restructuring is also demonstrated in chapter 5 via a spherical interface

placed in a time-reversed vortex field.

The marker tracking by itself is not restricted to any particular choice of interfacial

dynamics modeling; it may be used for sharp interface modeling such as ghost fluid

method. Such a sharp interface modeling avoids smearing of fluid properties and reduces









the spurious currents. So far, the only known work using marker based tracking with

sharp interface ghost fluid method is by Hao and Prosperetti (2004). They used indicator

function to compute the geometric properties of the interface for relatively simple

computations not exhibiting topology changes. The indicator function was used to assign

the fluid properties sharply as done for the temperature equation described in the

following chapter. For such a sharp interface method, several stability and accuracy

related practical issues caused by sub-grid undulations, interface reconstruction and

errors in curvature computation need to be examined closely.














CHAPTER 4
ADAPTIVE CARTESIAN GRID BASED COMPUTATION

The dynamic grid adaptation is based on the interface location and some feedback

from the flow solution. The grid is stored using an unstructured data format and flow

equations are solved using projection method. This chapter presents the grid

generation/adaptation method along with the algorithms for Navier-Stokes computation.

4.1 Data Structure

Figure 4-1 shows a face-based data-structure used for storage of an isotropically

refined grid. Since detailed information may be located in the work of Aftosmis (1997),

Ham et al. (2002), Singh et al. (2005), Singh and Shyy (2006), only a brief description of

the data-structure is represented here. The variable level in Figure 4-1 stores the

information about how many times a cell has been split from an initial state. The integers

i, j and k, together with the level information provide the size as well as the physical

location of a cell. Also, all the faces forming a cell are stored to provide cell-to-face

connectivity information. For every face, a single byte integer named orientation gives

the orientation of the face: faces normal to x-axis have orientation = 1, faces normal to y-

axis have orientation = 2 and faces normal to z-axis have orientation = 3. The variable

sideCelll is identity of a cell sharing this face in the direction of orientation and

sideCell2 is the cell on the other side of this face.

The size (dx, dy, dz) and center (xc, yc, zc) of a 3D cell with coordinates (level, i, j,

k) are computed using Equation(4.1). Other relevant information such as cell volume and

face area can be computed and stored in a very compact way. Since the volume of a cell









is depends on its level, it can be pre-computed and stored for up to the maximum allowed

level of refinement. Similarly the area of a face is a function of only its orientation and

the smallest of level of two side cells and hence it can be pre-computed and stored as

well.

(dx, dy, dz) = L 1 1
Nx 2zevel Ny 2evel Nz 2zee (4.1)
(xc, yc, zc) = ((i 0.5)dx, (j 0.5)dy, (j O.5)dz)

sideCell2
S sideCelll
t


Integer:: level
Cell Data Integer:: i, j, k
Integer:: cell-face list
Face Data Integer:: orientation
FInteger:: sideCelll, sideCell2


(a) (b)

Figure 4-1. Cartesian grid data-structure. (a) Cell and face data; (b) the convention for
two side cells for a face with orientation along x, y or z coordinates.

The grid generation process starts by uniformly dividing the entire domain into

prescribed number of cells in each coordinate direction. All the cells in this uniform grid

are referred to as base cells upon which adaptation is carried out. These base cells are the

largest cells possible in the computational domain and have their level coordinate set to

zero. To briefly explain the process, consider a 2-D domain in Figure 4-2(a) uniformly

split into four parts in x-direction and three parts in y-direction. Let the number of

uniform cells in x and y directions be Nx and Ny, respectively. In two-dimensions, each

cell has three integer coordinates (level, i, j). After the initial uniform splitting of the

computational domain, all the cells are marked as ground level cells i.e., level = 0. The

adaptation is carried out by splitting cells wherever desired. The final grid is obtained by


orientation

face










recursively splitting the cells to desired level of resolution. When a cell with coordinates

(level, i, j) in Figure 4-2 (b) is isotropically split, coordinates of the new cells shown in

Figure 4-2(c) are computed as:

level'= level +1

(i, ') (2i -1, 2j -1) (4.2)
(i ", j ") = (2i, 2j)

SLength=Lx, Nx=4



0-
II








(a)









(level, i, j) (level', i', j'evel', i",
( level', i", j"
(level', i', j" )
(b) (c)

Figure 4-2. Adaptive grid refinement. (a) A 4x4 base grid (level= 0), (b) coordinates of a
two-dimensional cell before splitting and (c) after splitting.

4.2 Grid Adaptation

Interface location and shape provides the geometry based adaptation criteria. In

addition to geometry, the flow-solution is also used as an indicator for adaptation. The

refinement biased algorithm for isotropic adaptation is similar to the work of Ham and et

al. (2002) i.e., when a cell flagged for refinement is to be broken; any cell blocking the









refinement is also broken. This makes the adaptation a refinement biased recursive

process. The coarsening however is not a recursive process; a cell is coarsened only if it

is permissible by the constraints due quality considerations and the data-structure.

Coarsening is achieved by combining four cells in 2D or eight cells in 3D only if all these

candidate cells are marked for coarsening. Following three restrictions for quality control

were imposed on the adaptation process:

1. Up to a prescribed number of cell layers (-10) across the interface are always at
maximum resolution to push errors due to non-uniform cells away from the
interface. Uniform grid around the interface also allows simple finite difference
solution of the Poisson equation for indicator function in chapter 3.

2. Cells sharing a face are not allowed to differ by more than one level of refinement
for smooth variation in grid size (Figure 4-3)

3. Corner cells are not allowed to differ by more than one level of refinement (Figure
4-3). This criterion is not a requirement when using collocated grid based
computation (Singh & Shyy, 2006). For the implemented staggered grid algorithm,
this is an essential requirement to facilitate a convenient programming
implementation.







Splitting not allowed




Figure 4-3. Cartesian cells sharing a face and the corner cells are not allowed to differ by
more than one level of refinement.

4.2.1 Interface Geometry based Adaptation

All the cells cut by the interface are flagged for refinement. To facilitate a smooth

transition in the cell size, up to five layers of cells around a flagged cell are also refined.

This process of flagging and splitting cells is carried out recursively until the prescribed









grid resolution is obtained. The Figure 4-4 shows the process of geometry based

adaptation starting from a base grid and refining up to six levels. For demonstration

purposes, only two additional layers of cells around an interface-cell were refined.












(a) (b)












(c) (d)

Figure 4-4. Geometry based adaptation. (a) Initial 5x5 base grid and grid with maximum
refinement level of (b) one, (c) three and (d) six.

4.2.2 Flow Solution based Adaptation

Cells away from interface (outside -10 cells across the interface) are adapted based

on the flow solution. The present implementation uses a curl based adaptation criteria

(Wang, 1998) that computes a parameter ; for each cell as shown in Equation(4.3). The

length scale / is estimated as the cubic root of cell-volume. The decision to refine or









coarsen a cell is made by comparing ;cel to the standard deviation (Equation(4.3)) using

the criteria in Equation(4.4).

=cell V&u|3=/2
1 (4.3)
Ncell =, cel

cel > y' -> Refine cell
(4.4)
celz < 0. c' > Coarsen cell

During the adaptation procedure, the Cartesian cell center values such as pressure,

temperature and face normal velocities need to be reconstructed for the newly created

cells and faces. The cell-centered values are reconstructed linearly while the face-normal

velocities are reconstructed using the divergence free reconstruction due to Balsara

(2001). Flow variable reconstruction during cell and face coarsening is performed simply

by averaging of the corresponding cell-centered or face-centered values.

4.3 Immersed Boundary Solution Procedure

The following seven steps present the overall structure of the computational

algorithm:

1. Interpolate marker velocities from background grid and due to phase change.
Compute the mass transfer rate ri on the markers.
Advect markers using 2nd Order Runge-Kutta for interpolated velocity and
forward Euler for the velocity due to mass transfer.

X* =Xn+U,
X =X + AtU

X X +AtU
(x )
n ** ( 'N
Xn + = Xn + X* At + 2p 1
2 21pp2 )

Estimate interface volume change due to mass transfer and explicitly correct
any volume error.
2. Check and restructure interface to maintain resolution.









3. Compute indicator function and smoothed fluid properties.
4. Compute surface tension forces and distribute as momentum source term.
5. Solve Navier-Stokes equation to advance flow field to next time level.
6. Check for geometry and solution based grid adaptation at prescribed intervals.
7. Check for interface reconstruction at prescribed intervals:
Perform reconstruction when required.
Explicitly correct any global phase volume error caused by reconstruction. The
magnitude of volume correction in a reconstructed interface is inversely
proportional to its volume so that lesser correction is applied to smaller entities.
If multiple interfaces are present, the reconstruction is selectively performed
only for the interfaces that request it based on the probe based criteria presented
in the previous section.

The first four steps and the step 7 have been described in the previous chapter. Of

the remaining two steps, the geometry and solution based adaptation criteria have already

been presented in the earlier sections. The following sections present the numerical

discretization and Navier-Stokes solution procedure.

The incompressible mass and momentum equations are solved using a projection

(Chorin, 1968; Kim & Moin, 1985) with staggered grid finite volume formulation. The

temperature equation is solved using collocated grid method. The pressure, temperature

and fluid properties are stored at the cell center and face-normal velocity is stored on

Cartesian cell faces (Figure 4-5(a)). The flow computation follows the following

sequence of steps:

Step 1: Solve for temperature field
The following temperature equation is solved using sharp interface method. The
spatial discretization of the terms is presented later. The fluid properties and
diffusion term are treated sharply. The convection and diffusion terms are
discretized using 2nd order Runge-Kutta and fully implicit method, respectively.

n+l n
Cn 1T T n+1/2 n+1
(pC) l-r" =-v.().uT) 2 +v.(KVT).l


Step 2: Velocity computation using projection method.
Predictor-step
Solve the momentum for an intermediate velocity field U* using Equation(4.5)
where all the know values such as surface tension source, gravitation, convection









and old time-step viscous term due to Crank-Nicholson method are lumped into
S". The temporal discretization of the convection term uses 2nd order Runge-
Kutta method. The pressure term is approximated using the old time pressure
field. Subsequently, remove the effect of pressure term by shifting the velocity
field back to obtain another intermediate velocity field U** using Equation(4.6).

AVp a, U=-VP".dA+f +S" (4.5)
At d4
dA
VP"
U**= U+At-- (4.6)
P

Corrector-step
Correct the predicted velocity field (**) using Equation(4.7). The pressure field
for this correction is computed by enforcing velocity-divergence condition and
solving the Poisson Equation(4.8) with conjugate gradient method. The
divergence of the new velocity field U"+' is zero if no phase change occurs else it
is computed using Equation(4.9) to account for the interfacial mass transfer
described in the previous chapter.

Vpn+l
U" = U** At-- (4.7)
p n+1

en+a 1 ndA U'' ndA V Un" dV2 (4.8)
Cell face Cell face cell


V-U= --- h VI (4.9)
A A P2)

4.3.1 Staggered Grid Arrangement and Approximations

A velocity constraint suggested by Losasso et al. (2004, 2005) is used to simplify

the spatial discretization. With this approach, the velocities and pressure gradients on fine

faces, face and face2 of a coarse cell shown in Figure 4-5(b) are constrained as U1 = U2

and Px, = Px2 i.e., face] and face2 have the same control volume. The face normal

gradients on these fine faces are approximated using Equation(4.10) so that same amount

of velocity correction is applied when using Equation(4.7). It was also shown by Losasso

et al. (2004, 2005) that applying the same pressure gradients on finer faces can produce









second order spatial accuracy for the computed pressure field. Although assigning the

same velocity to fine faces reduces the spatial accuracy near non-uniform cells, it is

accepted as a compromise for algorithmic simplicity especially when non-uniform cells

can be pushed away from the critical regions (e.g. interface) via suitable refinement

criteria.


PI = P2 = (P+ ) 2P (4.10)
2d


P,p,u


< d ->
(a) (b)

Figure 4-5. Staggered grid arrangement. (a) Spatial location of flow/fluid properties; (b)
single control volume used for fine faces of a coarse cell.

4.3.2 Spatial Discretization of Velocity Equation

The convection term of type V (OU) is represented in finite volume form integrated

over the control-volume faces as:

SV (U)dV = -U-.ndA= -U~ dA (4.11)
cv

The normal velocities U, are computed appropriately to maintain the divergence

free condition within each control volume i.e., _UfdA = 0. Using Figure 4-6(a) and the

velocity constraints, the control volume face-normal velocities are computed using face-

area weighted averaging that produces:


U, i

fae 1


Px2 2----U P
I Px2
face2
Il ___I1__







59


Uf = 0.5(U + U,)
Uf2 = 0.5(U + U3)
U,3 = 0.5(U, +U4) (4.12)
S(U5Areal + UArea2)
f4 = Areal + Area2

The momentum flux term q on control volume faces are computed using linear

interpolation of face-normal velocities. To conserve the momentum flux across non-

uniform grid cells, Figure 4-6(b) makes another approximation that computes the flux on

top control-volume face for velocity U1 and equally distributes to the corresponding

control volumes of faces for U2 and U3 at a T-node junction. The flux term q in Equation

(4.11) is computed using 2nd or 3rd order essentially non-oscillatory scheme (ENO) where

the grid is locally of uniform size. Since few layers of cells around interface are always at

uniform resolution, 3rd order ENO scheme is employed there. For non-uniform cells

(away from the interface) where uniform grid cells were not present to construct even a

2nd order ENO scheme, a simple central differencing was used for simplicity.


Ucf4
U5t U6


UI Ul f U2

U4 'I--Ucf3

I, Ucf2
I -U-,' U3
I I


I-- ---------- ----------------


f/2 f/2


I. -. ..
U,


4- Areal <- Area2 -
(a) (b)

Figure 4-6. Spatial discretization of convection term. (a) Control volume (CV) and face-
normal velocities (Ucf); (b) flux fat a T-node junction face (top of CV for U1)
is split equally (f/2) and passed to the control volumes of U2 and U3.









The momentum equation, with the known source terms lumped into (S,,Sv,Sw), can be

rewritten using Equation (4.13). The last term in Equation (4.13) comes from

V. (V u)for incompressible flows. The discretization of V t(Vu) requires

computation of control-volume face-normal velocity gradients and the terms due to

V j (V'u) require computation of mixed derivatives. Discretization of V. (Vu) follows

the strategy outlined by Figure 4-6 while the terms due to V (V'u) are added only

around the interface where the viscosity is spatially varying and it is easily discretized

there since the grid around the interface is uniform. Owing to divergence free nature of

the velocity and constant viscosity in most of the domain, the contribution due to terms of

V. jt(V'u) are neglected in cells away from the interface where the grid may be non-

uniform.

Du P + 8(u,v,w)
p- = --+S, +V-(uVu)+V- a
Dt ax ax
Dv aP (u,v,w)
p- --+ S, + V (Vv)+ V ,u (4.13)
Dt y
Dw OP + (u,v,w)
p- = --+ S,+V-(pVMw)+V -p ---
Dt 8z 8z

4.3.3 Spatial Discretization of Temperature Equation

Since the temperature is located at the cell-center, a collocated grid discretization

with auxiliary variable creation for non-uniform cells is used (Ham et al., 2002; Singh &

Shyy, 2006). Convection terms are discretized using 3rd order ENO (essentially non-

oscillatory) scheme around the interface. A 2nd order ENO scheme is used for cells away

from the interface if the required uniform grid stencil can be easily be found, else a









simple central differencing was used and was not found to produce any instabilities in the

computed results.

4.3.3.1 Sharp interface treatment

The fluid properties (pC, K) are used sharply based on the interface location. The

interface location, for this purpose, is approximated by indicator = 0.5 contour as was

done for interface reconstruction described in the previous chapter. The term V. (KVT)

is computed using central differencing for cells away from the interface. For interfacial

cells, the discretization is based on the sharp interface treatment of Gibou et al. (2002)

and Luo et al. (2005) using level-set method. The discretization in the interfacial cells can

be explained using Figure 4-7 where the diffusion term for cell i is computed using the

interface temperature as a boundary condition (Equation (4.14)).


S T K+1/2 -1/2 T-1



i terface


west i east i+1
Xx- I dx(414

& dx-delta



Figure 4-7. A one-dimensional interface between i and i+1.

4.3.3.2 Interface temperature treatment

Coming back to the temperature equation presented in Chapter 3, the role of the

source term Qs in Equation (3.3) is essentially to maintain the interface at saturation

temperature. The approach of Juric and Tryggvason (1998) treated this source term

implicitly using a Newton iteration technique. The work of Son et al. (1999) and Son









(2001) eliminated this source term by fixing the temperature inside the interface at

saturation temperature and solving only outside the interface. Later on, Morgan (2005)

suggested another approach by linearly correcting the temperature within one cell on each

side of interface. This linear correction approach is summarized as:

Step 1.
Solve temperature equation:
-+V uT =-(V .(KVT)) (4.15)
at pC

Step 2.
Correct temperature in cells immediately across the interface. Consider the
interface lying between cell 'i' and 'i+ 1' in the following figure. After solving the
temperature equation in step 1, the temperature at cell centers 'i' and 'i+ are
modified by linear interpolation of interface temperature (saturation temperature
Tsat) and neighboring cell-center temperatures as shown in Figure 4-8.

Corrected Ti
Ti -1-. /Corrected Ti+1
{T i-1- Tsat /
Ti.
T i+1 T i+2

I I I I
I I I I
Interface location
Interface temperature = Tsat

Figure 4-8. Linear interface temperature correction within one cell on each side of the
interface.

The present approach defines the interface at indicator = 0.5 and defines the

correction as a solution of V (KVT)= 0. This approach does not require explicit

knowledge of the minimum distance of a cell from the interface as required in the method

of Morgan (2005) using level-set tracking. This correction technique is convenient with

the immersed boundary method and can easily correct the temperature within any

prescribed thickness across the interface. This equation is referred to as the interface









temperature correction equation and it is solved sharply within a thickness of 1.2 grid-

cells on each side of the interface. Using Figure 4-9, this method will correct the

temperature of Ti and Ti+l while leaving Ti-1 and Ti+2 untouched. The interface saturation

temperature (Tsat) and temperatures of cells outside the correction zone are taken as

Dirichlet boundary condition for the correction equation. Since K is constant within each

domain, this equation reduces to V (VT)= 0 and it is enforced while solving the

temperature equation.

Interface T as BC


Ti-1 as BC

v.(KVT)=o v.(KVT)= T asBC
T Ti+2 as BC
Ti Ti+1


DD \ ~ I D


to be corrected

Figure 4-9. Presently employed interface temperature correction technique solves a
correction equation sharply on each side of the interface.

4.4 Summary

This chapter presented the adaptive Cartesian grid computation technique

employed in the current work. The discretization of various terms in the flow equation

and corresponding constraints on the grid were presented. The hydrodynamics was solved

using immersed boundary method while the temperature equation was treated sharply for

accurate interfacial mass-transfer computation.

Several improvements in the grid adaptation are possible. For example, a solution

based adaptation requires some measure of local solution error and a more general






64


criterion suitable for both low and high Reynolds number flows based on the truncation

or discretization error (Ham et al., 2002; Ferm & Lotstedt, 2003; Muzaferija & Gosman,

1997) can be considered. Flow solution procedure and numerical discretization also need

further research for in order to improve the spatial accuracy by eliminating restrictions

placed on the velocity and pressure field.














CHAPTER 5
RESULTS AND DISCUSSION

The tests are conducted in three categories: basic tests to establish the working of

adaptive grid flow solver and characterization of the interface handling techniques;

bubble/drop dynamics tests; phase change computations. The conservative interface

restructuring technique and reconstruction are tested for various aspects such as accuracy,

mass conservation. A spherical surface in a time reversed vortex field is used to test the

effectiveness of conservative restructuring. Interface reconstruction is tested for its

accuracy in terms of mass-conservation and solution disturbances it causes. The spurious

velocity currents and effect of reconstruction along with a set of multidimensional single

and multiple rising bubbles and binary drop collision are computed to establish the

accuracy of computational procedure. Computations of diffusion heat transfer between

concentric cylinders, natural convection, one dimensional phase change and a stationary

bubble growth are conducted to gauge the performance of the presently developed phase

change computation procedure.

5.1 Interface Restructuring: Interface in a Vortex Field

A spherical interface placed in a time-reversed vortex field (Figure 5-1) deforms

and should return back to the original spherical shape (Aulisa et al., 2004). The vortex

field with time period T= 4 was imposed on the Cartesian grid using Equation(5.1) and









the markers were advected by interpolating this velocity field from the Cartesian grid.

u(x, y, z) = cos(tt/T) sin2 (7cx) (sin(27cz)- sin(27cy))
v(x, y, z)= cos(7rt/T)sin2 (7y) (sin(2cx)- sin(27cz)) (5.1)
w(x, y, z) = cos(rt/T) sin2 (7cz)(sin(2cy) sin(27rx))

During the simulation period, markers are continuously added and deleted from the

interface to maintain the resolution. The interface at the end of one time period should

return back to the spherical shape with some errors due to restructuring and marker

advection using the interpolated velocity field from the Cartesian grid. The errors at the

end of one time period were estimated for the final surface area, volume and radius of

markers (distance from the expected sphere center) as shown in Equation(5.2). Table 5-1

presents the computed errors showing a better performance by the conservative

restructuring as compared to the non-conservative restructuring presented in Chapter 3.

Figure 5-1 shows the time history of the volume-error in the interface along with the grid

convergence studies. The conservative method performs considerably well showing much

smaller volume errors with better than quadratic convergence. Note that the volume

errors with conservative restructuring are caused solely by the integration of the marker

advection equation using the interpolated velocity field.

RMS_Radius_Error= ( nal-; ni,, )2/N


Area_Error= 100 1- Anl (5.2)

Volume Error = 100 1 Vfi ntal







67


Table 5-1.Error estimates using Equation(5.2) for interface in a time reversed vortex field
test comparing conservative and non-conservative restructuring.
Grid RMS Radius error % Area error % Volume error
(conservative, non- (conservative, non- (conservative, non-
conservative) conservative) conservative)
30x30x30 (0.2435, 1.0254)x10-2 (4.2898, -7.9391) (-0.5591, -14.0950)
60x60x60 (0.0091, 0.4212)x10-2 (1.5306, -2.4060) (-0.1178, -4.9882)
120x120x120 (0.0079, 0.1155)x10-2 (0.4169, -0.5978) (-0.0216, -1.0455)







St=0 t=T/4 t=T/2







t =3T/4 t= T
(a) (b)

Figure 5-1. Interface in a time reversed vortex field. (a) Initial position of the interface
with diameter 0.15 placed at (0.5, 0.75, 0.5) in a unit cube with its center at
(0.5, 0.5, 0.5); computed shape history for vortex field in Equation (5.1) with
period T= 4.


5-- Conservative
o0 101 Non-conservative
Conservative
SNon-conservative

0 10 -
> o E


S10



0 1 2 3 4 40 60 80 100 120
time Grid resolution
(a) (b)

Figure 5-2. Error estimate and grid convergence for interface in time reversed vortex
field. (a) Interface volume error on 30x30x30 grid; (b) final volume error on
30x30x30, 60x60x60 and 120x120x120 grids with dotted quadratic reference
rate.









5.2 Interface Reconstruction: Effect on Mass Conservation

The interface marker spacing and accuracy of the indicator-function computation

and hence the level-contour based reconstruction depends on the background grid

resolution. The aim is to measure interface volume-error caused by reconstruction i.e., the

difference between the original and reconstructed interface volumes. The tests were

performed using a spherical interface of radius r = 0.2 at the center of a domain with

dimensions (0.4, 0.4, 0.4). The interface was reconstructed with different background

grid resolution to confirm a quadratic convergence of the volume error (Figure 5-3)

observed by Shin and Juric (2002).

The volume errors produced by reconstruction are corrected explicitly by

perturbing the markers in local normal direction. It is important to assess the approximate

amount of perturbation a reconstruction may produce. Denoting the initial surface area

and volume of the interface as A and V, an approximate radial perturbation to correct a

volume defect AV may be estimated using Equation(5.3). This shows a second order

convergence behavior for the radial perturbation. To get further insight, this perturbation

should be measured in terms of the grid resolution. If the diameter of the interface is

resolved with N grid cells, the radial perturbation can be rewritten in terms of the grid

cell size A as shown in Equation(5.4). The last column of Table 5-2 shows the behavior

of this quantity and reveals that this perturbation is approximately one tenth of the grid

cell even for the coarsest grid used. This error too exhibits a quadratic convergence with

grid refinement. Considering the fact that fluid properties are smeared over nearly four to

five cells across interface making the immersed boundary method a first order accurate

method, the amount of radial disturbance caused by even moderately resolved grids (N >

20) is reasonably small.







69


Ar AV AV I A r (5
Ar --- (5.3)
A 4nr2 V 3


Ar -AV r AV (5.4)
V 3 V 6

10' -






S100
E 1


10-
I 1 I \
\
\
\ N
10-0


20 40 60 80 100
Number of cells per diameter


Figure 5-3. Grid convergence test for volume error due to reconstruction.

Table 5-2. Volume error for reconstruction of a spherical interface and approximate
radial perturbation required for correction.
Grid cell per-diameter (d/A) % Volume error Radial perturbation (Equation (5.4))
10 -7.5995 0.1266 A
20 -1.9887 0.0663 A
40 -0.4881 0.0325 A
80 -0.1233 0.0164 A


5.3 Single Phase Navier-Stokes Solver

A decaying vortex field and lid driven cavity tests are presented to establish the

accuracy of the hydrodynamic solver. The temperature equation is tested for the

discretization of the diffusion term using the diffusive heat transfer between concentric

cylinders maintained at constant temperatures. The specified temperature on the cylinder

geometries mimics the specified interface temperature and measures the error in the

diffusion term discretization.









5.3.1 Decaying Vortex Field

A vortex field (Kim & Choi, 2000) shown in Equation (5.5) was introduced in a

square domain of 1x1 with free-slip condition on all the walls. The center of the

computational domain is located at the origin of the coordinate system. The computed

decay of the velocity field at Reynolds number Re = 10 with (reference velocity and

length scales taken as 1.0) was compared with the theoretical value at time t = 0.3. Figure

5-4(a) shows the streamlines.

u(x, y, t) = cos(;x) sin(;ry)e -2 t Re
(5.5)
v(x, y, t) = sin(;x) cos(;ry)e -22tRe

The simulations were carried out on a series of adaptive Cartesian grids and

compared with uniform-grid computation with resolutions corresponding to the

maximum allowed resolution on the adaptive grid. The adaptive grids were generated by

recursively refining a 16x16 base grid and the computations were started with prescribed

maximum allowed refinement in the entire domain. The solution based adaptation was

performed at every 20 time step with the computational time step set at 0.001. An

adapted grid with two refinement levels at time t = 0.3 is shown in Figure 5-4(b).

The maximum error in numerically computed u-component of the velocity field

and the corresponding computation-time is presented in Table 5-3 with the time

normalized by the time taken for 16x16 uniform-grid computation. As seen in Figure 5-5,

uniform grid computations show quadratic convergence but adaptive grid shows an

average convergence rate of -1.5. Since it was observed that the maximum error occurs

along in the boundary cells which have non-uniform neighbor cells (Figure 5-4), the

computed error is a reasonable reflection of the accuracy of spatial discretization over







71


non-uniform Cartesian grids. Although the accuracy from uniform grid has been observed

to be superior, it is the computational cost that gives adaptive grid computations an edge.

Since most of the cells shown in adapted grids are of uniform sizes, these computations

may not offer the most attractive picture of computation-time savings that are expected to

be more pronounced for 3D computations.

















(a) (b)

Figure 5-4. Decaying vortex field. (a) Streamlines; (b) two refinement levels on a 16x16
base grid.

10-2


10-

a,
x
1 -


- -- Uniform
A Adaptive
- 2nd order reference


100 200 300
Square root of number of cells

Figure 5-5. Grid convergence test for decaying vortex field showing maximum velocity
error computed on uniform and adaptive grids.








Table 5-3. Max velocity error and computation time normalized with time taken for
16x16 uniform grid. The adapted grids are shown in Figure 5-4.
Max. grid Uniform grid computation Adaptive grid base computation
resolution Max. U error Normalized time cells Max. U error Normalized time
16x16 2.966x10-3 1.000 256 2.966x10-3 1.000
32x32 7.184x10-4 4.440 964 1.105x10-3 1.946
64x64 1.752x10-4 21.602 3616 4.111x10-4 10.725
128x128 4.344x105- 104.810 14548 1.611x10-4 90.64

5.3.2 Lid-driven Cavity

A lid driven cavity flow with Re = 100 was computed where the effect of grid

adaptation is more visible than the decaying vortex field in previous section. Figure 5-6

shows the streamlines of the flow computed on a 120x120 grid. Since the adaptation

criterion is based on curl of the velocity field, computed vorticity contours are presented

in Figure 5-7(a) and they are in accordance with the computations of Ghia et al. (1982).

--











Figure 5-6. Streamlines for lid-driven cavity simulation for Re = 100

The computations were carried out on three grids with maximum resolution of

30x30, 60x60 and 120x120 obtained by one, two and three refinement levels as shown in

Figure 5-7. Figure 5-8 shows the u-velocity profile along a vertical line at the center of

the domain for the uniform and adaptive grid computations showing no significant

difference among the computed results on adaptive and uniform grids. The exact values

of minimum u-velocity in Figure 5-8 are presented in Table 5-4. This table also provides









a comparison of the total number of grid cells and computational time for uniform and

adaptive grids. The second-last column of this table shows ratio of the total number of

cells from adaptive and uniform grid computations. The last column is the ratio of

computational time between adaptive and uniform grids indicating up to -70%

computational time savings on large grids.

Table 5-4. Lid driven cavity computation using uniform and adaptive grids. Comparison
of the number of cells, computation time and minimum u-velocity on the
vertical line through center.
Max. grid Uniform grid Adaptive grid
resolution #Cells Time Min u #Cells Time Min u Cell ratio Time ratio
30x30 900 1.000 -0.203 447 0.587 -0.204 0.497 0.587
60x60 3600 2.933 -0.211 1263 1.157 -0.210 0.351 0.394
120x120 14400 10.910 -0.212 4407 3.264 -0.212 0.306 0.299


(a) (b)


(c) (d)


Figure 5-7. Lid driven cavity. (a) Vorticity contour and final adaptive grids with
maximum refinement level of (b) one, (c) two and (d) three.











04
04
30x30
60x60
02 120x120 02 -
Uniform 120x120
0 C/ -- Adaptive 30x30
0 Adaptive 60x60
Adaptive 120x120
-02 0

-04 -04

-02 0 02 04 06 08 1 -02 0 02 04 06 08 1
U U
(a) (b)

Figure 5-8. The u velocity profile along a vertical line through the center of lid-driven
cavity using (a) uniform grids and (b) adaptive grids.

5.3.3 Diffusive Heat Transfer between Concentric Cylinders

The diffusion Equation(5.6) is solved for steady state solution with Dirichlet

boundary condition of constant temperature on the cylinder walls. The density, heat

capacity and conductivity are set to 1.0. The inner and outer surfaces of the cylinder are

maintained at temperature Ti = 2.0 and To = 1.0, respectively. With the inner and outer

radii denoted by Ri and Ro, the theoretical steady state temperature distribution T(R) in

circular coordinate system is shown in Equation(5.7).

The computations were performed on a 2D Cartesian grid. The geometry of

concentric cylinders (concentric circles in 2D) was modeled using an indicator function

as shown in Figure 5-9. The computations are performed where indicator >= 0.5 (space

between the concentric circles). Figure 5-10 shows the computed normalized temperature

distribution and grid convergence studies showing convergence rate between linear and

quadratic as expected based on the discretization of the diffusion terms near the

boundary.


pCa t= V (KvT) (5.6)












log(R/R,)

log(R,/R)


08

o
0
S06


04


02


0

-0
0.5 ._

a


0 I I I I I I II2 I I I I I


-04 -02 0
Radius

(b)


(5.7)


02 04


Figure 5-9. Indicator contour used for modeling the cylindrical computational domain. (a)
Computations are performed only within the domain with Indicator >= 0.5;
(b) cross-section of (a) showing outer and inner radii as Ro = 0.3 and Ri = 0.1,
respectively.

Inside of inner cylinder (R < Ri)


0005
0004
0003

0002


0.5

ID
QE
a
E
I-


0.001


Outside the outer cylinder (R > Ro)


------ L2 error
--6-- L infinity Error


Linear reference


\\
Quadratic '
reference


100
Grid Resolution
(b)


Figure 5-10. Computed temperature and error for heat transfer in a concentric cylinder.
(a) temperature distribution (T-Ti)/(To-Ti); (b) grid convergence study with
resolution 30 (30x30 grid), 60 (60x60 grid) and 120 (120x120 grid) showing
L2 and Linfinity error norms with reference quadratic and linear lines.


T(R) -T

T-T


150 200


I I I . I . I . I


Ns ''









5.3.4 Natural Convection in a Square Cavity

The temperature equation solver was tested for a natural convection case with

Prandtl number (Pr) = 0.71 and Rayleigh number (Ra) = 1.0x105 using Boussinesq

approximation. The computations were performed in a 2D unit square (L = 1) with top

left and right walls at hot and cold temperatures (difference = AT), respectively. The top

and bottom walls were assigned adiabatic boundary conditions. Results in Table 5-5

show good agreement with the computations by Davis and Jones (1983) for the minimum

and maximum Nusselt numbers on the left wall (Nux=o = -OT/8x with temperature and

lengths non-dimensionalized with AT and L) and vmax as the maximum vertical velocity

on a horizontal line through the center of the domain. Figure 5-11 shows the computed

streamlines and temperature isotherms.


(a)
Figure 5-11. Natural convection with Pr =
the temperature isotherms.

Table 5-5. Natural convection computation
Measured parameter Grid resolution
60x60 120:
umax 68.29 68.5
Nux=o min. 0.721 0.72
Nu .. ima\ 7.952 7.78


(b)
0.71 and Ra 1.0x105. (a) Streamlines and (b)


Davis and Jones (1983)

68.59
0.729
7.717









5.4 Spurious Velocity Currents

Spurious velocities or parasitic currents (Figure 5-12(a)) are unphysical velocity

fields created due to numerical errors causing imbalance of interfacial stresses. These

currents are best seen for static bubble simulations where the theoretical velocity is

identically zero everywhere and the pressure drop across interface balances the surface

tension forces dictated by the Young-Laplace Equation(5.8) (Figure 5-12(b)). There are

several reasons ranging from the accuracy of surface tension computation to errors in

computation of pressure and grid anisotropy effects. A further detailed insight into the

origins and nature of these currents can be found in the work of Torres and Brackbill

(2000) and Harvie et al. (2005). Since the magnitude of these currents for continuous

interface methods are found to be directly proportional to the surface tension and

inversely proportional to the viscosity (Lafaurie, 1994), they have the potential to

destabilize a computation and restrict the accuracy and the range of computational flow-

parameters (fluid properties, bubble sizes) that can be successfully simulated.

The expected magnitude of non-dimensional spurious velocities (Capillary number)

based on the observations of Shin and Juric (2002) and Sousa et al. (2004) is of the order

of 10-4 and it has been hypothesized by Lafaurie (1994) that this magnitude remains

largely unaffected by the Laplace number. To offer a perspective, the volume of fluid

computations of Lafaurie(1994) observed spurious current of the order of 10-2 that limited

their simulations to Laplace number (La) < 10 while the current code was tested up to

Laplace number 12000 and found to be stable. The properties of the fluid outside the

interface and bubble diameter were taken as relevant reference scales.

A spherical bubble of diameter d = 0.2 was placed in a cubic domain of size (3d,

3d, 3d) with slip-wall boundary condition. The surface tension a was set to 0.1 producing







78


a theoretical pressure jump of 2.0. The computations were performed on 25x25x25,

50x50x50 and a 100x100x100 grid. All the time axes in the plots for the spurious velocity

computations are non-dimensionalized with the capillary time scale of d[l/o.

Ap= CK (5.8)

2 -

o u : : : : : :: ,\ p u, fo/ t.t a ju o ..
,1 5


-. ..- /" -
05






A Laplace number = 250 case was performed on 25x25x25 and 50x5x5 grids.
.. i .. .... ......- 03 -02 -01 0 01 02 03
(a) (b)

Figure 5-12. Static bubble computation. (a) Spurious currents of magnitude Ca 10-4, (b)
computed pressure for theoretical jump of 2.0.

5.4.1 Effect of Fluid Property Jump and Grid Resolution

A Laplace number = 250 case was performed on 25x25x25 and 50x50x50 grids.

The magnitude of spurious current reduced by an order of magnitude on finer grid as

shown in Figure 5-13(a) but was found to exhibit no concrete improvements on finer

grids. Since the pressure drop across the interface is captured over a few cells, an error in

computed pressure drop was measured with an error norm in Equation(5.10) used by

Brackibill et al. (1992) where Nin is the number of cells inside the interface. The pressure

drop Api in Equation (5.10) is the difference in the pressure computed at the cell center

and average pressure outside the interface (pou,) as shown in Equation (5.9) for

computing the average pressure drop Apum. The observed convergence rate of -0.3 by

Brackbill et al. (1992) is supported by the computed data in Table 5-6 and Figure 5-13(b).







79


The simulated Laplace numbers in the range of 250 to 12000 were found to have no

appreciable effect on the Capillary number as shown in Table 5-7. High density and

viscosity ratios (1000 and 100 respectively) were also tested and found not to cause very

serious degradation as shown Table 5-8 and Table 5-9.

1 N, 1 Now
Ap,, = n,,, Pont = --- j (5.9)
",n =-1 Nout =l1


Pressure drop error = / (5.10)
N Ap exact


S(25, E,)
0- C
0S --1 25x25x25 C
S50x50x50 o




.- -I (E2/E

I I I I

40 60 80 100, E)
', I, _, I, _, 40 60 80 100
10 20 30 40 Grid resolution
(a) (b)
Figure 5-13. Grid convergence tests for static bubble. (a) Capillary number evolution, (b)
error in pressure drop (Equation(5.10)) using data in Table 5-6.

Table 5-6. Pressure drop for an inviscid fluids with density ratio =1.0. Data was taken
after three non-dimensional time step of 5.0x10-3.
Grid resolution Pressure drop error (Equation (5.10))
25x25x25 0.197
50x50x50 0.151
100x100x100 0.117


Table 5-7. Effect of Laplace number on spurious velocity.
Laplace number (La = opd d /2 ) Capillary number (Ca = Um,,/cr)
250 5.32X10-4
12000 5.35X10-4
12000 5.35x1 0-4









Table 5-8. Effect of density ratio on pressure drop and spurious current on 50x50x50 grid
(viscosity ratio = 1.0, La = 250)
Density ratio = p,/p2 Pressure drop= Pm Capillary number= Umxp/a
/APexact


1 0.90 5.32x10-4
10 0.89 5.27x10-4
100 0.87 5.15x10-4
1000 0.87 6.35x10-4

Table 5-9. Effect of viscosity ratio on pressure drop and spurious current on 50x50x50
grid (density ratio = 1000, La = 250)
[, Ap Capillary number = UmVl/
Viscosity ratio = Pressure drop=p num eact
[1.2 A }exact
10 0.87 6.52x10-4
100 0.86 6.61x10-4

5.4.2 Effect of Interface Reconstruction

The current level-contour based interface reconstruction is expected to have the

same characteristics as shown in the tests of Shin and Juric (2002). The effect of

reconstruction on mass conservation showing quadratic convergence with grid refinement

has already been established. Here the impact of reconstruction on spurious velocities is

assessed.

The tests were performed for Laplace number 2500 on a 50x50x50 grid. Temporary

spikes in spurious velocities were observed with every reconstruction (Figure 5-14(a)) as

also seen by Shin and Juric (2002). These spikes tend die down in time unless the

reconstruction is performed too frequently depending on the time scale of the problem.

Figure 5-14(a) shows that if reconstruction is performed every 10 Capillary time scale,

the temporary increase in spurious velocity tends to die down to the profile exhibited if

no reconstruction were performed. However, if the reconstruction is performed at every

Capillary time scale, the spurious currents are an order of magnitude larger and continue

to grow. Further observation can be made from Figure 5-14 (b) which shows the interface







81


normal velocity at a selected marker point. It shows that the oscillatory behavior caused

by the surface tension effects quickly die down if no reconstruction is performed during

the simulation. With reconstruction every 10 Capillary time scale, the amplitude of

velocity does not decay but remains stable with mean value of nearly zero. This is to be

expected based on the Laplace number of 2500 representing viscous time scales 2500

times slower than the Capillary time scales. But the reconstruction every Capillary time

scale clearly shows a completely different mean value that is several orders of magnitude

higher. These observations further strengthen the belief that reconstruction should be

performed as sparingly as possible.


o- -- ------- .......
. No r u rf l l T
-- No reconstruction
S-- Every Capillary time scale No reconstruction
E Every 10 Capillary time scale Every Capillary time scale
S t ---- Every 10 Capillary time scale
- Every 10 Capillary time scale
0 G


CI I


0 50 100 150 000
time time
(a) (b)

Figure 5-14. Effect of reconstruction frequency on spurious velocities. (a) Capillary
number evolution, (b) normal velocity at a selected marker point.

5.5 Buoyancy Driven Single Rising Bubble

A set of rising bubbles are presented next to establish the accuracy of the present

multidimensional immersed boundary procedure and implementation for handling

complex interfacial flows.









5.5.1 Grid Convergence Test

A single bubble of diameter d is placed in a computational domain of size and

orientation indicated in Figure 5-15(a). The bubble rises in the positive Y-coordinate

direction under the influence of gravity. The non-dimensional parameters used in the

simulations are Reynolds number (Re), Morton number (M), Eotvos number (Eo),

density and viscosity ratios. The properties of fluid outside the interface and bubble

diameter are taken as the required scales.











-s y







3d3d
3d
(a) (b)

Figure 5-15. Computational set-up for a single rising bubble. (a) Computational domain,
(b) initial adaptive grid with 3 levels of refinement over 7x21x7 base grid.

Simulations of Eo = 1.0 and M = 1.0x103 representing a high surface tension and

low Reynolds number case were conducted. To test the grid convergence, simulations

were conducted on two grids with three levels of refinement over a 7x21x7 base grid and

a 10x30x10 base grid. The 7x21x7 base grid produces maximum resolution equivalent to