UFDC Home  myUFDC Home  Help 



Full Text  
THREEDIMENSIONAL MARKERBASED 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 EulerianLagrangian ............................ ............. .. ............. 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 odim ensional interface.................................. ............... 42 3.3.2.2 Threedim 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 avierStokes Solver ........................................ .....................69 5.3.1 D ecaying V ortex Field ........................................ ......... ............... 70 5.3.2 Liddriven 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 eadon C oalescence ........................................ .......................... 87 5.6.2 O ffaxis C oalescence................................................... ............... ............90 5.7 B inary D rop C ollision................................................ ............................. 91 5.7 Phase Change Com putation.......................................... .......................... 97 5.7.1 1D 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 11 Some industrial and natural occurrences of twophase flows. ............. ................2 12 Key marker based 3D tracking computations in literature............... ...................5 21 Summary of some representative methods and related issues. .............................14 22 Summary of recent advances in interface tracking. .......................................15 51 Error estimates using Equation(5.2) for interface in a time reversed vortex field test comparing conservative and nonconservative restructuring.................................67 52 Volume error for reconstruction of a spherical interface and approximate radial perturbation required for correction. ............................................. ............... 69 53 Max velocity error and computation time normalized with time taken for 16x16 uniform grid. The adapted grids are shown in Figure 54....................................72 54 Lid driven cavity computation using uniform and adaptive grids .........................73 55 N natural convection com putation. ........................................ ......................... 76 56 Pressure drop for an inviscid fluids with density ratio =1.0..................................79 57 Effect of Laplace number on spurious velocity. .....................................................79 58 Effect of density ratio on pressure drop and spurious current on 50x50x50 grid (viscosity ratio = 1.0, La = 250) .................................................................. 80 59 Effect of viscosity ratio on pressure drop and spurious current on 50x50x50 grid (density ratio = 1000, La = 250)........................................................ ............ 80 510 Grid convergence test for Eo = 1.0, M= 1.0x103. The density and viscosity ratios are set to 20.0. ........................................................................ 83 511 Computation of rising bubbles with different terminal shape regimes ...................87 512 Computational parameter for rising bubblecoalescence tests .............................88 513 Properties of tetradecane and nitrogen ...... ..............................93 514 Binary drop collision computation with density ratio = 666.081 and viscosity ratio = 1 7 9 .2 8 ....................................................................... 9 4 515 Computational parameters for onedimensional phase change computation .........98 LIST OF FIGURES Figure p 11 Illustration of a typical multiphase flow ....................................... ............... 1 12 M ixed EulerianLagrangian tracking. ............................................. ............... 5 21 Illustration of Lagrangian interface tracking........... ..... .................10 22 Schematic of Eulerian tracking using level set and volume of fluid...................... 11 23 Interface in a time reversed vortex field using levelset method. ..........................12 24 A conceptual interface reconstruction using piecewise linear segments in the volum eoffluid m ethod. ........................................... ......................................... 13 25 A surgical alteration of the interface data to accomplish topology change .............15 26 A threedimensional marker based film boiling computation ..............................15 27 Illustration of conceptual differences between markerbased tracking with and without connectivity inform ation. ........................................ ........................ 17 28 Surface tension, pressure and normal component of shear stress acting on an in terfa c e ............................................................................................... 1 8 29 Sharp and continuous interface methods............................................ .............21 210 The effect of fluid property jumps on the computation time for a rising bubble.....23 211 A cellbycell and blockbyblock grid adaptation. ............................................25 212 A treebased adaptive Cartesian griddata. ................................... ............... 25 31 Schematic of the immersed boundary method. ....................................................28 32 Probe based technique for temperature gradient computation on the interface.......30 33 A triangulated interface and the datastructure. ........................................... ........... 31 34 Indicator function and the material flags for cells completely outside or inside the in te rfa c e .......................................................................... 3 3 35 Momentum source term due to surface tension force on a circular interface. .........35 36 Computation of the unit normal and tangent vectors on interface triangles. ...........35 37 Curvature computation for a unit circle using Equation(3.18) and cubic spline with (a) tw elve and (b) sixty points ....................................................... ...................36 38 Interface restructuring approach .........................................................................38 39 Conservative marker deletion for twodimensional interfaces..............................39 310 Conservative marker deletion for threedimensional interfaces ............................39 311 Levelcontour based interface reconstruction. ................................. ...............40 312 Interface reconstruction criterion uses two probes from the interface along the normal vector and bilinearly interpolates the indicator values Indl and Ind2........41 313 Two dimensional reconstructed interface edges (el, e2 and e3) and corresponding connectivity data. ................................................... ................. 42 314 Illustration of 2D cellbycell interface reconstruction. .........................................43 315 A 2D cell with four markers creating two locally disconnected edges ..................44 316 A reconstructed interface segment in a cell. ................................. .................45 317 Some common reconstructed interface segments and integer vertexflags. ...........45 318 Difficulties due to degeneracy of reconstructed interface data and implemented re m e d y ........................................................................... 4 6 319 Orientation setting of reconstructed interface. ............. ................. ..................... 47 320 Examples of 3D interface reconstruction. ...................................... ..................47 41 Cartesian grid datastructure. .............................................................................. 51 42 A daptive grid refinem ent. ............................................................. .....................52 43 Cartesian cells sharing a face and the corner cells are not allowed to differ by more than one level of refinement.........................................................................53 44 Geom etry based adaptation. ............................................. ............................ 54 45 Staggered grid arrangem ent. ............................................ ............................ 58 46 Spatial discretization of convection term ..................................... ...............59 47 A onedimensional interface between i and i+. .............................................. 61 48 Linear interface temperature correction within one cell on each side of the in te rfa c e .......................................................................... 6 2 49 Presently employed interface temperature correction technique solves a correction equation sharply on each side of the interface. .............................. ......... ...... .63 51 Interface in a time reversed vortex field ................. ............. ....... ........... 67 52 Error estimate and grid convergence for interface in time reversed vortex field.....67 53 Grid convergence test for volume error due to reconstruction..............................69 54 D ecaying vortex fi eld. ...................................................................... ...................71 55 Grid convergence test for decaying vortex field showing maximum velocity error computed on uniform and adaptive grids ...................................... ............... 71 56 Streamlines for liddriven cavity simulation for Re = 100....................................72 57 L id driven cavity. ............................................. ................... .... ...... 73 58 The u velocity profile along a vertical line through the center of liddriven cavity using (a) uniform grids and (b) adaptive grids..................................................... 74 59 Indicator contour used for modeling the cylindrical computational domain. .........75 510 Computed temperature and error for heat transfer in a concentric cylinder. ...........75 511 Natural convection with Pr = 0.71 and Ra = 1.0x105. ...........................................76 512 Static bubble com putation ............................................... ............................. 78 513 Grid convergence tests for static bubble. ..................................... ............... 79 514 Effect of reconstruction frequency on spurious velocities ............. ...............81 515 Computational setup for a single rising bubble. .............................................. 82 516 A spect ratio com putation .............................................................. .....................83 517 Bubble rise with Eo = 1.0, M= 1.0x103, density and viscosity ratio of 20.............84 518 Computational domain for cases in Table 511 .......................................... ........... 85 519 Computed shapes from Table 511 .........................................................................86 520 Computed case forM = 0.971, Eo = 97.1. .................................... .................86 521 Computational setup for headon coalescence test. ..........................................88 522 H eadon coalescence. ...... ........................... ........................................... 89 523 Grid and streamlines at t = 0.125 sec for headon coalescence..............................89 524 Instance of interface reconstruction during headon coalescence of two rising bubbles at t = 0.116 s .................................... ................ ................... 90 525 Experimental photographs of bubble coalescence taken at 0.03s time intervals. ....90 526 O ffaxis coalescence. ...................................................................... ................... 9 1 527 Binary drop collision regimes marking the location of computed parameters with so lid c irc le s ................................................................................................ ..... 9 3 528 Impact parameter B is defined as B=hd. ...................................... ............... 93 529 Binary drop collision Case 1 from Table 514..................................................100 530 A cross section of the computational grid and the time history of grid datasize for Case 1 from Table 514........................................ .......... .... ................. 101 531 Binary drop collision Case 2 from Table 514..................................................102 532 A cross section of the computational domain showing the interface and the velocity vectors for Case 2 from Table 514.................................... ........................ 102 533 A cross section of the computational grid and the time history of grid datasize for Case 2 from Table 514........................................ .......... .... ................. 103 534 Binary drop collision Case 3 from Table 514..................................................104 535 A cross section of the computational grid and the time history of grid datasize for Case 3 from Table 514........................................ .......... .... ................. 105 536 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 514..106 537 Binary drop collision Case 4 from Table 514..................................................107 538 A cross section of the computational grid and the time history of grid datasize for Case 4 from Table 514........................................ .......... .... ................. 108 539 A cross section of the computational domain showing the interface and the velocity vectors for Case 4 from Table 514.................................... ........................ 109 540 One dimensional phase change setup ...................................... ........................... 110 541 Comparison of computed and theoretical solution for Case 1 in Table 515.........110 542 Onedimensional phase change computation for Case 2 in Table 515 showing interface location on the left and temperature distribution at t = 0.1 on the right.. 111 543 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 Diracdelta and discretized Diracdelta 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 THREEDIMENSIONAL MARKERBASED 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, multilevel, threedimensional 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 markerspacing, and a connectivity preserving level contourbased 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 interfacebased 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 11 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 11. 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 liquidliquid, gasliquid, liquid or gas with solid. Table 11 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 11 is a prohibitively difficult task, the scope of the present work is controlled by focusing only on incompressible liquidliquid or gasliquid flows for bubble and drop dynamics problems. Table 11. Some industrial and natural occurrences of twophase flows. Phase Phase configuration Occurrence combination GasSolid Dispersion of solid Aerosol, particulate matter pollution; filters and particles in gas particlecollection devices; Pneumatic transport of particles; gasfluidized beds LiquidSolid Dispersion of Cells and particles in blood; hydraulic transport of particles in liquid particles; liquidfluidized beds GasLiquid 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 Liquidliquid Dispersion of liquid Emulsion and creams; breakup 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 oilwater 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 levelset; Welch and Wilson (2000) using volumeoffluid; 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 computationsystem. 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: * Multifluid 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 NavierStokes computation. The interface is tracked on a stationary Cartesian grid using time dependent triangulated surfaces (Figure 12). 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 EulerianLagrangian type. Although this approach has several advantages due to explicit knowledge of the interface shape/location for computing the geometric properties, the difficulties in datahandling have deterred a wider usage. Table 12 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 triangulatedinterface 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 levelcontour 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 12. Mixed EulerianLagrangian tracking. (a) A two dimensional interface using connected markers on a background grid and (b) a triangulated interface. Table 12. 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 datastructure 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 nonuniform nature of the grid makes staggered grid relatively more difficult to implement than the corresponding collocated grid (Singh & Shyy, 2006), but provides a compact pressurevelocity 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 coarsefine gridinterfaces 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 levelset 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 levelset 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 connectivityfree 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 masstransfer computation during phase change. Strategies to incorporate levelset 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 NavierStokes computation algorithm for velocity and pressure computation. c. Development of a sharp interface, collocated grid computation procedure for temperature equation in heattransfer/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 NavierStokes 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 EulerianLagrangian. 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 stateoftheart developments. 2.1.1 Lagrangian Lagrangian method discretizes the computational domain with the interface as a boundary for creating a body fitted grid (Figure 21). 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 21. Illustration of Lagrangian interface tracking. 2.1.2 Eulerian Eulerian methods use a scalar function to represent the interface implicitly. The levelset (Osher & Fedkiw, 2001, 2003) and volumeoffluid (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 levelset method and volume fraction in the volumeoffluid method (Figure 22) 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 22. 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 EulerianLagrangian 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 12). 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 21 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, levelset method (LS) suffers from erroneous mass loss that may become unacceptably large and detrimental. Such a massloss 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 timereversed vortex field. Figure 23(a) shows the sever loss in volume of a spherical interface in such a field using levelset method for interface tracking. This drawback was ameliorated by Enright et al. (2002) using a particle levelset method that uses Lagrangian markers around the interface to avoid erroneous breakup/deletion of the levelset characteristics (Figure 23(b)). "p"i(^ W% c2A ......  fei8 1' 31% 641k .me.M . (a) (b) Figure 23. Interface in a time reversed vortex field using levelset 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 levelset method, the volumeoffluid tracks the volume (mass) explicitly and hence does not suffer from massloss problems. However, volumeoffluid methods have been observed to produce spurious flotsam ('floating wreckage') and jetsam ('jettisoned goods') in the regions of high interfacecurvature underresolved 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 24). 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 24. A conceptual interface reconstruction using piecewise linear segments in the volumeoffluid method. The mixed EulerianLagrangian 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 EulerianLagrangian 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 twodimensional interfaces (Figure 25), 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 26). 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 22 presents some key developments by combining various individual techniques to exploiting their attractive features. Table 21. 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 geometrycomputation. Severe mass loss/gain in under (Levelset) 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 datastructure and Eulerian subgrid levels, bookkeeping for 3D interfaces. Lagrangian Flow computations on Eulerian grid. Difficult for topological change. Table 22. Summary of recent advances in interface tracking. Base method Levelset (LS) Volume of fluid (VOF) LS and VOF Marker tracking Improved method Particle levelset (Enright et al., 2002) Hybrid markerVOF (Aulisa et al., 2003, 2004) Coupled LS VOF (Sussman, 2003) The pointset 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 levelset interface. * Employs makrer based interface to avoid interface reconstruction for geometricproperty computation. * Improved performance in underresolved areas. * Excellent massconservation due to VOF. * Easy geometrycomputation using Levelset. * A set of unconnected. * Lack of connectivityinformation makes it easier to handle topological changes. * Geometrycomputations require splinefitting. * A set of unconnected triangles. * Possible to compute geometry (curvature, normal) without any surface fitting or interpolation. * Levelcontour reconstruction for topology change. 1 6 Figure 25. A surgical alteration of the interface data to accomplish topology change. t~Ep~ a Figure 26. A threedimensional marker based film boiling computation. (a) A snapshot of the computational domain and (b) a triangulated interface before and after pinchoff 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 datahandling, 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 levelcontour based reconstruction with connectivityfree 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 markerspacing on the interface. As seen in Figure 27(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 connectivityfree 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 27(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 27. Illustration of conceptual differences between markerbased 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 fluidproperty 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 28 shows the pressure (P) and the normal shear stress components (n.rn) 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 nT n Fluid 2 Figure 28. 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 cutcell method of Ye et al. (1999, 2001). This method cuts the background Cartesian grid with the interface and creates subdomains (Figure 29(a)) that act like body fitted grids with arbitrary shaped cells on the interface. For stability reasons, the cutcells are merged with suitable neighbors to remove arbitrarily small cells. The governing equations are solved separately in each subdomain 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 ghostfluid techniques use the jump conditions to construct boundary conditions for flow variables on the interface. As an example, the ghostfluid technique solves the governing equations in each phase separately by creating ghostcell 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 cutcells. 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 29(b) is usually three to four cells in thickness. The continuous interface method has been widely used with mixed EulerianLagrangian tracking (i.e., the immersed boundary method) (Peskin, 1977; Tryggvason et al. 2001; Francois & Shyy, 2003) as well as Eulerian tracking such as levelset and volumeoffluid. 2.2.3 Challenges and Recent Advances The iterative procedure of coupling the phases across the interface and the complex procedure of cellcutting 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 levelset based interface representation. Although successful demonstration of ghost fluid type approach with EulerianLagrangian 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 YoungLaplace 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 512(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 nonzero velocity field termed as parasitic currents or spurious velocities (Figure 512(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 104 (in terms on Capillary number Ca) (Tryggvason, 2001). The first order accurate ghost fluid technique could also keep these velocities to the order of 106 or lower (Kang et al., 2000). Several attempts using improved curvature/surfacetension 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 29. Sharp and continuous interface methods. (a) Sharp interface cutcell 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 heartvalves; 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 microscale 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 flowgoverning 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 computationtime. 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 buoyancydriven rising bubble using the immersed boundary method. The Reynolds number, Weber number and Froude number were set to 100, 4 and 1.0. A Wcycle multigrid technique was used to solve the pressure Poisson equation. The horizontal axis in Figure 210 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 BDI 0.4 E!! 1 0.4 X0 0.35  1 """' Il'O20.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 210 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 datasize. 2.3.1 Cartesian Grid Data Structures The approach for generating adaptive Cartesian grids is divided into two categories: a cellbycell adaptation (Figure 211(a)); a blockbyblock adaptation (Figure 211(b)). Adaptive mesh refinement known as AMR is an example of a blockbyblock adaptation (Berger & Oliger, 1984; Sussman, 1999, 2005) that interprets the grid as a union of uniformly refined gridblocks (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 datasize 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 datasize, cellbycell adaptation is a more flexible approach. A cellbycell adaptation uses either a treebased 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 treebased datastructure (Figure 212) offers a compact storage system. All the connectivity information is implicitly contained in the location of cells in the treestructure. One of the most common queries in a numerical simulation is the cellcell 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 1020% of the CPU time for CFD solvers is dedicated to tree traversals answering gridconnectivity queries. (a) (b) Figure 211. A cellbycell and blockbyblock grid adaptation. (a) A cellbycell adaptation in the present work; (b) schematic of a blockbyblock 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 212. A treebased adaptive Cartesian griddata. (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 directionbased refinement strategy. However, the resulting datastructure 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 cellbycell 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 NavierStokes 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 31(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 31(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 masstransfer 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 31. 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 32 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 (KVTn = 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 masstransfer 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 nondivergencefree 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 32. Probe based technique for temperature gradient computation on the interface. 3.1.2 Material Property Smoothing Figure 33 shows a common finiteelement type datastructure used to store the triangulated interface. It assigns a unique identification (integer number) to each node (marker) and stores its coordinates along with the trianglenodes. 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 nodel node2 node3 1 2 ntria Figure 33. A triangulated interface and the datastructure. 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(xxint) is a Diracdelta function nonzero 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 (45) 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 34). 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(332h+ 1+4h(1h) :ifh e[0,l) 8 (3.14) h= 1=<(52h 7+4h(3h)) :ifhe(1,2] 0 :else A raytracing algorithm popular in computergraphics and computationalgeometry 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 computergraphics rendering is used. Unlike the raytracing algorithm, the painter's algorithm does not require expensive computation of threedimensional 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 34. 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 35 shows the smoothed nature of the momentum source for a twodimensional 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 twodimensional interfaces (Francois & Shyy, 2002, 2003; Ye et al., 2001) and surface fitting for threedimensional 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 (AlRawahi & 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 lineintegral 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 AlRawahi and Tryggvason (2004) shown in Figure 36. To ensure conservation (i.e., the net surface tension force on a closed interface is zero), the crossproduct 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 37 comparing curvature of a unit circle using present method and a cubicspline 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 AlRawahi 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 35. Momentum source term due to surface tension force on a circular interface. (a) The xcomponent, (b) ycomponent and (c) the xcomponent along a horizontal line though the center of circle. tl1t2 Stl t21 t2 Ip1l pl p3 \p3 Sp3p2 triangle tl p3 p2I t3 p2 p Iv2 11 p2 triangle => pl, p2, p3 Figure 36. 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 37. 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)dxdydzm + P2 n (3.20) l (c ell) 2pI 2 ) For simplicity, the marker movement due to background fluid velocity is solved using 2nd Order RungeKutta 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 phasevolume 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 (107%). This correction may be performed either at every timestep 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 masstransfer 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 Cartesiancellsize) will appear as a 'hole' on the interface. On the other hand, excessively small edges can produce unphysical subgrid 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 38) do not preserve the interface volume. More precisely, it is the marker deletion that is nonconservative. Such nonconservative behaviors are demonstrated by Sousa et al. (2004) who encountered up 5% volume errors due to interface restructuring. (a) Pi (b)  (b) Figure 38. Interface restructuring approach. (a) Long edges are split at the midpoint and (b) a small edge pp2 is collapsed to its midpointp3. The presently employed restructuring approach is same as in Figure 38 with an added correction step to the marker deletion procedure in order to preserve the interface volume. This correctionalgorithm 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 subgrid undulations on the interface. The underlying procedure is explained using a 2D example shown in Figure 39 that shows a marker removal by collapsing the edge formed by point p2 and p4 to a point p3 resulting in a net phasearea loss. To correct this error, first a reference pointpref atp3n is selected (Figure 39(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 39(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 pn and compute a reference volume by adding the tetrahedralvolumes 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 310(b). 4. Relocate the point p to p + n (Q0/Q) for preserving the local phase volume Q. Figure 39. Conservative marker deletion for twodimensional interfaces. (a) Old interface pip2P4P5 and new interface pip3p5 with nonconservative 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 310. Conservative marker deletion for threedimensional interfaces. (a) A reference pointpref is defined atpn (p is edge midpoint) 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 gridcell 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 311 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 312 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 cellsize. 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 312. Interface reconstruction criterion uses two probes from the interface along the normal vector and bilinearly 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 celledges are connected to define interface segments within the Cartesian cells (Figure 313). The following sections present the overall reconstruction procedure using a two dimensional interface and the ideas are extended to threedimensional 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 313. Two dimensional reconstructed interface edges (el, e2 and e3) and corresponding connectivity data. 3.3.2.1 Twodimensional interface The indicator function values computed at cell centers are linearly interpolated to Cartesian cellvertices. 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 314(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 314(b) all same physical location i.e., indicator = 0.5 contour is at the Cartesiancell 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 314. Illustration of 2D cellbycell interface reconstruction. (a) Markers on edges with opposite vertexflag 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 nonconnected 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 cellcenter. 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 315. 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 Threedimensional interface The basic approach for threedimensional 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 316 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 317 shows some typical reconstructed polygons with corresponding vertexflags. The last two cases in Figure 317 contain six markers each but one has a single sixsided 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 316. 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 317. Some common reconstructed interface segments and integer vertexflags. 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 318(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 318(b)). Similar to 2D case, whenever required the cells are locally broken into tetrahedrons (Figure 318(c)) by connecting the cell center to cellvertices 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 318. 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 319(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 319(b)). The triangles on the entire interface can be correctly oriented starting with proper 47 orientation of any single triangle. Figure 320 shows examples of reconstruction creating a rupture and merger of two interfaces. VI nI r P Lr:: q   q q tria => p, q, r nghtria => r, q, s (a) (b) Figure 319. 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 crosssection (a) . . . 3D reconstructed Figure 320. 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 datamanagement. 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 timereversed 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 subgrid 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 NavierStokes computation. 4.1 Data Structure Figure 41 shows a facebased datastructure 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 datastructure is represented here. The variable level in Figure 41 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 celltoface connectivity information. For every face, a single byte integer named orientation gives the orientation of the face: faces normal to xaxis have orientation = 1, faces normal to y axis have orientation = 2 and faces normal to zaxis 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 precomputed 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 precomputed 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:: cellface list Face Data Integer:: orientation FInteger:: sideCelll, sideCell2 (a) (b) Figure 41. Cartesian grid datastructure. (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 2D domain in Figure 42(a) uniformly split into four parts in xdirection and three parts in ydirection. Let the number of uniform cells in x and y directions be Nx and Ny, respectively. In twodimensions, 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 42 (b) is isotropically split, coordinates of the new cells shown in Figure 42(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 42. Adaptive grid refinement. (a) A 4x4 base grid (level= 0), (b) coordinates of a twodimensional 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 flowsolution 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 datastructure. 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 nonuniform 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 43) 3. Corner cells are not allowed to differ by more than one level of refinement (Figure 43). 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 43. 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 44 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 44. 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 cellvolume. 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&u3=/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 cellcentered values are reconstructed linearly while the facenormal 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 cellcentered or facecentered 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 RungeKutta 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 NavierStokes 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 NavierStokes 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 facenormal velocity is stored on Cartesian cell faces (Figure 45(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 RungeKutta and fully implicit method, respectively. n+l n Cn 1T T n+1/2 n+1 (pC) lr" =v.().uT) 2 +v.(KVT).l Step 2: Velocity computation using projection method. Predictorstep 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 timestep viscous term due to CrankNicholson 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 Correctorstep Correct the predicted velocity field (**) using Equation(4.7). The pressure field for this correction is computed by enforcing velocitydivergence 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 VU=  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 45(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 nonuniform cells, it is accepted as a compromise for algorithmic simplicity especially when nonuniform 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 45. 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 controlvolume 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 46(a) and the velocity constraints, the control volume facenormal velocities are computed using face area weighted averaging that produces: U, i fae 1 Px2 2U 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 facenormal velocities. To conserve the momentum flux across non uniform grid cells, Figure 46(b) makes another approximation that computes the flux on top controlvolume face for velocity U1 and equally distributes to the corresponding control volumes of faces for U2 and U3 at a Tnode junction. The flux term q in Equation (4.11) is computed using 2nd or 3rd order essentially nonoscillatory 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 nonuniform 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 'IUcf3 I, Ucf2 I U,' U3 I I I   f/2 f/2 I. . .. U, 4 Areal < Area2  (a) (b) Figure 46. Spatial discretization of convection term. (a) Control volume (CV) and face normal velocities (Ucf); (b) flux fat a Tnode 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 controlvolume facenormal 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 46 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 cellcenter, a collocated grid discretization with auxiliary variable creation for nonuniform 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 levelset method. The discretization in the interfacial cells can be explained using Figure 47 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 T1 i terface west i east i+1 Xx I dx(414 & dxdelta Figure 47. A onedimensional 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 cellcenter temperatures as shown in Figure 48. Corrected Ti Ti 1. /Corrected Ti+1 {T i1 Tsat / Ti. T i+1 T i+2 I I I I I I I I Interface location Interface temperature = Tsat Figure 48. 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 levelset 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 49, this method will correct the temperature of Ti and Ti+l while leaving Ti1 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 Ti1 as BC v.(KVT)=o v.(KVT)= T asBC T Ti+2 as BC Ti Ti+1 DD \ ~ I D to be corrected Figure 49. 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 masstransfer 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 massconservation 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 timereversed vortex field (Figure 51) 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 51 presents the computed errors showing a better performance by the conservative restructuring as compared to the nonconservative restructuring presented in Chapter 3. Figure 51 shows the time history of the volumeerror 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 51.Error estimates using Equation(5.2) for interface in a time reversed vortex field test comparing conservative and nonconservative restructuring. Grid RMS Radius error % Area error % Volume error (conservative, non (conservative, non (conservative, non conservative) conservative) conservative) 30x30x30 (0.2435, 1.0254)x102 (4.2898, 7.9391) (0.5591, 14.0950) 60x60x60 (0.0091, 0.4212)x102 (1.5306, 2.4060) (0.1178, 4.9882) 120x120x120 (0.0079, 0.1155)x102 (0.4169, 0.5978) (0.0216, 1.0455) St=0 t=T/4 t=T/2 t =3T/4 t= T (a) (b) Figure 51. 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 Nonconservative Conservative SNonconservative 0 10  > o E S10 0 1 2 3 4 40 60 80 100 120 time Grid resolution (a) (b) Figure 52. 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 indicatorfunction computation and hence the levelcontour based reconstruction depends on the background grid resolution. The aim is to measure interface volumeerror 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 53) 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 52 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 100 20 40 60 80 100 Number of cells per diameter Figure 53. Grid convergence test for volume error due to reconstruction. Table 52. Volume error for reconstruction of a spherical interface and approximate radial perturbation required for correction. Grid cell perdiameter (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 NavierStokes 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 freeslip 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 54(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 uniformgrid 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 54(b). The maximum error in numerically computed ucomponent of the velocity field and the corresponding computationtime is presented in Table 53 with the time normalized by the time taken for 16x16 uniformgrid computation. As seen in Figure 55, 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 nonuniform neighbor cells (Figure 54), the computed error is a reasonable reflection of the accuracy of spatial discretization over 71 nonuniform 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 computationtime savings that are expected to be more pronounced for 3D computations. (a) (b) Figure 54. Decaying vortex field. (a) Streamlines; (b) two refinement levels on a 16x16 base grid. 102 10 a, x 1    Uniform A Adaptive  2nd order reference 100 200 300 Square root of number of cells Figure 55. Grid convergence test for decaying vortex field showing maximum velocity error computed on uniform and adaptive grids. Table 53. Max velocity error and computation time normalized with time taken for 16x16 uniform grid. The adapted grids are shown in Figure 54. Max. grid Uniform grid computation Adaptive grid base computation resolution Max. U error Normalized time cells Max. U error Normalized time 16x16 2.966x103 1.000 256 2.966x103 1.000 32x32 7.184x104 4.440 964 1.105x103 1.946 64x64 1.752x104 21.602 3616 4.111x104 10.725 128x128 4.344x105 104.810 14548 1.611x104 90.64 5.3.2 Liddriven 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 56 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 57(a) and they are in accordance with the computations of Ghia et al. (1982).  Figure 56. Streamlines for liddriven 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 57. Figure 58 shows the uvelocity 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 uvelocity in Figure 58 are presented in Table 54. This table also provides a comparison of the total number of grid cells and computational time for uniform and adaptive grids. The secondlast 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 54. Lid driven cavity computation using uniform and adaptive grids. Comparison of the number of cells, computation time and minimum uvelocity 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 57. 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 58. The u velocity profile along a vertical line through the center of liddriven 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 59. The computations are performed where indicator >= 0.5 (space between the concentric circles). Figure 510 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 59. Indicator contour used for modeling the cylindrical computational domain. (a) Computations are performed only within the domain with Indicator >= 0.5; (b) crosssection 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 510. Computed temperature and error for heat transfer in a concentric cylinder. (a) temperature distribution (TTi)/(ToTi); (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 TT 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 55 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 nondimensionalized with AT and L) and vmax as the maximum vertical velocity on a horizontal line through the center of the domain. Figure 511 shows the computed streamlines and temperature isotherms. (a) Figure 511. Natural convection with Pr = the temperature isotherms. Table 55. 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 512(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 YoungLaplace Equation(5.8) (Figure 512(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 nondimensional spurious velocities (Capillary number) based on the observations of Shin and Juric (2002) and Sousa et al. (2004) is of the order of 104 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 102 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 slipwall 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 nondimensionalized 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 512. Static bubble computation. (a) Spurious currents of magnitude Ca 104, (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 513(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 56 and Figure 513(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 57. High density and viscosity ratios (1000 and 100 respectively) were also tested and found not to cause very serious degradation as shown Table 58 and Table 59. 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 513. Grid convergence tests for static bubble. (a) Capillary number evolution, (b) error in pressure drop (Equation(5.10)) using data in Table 56. Table 56. Pressure drop for an inviscid fluids with density ratio =1.0. Data was taken after three nondimensional time step of 5.0x103. Grid resolution Pressure drop error (Equation (5.10)) 25x25x25 0.197 50x50x50 0.151 100x100x100 0.117 Table 57. Effect of Laplace number on spurious velocity. Laplace number (La = opd d /2 ) Capillary number (Ca = Um,,/cr) 250 5.32X104 12000 5.35X104 12000 5.35x1 04 Table 58. 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.32x104 10 0.89 5.27x104 100 0.87 5.15x104 1000 0.87 6.35x104 Table 59. 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.52x104 100 0.86 6.61x104 5.4.2 Effect of Interface Reconstruction The current levelcontour 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 514(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 514(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 514 (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 514. 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 515(a). The bubble rises in the positive Ycoordinate direction under the influence of gravity. The nondimensional 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 515. Computational setup 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 