Citation |

- Permanent Link:
- http://ufdc.ufl.edu/AA00029939/00001
## Material Information- Title:
- Computational study of leukocyte rheology based on a multilayer fluid model
- Creator:
- Kan, Heng-Chuan, 1967-
- Publication Date:
- 1997
- Language:
- English
- Physical Description:
- vi, 117 leaves : ill. ; 29 cm.
## Subjects- Subjects / Keywords:
- Cytoplasm ( jstor )
Deformation ( jstor ) Eggshells ( jstor ) Flow distribution ( jstor ) Leukocytes ( jstor ) Momentum ( jstor ) Neutrophils ( jstor ) Newtonianism ( jstor ) Velocity ( jstor ) Viscosity ( jstor ) Aerospace Engineering, Mechanics and Engineering Science thesis, Ph. D ( lcsh ) Dissertations, Academic -- Aerospace Engineering, Mechanics and Engineering Science -- UF ( lcsh ) - Genre:
- bibliography ( marcgt )
non-fiction ( marcgt )
## Notes- Thesis:
- Thesis (Ph. D.)--University of Florida, 1997.
- Bibliography:
- Includes bibliographical references (leaves 109-116).
- General Note:
- Typescript.
- General Note:
- Vita.
- Statement of Responsibility:
- by Heng-Chuan Kan.
## Record Information- Source Institution:
- University of Florida
- Holding Location:
- University of Florida
- Rights Management:
- The University of Florida George A. Smathers Libraries respect the intellectual property rights of others and do not claim any copyright interest in this item. This item may be protected by copyright but is made available here under a claim of fair use (17 U.S.C. Â§107) for non-profit research and educational purposes. Users of this work have responsibility for determining copyright status prior to reusing, publishing or reproducing this item for purposes other than what is allowed by fair use or other copyright exemptions. Any reuse of this item in excess of fair use or other copyright exemptions requires permission of the copyright holder. The Smathers Libraries would like to learn more about this item and invite individuals or organizations to contact the RDS coordinator (ufdissertations@uflib.ufl.edu) with any additional information they can provide.
- Resource Identifier:
- 027599631 ( ALEPH )
37266911 ( OCLC )
## UFDC Membership |

Downloads |

## This item has the following downloads: |

Full Text |

COMPUTATIONAL STUDY OF LEUKOCYTE RHEOLOGY BASED ON A MULTILAYER FLUID MODEL By HENG-CHUAN KAN 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 1997 ACKNOWLEDGMENTS I would like to express my sincere thanks toward my advisor, Dr. Roger Tran-Son-Tay, for his expert guidance, encouragement, and financial support during the course of my graduate studies. He has been a constant source of inspiration for me both professionally and philosophically. The same gratitude is also given to my cochair, Dr. Wei Shyy, for his understanding and help in so many ways. Very special thanks go to Dr. H. S. Udaykumar for his tremendous help and enthusiasm in this work. In spite of spending the countless and frustrated hours in front of a variety of computers with me during the early debugging process, his patience and belief kept this work going. I would like to thank the members of my supervisory committee, Dr. Renwei Mei, Dr. Richard B. Dickinson and Dr. Nicolae D. Cristescu, for their fruitful advices. I wish to thank my colleagues at the Laboratory of Cellular Mechanics and Biorheology for their direct and indirect contributions and entertaining discussions. In particular, I thank Emmanuelle, David, Neal, Jie, Philou, and Philippe for their friendship and assistance. Words are not enough to describe my sincere gratitude to my parents, Pao-Hsiang Wang and Che-Chang Kan, for their persistent support and love. Their emphasis on higher education for their children pave the opportunity for me to complete this work. I also wish to thank my brothers, Po-Wen Kan and Feng-Jung Kan, for their continuous encouragement throughout these years. Finally, I reserve my deepest appreciation to my wife, Min-Hui. I will be forever indebted to her endless sacrifices and love. To her, I dedicate this dissertation. TABLE OF CONTENTS ACKNOWLEDGMENTS .......................................... ii A BSTRA CT ... ................................................... v 1. INTRODUCTION AND BACKGROUND............................ I 1.1. M otivation and Objectives .................................... 1 1.2. Importance of Compound Drops Dynamics ....................... 2 1.3. Significance of Leukocyte Rheology ............................. 2 1.4. Observed Rheological Behavior of Leukocytes .................... 4 1.5. Current Status of Leukocyte Modeling ........................... 6 1.6. Summary ................................................. 9 2 THEORETICAL ANALYSIS ..................................... 12 2.1. Background ... ............................................ 12 2.2. ELAFINT: A Mixed Eulerian-Lagrangian Solution Algorithm ........ 15 2.2.1. Characteristics of the Numerical Methodology ............... 15 2.2.2. Interface Tracking Procedure ............................ 16 2.2.3. Formulation of the Pressure-Based Flow Solver ............. .18 2.2.4. Discretization Procedure of the Flow Field Equations ......... .20 2.3. M omentum Interpolation ..................................... 23 2.4. Interactions between Interfacial Physics and Flow Field Variables ..... 28 2.4.1. Cut-Cell Procedures for Solid-Liquid Boundaries ............ 28 2.4.2. Immersed Boundary Technique for Fluid-Fluid Interfaces ...... 31 2.5. Mass Conservation in Individual Phases .......................... 36 2.6. Summary ........... ................ ................... 39 3 ASSESSMENT OF COMPUTATIONAL TECHNIQUE ................. 40 3.1. Flows With Fixed Arbitrary-Shaped Internal Boundaries ............ 40 3.2. Moving Boundary Flow Problems ............................... 47 3.2.1. Deformation and Recovery of Viscous Drops in Extensional Flows 47 3.2.2. Deformation of Drops in Extensional Flows with Inertia Effects. 58 3.2.3. Deformation of Drops in Constricted Tubes ................. 65 3.3. Summary ................................................. 71 4 DYNAMICS OF COMPOUND DROPS ............................ 72 4.1. Problem Statement and Mathematical Formulation ................. 72 4.1.1. Computational Setup ................................... 76 4.1.2. Selection of Model Parameters ........................... 77 4.2. Deformation Analysis ....................................... 77 4.2.1. Steady-State Configuration ( Ca < Cacr ) ................... 77 4.2.2. Unsteady Response ( Ca > Cacr ) .......................... 79 4.3. Recovery Analysis .......................................... 81 4.3.1. Effect of the Capillary Number ........................... 81 4.3.2. Effect of the Core Viscosity .............................. 84 4.3.3. Comparison of Simple and Compound Drops ................ 86 4.4. Summary .................................................. 88 5 LEUKOCYTE MODELING ..................................... 89 5.1. Choice of Model Parameters ................................... 89 5.2. Effect of a Viscous Cytoplasmic Layer ........................... 90 5.3. Effects of Nucleus Deformation and Time Scales ................... 93 5.3.1. Incompatible Nucleus and Cytoplasm Time Scales ............ 93 5.3.2. Compatible Nucleus and Cytoplasm Time Scales ............. 96 5.4. Discussion And Implications to Leukocyte Rheology ............... 102 6 SUMMARY AND FUTURE RESEARCH ............................ 105 6.1. Conclusion ........................ ........................ 105 6.2. Future Research ......................................... 106 6.2.1. Experimental Work ................................... 107 6.2.2. Numerical Algorithm .................................. 107 REFERENCES ............................................... 109 BIOGRAPHICAL SKETCH ........................................... 117 Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirement for the Degree of Doctor of Philosophy COMPUTATIONAL STUDY OF LEUKOCYTE RHEOLOGY BASED ON A MULTILAYER FLUID MODEL By Heng-Chuan Kan May 1997 Chairperson: Dr. Roger Tran-Son-Tay Major Department: Aerospace Engineering, Mechanics and Engineering Science Knowledge of the rheological properties of leukocytes (white blood cells) is essential not only for the comprehension of microcirculatory flow dynamics, but also for the understanding of their functions and behavior in health and disease. These properties, together with the leukocyte's structural characteristics, determine the deformability of the cell, especially during large deformations, such as those involved in their release from the bone marrow or extravasation into the interstitium. Existing viscoelastic solid, Newtonian and nonNewtonian liquid drop models used to calculate the cell rheological properties do not adequately explain reported experimental observation. This is because such models do not take into account the morphology of the cell and its multi-layer structure. In this work, the deformation and recovery of a compound liquid drop (three-layer fluid model) are investigated by means of a mixed Eulerian-Lagrangian computational methodology which allows large viscosity and capillarity differences between layers. Analysis of the deformation process confirms published results and clearly indicates that the history of the deformation is stored in the compound drop in terms of the shape of its components. The recovery process is found to provide information on the hydrodynamics of compound drops that is not supplied by the deformation analysis. It is discovered that a compound drop does not recover like a single phase Newtonian liquid drop except when the core is sufficiently deformed and its material properties are such that the core deformation/recovery time scale is compatible with that of the shell layer. A difference in the core and shell layer time scales causes an initial rapid recoil of the drop during which the shell fluid is the sole participant in the hydrodynamics, followed by a slower relaxation period during which the core and shell layer couple. The large deformation and subsequent recovery processes of compound drops also provide new insights into the complicated behavior of leukocytes. This three-layer fluid model describes the morphology of leukocytes better than existing ones. It consists of an outer interface (cell membrane), containing a shell layer (cytoplasm), and a core (nucleus). Based on the compound drop dynamics, it is concluded that (1) the nucleus plays an essential role in the cell response, and (2) the energy stored in the deformed interfaces can be important. The present investigation indicates that unless the nucleus and its deformation are included in the analyses, neither Newtonian nor non-Newtonian drop models for leukocytes can yield a reliable picture of the hydrodynamics of such cells. CHAPTER 1 INTRODUCTION AND BACKGROUND 1.1. Motivation and Objectives Knowledge of the rheological properties of leukocytes (white blood cells) is essential not only for the comprehension of microcirculatory flow dynamics, but also for the understanding of their functions and behavior in health and disease. These properties, together with the leukocyte's structural characteristics, determine the deformability of the cell, especially during large deformations, such as those involved in their release from the bone marrow or extravasation into the interstitium. This deformability is also essential in the leukocyte's response to disease/infection in humans because it determines the ability of the cell to flow and deform in capillaries and/or migrate in tissue. Therefore, obtaining the correct mechanical description for the leukocytes has been a focus of research in the recent years. Existing viscoelastic solid, Newtonian and non-Newtonian liquid drop models used to calculate the cell rheological properties do not adequately explain reported experimental observation. This is because such models do not take into account the morphology of the cell and its multilayer structure. The three-layer (or compound drop) model has been suggested to be a good model for the leukocytes (Skalak et al., 1990; Dong et al., 1991; Hochmuth et al., 1993). This model gives a better representation of the morphology of the leukocyte and its interior over the existing ones. It consists of a cell membrane (outer interface), containing a cytoplasmic layer (shell), and a nucleus (core and inner interface). However, the three-layer model has not hitherto been approached analytically due to the limitation of available tools. The objectives of this work are: (1) to develop a general computational methodology for investigating the rheological behavior of multilayer drop systems, (2) to test the validity and accuracy of the algorithm, (3) to study the hydrodynamics of compound drops, and (4) to provide a better understanding of leukocyte rheology. 1.2. Importance of Compound Drops Dynamics The deformation and subsequent relaxation of a compound liquid drop is a fundamental fluid mechanics problem of medical importance. The underlying idea that makes the application of compound drops attractive and has spawned several new technologies (Jackson and Lee, 1991; Powell and Mahalingam, 1992; Batich and Vaghefi, 1994; Goel et al., 1994; Martin and Parthasarathy, 1995; Hassan et al., 1996; Thies, 1996) is that a core of active ingredient (whether solid, liquid or gas) can be protected by enclosing it within a shell of a second material. The active ingredient, which must not of course attack the chosen encapsulating material, is released at a determined rate, time and location to achieve certain localized effects. These include controlled release of drugs (Tai and Sun, 1993), flavors (Jackson and Lee, 1991; Gorski, 1994), enzymes (Hassan et al., 1996) or latent heat (Goel et al., 1994; Mulligan et al., 1996). While many aspects of compound drops have been studied in the literature (Johnson and Sadhal, 1985), including in particular deformation in an extensional flow (Stone and Leal, 1990) and stability (Landman, 1985; Oguz and Sadhal, 1987), the range of conditions under which compound drop dynamics has been investigated does not span the parameters of interest in this work. 1.3. Significance of Leukocyte Rheology The persistence and widespread occurrence of diseases associated with defects in the circulation of blood has largely provided the impetus for the study of blood rheology. Human blood is a suspension of erythrocytes (red blood cells), leukocytes (white blood cells) and platelets in plasma. Generally speaking, blood is considered as a Newtonian fluid in large vessels at shear rates greater than 100 s-1. When the shear rate is less than 10 s-1, blood behaves like a Casson fluid (Cokelet et al., 1963). However, for blood vessels less than 500 [tm in diameter, blood cannot be considered as a continuous fluid. Blood must be treated as a two phase fluid (core and plasma layers) in vessels with diameters between 30 to 500 tm. Several effects (e. g., Fahraeus, and Fahraeus-Lindqvist) occur in that diameter range and become stronger as the vessel diameter decreases (Fung, 1993). For blood vessels less than 20- 30 [tm in diameter, the presence of individual cells cannot be neglected. The flow is creeping and the rheological properties of individual blood cells have to be taken into account. The rheological behavior of erythrocytes has been intensively studied and found to be adequately modeled by a Newtonian fluid enclosed in a viscoelastic membrane (Evans and Hochmuth, 1976; Evans and Skalak, 1980). On the other hand, the search for a suitable model to describe the rheological behavior of leukocytes continues. Rheological studies of leukocytes have been neglected because they occupy less than 1% of the total volume content of blood. However, since leukocytes have found to have a larger volume and a lower deformability than red blood cells, the number of studies has steadily increased. Leukocytes, because of their relative size and rigidity, play a significant role on the overall circulatory flow. They can affect in particular the apparent viscosity of blood and therefore the pressure drops encountered across vessels and capillaries (Pries et al., 1994). In addition, the ability of the leukocytes to flow and deform in capillaries and/or migrate in tissue is essential in the response to disease/infection in humans. Therefore, knowledge of the rheological properties of the leukocytes is critical for understanding microcirculatory flow dynamics and its effects in health and disease. Given the microscopic size of the blood cells (about 8 - 10 [tm in diameter), experimental efforts face obvious hurdles. However, micro-manipulation techniques have been successfully developed not only for deforming but also for evaluating mechanical and chemical effects on blood cells (Schmid-Schonbein et al., 1981; Evans and Kukan, 1984; Needham and Hochmuth, 1990; Usami et al., 1992; Skierczynski et al., 1993; Tran-Son-Tay et al., 1994a). One common way of determining the rheological properties of the leukocyte has been to deform the cell by aspiration into a micropipet (3 -6 p.m in diameter) and then to eject the cell and observe the recovery characteristics. The recovery time and path then provide an estimate of the "apparent viscosity" of the cell. The apparent viscosity is an average rheological property of the whole cell, including all the inclusions. 1.4. Observed Rheological Behavior of Leukocytes The most populous type of leukocytes in human blood is the neutrophil. The neutrophil contains a relatively small segmented nucleus, occupying about 21% of its volume, and a cytoplasm that is rich in granules. In contrast, in a lymphocyte, another type of leukocyte, the nucleus occupies approximately 44% of the cells (Schmid-Schonbein et al., 1980). A neutrophil is spherical in a passive, non-reactive state. The membrane of the neutrophil is highly ruffled, giving it an excess surface area (estimated at 110% - 120% of the undeformed sphere), and permits deformation at constant volume. In spite of the intensive research effort, to date, a variety of behaviors of the leukocyte cannot be explained clearly. For example, when a neutrophil is aspirated into a micropipet it flows into the pipet when the aspiration pressure exceeds a value called the critical suction pressure. Neutrophils rapidly aspirated into a micropipet with a diameter smaller than the cell diameter enter at an initial velocity much faster than the steady-state velocity. Also, as the rate of deformation increases, the cell viscosity is found to decrease, which has been speculated to be due to a shear-thinning effect (Needham and Hochmuth, 1990; Tsai et al., 1993). Needham and Hochmuth (1990) observe that after the cell enters into the pipet, the subsequent deformation occurs at a constant viscosity, independent of the aspiration pressure. However, as the excess aspiration pressure increases, the apparent viscosity of the cell decreases. Furthermore, following small but rapid deformations (in less than 0.3 seconds) of the cell, a rapid initial recovery (around 0.3 seconds) is obtained, followed by a slower phase (Dong et al., 1991). On the other hand, when the neutrophil is aspirated and held inside the pipet for a short time period of 5 seconds or less, and then expelled, the cell recovers faster Cell Apparent Viscosity Authors Experiments Operating conditions g(Poise) Dong et al. (1988) Evans and Yeung (1989) Needham and Hochmuth (1990) Dong et al. (1991) Tran-Son-Tay et al. (1991) Dong and Skalak (1992) Hochmuth et al. (1993) Tsai et al. (1993) Waugh and Tsai (1994) Aspiration Recovery Aspiration Small deformation Large deformation High deformation rate Aspiration Low deformation rate Aspiration High deformation rate Recovery Large deformation Large deformation Recovery (short holding time) Large deformation (long holding time) Aspiration Large deformation Recovery Small deformation Aspiration Low aspiration pressure Aspiration High aspiration pressure Aspiration High deformation rate Table I. A list of the published values of neutrophil cytoplasmic viscosity based on a fluid model under different experimental conditions. than when held for a longer time with an initial rapid recoil (Tran-Son-Tay et al., 1991). Finally, neutrophils that are only slightly deformed in a pipet before ejection, yield much smaller values of apparent viscosity (Hochmuth et al., 1993) than cells that are significantly deformed. This is analog to the finding that the cell appears to be more viscous when it flows into smaller diameter pipets (Evans and Yeung, 1989). It appears that the rate of deformation, the extent of deformation, the excess aspiration pressure and the holding time all affect the cell recovery characteristics, resulting in different apparent viscosities. Table I summarizes the apparent viscosity values for neutrophils re- 300 1000-2000 1000 2000 200-2000 800-1500 1000-2000 50-600 600 5000 500 1000-1500 ported by various groups. These results obtained from the aspiration (flow into a pipet) and recovery experiments. Therefore, these curious properties of the cell have yet to be fully explained by existing theoretical models based on the Newtonian or the non-Newtonian framework. 1.5. Current Status of Leukocyte Modeling Much effort has been done over the years on obtaining the correct rheological model for leukocytes. In earlier studies, Bagge et al. (1977) show that granulocytes deform as they flow down a tapered pipet. Transit times down the tube were on the order of seconds, for pressure drops of magnitudes found in-vivo. When expelled from the pipet, the granulocytes recover their original spherical shape in times ranging from 20 to 80 seconds. In the work of Schmid-Schonbein et al. (1981), the leukocyte is modeled as a simple viscoelastic solid with a Maxwell element in parallel with an elastic element, as shown in Figure 1.1. kl k2 Figure 1.1. Illustration of a standard viscoelastic solid model for the leukocyte. The model constants t and k's represent the viscous and elastic coefficients. The elastic element is thought to impart shape memory, which enables the cell to recoil to the resting shape. However, Evans and Kukan (1984) show that the deformation of a leukocyte into a pipet is a continuous flow process with no approach to a static deformation limit. It is therefore suggested that the leukocyte be modeled as a Newtonian liquid drop with a Figure 1.2. Illustration of a two-layer Newtonian liquid model with constant cortical tension for the leukocyte. y represents the isotropic surface tension coefficient. prestressed cortical tension (Yeung and Evans, 1989), as shown in Figure 1.2, rather than a viscoelastic solid. The shape memory then would be induced in much the same way as in drops, by ascribing a cortical tension (corresponding to surface tension in drops) to the membrane enclosing the cytoplasm. Based on this representation, Yeung and Evans (1989) simulate the flow of a neutrophil into a pipet by imposing a constant "ring" force at the pipet tip and a constant orifice pressure inside the pipet to satisfy the free slip condition. Although the liquid-like behavior of neutrophils is widely accepted, the precise nature of this liquid is not established. In their experiments, Needham and Hochmuth (1990) observed that the cell exhibits a non-Newtonian shear-thinning behavior when flowing into a pipet at high rates of deformation produced by large aspiration pressures. The non-Newtonian nature of the cell is also suggested by its contents, for example actin filaments, microtubules and other suspended particles, even though, their role may increase the cell viscosity only. It is clear from experimental observations, as summarized in the previous section, that the apparent viscosity changes with rate and extent of deformation. Thus, the behavior of neutrophils cannot be described in the same fashion as simple Newtonian liquid drops. These types of behaviors have been thought to indicate a "fading elastic memory" reminiscent of Maxwell liquids. Dong et al. (1988) use a Maxwell liquid model with a constant corti- Figure 1.3. Illustration of a two-layer Maxwell liquid model with constant cortical tension for the leukocyte. cal tension, as illustrated in Figure 1.3, to analyze the behavior of neutrophils flowing continuously into a pipet. The deviatoric stress t is obtained in such a model from: Tij -- J = 2 flij (1.1) where t and k are model constants controlling the viscous and elastic properties of the Maxwell liquid. The strain rate Tl is obtained from a vi avj i 2 x axj i (1.2) where v is the fluid velocity and x represents spatial coordinate. However, experimental data could be correlated with their model only if the material properties of the Maxwell liquid were continuously varied as the deformation proceeded. Hence, their model cannot be considered to yield a fundamental description of the leukocyte in terms of a Maxwell fluid. Another attempt to obtain a non-Newtonian description of the leukocyte is made by Tsai et al. (1993) through a power-law representation of the cell apparent viscosity. The expression used by these authors is of the form: t = ILO( Y-a)-b = ( (1.3) where R0 is the viscosity at a reference shear rate io. The constant b is obtained empirically. The power-law model essentially tries to link the shear rate dependent with the cytoplasmic viscosity to explain the shear-thinning behavior observed in the micropipet experiments. Although the power-law fluid model could be matched, for suitable initial conditions, to a fairly wide region of the deformation path of leukocytes, the rapid recoil region could not be predicted. 1.6. Summary Numerical study of leukocyte modeling has relied on methods based on semi-analytical approaches (Schmid-Schonbein et al., 1981; Tran-Son-Tay et al., 1991) or finite-elementbased structural analysis (Dong et al., 1988; Dong and Skalak, 1992). It is our opinion that the full scope of the numerical techniques available to the computational fluid dynamics community has not been brought to bear on the study of leukocyte modeling. It is clear from the review of current status of leukocyte modeling in the previous section that the existing viscoelastic solid, Newtonian and non-Newtonian liquid drop models fail to explain the complex rheological behavior of leukocytes. This is because such models do not take into account the internal structure of the cells. As Hochmuth et al. (1993) point out in the conclusion to their review of the state of affairs in leukocyte modeling, it seems that the correct model of the leukocyte can only be obtained when it is true to the morphology of the cell and its contents. However, since the internal structure of the cell is extremely complex, how can a morphologically consistent model be arrived at? One strategy, as suggested in the three-layer representation (Skalak et al., 1990; Dong et al., 1991; Hochmuth et al., 1993) shown in Figure 1.4, is to simplify the situation by factoring in the most important physical components, namely the membrane cortex, the cytoplasm and the nucleus and to temporarily ignore the contributions of the smaller intracellular particles including the microtubules and filaments. This is a good first order approximation since the cytoplasm is mainly composed of water. The cytoplasmic elements are assumed to only produce a higher cytoplasmic viscosity. layer (membrane) Inner layer (nucleus) Middle layer (cytoplasm) Figure 1.4. schematic of a three-layer or compound drop model for the leukocyte. In the three-layer model or compound drop under study here, the outer layer is the cortical region surrounding the cytoplasm. The cortex is comprised of an actin network. Its thickness is estimated to be of the order of 0.1 microns (Zhelev and Hochmuth, 1994) compared to the cell diameter of about 8 microns. Thus, the cortical layer can be considered to be a thin shell with constant tension, estimated at 0.024 to 0.038 mN/m (Evans and Yeung, 1989; Needham and Hochmuth, 1990; Zhelev and Hochmuth, 1994). The second layer, namely cytoplasm, is of a more complex nature and under the three-layer viewpoint, some scenarios have been advanced. It could be a two-fluid mixture, with a highly viscous liquid intermixed with a less viscous liquid. Osmotic stresses may exist between these two liquids. When the cell is deformed, the less viscous liquid would flow preferentially. The resulting osmotic stresses would tend to restore the equilibrium state between the two fluids. With regard to the third layer, the nucleus, researchers (Skalak et al, 1990; Dong et al., 1991; Hochmuth et al., 1993) suggest that it can be modeled as a highly viscous, solid-like component which retains its shape when the cell is deformed. Therefore, the leukocyte is modeled as a thin cortical shell with a persistent isotropic surface tension (outer surface) surrounding a thick layer of cytoplasm (shell) with a nucleus inside (core and inner interface). This representation is also known as a compound drop model. The layout of this thesis is as follows. In Chapter 2, the details of the salient features of the computational methodology employed are given. In Chapter 3, several well-defined and challenging flow problems are tested to assess the capabilities of the current numerical technique in terms of accuracy, robustness and flexibilities. The large deformation and recovery behavior of viscous compound drops are given in Chapter 4. This phase of the study elucidates the hydrodynamics of compound drops in relation to simple Newtonian drops and provides a basis for the development of a three-layer leukocyte model. In Chapter 5, the hydrodynamics of the three-layer model and its implications for the rheological behavior of the leukocytes are discussed. Finally, Chapter 6 summarizes the findings of this study and provides a list of suggestions for future research. CHAPTER 2 THEORETICAL ANALYSIS 2.1. Background The challenges posed by fluid flow problems involving moving boundaries have resulted in a variety of computational techniques directed toward their solution (Shyy et al., 1996). The moving boundaries may represent material discontinuities in the flow field, and their dynamics is coupled with that of the flow. Numerical methods applied to such problems are required to follow the evolution of the often complex shapes assumed by the moving boundaries. To tackle these evolving interfaces, as a natural extension of adaptive grid methods for stationary complex geometries, algorithms based on body-fitted coordinates have been applied to moving boundary problems (Kang and Leal, 1987; Glimm et al., 1988). A boundary-conforming grid arrangement is convenient and accurate; the discontinuity across the interface is always maintained and interfacial conditions are applied at the exact locations without any smearing or redistribution. However, when the interface deformation becomes large, it is difficult to adequately resolve the geometrical complexities while maintaining desired mesh control. Mesh skewness and stretching may impact negatively on accuracy and efficiency of the flow solver. The boundary-fitted grid method presents further difficulties when interfaces merge or fragment. In simulating multi-layer fluid phenomena involving highly deformed interface shapes, purely Eulerian methods (Sethian, 1990; Gunstensen et al., 1991; Kothe and Mjolsness, 1992; Grunau et al., 1993; Sussman et al., 1994) have come to be widely used. In such methods, a fixed Cartesian grid is usually used, so problems associated with grid generation are circumvented. Eulerian methods based on the Volume of Fluid (Hirt and Nichols, 1981), Level-Set (Sethian, 1990) or Lattice Boltzmann Equation (Grunau et al., 1993) have been developed to handle free-surface flows, such as wave-breaking, sloshing and splashing type problems, and multi-phase fluid flows in complex geometries. In purely Eulerian methods, the details of the interface shape are not explicitly needed, since the interface is deduced based on a computed field variable. The interface information can be obtained from an appropriate computed field variable in order to describe the details of the interface shape and curvature. An accurate estimation of the interface curvature is very important when capillarity plays a significant role in the interfacial physics. In tracking highly distorted interfaces with good accuracy, a combination of the strengths of Eulerian and Lagrangian methods has been shown to be useful. This mixed Eulerian-Lagrangian approach has been used in recent years for problems involving fluid-fluid (Unverdi and Tryggvason, 1992; Nobari et al., 1996) as well as solid-liquid interactions (Udaykumar and Shyy, 1996; Juric and Tryggvason, 1996). These methods employ fixed structured grids, and afford simplicity and availability of well-tested flow solvers. The only moving component is the interface which is a lower-dimensional surface overlying the fixed grid and is explicitly tracked. One limitation of the methods based on fixed grids is that since they have traditionally employed Cartesian grids, complex flow geometries are difficult to incorporate. A complex geometry facility is required for studying the interaction of the moving boundary with flows in arbitrary domains. The representation of a solid boundary on a fixed Cartesian grid layout calls for special treatment, since the control volumes through which the irregular boundary passes become fragmented into solid and liquid regions. This situation is illustrated in Figure 2.1. Quirk (1992) details some of the procedures involved in dealing with cell fragments or cut cells in the framework of compressible, inviscid flows around stationary obstacles. However, the calculation of fluxes is not presented in detail in that work. Also in connection with compressible, inviscid flows, Pember et al. (1995) use local refinement and a redistribution procedure to enforce conservation in the vicinity of the solid-fluid boundary. Arbitrary-shaped boundaries running through a Cartesian mesh have also been dealt with using Interfacial cell 0 o o 0' Solid-liquid interface 0 0 0 0 W P E o oS Bulk cell Control point (i,j) 1 2 Figure 2.1. Illustrations of (a) the situation when a solid object in represented on a Cartesian grid, (b) a typical control volume encountered in the vicinity of the solid-liquid interface (Interfacial cell). finite element methods (Young et al., 1991). In Zeeuw and Powell (1990) and Bayyuk et al. (1993), a Cartesian grid is employed to track the motion of solid objects through an inviscid compressible fluid. Goldstein et al. (1993, 1995) employ a modified immersed boundary technique (Peskin, 1977) to impose no-slip boundary conditions over arbitrary-shaped solid boundaries. This facility is useful for complex geometries, an important consideration for the spectral methods in use therein. Turbulent flows around solid obstacles have been calculated using this approach and the results have been shown to compare favorably with experiment. The method, as used in these works, treats the solid boundary as a series of steps conforming with the Cartesian grid layout, i.e., in piecewise constant fashion. In contrast, a piecewise linear representation of the free surface is used by Miyata (1986) for the simulation of wave breaking; the calculations are performed on a Cartesian grid using a finite-difference technique. In a finite volume setting, Udaykumar and Shyy (1995a) use a cut-cell technique to compute incompressible, viscous flows in arbitrary geometries. In that work, it is demonstrated that by means of a suitably defined consistent mosaic of control volumes, it is possible to obtain accurate solutions for flows on Cartesian meshes. However, the staggered variable layout in Udaykumar and Shyy (1 995a) renders the measures adopted to handle the internal boundaries rather tedious. In this work, we present a fixed grid methodology for the simulation of flow through complex geometries with fluid-fluid moving boundaries, considerably simplified by the choice of a collocated variable layout (Lien and Leschziner, 1994). The performance of the flow solver in terms of stability and accuracy is thoroughly tested for various flow conditions. 2.2. ELAFINT: A Mixed Eulerian-Lagrangian Solution Algorithm 2.2.1. Characteristics of the Numerical Methodology There are two distinct aspects to the computational framework, called ELAFINT, which stands for Eulerian-Lagrangian Algorithm for INterface Tracking, reported in Udaykumar et al. (1995a, 1996). These are: (1) A pressure-based flow solver to simulate the incompressible fluid flow over a wide range of Reynolds numbers. (2) A procedure to deal with the interaction of fixed arbitrary geometries and evolving moving boundaries with the flow. The first aspect is straightforward and approached in a conventional manner using the robust pressure-based solution procedure detailed by Patankar (1980) and in the context of more modem developments by Shyy (1994). The scheme to track the moving fronts and describe complex geometries is summarized below, following which the flow solver and its discretization procedure are described. 2.2.2. Interface Tracking Procedure The interfaces are described by means of markers indexed sequentially and separated into objects. There is no restriction on the number of objects or on the end conditions that can be imposed on these curves. As shown in Figure 2.2, the objects can be open or closed. Furthermore, they can be fixed or moving. They can enclose different phases or materials with different properties. Once the interface description is obtained in terms of markers and their positions, the shape characteristics such as normals and curvatures can be extracted. The convention adopted by the algorithm is to consider the normal to the interface as pointing from phase II to phase I, as shown in Figure 2.2, such that when one traverses the interface along the marker sequence, phase II lies to the right. In order to handle multiplej+1 j+2 PHASE II P n PHASE I i i+l i+2 PHASE II Figure 2.2. Illustration of objects, markers and normal convention in the description of interfaces. valued interfaces, the interface tracking algorithm employs piecewise polynomial fits to compute the necessary shape derivatives at each marker location. Facility exists to construct quadratic or cubic-spline functions xi(s) and yi(s). The polynomial fits are obtained by parametrizing the curve as a function of arc length s so that, for the ith marker Xi = xi(s), yi = Yi(s) (2.1) Then, the normal and curvature are obtained by taking appropriate derivatives of these piecewise polynomials. Thus, the components of the normal are, Ys xsnx (xs2 + ys2)1/2 ' (Xs2 + ys2)1/2 (2.2) The curvature x is, for the two-dimensional case: S^ ysXss - YssXs x= V-n= (ys2 + xs2)3/2 (2.3) and the axisymmetric case: Xs + YsXss - YssXs x = + Y(Ys2 + Xs2)3/2 (Ys2 + Xs2)3/2 (2.4) Interfacial markers are advected to their new positions by a Lagrangian translation. Thus, x 1 xk u (2.5) yek+l -- y f + At ve(25 Yk+ =Yek + At v (2.6) where the subscript f indicates the interfacial marker number and the superscript k implies the time level. It is demonstrated (Udaykumar and Shyy, 1995a, b; Shyy et al., 1996) that the interface tracking procedure can handle extreme deformation, including topological changes such as repeated mergers and breakups. It is also established that the interface tracking facility is versatile and accurate. Once the objects are described on the computational domain, the interaction with the flow field needs to be established. This is mainly done by providing the following information: (1) The cell in which each interface marker lies. (2) The markers contained within each interfacial cell (i.e., control volume through which the interface passes). Since a Cartesian grid is employed, establishing these relationships is a straightforward task. Based on the information mentioned above, the flow solver can incorporate either changes in the control volume definitions (as in the cut-cell approach) or appropriate source term contributions (as in the immersed boundary technique). These two approaches enable the effects of the internal stationary as well as moving boundaries to be communicated with the flow field. In turn, based on the flow field developed, the motion of the interface can be determined. In Udaykumar and Shyy (1995b), moving as well as stationary solid-liquid boundaries are handled by the cut-cell technique. In this work, the cut-cell technique is retained for the stationary solid-liquid boundaries, while the immersed boundary technique (Peskin, 1977; Fauci and Peskin, 1988; Unverdi and Tryggvason, 1992) is used to treat moving fluid-fluid boundaries. This, along with the use of a collocated variable arrangement, enhanced the robustness and simplicity of the algorithm. Extensive testing of the algorithm has been performed for a wide variety of flows and the results thereof are presented elsewhere (Shyy et al., 1996; Udaykumar et al., 1996). 2.2.3. Formulation of the Pressure-Based Flow Solver Consider the two-dimensional planar or axisymmetric, incompressible Navier-Stokes flow. The governing equations in conservation forms are: Continuity equation: -+ a(gu) + y 1 a yv) = 0 aBt ax ym a (2.7) Momentum equations: a(u) + - (guu) + (ymgVU) at ax y ay ap + Ou I a--M au. =- + (Y _) + y ' + S x ax x y (2.8) a(Qv) a (v) + -(Quv) + '1 ( vv) at ax - y my ap my a av 1 a m av = - - + ( ) + (Y mIpT) + SV (2.9) Energy equation: a(Q7) aa aTm1 aanaaT + (0uT) + 1 8 m - T (U) 1 . Y T) at axym ay a v) axax yma yay (2.10) where m=O for planar and 1 for axisymmetric case. In the above, g is the fluid density, u and v are the fluid velocity components, T is the temperature, p is the pressure and j is the coefficient of viscosity and Su and S, are the source terms in x- and y-momentum equations which includes, for example, the buoyancy term under the Boussinesq approximation and the surface tension contributions discussed later. The boundary conditions to be applied on the interface for cases involving no mass exchange between phases are the continuity condition and normal stress balance: (Vn) = (Vn)2 = (Vn)1 (2.11) P2P PI + 2 2- 1U I1 (2.12) where the subscripts 1 and 2 represent the two adjoining fluids and I represents the interface, Vn is the normal velocity component, y is the surface tension and X is the curvature of the interface, with the normal pointing from fluid 2 into fluid 1, as shown in Figure 2.3. normal to interface Underlying fixed . Cartesian grid nx, Li uid normal to solid surface Liquid 1 Figure 2.3. Ilustration of the interaction of the solid-liquid and liquid-liquid boundaries with Cartesian grid. The solid-liquid is treated using cut-cell technique, while the liquid-liquid interface is modeled with the immersed boundary technique. 2.2.4. Discretization Procedure of the Flow Field Equations The discretization of the governing equations is carried out using a control volume formulation, which with particular reference to two-dimensional planar flow. Similar treatment applies for the axisymmetric case. Consider the conservation law for the variable 0, defined to be (a) 1 for the continuity equation (b) u and v for the x- and y- momentum equations and (c) specific enthalpy (h= Cp, T) for the energy equation. In the general case, we have at + V - (i) = (V - IP) + S (2.13) The control volumes in the vicinity of a solid boundary passing through the grid are irregularly shaped, and in the general case can assume a five-sided shape as shown in Figure 2.1(b). In order to evaluate the fluxes through the cell faces, we integrate Eq. (2.13) over the control volume and employ the divergence theorem for two-dimensional problems to obtain dA + (9 i) * Rdl = (IV#) * Edl + SdA A I 1 A (2.14) where the second term on the left hand side is now the line integral of the outward normal convective fluxes and the first term on the right side is the line integral of the outward normal diffusion fluxes through the faces of the control volume. We now proceed to discretize each of the terms in Eq. (2.14) for the control volume shown. Thus Qk+ l0k+l I k k ag = Aij A f t at A (2.15) where the superscript k, k+1 indicate the time levels and 6t is the time step size. Ac is the area of the irregular control volume. The time derivative is discretized in the backward Euler form leading to an implicit-in-time solution procedure which relaxes the time step constraints in comparison with explicit schemes. Next, 5 (o)*dl = (Qeu~pe)k+ dle e=1 (2.16) is the summation of the convective normal effluxes through each cell face of the control volume. The superscript k+1 indicates the implicit nature of the scheme. Following Shyy et al. (1996) and Udaykumar et al. (1996), the diffusion fluxes are computed as 1- 5 (IPV) - dl = (Fe,( ) )k+dle Se 1 (2.17) Let us denote the source term as follows, SdA = A (2.18) Here S contains contributions from the pressure and body force terms. Substituting the above Eqs. (2.15-2.17) in Eq. (2.14) one obtains the discretized form as qktk+1 k+1 CJk 5 5 A v6t + (9uf )k+Idle = (F ) dle e = 1 e= 1 (2.19) In particular the discrete form of the continuity equation can be written as ek+ I .k ) 5 i Ac, + (un)edle = 0 6t e = 1 (2.20) Now multiplying Eq. (2.20) by the value ij and subtracting from Eq. (2.19) gives Okktl + 0.)Ac 5 i --idt- AY + Que(o - i)k+dl = Z(F(n)e)k+ldle + 6t Ac+In e=i e=I (2.21) In the notation of the pressure-based methodology (Patankar, 1980; Shyy, 1994), we cast the final discretized form as Apop = aAp/N + asos + agfE + aw + b (2.22) where b = S+ rTn)dl', - Qlundli5 - P) 1 (2.23) Thus, a five-point stencil is adopted, and the value of the variable q at point P is dependent on its four adjoining neighbors. The form of the coefficients aEW,N,S depends on the convection scheme employed. In this work, we employ the second-order central difference scheme for the convection term. Diffusion and pressure terms are also discretized using central differences. When cast in the matrix form the above equation reads [C] [0] = [B] (2.24) where [C] is a pentadiagonal coefficient matrix, P is the solution vector and B is a source vector. The procedure usually adopted is to solve the system of equations in iterative fashion employing a line SOR method which calls for a tridiagonal matrix solution procedure. A pressure-correction procedure provides the pressure field consistent with mass conservation. In the pressure-based methodology, it is common practice to use a staggered layout of variables. However, here we use a collocated arrangement of the dependent variables. This calls for special measures, called momentum interpolation (Rhie and Chow, 1983; Lien and Leschziner, 1994), when estimating cell face velocities. The details of momentum interpolation will be discussed next. 2.3. Momentum Interpolation One of the key issues in using a pressure-based solver is the grid arrangement. The staggered grid layout has been preferred for incompressible flow computations for many decades due to its many advantages (Patankar, 1980; Shyy, 1994). However, the collocated arrangement has desirable attributes in terms of simplicity and economy, has come into common use. For example, in the present work, in dealing with the internal boundaries, whether by the cut-cell or the immersed boundary technique, the use of a collocated variable layout leads to the treatment of a single set of control volumes for all the dependent variables (namely, u, v, p and T). In the case of a staggered system, three sets of control volumes, for u, v and the scalars need to be considered and the data structures associated with the cut-cell as well as immersed boundary technique need to be computed and stored. A particular instance is the identification of the phase of the control volume. When the staggered layout is used the phase in which the control points for the u, v and scalar variables lie needs to be determined separately, while in the collocated layout only one set of control volumes needs to be handled. This leads to a considerable savings both in terms of the programming effort and also in the computational effort expended in obtaining the information regarding the interaction of interfaces with the flow field. These two variable arrangements are illustrated in Figure 2.4. One potential difficulty of the collocated grid is the occurrence of unrealistic oscillations in the pressure and velocity fields, leading to numerical instabilities for certain flow problems (Patankar, 1980). This is due to the fact that the collocated grid arrangement leads to decoupling of pressure and velocity fields when the central difference operator is employed in discretizing gradients in the flow field. However, with the incorporation of momentum interpolation, introduced by Rhie and Chow (1983), it is possible to prevent the onset of checkboard phenomena resulting from the decoupling between the pressure and velocity fields. Lien and Leschziner (1994) have applied this measure to obtain solutions for a wide range of flow problems, including shocked flows. In our implementation, we use the version of"momentum interpolation", recommended by Lien and Leschziner (1994) to evaluate the velocities at the faces of the control volume. Essentially, momentum interpolation introduces a third-order dissipation in estimating the cell face velocity. This can be seen as follows. The discretized momentum equations for up and uE based on finite volume formulation, presented for a general transport variable 9 in Eq. (2.13) are given in the one-dimensional case by Apup = anbUnb + Sp - (Pe - Pw)p nb=E,W (2.25) AEUE = anbUnb + SE - (Pe -PW)E nb= EE,P (2.26) where subscript nb implies neighboring points in the stencil, Sp and SE represent the source terms, A = anb and AE = anb. The relative positions of the points subnb=E,W nb = EE,P scripted P, E, W, EE and WW are shown in the Figure 2.5. u - control volume p - control volume A v - control volume u, v and p - control volume -I > 0> 0 > 0 A T > 0 I 0 A -4 > 0 >- > 0 A Figure 2.4. Schematics of (a) staggered grid layout. (b) collocated grid layout. - - . - I - ww w e ee uw Up Ue L WW W P E EE "wT "wT %-"lTB TEE Figure 2.5. Schematic of 1-D collocated grid layout. Equation (2.25) and Eq. (2.26) can be rewritten as Hp 1 up - Hp A(Pe - Pw)p Ap AP (2.27) HE 1 UE _ E 1 (Pe - PW)E AE AE (2.28) where Hp = anbunb + Sp and HE anbUnb + SE. In performing the internb= EW nb = EE,P polation, it is assumed that He 1 (Hp HE Ae 2 AP AE) (2.29) Then, substituting Eq. (2.27) and Eq. (2.28) into Eq. (2.39), the interpolation for the face velocity Ue can be expressed as Ue = up - I(Pw - Pe)p + UE A(Pw - Pe)E + P - PE) P AE AE (2.30) or alternatively, Ue = (Up + UE) + {!P - PE) - 1 (pw - Pe)P + Pw - Pe)E (2.31) In Eq. (2.31) the face velocity of control volume is directly linked to the two adjacent cell-center pressures through the term I above. This ensures that strong coupling between the pressure and velocity fields is retained. The second term, indicated as II, is the pressure smoothing term, which eliminates the non-physical wavy pressure or velocity fields. In fact, it can be shown that Eq. (2.31) can be recast, following some manipulation (Lien and Leschziner 1994), as Ue = "(Up + UE) + oX3 (2.32) This brings to light the fact that momentum interpolation adds a third-order term to the cell face velocity, relative to straightforward averaging of adjacent cell-center velocities. Thus spurious pressure oscillations are eliminated. Typically, iterative procedures are used to solve the discretized momentum equations, and under-relaxation can be implemented in Eq. (2.27) in the form up = a up + (1 - a) up (2.33) where a is the under-relaxation factor, up is the current iterative value, u* is the value of up at the previous iteration, while up is the value obtained after relaxation. Similarly, UE = a uE + (1 - a) uE (2.34) Therefore, the cell face velocity, after including the relaxation factor is Ue = a ue + (1 -a) ue = a b(Up + uE) + _(PP - PE) - (Pw - Pe)P + (Pw - Pe)E 1 EA 2[1Ap" AEE] + (1 - a) Ue (2.35) The momentum interpolation as obtained from Eq. (2.31) yields flow fields independent of the relaxation factor, thus eliminating the relaxation factor dependency detected in early versions of momentum interpolation (Majumdar, 1988). From Eqs. (2.33-2.34), Up and uE can be rewritten as up= a a u* (2.36) uE 1-a u* UE a a u (2.37) Thus, Eq. (2.35) becomes ï¿½* 1 *- ** [ * l1 * *] ue = j~U(* + u *) + (1 - ) e - (Up + uE) + a { (PP - PE) - (Pw - Pe)p + L(Pw -Pe)E] 2Ap AE (2.38) Equation (2.38) is the form of the momentum interpolation suitable for steady flow problems. The extension of Eq. (2.38) for unsteady flow problems will include the time-dependent term (9Axuo/A t) where superscript o indicates the previous time value. The appropriate form is S I Ap ** AE ** * 1(Ap * + AE * UAe = P E l ) e 2 Ae UP + Ae UE a_ oï¿½Ax o I Goï¿½Ax o oï¿½Ax uo] + Ef"XeuO [() )uO + u' ) +Ae At ]ee 2 t *u+ At ]EE ee P2 + 7(Pp - PE) - 4tPw - Pe)p + Pw - Pe)Ei}A (2.39) Although Eq. (2.39) applies to a I-D, uniform grid situation, it can be easily extended to higher dimensions and non-uniform grids. 2.4. Interactions between Interfacial Physics and Flow Field Variables 2.4.1. Cut-Cell Procedures for Solid-Liquid Boundaries In accounting for the presence of the internal solid boundaries which define complex geometries on Cartesian grids, the cut-cell approach developed by Udaykumar and Shyy (1995a) is applied. In this approach, when the interface passes through the grid, the scenario shown in Figure 2.1 is encountered. The control volume designated by P is shown to be cut by the solid surface passing through the grid. The cut-cell procedure rearranges the control volumes in the vicinity of the interface to form a fully consistent mosaic of cells. This is obtained by redefining the interfacial cells (i.e., the control volumes in the vicinity of the interface), so that the control volume fragments lying in each phase are absorbed into appropriate partner cells (Udaykumar and Shyy, 1996). Figure 2.6. Illustration of the situation resulting when a complex solid geometry is treated on a fixed Cartesian grid. The cut interfacial cells are reassembled to provide a closed, consistent mosaic to maintain flux conservation. The assembly of cells is performed by absorbing fragments of the control volumes into appropriate cells to reconfigure the cells in the vicinity of the boundary as shown by the dotted lines. Also shown in the normal probe for obtaining gradients at the interface. An example is shown in Figure 2.6, where the dotted lines indicate the newly defined control volumes. The fluxes through these redefined control volume faces are thereby evaluated in such a way as to maintain conservation and consistency of the fluxes across each control volume face. For example, the flux of a conserved variable is computed for a cell that is abutted by two neighbors to the left (as in the case of the hatched cell in Figure 2.6), such that the flux through the west face of the hatched cell is the sum of the fluxes through the east faces of the abutting cells. This procedure obviates any smearing of information at the interface since the boundary condition is imposed at the exact location of the internal boundary in each control volume. This is in contrast to the means employed by Goldstein et al. (1993) in imposing the no-slip condition. In that work, the no-slip boundary condition is imposed via an external concentrated force field which is computed based on the developing flow field in the vicinity of the interface. Thus, a feedback control approach is used to compute the shear stresses necessary to obtain a no-slip condition. Information regarding the no-slip surface is redistributed over the finite bound of the 6 - function used in the immersed boundary representation of the solid surface. In this work, for the evaluation of the source term b in Eq. (2.23), the diffusion flux from the interface is obtained by describing the field O(x,y) for any transport variable such as velocity or temperature. This is done in each phase separately around the interface by means of a biquadratic function and the gradient normal to the interface required in obtaining the diffusion flux is obtained by extending a normal probe from the interface into each phase. The length of the normal is one grid spacing h. This situation is illustrated for the liquid phase in Figure 2.6. The value ol+dn at the end of the probe is obtained from the biquadratic function representation in the required phase as O(x,y) = ax2+by2+cxy+dx+ey+f. The coefficients are obtained by inverting the 6x6 matrix obtained by the known values of at six suitably chosen points in the required phase around the end of the probe. The inversion of the matrix is performed by LU decomposition. Thus, the required normal gradient of the field O(x,y) is computed as follows: o_ 021+dn I an dn (2.40) where dn is the normal probe length. The pressure at the solid boundary is obtained by computing the pressure in a similar manner to the above at the reference point and extending the BP same to the interface, which implies that = 0 at the solid boundaries. In the case of the an stationary solid boundaries in this work, the convective interfacial fluxes do not contribute to the momentum. However, for moving solid boundaries, the treatment is given in Udaykumar et al. (1996). The procedures for the treatment of the internal solid boundaries lead to stable, accurate computations of flows through complex geometries. 2.4.2. Immersed Boundary Technique for Fluid-Fluid Interfaces When the fluid-fluid interface (i.e., the moving boundary) lies on a fixed grid, some means of communication between the interface and the Eulerian framework for the field equations needs to be devised. One such technique is the cut-cell method described above and in more detail in Udaykumar et al. (1996). In this work, we adopt the immersed boundary technique originated by Peskin (1977), and employed extensively by Tryggvason and coworkers (Unverdi and Tryggvason, 1992; Juric and Tryggvason, 1996; Nobari et al., 1996) for a variety of two-phase flow problems. In the immersed boundary technique, the presence of the interface and its effect on the flow field is established by means of source terms in the governing equations, i.e., the information from the Lagarangian interface component is redistributed to the Eulerian framework. For the particular case of two-fluid interfaces under constant interfacial tension, the source term represents the normal stress contribution from the interface. The momentum jump conditions, given by Eqs. (2.11-2.12) at the interface, is incorporated in the x- and ymomentum equations. Specifically, the contributions to the normal stress balance from the capillarity term in Eq. (2.12) is redistributed to the surrounding grid points; it is supplied to the source term Su and S, in Eqs. (2.8-2.9). This is done with the aid of a 6 function of bound 4h as illustrated in Figure 2.7 (where h is the grid spacing). This procedure leads to the conservative redistribution of the source terms to the grid, due to the properties of the delta func- / . Transition area I I I I I / :II I I D(x) I I I I I I I I I I I I I I I I I I I I'* ' ' I I -1 I I I 4h Figure 2.7. Illustration of the form of the delta function in the vicinity of the interface. h is the grid spacing. The bound of the delta function is shown. tion. The discrete delta function at each point in the flow field assumes the form originally proposed by Peskin (1977): 2 ({ 1( + cosh(xi - (xi)e) if I - (x-)el ER 0, if I - (x-)el R (2.41) (2.41) where Y and (x-)e are the locations of points in the flow field and on the interface, respectively, ( represents the interfacial marker index. R is set equal to 2h in our computations. Each marker contributes to the momentum source term in its proximal control points, i.e., those lying in the circle of radius 2h centered at the marker. The control points experiencing the influence of the capillary forces are marked in Figure 2.8 by closed dots. The force contribution from each interfacial marker to the momentum equations is at any point P given by markers contributing to S. radius of circle each marker = radius of circle = 2h around point P Figure 2.8. Illustration of the form of the delta function in the vicinity of the interface. h is the grid spacing. The bound of the delta function is shown. fpf -= Y"efleASf (2.42) where A Se is the arc length corresponding to the (th marker. Thus, the source term in the momentum equations at a given point in the flow field depends on the proximity of that point to the fluid-fluid interface. The magnitude of the interfacial forces experienced is dictated by the delta function and is the sum of the contributions of all the interfacial markers lying within the distance of 2h from the point. Therefore, the normal stress contribution to the momentum equation is given by the expression, Fp(x- = (- (x)fpe (2.43) This procedure of redistribution of interfacial stresses as sources to the field equations leads to the smoothing of the discontinuity at the interface. The stability and robustness of this method is enhanced by such smoothing. The continuity condition at the interface for material interfaces, Eq. (2.11), is enforced by extracting the normal velocity of the interfacial markers from the grid information by means of a delta function. The 6 function form is j (2.44) where the index corresponds to the proximal (in the region given by R=2h) grid points. The interface is advected with this velocity. Another feature that needs to be mentioned is the handling of property jumps across the interface. In other approaches using a fixed grid (Kothe et al., 1992; Sussman et al., 1994) some field variable is employed to smooth the property jumps across the interface. For example, a Heaviside function can be employed such that: 1, q5< -a 0, 0> +a S+ sin otherwise / (2.45) where a is the bound region around the interface which is equal to 2h in our computations. Usually this bound is taken to be proportional to the grid size. Based on the Heaviside function, the material properties are assigned as : P =1 (,02 - Pi)G(O)a (2.46) where P is any material property such as viscosity and density. The subscripts 2 and 1 indicate the two adjoining phases, in line with our convention for the normal. In the RIPPLE approach (Kothe et al., 1992), which is based on a VOF representation of the phase field, q is taken as a color function dependent on the fluid fraction field, while in the level-set algorithm the t value corresponds to a distance function. In these methods the field variable 95 is an essential part of the interface tracking procedure. In the immersed boundary technique, Unverdi and Tryggvason (1992) employ an indicator function which is obtained by solving a Poisson equation for the entire domain. The indicator function has the characteristics of the Heaviside function but has to be computed every time the interface moves. Thus, the calculation of an indicator function involves the solution of an additional field equation and does not fully exploit the available information from the explicitly tracked interface. In our algorithm we have exploited the Heaviside representation to achieve the property smoothing, based on explicit knowledge of the interface and the information indicated in Section 2.2.1. Once the marker locations with respect to the grid are known, the property jumps can be assigned based on the Heaviside function. In the present case, the Heaviside function takes the form: .(x - (xi)e + -sin-(xi - (xi)e) if 1Y - (x)el E a 1, sgn(- - (x)e)I - (x)e1 < - a 0, sgn(- - (x)e)If - (x-)el > + a (2.47) (2.47) where xi represents the coordinate of the control point where the Heaviside function is required. The smoothing of the material properties is performed with Eq. (2.47) above. The material properties in the bulk of the phases are assigned first. Thereafter the transition region between any two phases is smoothed by defining the Heaviside function locally in this region. To accomplish this, one traverses along the interface and modifies the value of the material property according to Eq. (2.47), in the cells adjacent to the interface and lying within the bound a of the interface. These operations are confined to the region close to the interface only, hence reducing the computational burden. This procedure for smoothing material discontinuities holds when there are multiple interfaces, including cases in which interfaces may enclose other interfaces. It is noted that computational difficulties such as spurious oscillations and numerical instabilities are often encountered when the property gradient across the interface is steep. By applying the Heaviside function to smooth the discontinuities across the interface, the choice of bound region a needs to be made with care. A large a can cause the solution to be overly smeared and to substantially lose physical fidelity while a small a will be insufficient to stabilize the numerical solution. Another issue is the selection of the smoothing function. Instead of choosing the smoothing function which is linearly proportional to the mesh size, a higher order smoothing function can be used (Grunau et al., 1993). 2.5. Mass Conservation in Individual Phases In the scheme presented above, the interaction of the interface tracking algorithm as provided by the immersed boundary technique, with the flow solver does not enforce conservation of mass in the phase enclosed by the interface. This is because the divergence-free velocity field is enforced in each control volume by the flow solver, the velocity employed to evolve the interface is interpolated from such a velocity field. This interpolated velocity field may violate the divergence-free condition, resulting in a mass leak from the phase enclosed by the interface. Peskin and Printz (1993) suggest a method based on the modification of the divergence operator in order to enforce mass conservation. In this work, we choose to achieve the same result by enforcing mass conservation in our iterative solution process so that the correct level of hydrostatic pressure is induced in the enclosed phase. This is done in the following way. In order to satisfy the mass constraint we require, in each enclosed phase: mbubble = pI ds = constant 6 (2.48) where Hy is ;ry2 for axisymmetric and y for planar case, g is the density of the fluid inside the enclosed mass, such as a drop, s is the arc length and snax is the total arc length of the curve representing the boundary of the enclosed mass. In Eq. (2.48), a bubble is used as the example; obviously, the same procedure is applicable to other cases, such as a drop. Since the pressure correction formulation only provides an estimate of the pressure gradients necessary to satisfy the continuity equation cellwise (i.e., across control volume faces), global conservation information is not provided by this means. Thus, the continuity constraint as expressed by the pressure correction equation establishes the hydrodynamic part of the pressure inside the enclosed mass consistent with the flow field, but provides no information regarding the hydrostatic component. The dynamic form of the Young-Laplace equation at the interface establishes the normal stress balance as: A (Pd+ Pc) + A (viscous terms) = y(XI + x2) (2.49) where A represents the jump across the interface and the viscous terms take the appropriate forms in planar and axisymmetric cases, as given in Eqs. (2.8-2.9). In the above equation, subscript d indicates the hydrodynamic and c the hydrostatic parts of the pressure. The hydrostatic pressure must be such that it satisfies the mass constraint. Kang and Leal (1987) give the appropriate expression for pc term as a time-dependent parameter C(t). This constant arises when the pressure at the interface is obtained by integrating the component of the momentum equations in the direction tangential to the interface in each phase. This is easier to see in the case where the flow on each side of the drop is obtained separately and then matched at the interface through the jump condition Eq. (2.12), as is the practice in the boundary-fitted formulation employed by Kang and Leal (1987). In related work, the practice of rescaling the interface shape periodically to maintain the enclosed mass constant is adopted (Stone and Leal, 1989a). In the work presented by Nobari et al. (1996), it is mentioned that the mass of drops is maintained, but the mechanism by which this is established by their solver is not clear. In recent work on the dynamics of menisci in crystal growth, Shyy and Rao (1995), with the pressure-correction methodology used here, maintain the volume of the meniscus by employing an iterative procedure to obtain the proper value of the hydrostatic pressure, consistent with Eq. (2.49) as well as mass conservation. However, the computations use the boundary-fitted technique and the implementation of the correction of mass conservation need to be modified in the present context. A method for correcting the hydrostatic pressure in the enclosed mass needs to be devised so that it does not lose or gain mass. To achieve this goal, we devise an indirect method of establishing mass conservation similar to that adopted by Kang and Leal (1987). The desired result in enforcing mass conservation is to establish an additive constant to the pressure inside the enclosed mass. The stress balance equation in Eq. (2.49) then becomes: A (Pd + Pc + 6Pc) + A (viscous terms) = y(xI + X2 + 6X) (2.50) The correction to the pressure, 6pc, is responsible for changes in the drop shape represented by 6x so that upon integrating as in Eq. (2.49) over this new shape, the volume is conserved. In our procedure, during the iterative process of establishing the flow field, we impose the mass constraint on the interface motion by establishing the shape correction 6x and inducing the pressure correct increment 6Pc. At each iteration, this shape correction is performed after the interface has been advected to its new position based on the kinematic condition. Thus the fresh iterate of the x-coordinate of interfacial marker position is obtained by (XeNk+ 1) = X k + 6tVnnx (XeN,k+ )c = (Xe N,k+ I)p + 6nnx (2.51) where the subscripts p and c represent the predicted and corrected values, the subscript f is the interfacial marker index, the superscripts N and k represent the iteration number and time step, respectively. Similar expressions are obtained for the y-coordinate. The shape correction in the normal direction 6n is performed in an iterative fashion by means of a bisection method, so that the mass of the drop as computed in Eq. (2.48) remains a constant. The value of 6,n is obtained from the knowledge of the arc length measure of the interface and the estimate of the mass loss. The bisection method used to obtain the value of 6,n converges in a few iteration (typically < 5). The entire mass conservation measure is enclosed in the outer iteration loop for the entire solution procedure, including interface movement, flow field computation and mass constraint imposition. The procedure is then required to achieve a preset tolerance for convergence levels of the flow field, interface shape and enclosed mass. It should be emphasized that the present procedure is of a predictor-corrector type, and, at convergence should yield a drop shape independent of the manner 6Pc is imposed. 2.6. Summary In this chapter we have described the fundamental elements of the computational methodology called ELAFINT. The numerical framework developed here have several improvements over the previous version of the code that can be summarized as follows: (1). A collocated variable approach is adopted. The stability and accuracy of the algorithm in computing flows of comparable complexity to those handled by the staggered grid approach was established. The collocated layout considerably simplifies the algorithm. (2). The cut cell technique is facilitated by the collocated grid layout, and is rendered less tedious. While previously the cut cell technique had to contend with three types of control volumes, in the present case only one control volume is used. (3). The fluid-fluid interface is modeled using the immersed boundary technique. The cutcell approach was used in previous work on solidification for the moving boundary calculations. For fluid-fluid interfaces however, this method proves to be sensitive to noise and difficult to stabilize. This is because the cut-cell method employs no dissipation, so that the boundary condition application is exact. In other moving boundary treatments based on a fixed grid, some smearing is added near the discontinuity. This is the case in the immersed boundary technique (Unverdi and Tryggvason, 1992), the VOF method (Kothe et al., 1992) and the level-set method (Sussman et al., 1992). The choice of the immersed boundary technique is made due to its compatibility with explicit interface tracking where ambiguity in regard to interface location is avoid. CHAPTER 3 ASSESSMENT OF COMPUTATIONAL TECHNIQUE In the preceding chapter, the elements of the numerical framework developed for handling moving boundary problems involving fluid-solid interaction as well as multi-fluid flows are described. In this chapter, the focus is to illustrate the accuracy and versatility afforded by the combination of the cut-cell method and the immersed boundary technique. Towards this objective, we present several designated cases to evaluate the performance of the cut-cell method in the current collocated-variable framework. This cut-cell approach was tested in previous work (Shyy et al., 1996; Udaykumar et al., 1996) for a staggered variable arrangement. Results given by the immersed boundary technique in the present finitevolume, pressure-based framework are compared to known solutions of viscous drop dynamics. For this particular application, reliable experimental and numerical results exist, at least for highly viscous flow, with which we compare our results. Finally, we put the cut-cell and immersed boundary representations together to enable the simulation of drop dynamics in a constricted tube. This last aspect demonstrates the ability of the present method to handle multi-fluid flows in arbitrary geometries under viscosity- or inertia-dominated flow conditions. 3.1. Flows With Fixed Arbitrary-Shaped Internal Boundaries One objective of this work is to develop a facility for computing complex geometry flows on Cartesian grids. This feature will remove restrictions regarding the flow domains in which the moving boundary simulations can be performed, without having to resort to body-fitted grids or unstructured grid layouts. We first compare solutions for a stationary interface with those from well-tested solution procedures using body-fitted coordinates (Shyy, 1994). This is one of several exercises 0.9 08 07 0 0.1 0'2 0.3 0 4 05 06 0.7 08 0.9 Figure 3.1. Comparison of a deformed cavity at Re= 1000, for results using the cut-cell technique on the collocated grid with the solution on a boundary-fitted staggered grid. (a) Streamline contour for the collocated grid. (b) u-velocity along vertical centerline for the present collocated grid case (solid line) and the staggered grid case (open circles). (c) v-velocity comparison along vertical centerline. 77) performed to validate the conservation and consistency characteristics of the cut-cell method used to treat the complex geometry. The interface is sufficiently deformed that all the possible types of cut cells are encountered. The square cavity is a frequently adopted test bed for numerical experiments for incompressible flows due to available benchmarked solutions for comparisons. We deform the base of the cavity, the amplitude being 10% of the base. A 121 x 121 Cartesian grid is employed. We first present the results for a driven cavity flow, where the top wall of the cavity is pulled at velocity U=I corresponding to a Reynolds number of 1000. The results from the present method are compared with a previously benchmarked body-fitted staggered grid formulation (Thakur et al., 1996) with the same grid size. The results are shown in Figure 3.1. The streamline pattern is shown in Figure 3.1(a). A quantitative comparison can be seen in Figures 3.1(b) and (c) where the centerline velocities are plotted. As seen in this figure, the centerline profiles obtained from the cut-cell technique using the collocated grid (solid line) agree well with the body-fitted staggered grid calculations (open circles). We next demonstrate the ability of the complex geometry treatment to simulate more complex flow situations. In Figure 3.2, the streamlines and pressure contours for a driven cavity with top and bottom moving walls, and enclosing a semicircular cavity at the bottom wall are shown. The smaller cavity is indicated with dotted lines in Figure 3.2, which represents a solid boundary with zero internal thickness. This configuration gives an indication of the versatility of the present algorithm. Figures 3.2(a) and (b) show the streamlines and the pressure contours in the two cavities, respectively. As can be seen from the figure, for the Reynolds number of 1000, the multiple recirculation zones, including the one in the smaller cavity are well captured. It would be a significant task to fit a good boundary-fitted grid to this configuration. An added advantage that may be realized is that the grid can be refined easily in the present Cartesian framework without consideration of the solid boundary shape or the number of solid boundaries present. (a) U= I 1.0 0.9 0.8 0.7 0.6 y 0.s5 0.4 0.3 0.2 0.1 0.0 0.0 0.2 0.4 0.6 0.8 1.0 U= 1 x (b) 1.0 0.9 0.8 0.7 0.6 Y 0.5y o 0.4 0.3 0.2 0.1 0.0 0.0 0.2 0.4 0.6 0.8 1.0 x Figure 3.2. Double-Driven cavity flow with obstruction. Re= 1000, 121 x 121 grid. Dotted line shows the shape of the enclosed cavity. (a) Contour plot for streamline. (b) For pressure. Figure 3.3 shows the results of a natural convection flow in a square cavity with multiple objects in the flow field. The results are for a flow with a Prandtl number of 1 and Rayleigh number of 106. The sidewalls of the cavity are heated, the non-dimensionalized temperature of the left wall maintained at 0 and that of the right wall at 1. The top and bottom walls are adiabatic. Each of the cylinders is maintained at a constant temperature of 0. In Figure 3.3(a) we show the isotherms and in Figure 3.3(b) the velocity vectors. The skew symmetry of the temperature field, obtained in such a cavity without obstacles, is broken by the placement of the five shapes. In particular, the temperature gradients in the right bottom corner appear to be enhanced by the presence of the cylinder in that location while those in the lower left corner are diminished. This would obviously have a significant impact on the heat transfer within the cavity shown. This type of a facility would be useful in computing heat transfer around complex geometries, and can be combined with phase change as demonstrated in Udaykumar et al. (1996). No restriction exists on the number or shape of the obstructions that can be inserted in the domain. Since all the measures dealing with the interior solid boundaries are handled by means of one-dimensional arrays, the expense of incorporating internal boundaries varies linearly with the number of objects. No grid generation measures are required when a new object is introduced, thus objects can be added or deleted without affecting the existing geometry globally. Finally, a transient calculation is performed in an open channel configuration. The case shown in Figure 3.4 demonstrates the usefulness of this procedure for the time-dependent computation. Here, the inlet is at the left of the domain. The Reynolds number in this case is 100. A 121 x81 grid is used. The development of the recirculation zone behind the cylinder is shown until a steady-state is reached at a non-dimensional time t=3.5. In Figure 3.4(e), we show the streamlines obtained for this geometry from a steadystate calculation on a well-tested boundary-fitted staggered grid (Thakur et al., 1996). The flow fields obtained from the two methods are identical. In particular, the recirculation zone (a) __1.0 T=0 0.8 T=0 0.6 y 0.4 T= 1 0.2 0.0 . . I 0.0 0.2 0.4 0.6 0.8 1.0 Adiabatic wall x (b) 1.0 0 .8 \" I 'l J l I . . . . . . ..0.6 / 1/ y1 ---------1\\\ .... ----------0. \\\\\\' \ -' \\\l \\\N \"... ... .. '' /\ 02 - % \\\ " . - / 0.0 T 7 -_ 7 0.0 0.2 0.4 0.6 0.8 1.0 x Figure 3.3. Natural convection flow in a cavity. Pr= 1, Ra= 106, 121 x 121 grid. (a) Contour plot for temperature. (b) Velocity vectors. (a) t*=0.5 Inlet --- --- -_ ........ (b) (c).. (d). Cylinder Cylinder t *=1.5 t*=2.5 Exit t*=3.5 (e) Figure 3.4. Open channel flow, transient calculation. Re=100, 121x81 grid. (a) to (d) Streamline plot showing development of the steady recirculating flow behind the cylindrical obstruction computed using the cut-cell method. (e) Steady-state streamlines from the body-fitted calculation. in each case extends up to x=3.61 in the rear of the cylindrical obstruction. The slight differences in the streamlines between (d) and (e) are due to the differences in data representation on the Cartesian and curvilinear grids which is interpreted differently by the contouring fa- _- -. - ------ cility. It is also noted that the cylindrical obstruction has been shown in this figure by means of a phase fraction contour. The step-like nature of the obstacle is a limitation of the contouring procedure. As far as the solver itself is concerned, the surface is represented smoothly by means of the cut-cell technique. 3.2. Moving Boundary Flow Problems 3.2.1. Deformation and Recovery of Viscous Drops in Extensional Flows We now present results for the drop deformation and recovery studies. Consider the motion of an incompressible, neutrally buoyant drop in a straining flow described by the Stokes equations. The non-dimensionalization of the governing equations is performed using the reference quantities as described below. We define the viscosity ratio between drop and suspending fluids asA = Pitl/o. The characteristic velocity, pressure, viscosity and length are the following: u, = y/Uc, Pc = (Ic ic / fc), c = uo is the viscosity of the suspending fluid, Ec = a = the radius of undeformed circular shape and y is the surface tension. Then, the non-dimensionalized governing equations for creeping flows become Continuity equation: V. = 0 (3.1) Momentum equation for suspending fluid: 0= - Vp + V2 + xid( - (,))ds J (3.2) for the drop: 0 = - Vp + AV2 + x( - (e))dS J (3.3) In the above equations, V and p represent the velocity and pressure, respectively, hr is the unit normal vector, defined as positive when pointing outward from the drop phase to the suspending fluid phase, x = V - f represents the local mean curvature of the interface, as defined in Chapter 2.2.2; 6 denotes the delta function, Y and (x)e are the location of points in the flow field and on the interface. The integration term shown in Eqs. (3.2) and (3.3) represents the Eulerian contributions of the normal stress discontinuity due to the capillary force evaluated on the interface. This force is transmitted from the interface to the fluid by use of the Immersed Boundary Technique originated by Peskin (1977). The detailed description of the numerical implementation of the force term is given in Chapter 2.4.2. The far field condition for planar flow is 2 0 xi [0 __C 0 -2 0 7 Ca 0 2 0 =[0 (3.4) and for axisymmetric flow, 2 0 O' X V a 0 -1 0 2 0 0 -1] (3.5) where the capillary number Ca which represents the ratio of viscous force to surface tension is given by (Po Ga/ y). G is the strain rate of the imposed flow. The discretization of the above equations is described in Chapter 2.2.4. The flow condition and grid layout is indicated in Figure 3.5. The grid distribution is nonuniform; in the far field, i.e., away from the inner region shown, the grids are stretched linearly toward the boundaries of the domain, while in the inner region a uniform fine grid is used. This yields better resolution in the desired region of the flow field, i.e., in the vicinity of the interface. The flow field is imposed as shown, thereby exploiting the symmetry of the problem. The incoming flow at the top and outgoing flow at the right are held constant as specified by Eq. (3.4) and Eq. (3.5) for planar and axisymmetric cases, respectively. The distribution of markers along the interface is maintained at a constant arc length spacing. The arc length spacing used in all computations presented here is 0.8h, which was found by Inomng flow Regio of coarse grid Incoming flow Axis of symmetry Region of fine grid Outgoing flow Y V S~Axis of symmetry Sx Location of droplet Figure 3.5. Illustration of domain of computation and boundary conditions for drop deformation simulation. experimentation to provide adequate resolution and accuracy. The redistribution along the interfacial curve is performed at every time step to prevent undesirable depletion or accumulation of markers. The dimensions L and W are, respectively, defined as the deformed extents of the drop in the x- and y- directions. At the outset, a suitable grid is determined for the moving boundary simulations. A grid independence study was first performed to obtain the steady-state shapes of the drop in extensional flow at a capillary number of 0.1, which is slightly below the critical capillary number. In Figure 3.6, we compare the results obtained from calculations on an 81x81 and 161x161 grid with and without enforcing mass conservation condition. As can be seen, the results from the two grids with mass conservation imposed are in close agreement. We 1.5 1 0.5 -0.5 -1 -1.5 -1.5 -1 -0.5 0 0.5 1 1.5 Figure 3.6. Grid independence test for Ca=0. 1, X= 1. The final steady-state shapes are shown for (a) ..... : 8 1x8 1 grid without mass conservation (b) - -: 161xl61 grid without mass conservation (c) - : 81 x8 1 grid with mass conservation (d) -- : 161x161 grid with mass conservation. henceforth retain the same grid density as that corresponding to the 81 x8 1 grid in the present case. The number of grid points will therefore scale with the domain size. The moving boundary capability using the immersed boundary technique is first tested against the results of Stone and Leal (1989a, b) for steady axisymmetric drop shapes under an uniaxial extensional flow. For a range of capillary numbers below the critical capillary number Car the drop deformation culminates in a steady-state shape. Stone and Leal (1 989a) have proposed that the deformation parameter D=(L-W)/(L+W) is an indicator of the final bubble shape at steady-state, where L is the major radius and W the minor radius of the ellipsoidal drop as shown in Figure 3.5. This quantity therefore provides an estimate of the aspect ratio of the steady-state shape. In Figure 3.7, we compare the D vs Ca curves obtained from our calculations with those of Stone and Leal (1989a,b) for viscosity ratio A=1. The grid in this case was 61x61 in the inner region of dimensions 3x3 units. The arc length spacing of the interfacial points, following testing to determine accuracy, was fixed to be 0.8h in this and following calculations. The number of points on the interface was 45 at steady state. As can be seen from the figure, our results are in good agreement with those of Stone and Leal (1989a). Furthermore, the critical capillary number of Cacr =0.12 obtained in the present case also corresponds to the value in Stone and Leal (1989a). Figures 3.7(b) and (c) show the flow field and pressure contours in the case of Ca=O. 1. It is clear from the velocity vector plot that the flow at steady-state is tangent to the surface of the drop and the single-domain treatment of the interfacial forces yields the desired shape and flow field by the balance of normal stresses in Eq. (2.12). The large pressure gradient across the interface is also seen clearly in this figure. We next study the unsteady deformation and recovery characteristics of drops under imposed straining flow with Ca > Cacr. The inner region of fine grid in this case is 6x6 units. The domain of computation is 8x8 units in size. In the inner region 121 x 121 grids are used, while in the coarse grid region 20 points are used along each direction and the grid is linearly stretched. In Figure 3.8(a), we show the shapes of the drop at different times for the Capillary number of 0.14. In this case, the drop deforms to a maximum extent of 5.4 times the initial radius at non-dimensional time t=35. The number of points on the interface increases from 27 initially to 140 at the final stage. To study the recovery of the drop, we turn off the flow at the stages corresponding to the shapes at t=28 and t=35 and track the recovery path in the two cases. The results are presented in Figures 3.8(b) and (c). As can be seen, for the smaller deformation case the interface recovers back to the spherical drop shape, while for the larger deformation case, the recovery process results in the impending breakup stage shown. Although the further development of the drop after this stage can be handled by the code, we stop the computations (a) 0.41 03h D 0.2 1 35 0 0.2 0.4 0.06 0.08 0.1 0.12 0.14 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0n0 05 10 nn05 1.0 C 2746 0 B 2535 A 2-113 - 9 1.888 08 a8 1.680 7 1269 6 0.889 5 0.636 06 4 0213 3 0.045 2 -.0.026 0.4 1 -0.209 D.4. 02 0 0 05 1.0 15 Figure 3.7. (a) Comparison of D=(L-W)/(L+W) vs. Ca obtained from the present method (solid line) with the results (dashed line) of Stone and Leal (1989a). The critical capillary number is approximately 0.12 in both studies for X= 1. (b) Velocity vectors at steady-state for Ca=0.1. (c) Pressure contours at steady-state for Ca=0.1. . . . . . . . .N . . . - - o- -N - - - - -. . . -. - -. - -.--. - - -- - - --. . . . . . . . . . . .. N. . . . . . N N 'oN " . N'- N-,'. -.- - - - - -.\ - - - -"- (a) t = 0, 7, 21, 28, 35 (b) t = O, 5, 10, 25 (c) t = 0, 10, 27, 31 Figure 3.8. Deformation and relaxation of the drop for capillary number of 0.14 (critical capillary number is 0.12). (a) Deformation stage, shapes shown at t=0, 7, 21,28 and 35. At t=35, the deformation is 5.4 times the original undeformed radius. (b) Recovery of the drop from shape shown at t=28 in (a). (c) Recovery of the drop from shape shown at t=35 in (a). at this stage. From Figure 3.8(b), it is seen that the drop returns to a sphere of radius 1, which demonstrates that the mass conservation procedure implemented holds well in the course of the large deformations experienced by the drop. It may be surmised from Figure 3.8(c) that a small satellite drop is formed along with two larger drops in the recovery from the large deformation. This behavior of the drops with extent of deformation for the same capillary number is in qualitative agreement with the results of Stone and Leal (1989b). In Figure 3.9, we show the flow fields resulting from the computations. In Figure 3.9(a), the velocity vector is at the top and the pressure contour is at the bottom for a drop deformed at t=35 and Ca=0.14. As can be seen, the flow in this case tends to draw the drop out into a filament. Figure 3.9(b) shows the velocity vectors for the recovered shape after release from a small deformation. The velocity magnitudes are very low (O(10-4)) in this nearly quiescent state of the flow field. The corresponding pressure contours are shown underneath. It can be seen that the pressure levels inside the drop are uniform and equal to approximately 1.96 in value, which is close to the expected static value for the given initial drop shape. In Figure 3.9(c), the velocity vectors and corresponding pressure contours inside a drop corresponding to the recovery case shown at t=31 in Figure 3.8(c) are shown. The impending breakup of the drop is evident from the large velocity and pressure values observed in this figure. The capillary induced flow which results in the breakup of the drop is clearly seen in regions of large negative curvature in Figure 3.9(c). Due to the lack of resolution in the region of breakup, the values in this figure are not of significance. It is noteworthy that the collocated grid solution framework holds up in the presence of the large pressure gradient seen in this case (maximum pressure value=24, minimum=-0.8). Figure 3.10 presents the recovery characteristics of viscous drops in uniaxial extensional flow for the viscosity ratios A =1 and 10. In Figure 3.10(a), we plot the deformation of a drop for the case of Ca=0.13, which is slightly above the critical capillary number, and A=1, i.e., the viscosities of the drop and the suspending fluid are the same. The drop is deformed to a certain L and a step change in the flow to Ca=0.065 is made. The step change in the imposed flow strength is made at three different points along the deformation curve. The shapes along the deformation and recovery curves are also shown in the inset boxes. These curves seek to demonstrate the accuracy of the present calculations for long periods and large magnitudes of drop deformation. The results here agree very well with those of Stone and Leal (1989b). In particular, the fate of the recovering drops from three different Figure 3.9. Velocity vectors and corresponding pressure contours plots. (a) deformed drop at the time t=35 in Figure 3.8(a); (b) drop at a relaxation time t=25 in Figure 3.8(b); (c) drop at a relaxation time t=31 in Figure 3.8(c). on- -o -----------s------- :i i ii , , ! ! i- --- - -- --- --- -- -- -- -- - - i i - - - - - - - - - - - - . . iiiii------- - ... . . initial deformations, following a step change in the flow strength match with the predictions in Stone and Leal (1989b). The initial large deformation of the drop due to stretching is evident in the inset marked I. When the released drop takes path II, i.e., for a L value before release of 2.79, the drop recovers to a steady-state shape corresponding to Ca=O.065. For a value of L=3.07, no recovery is possible. In this case, as can be seen in the recovery curve, the drop initially shrinks. However, during this process a waist develops at the equator of the drop and results in the extension of the drop as the flow induced by the capillary forces accompanying this negative curvature drives the fluid outward, as seen in the extreme case in Figure 3.9(c). This situation is unstable and leads to the further development of the waist and subsequent breakup of the drop. Thus, drop breakup in such flows is dependent on the competition between the imposed flow fields and the flow induced by capillarity. Therefore, although a waist exists in the deformation curve marked I, in this situation the flows due to stretching and capillarity reinforce each other, while in the recovery phase they oppose each other. The result of the recovery process depends on the strength of the destabilizing flow induced by the capillary forces acting at the waist, which depends on the curvature at the waist. In Figure 3.10(b) for example, it is seen that the drop in case IV, which represents L=3.24 before release, rapidly heads towards breakup. This is in contrast to the other cases, where there is an initial tendency to respond to the imposed flow field before capillarity takes over. It is therefore clear that questions regarding recovery, number of drops resulting from breakup, sizes of drops and so on can only be answered in the context of the flow fields imposed, relative properties of drops and suspending fluids and extent of deformation. The rates of deformation obtained from our calculations agree very well with those in Stone and Leal (1989b). Some differences in times of deformation exist between the two calculations in the initial phase of the deformation. These may be due to the differences in initial conditions or capillary numbers chosen for the deformation. A notable point is that the time scale of deformation of drops in viscous extensional flows is expected to be propor- (a) 1o' IVV 1 ZIi ----D C) D ID DIIII i0 100 0 50 100 150 t (b) 1 o' IV I IV Dcz-D Iv CD C@D i D I Ill III 100 I~: I II 0 100 200 300 400 500 600 700 t Figure 3.10. Recovery curves for viscous drop subjected to uniaxial extensional flow as a test of the accuracy of the method for large deformations and property jumps. (a) The recovery curve for X=l. (b) Recovery curve for X=10. tional to (I+A) (Stone and Leal 1989a). A comparison between Figures 3.10(a) and (b) reveals the validity of this scaling. 3.2.2. Deformation of Drops in Extensional Flows with Inertia Effects In order to demonstrate the ability of the algorithm to deal with convective drop dynamics, we investigate the effects of an imposed uniaxial extensional flow. The scales chosen for this problem are the same as the ones for the viscous drop presented in the previous section. However, in this case, with the inclusion of the convection terms, the non-dimensionalized governing equations are Continuity equation: V = 0 (3.6) Momentum equation for suspending fluid: + V(VV) = - Vp + eV2V + e x@(- - (e))ds Ot R, W (3.7) for the drop: av +V(V)= -Vp+-V2V+ 1-, X k( (e))ds at Re W( ())ds (3.8) Under convective conditions the relevant non-dimensional parameters that arise are the Reynolds number (Re=2Qc(Gec)ec/lc, based on the drop radius) and Weber number (We=29c(Gfc)2fcy). The former provides a measure of the relative dominance of inertia to viscous effects, while the latter is the ratio of inertia to capillary forces. In the results presented the imposed flow corresponds to a Reynolds number of 10. We attempt to obtain the critical Weber number for a drop under such an imposed linear flow condition. Kang and Leal (1987) have obtained the critical Weber number of a bubble at this Reynolds number as lying between 0.9 and 1.0. In that work, the density and viscosity inside the bubble was considered to be negligible so that the pressure inside the bubble remains a constant, and the interior of the bubble is passive, in that it has no effect on the deformation characteristics. Under this assumption they fit a boundary-conforming grid to the exterior of the bubble and obtained the transient deformation characteristics of the bubble surface. This treatment is a limitation imposed by the boundary-fitted grid approach. If calculations were to be performed inside the bubble also, grid-generation would become complicated. In our framework this limitation does not exist. On the other hand, it is not a simple matter to render the flow field inside the bubble entirely passive, since the present method is a single-domain method. This rendering of a part of the domain passive was accomplished in previous calculations involving stationary and moving solid boundaries, where there is no flow in the solid region. In the case of two-fluid interfaces this is not as straightforward to do due to the fact that the interface velocity is governed by the kinematic condition at the boundary, and the boundary velocity is obtained from a single-domain formulation, Eq. (2.44). Thus information regarding the velocity field is required from both sides of the interface, which means that the flow field needs to be computed in both phases. Therefore, in performing our calculations in this section we will only be able to obtain a qualitative comparison with Kang and Leal (1987), and we will support our result with expectations based on the physics of the problem. The flow configuration remains the same as shown in Figure 3.5. The density and viscosity ratios used are fP = Qi=/o = 0.02 and A = pilPo = 0.02, respectively. In Kang and Leal (1987), it was shown that the deformation of the drop, in particular the availability of a steady-state solution for a given Weber number was dependent on the initial state of the drop in the flow field. Although for certain Weber numbers close to but below critical, gradually increasing the flow strength from a lower Weber number yielded steady-state solutions, imposing a step change in the flow field could drive the drop toward unsteady, continuous stretching. Such behavior was noticed in our simulations also. Therefore, for the cases presented here, we obtained the flow at a given We, under a particular combination of density and viscosity, starting from the flow field corresponding to a value of We 0.1 lower. The curves in Figure 3.11(a) show the variation of D with We for the im- posed flow at Re=O10. To start this curves, we first obtained the flow field at We=0.5 by imposing a flow field instantaneously. This led to an oscillation of the drop due to the competition between inertia and capillarity effects. This type of behavior is also reported in simulations by Mashayek and Ashgriz (1996). For example, P=0.02 and A=0.02, for upon imposition of the flow field the drop at We=0.5 stretches rapidly from its initial circular shape to a value along the x-axis of L=1.4, and subsequently relaxes in an oscillatory manner to asteady-state value of L=1.17. Stepping up the Weber numbers in increments of 0.1 thereafter yields the curve shown in Figure 3.11(a). For3=0.02 andA=0.02, and at a Weber number of 1.2, a steady-state did not exist. If these fluid properties are changed, then the critical Weber number changes accordingly. A closer examination of the case of P=0.02 and 2=0.02 can be seen from the curve of L vs t shown in Figure 3.11(b). There the deformation histories for We=0. 7, 1.1 and 1.2 are plotted. As can be seen, these curves show the oscillatory relaxations mentioned above. However, in the case of We= 1.2, the oscillations carry the interface away from a steady-state and instability results. Note that for the case of a deforming bubble, Kang and Leal (1987) obtain 0.9< Wecr <1.0 while our results appear to indicate 1.1< Wecr <1.2. This discrepancy in the critical value of Weber numbers is somewhat puzzling. However, as pointed out in Ryskin and Leal (1984), the critical Weber numbers for bubbles (3 -+ 0, A -+ 0) and drops are expected to be different. This is due to the non-negligible density and viscosity of the fluid in the drops. In our computations, we expand the parameter range of density and viscosity combinations. To gain more insight, the following description is given. Considering the normal stress balance at the interface represented by Eq. (2.12), the pressure on each side of the interface needs to be obtained. This is achieved if one defines a system of coordinates at the interface, such that r and & are coordinates oriented along the tangential and normal directions to the interface. Now, the pressure in each phase at the interface may be obtained by solving the momentum equation at the interface, which may be written in the form (Ryskin and Leal 1984) 0.3 F 0.1 F 0.4 0.5 0.6 0.7 0.8 We 0.9 1 Figure 3.11. Deformation of a bubble under uniaxial extensional flow at Re= 10. (a) Curve shows the deformation parameter D vs Weber number We. The critical Weber number for this case lies between 1.1 and 1.2. (b) L vs t plot for three different Weber numbers (We=0.7, 1,1, 1.2). -o- P3=0.02, k=0.02 -*-. 3=0.02, X=0.002 -+-:. 3=.1, X=0.01 .-x-.: "=1, k=0.01 x / , - -~. - V(p+ ) = ( xo )-2V x (3.9) where wto is the vorticity. Integrating this expression along the ril-direction and considering the idealized case of inviscid, irrotational flow, one obtains: Pi/o = (Pc)i/o - 1/o(V,2)i/o (3.10) where the subscript i/o implies that the expression applies to both the inner and outer phases of the interface and V, is the r (tangential)-component of the velocity at the interface. The component of pressure pc is hydrostatic. The difference between the inner and outer hydrostatic pressure components is adjusted as described previously for enforcing mass conservation. Under this idealization, the dynamic pressure plays the dominant role in the normal stress balance at the interface. This is an important term that distinguishes the behavior of bubbles from the inviscid drops. In the case of the bubble, the dynamic contribution to the pressure from the inner phase is negligible and the stability of the bubble is influenced only by the outer dynamic contribution. For the case of an inviscid drop with the densities /=1, the kinematic condition implies (V4)i = (V)o. Therefore, under this condition, when the pressures at the interface as given by Eq. (3.10) are inserted into Eq. (2.12), the dynamic contributions to the normal stress balance exactly cancel. In inertia dominated flows it is the dynamic contribution that dominates in the stress balance. Therefore, in this situation the behavior of the inviscid drop will be vastly different from that of the bubble where the inertia effect is negligible, and different densities of the drop will cause different critical Weber numbers. For the case of a drop, since the dynamic contributions from the inner surface partially counteracts the dynamic pressure contribution on the outer phase, the critical Weber number may be expected to be higher than that for a bubble. In the simulations presented here, the ratios for density and viscosity are f=O.02 and A=O.02. To see what the flow field characteristics are in this case, we plot in Figure 3.12(a) the velocity vectors for a We=1.2, i.e., the critical value. The pressure contours and profiles are shown in Figures 3.12(b) and (c), re- 63 (a) (b) j . . . .. .. ..- . 98 (c) (d) 22 11 P u. 00 p u -1 1 -2 -1 0 20 40 60 80 0 20 40 60 Figure 3.12. flow field characteristics for the case of drop in extensional flow at Re=10, We=1.2, 3=0.02, X=0.02. (a) Velocity vectors in the region containing the drop. (b) Pressure contours in the same region. (c) Pressure profiles at 8 different sections along the x-axis starting from i=4 to i=95 (drop extends from i= 1 to i=98). Arrow indicates increasing i. (d) u-velocity profiles at the same stations as in (c). spectively, at 8 equally-spaced sections, starting from i=4 to 95, along the x-axis in the region occupied by the drop that extends from i=1 to 98. The arrow in Figure 3.12(c) shows the direction of increasing i. As seen in Figures 3.12(b) and (c), there is ajump in the pressure at the interface location. It is seen that although the pressure at each x-station inside the drop is practically constant, there is a variation of the pressure in the drop in the tangential direction. In fact, the pressure rises somewhat as one goes from the top of the drop toward the side. Although this rise in pressure is not as high as that on the outer phase as can be clearly shown in Figure 3.12(c), this may lead to the enhanced stability in the present case when compared to a bubble where the pressure inside is a constant and is solely hydrostatic. Note that in the interior of the drop the pressure along the interface rises from the value of about 1.7 to 1.95 (6p=0.25), while outside the drop it decreases from -0.8 to about -2.3 (bp=-1.5). Therefore, the tangential increment of the pressure inside the drop is non-negligible compared to that outside. In Figure 3.12(d) we present the u-velocity profiles at the same sections as the p-profiles in Figure 3.12(c). The arrow in Figure 3.12(d) also indicates the direction of increasing i. As can be seen large velocity gradients exist in the interior of the drop due to the strongly recirculating flow present there. However, the viscous stresses in the drop, in relation to that outside, may not significantly influence the drop stability via Eq. (3.12) in view of the inertia-dominated flow and due to the viscosity inside the drop being 0.02 times the outer viscosity. It is not straightforward to quantify the individual contributions of the viscous stresses and pressure terms in the balance expression given by Eq. (2.12) due to the non-linearity of the flow field. At the present time, based on the flow quantities displayed in these figures, it is proposed that the discrepancy in the Wecr value is due to the fact that even for the fairly large density and viscosity jumps considered here, the dynamics internal to the drop is nonnegligible. Further investigation of these aspects is required and is in progress. Suffice it to say that, as pointed out by Ryskin and Leal (1984), the results from studies of the behavior of bubbles do not carry over directly to inviscid drops, except at negligible Reynolds num- bers; for the imposed flow at small Reynolds numbers, the pressure inside the bubble and inviscid drop are constant, since there is no dynamic contribution to pressure. 3.2.3. Deformation of Drops in Constricted Tubes In the above sections we have demonstrated the ability of the algorithm to 1) handle complex flows in the arbitrary domains represented on a Cartesian grid, 2) to accurately obtain the dynamics of interfaces under the balance of inertia, viscous and capillary forces. In this section we tie these aspects together in studying the deformation of drops in viscosity and inertia dominated flows through a constricted tube. The interaction of drops with the geometry in which they are constrained to traverse has received some attention from experimental (Olbricht and Kung, 1992; Olbricht and Leal 1993) and numerical viewpoints (Tsai and Miksis, 1994; Manga, 1996a, b). Since the major focus has been on the behavior of multiphase flow through porous media with applications in oil recovery, such flows have been restricted to the creeping flow regime. The boundary integral technique is applied to simulate such a problem (Tsai and Miksis, 1994). In recent work, Manga (1996a, b) applied the same technique to study the behavior of drops in geometries such as a driven cavity and bifurcating channel, again restricted to Stokes flow. To our knowledge, numerical work on the behavior of drops in fixed geometries with significant inertia effects has been restricted to the study of drops in straight capillaries (Lee et al., 1995). While the boundary integral technique is unable to handle flows with convection, purely Eulerian methods have to contend with limitations imposed by the Cartesian grid on which they are based. The present algorithm does not suffer from these limitations and thus will be employed below to investigate the flow of drops in a constricted tube at non-vanishing Reynolds numbers. Figure 3.13 shows the results of the simulation of drop motion through a constricted tube at small Reynolds number. For this case, we have employed the same conditions as one of the cases in Tsai and Miksis (1994). The constriction in the tube is a sine-wave with ampli- 0 2 4 6 8 10 Figure 3.13. Flow of a drop through a constricted tube for the case of Re=0.0, Ca=0.1, 13=1, X=0.01, a=0.9. (a) The interface shapes shown at non-dimensional times of 0.0, 0.25, 0.5, 0.75, 1.0, 1.25, 1.5, 1.75, 2.0. (b) ,(c) and (d) Velocity vectors for the times 0.25, 1.25 and 1.5 respectively. - - - iiiii 7 7- . :-1 N N - - - tude 0.6R and wavelength 2.0 units, where R is the radius of the tube. For this case Ca=0.1, the density ratio P = Qi/Qo = 1, the viscosity ratio A = pilo = 0.01, the ratio of the undeformed drop radius to the tube radius, a=0.9. As in the case of Tsai and Miksis (1994), the circular drop is placed initially such that its center is at x=4.0. The initial shape of the drop is shown by the dotted lines in Figure 3.13(a). The evolution of the drop is shown at equal intervals (6t=0.25) of time. The interface shapes at these instants agree remarkably well with the results of Tsai and Miksis (1994) using the boundary integral technique. As can be seen, the drop deforms essentially in conformity with the shape of the constriction for this case. In Figures 3.13(b) to (d) we show the velocity vectors for this case at instants before, during and after passage of the drop through the constriction. The effect of capillarity can be seen in these figures from the turning of the velocity vectors in the region where the flow field is not strong, such as at the regions close to the wall. Also to be noticed is that at the constriction, the front of the bubble is elongated whereas the rear is flattened by the acceleration of the flow at that location. As an additional element in the physics regarding to how the inclusion of inertia alters the dynamics of the drop, we simulate the flow of the same drop in the constricted tube, but under convective flow conditions. If the convection is not negligible, the boundary integral technique used by Tsai and Miksis (1994) cannot be applied due to the non-linear field equations. The Reynolds number based on the tube radius is 50. For the case shown in Figure 3.14, Ca=0. 1, = 1, A=0.01 and a=0.9. The response of the drop in this situation is significantly different. Due to the presence of the recirculation zone aft of the constriction the picture changes drastically. For the area ratio and Reynolds number considered, there appears to be an unsteady motion in the region following the constriction. A region of recirculating flow is formed and travels downstream of the throat. This region of the flow has a profound effect on the interface shape. The interface gets stretched into an elongated shape downstream of the constriction due to the higher flow speeds there on account of the rapid area change. Since the material boundary is convected with the flow the elongated drop shape 68 (a) o -ï¿½ \ I I I 0 2 4 6 8 10 12 (b) " = = = - " . ...:" . . . . . . . ...= = = (c) (d) (e) Figure 3.14. Entry of a drop into a constricted tube for Re=50, Ca=0. 1, 3=1, =OO0.01, a=0.9. (a) drop shapes at t=0, 0.25,0.625,0.875, 1.0 and 1.12. (b)-(f) Velocity vectors at t=0.25, 0.625, 0.875, 1.0 and 1.12, respectively. is formed. As the drop develops however, the shedding recirculating flow region has an effect on the interface and draws the front end of the interface into the shape seen in Figure 3.14(a). In this case, as the interface deforms, the rear of the drop is pushed toward the constriction and flattened. Ultimately, the drop collides with the wall, at which point we terminate the calculations due to uncertainty in the subsequent behavior. One possibility is that the drop will eventually pass through the constriction, as the film of liquid squeezed between the drop and the wall of the tube may act as a lubricating layer, facilitating the motion of the interface relative to the wall. This behavior cannot be captured due to the resolution used by the grid. Alternatively, the collision with the wall may result in rupture of the drop surface leading to breakup. Since the impending events are somewhat nebulous, we present results only up to the stage where the collision is detected by the algorithm. As a contrasting case, we perform the flow simulation with convective effects on a smaller drop, for which a=0.75. All other conditions are identical to the previous case. From Figure 3.15, it is seen that the dynamics is very similar to the case of the larger drop, except that eventually the drop passes through the constriction. Since the calculation was carried further than in the previous case, it is seen that the translating recirculating flow structure downstream of the constriction continues to deform the drop resulting in the flattened structure of the fore end of the drop. For the given Ca=0. 1, the inertia effects appear to dominate over the capillary forces as there is hardly any difference in the behavior of the smaller drop in comparison to the larger drop. Due to the larger surface curvatures of the smaller drop, in the case of sufficiently small Capillary numbers, it is expected that the drop size will have an effect on its interaction with the flow field. Also, due to the unsteady nature of the flow field, the precise shape of the drop downstream of the throat will depend on the relative phase of shedding in that region and drop entry aft of the throat. Also, from preliminary studies for a lower Reynolds number, where shedding is absent, the behavior of the drop is noticed 70 (a) I ï¿½f I I 0 2 4 6 8 10 12 (b) -~~ ~ ~ ~ -~j~~ - ----(d) (e) (f) (g) Figure 3.15. Entry of a drop into a constricted tube for Re=50, Ca=0., =13=, 1 =0.01, a=0.75. (a) Shapes of the drop at t=O.0, 0.25, 0.5, 0.625, 0.875, 1.25, 1.375, 1.5. (b)-(g) Velocity vectors at t=0.25, 0.5, 0.625, 0.875, 1.25, 1.5, respectively. to be significantly different from that presented here. These studies are currently being investigated but are not part of the scope of the present thesis. 3.3. Summary In this chapter, we have focused on the numerical aspects of the algorithm and have tested several well-defined and challenging flow problems. The performance of the method has been demonstrated. The applications include: (1). Flows in complex 2-dimensional geometries under significant convection, which obviates the need for grid-generation. In this section, we have compared results where possible to those given by a boundary-fitted code and shown good comparison. (2). The results of the deformation and recovery of viscous drops, subjected to an uniaxial extensional flow, show good agreements with those of Stone and Leal (1989a, b). Thus, it is established that the viscous-capillary balance is well captured by the numerical technique. (3). The deformation of drops under uniaxial extension flow under convective conditions has been studied. In contrast with the work by Kang and Leal (1987), here we consider the motion of the fluid inside the drop. It is interesting to note that the estimated critical Weber number for the present simulation is somewhat higher than that predicted by Kang and Leal (1987) for a bubble with no internal dynamics. As discussed in Reskin and Leal (1984), this discrepancy in predicting critical Weber numbers is expected. In line with Mashayek and Ashgriz (1996), the relaxation oscillations of the drop under the competition between inertia and capillary forces is observed. (4). Finally, the various aspects of the algorithm tested above are grouped together in order to study the behavior of drops in constricted tubes. For low Reynolds number flows, we show that the results of Tsai and Miksis (1994) are exactly reproduced. But we go beyond this creeping flow regime, and find some interesting behaviors of the drop in the region downstream to the throat at Re=50. CHAPTER 4 DYNAMICS OF COMPOUND DROPS The focus of this chapter is to analyze the deformation and, in particular, subsequent recovery of 3-layer or compound drops in order to shine some light on our current understanding of the leukocyte rheology. Recovery experiments have been extensively used in recent years in order to characterize the rheological properties of leukocytes (white blood cells). In these experiments, a cell of about 8 [tm in diameter is aspirated into and expelled from a 3 -7 [im diameter glass pipet. The advantages of the recovery experiments over those involving the flow of a cell into a pipet are that they are easier to perform and analyze, and that the potential adhesion of the cell to the pipet is not an issue. Using micropipet techniques, mechanical properties of leukocytes have been evaluated with Newtonian and non-Newtonian drop models. However, reported experiments indicate that these existing models fail to explain the complex rheological behavior of leukocytes because they do not consider the internal structure of the cells. It has been suggested that a more adequate model for describing the cell's behavior is the compound liquid drop model (Hochmuth et al., 1993). 4.1. Problem Statement and Mathematical Formulation We consider the dynamics of a compound drop immersed in a uniaxial extensional flow and its subsequent recovery when the imposed flow is removed. The physical situation corresponds to three incompressible Newtonian fluid layers of density Qi and viscosity ti, which occupy the regions Qi (i=l, 2, 3). The regions QD , Q2 and Q3 represent, respectively, the suspending fluid, shell and core as shown in Figure 4.1(a). The surface tension coefficients 012 and 023 at the two interfaces F12 and F23 are assumed to be constants. The dimensions R2 and R3 are the undeformed radii of the outer and inner drops, respectively. (a) n 2 r2 Q3 R3 223 (b) (' Axis of symmetry Incoming flow coarse grid Axis of symmetry fine grid Outgoing flow Y \1.- - L3 x , L2 - Axis of symmetry Figure 4.1. Illustration of the parameters in the computation of the compound drop deformation and recovery. (a) The definition of model parameters. (b) Computational set-up and boundary conditions. The choices of characteristic velocity, length and pressure are as follows: CY- R2P (J I R2, p It where = tI is the viscosity of the suspending fluid, C ~ 012 is the surface tension of the cytoplasm-suspending fluid interface. Then, the dimensionless governing equations for axi- symmetric creeping flows are: (1) Continuity equations in region Qi (i=1, 2, 3): V - Vi = 0 (4.1) (2) Momentum equation for suspending fluid (region 01): 0 = - Vp + V2V + f12 (4.2) and for each of the enclosed fluids, namely the shell (Q2): 0 = -Vp + X3V2V + f12 + f23 (4.3) and the core (Q3): 0= -Vp + X3V2V + f23 (4.4) In the above equations, k2 = 02/t1 is the viscosity ratio of shell fluid to suspending fluid, X3 = 03/0l is the ratio of core viscosity to suspending fluid viscosity, Vi and Pi represent the velocity and pressure, respectively, in region Qi. The forces f12 and f23 denote the Eulerian contributions of the normal stress discontinuity due to the capillary forces evaluated on the interfaces F12 and F23. These forces are transmitted from the interfaces to the fluid by use of the Immersed Boundary Technique (Peskin, 1977; Unverdi and Tryggvason, 1992). The form of the force contribution is given by: f12 = fY22n28(x - (xe))dF12 f23= f y3x3n3(x - (xe))dF23 J J (4.5) where ni is the unit normal vector, defined as positive when pointing outward from each region i, as shown in Figure 4.1(a), xi = V * ni represents the local mean curvature of the interface; x and xe are the location of points in the flow field and on the interface, 6 denotes the delta function, Y2 = 012/0 and Y3 = 023/0 are the dimensionless surface tensions for the shell and the core, respectively. The detailed description of the numerical implementation of these force terms is given in Chapter 2.4.2. For the imposed axisymmetric uniaxial extensional flow, the far field condition is V =ER=G 0-1 0 2 0 0 - (4.6) where E is the strain-rate tensor, R is the position vector and G represents the shear rate of the imposed flow. Due to the choice of the velocity scale, the nondimensional boundary condition at infinity becomes V . = a [ 02 0 i0 a V =C 0 -1 0 S 0 - (4.7) where the capillary number Ca provides a measure of the relative dominance of viscous forces to surface tensions and is given by (tIGR2/Ol2). The boundary conditions to be satisfied on each interface for cases involving no mass exchange across it are the continuity condition and the stress balances: (V)rl2 = (V)l = (V)2 (4.8) (V),23 = (V)2 = (V)3 (4.9) n2 T1 - X2n2 T2 = y2n2(V, n2) (4.10) n3 T2 - X33 = Y3n3(V2 1 n3) (4.11) where T is the stress tensor, defined as follows. T = - pl + [ (VV) + (VV)T] (4.12) It is noted that Eqs. (4.10-4.11) are implicit in the immerse boundary representation of the capillary force terms, f12 and f23, in Eqs. (4.2-4.4). The advection of the interface is performed with the kinematic condition dX - ((V), n n dt ' J i (4.13) where X represents a marker point on the interface. Considerations regarding accuracy and conservation of mass in each region of the compound drop are addressed in Chapter 2.5. 4.1.1. Computational Setup The choice of an extensional flow to pursue our study is motivated by the need to employ a simple linear flow field which does not add complexity to the physics. In addition, the behavior of single phase drops (henceforth simple drops) have been extensively studied from the point of view of deformation as well as recovery under imposed extensional flows (Stone et al., 1986; Stone and Leal, 1989a, b), and thus provide results for comparison. Linear uniaxial extensional flow is distinguished from simple shear by the absence of vorticity. No net translation of the center of mass of the drop occurs and the symmetries dictate that no tank-treading type motion is possible. The effects of this straining flow on a compound drop are investigated by solving the Stokes equations. The flow conditions and computational setup are indicated in Figure 4.1(b). The symmetries of the flow are exploited to mitigate the computational effort and only the first quadrant is treated. The distribution of markers along the interface is maintained at a constant arclength spacing. The redistribution along the interfacial curve is performed at every time step to prevent undesirable depletion or accumulation of markers. The arclength spacing used in all computations presented here is 0.8h, which was found by experimentation to provide adequate resolution and accuracy. The number of markers on the interface varies as the interface stretches. The grid distribution is nonuniform; in the far field, i.e., away from the inner region shown, the grid spacing increases linearly toward the boundaries of the domain, while in the inner region a uniform fine grid is used. This yields better resolution in the desired region of the flow field, i.e., in the vicinity of the interface. The flow field is imposed as shown, thereby exploiting the symmetry of the problem. The incoming flow at the top and outgoing flow at the right are held constant as specified by Eq. (4.7) for the axisymmetric case. The dimensions L2 and L3 are, respectively, defined as the final deformed extent of the shell and core in the x-direction just before recovery begins. 4.1.2. Selection of Model Parameters Even for Newtonian fluids, the compound drop model is represented by five parameters (X2, X3, Y2, Y3 and R3) as opposed to two (k2 and Y2) for a simple drop. Apart from the five model constants, an additional parameter that enters in this study is the capillary number of the imposed flow. One factor that can be fixed is the value of R3. This is based on experimental observation that the nucleus occupies a volume of 21% of the neutrophil. This corresponds to a spherical nucleus of nondimensional radius, R3=0.585. Since the values of the other parameters and the nature of the nucleus and shell are not settled (Dong and Skalak, 1992; Hochmuth et al., 1993; Waugh and Tsai, 1994), we assign these parameters by drawing on experimental observations. It is generally agreed upon that the nucleus is highly viscous compared to the shell. Therefore, we will examine the effect of introducing a viscous core in a Newtonian liquid drop and identify the flow and deformation conditions that are meaningful to explore. 4.2. Deformation Analysis 4.2.1. Steady-State Configuration ( Ca < Cac ) We first examine the steady-state shapes of the simple drop against the results of Stone and Leal (1989a) with X2=1 and Y2=l. In the case of a simple drop, for a range of capillary numbers below the critical capillary number, Cacr, the droplet deformation culminates in a steady-state shape. The quantity D=I(L-W)/(L+W)I is a sensitive indicator of the final drop shape at steady-state, where L is the non-dimensional major radius and W the minor radius of the ellipsoidal drop as shown in Figure 4.1(b). This quantity therefore provides an estimate of the aspect ratio of the steady-state shape. 78 (a) 0.40 IV 0.35 IV 0.30 shell 0.25 -I/ D 0.20 0.15 -00 II 0.10 / core 0.05 0.00 .. . . 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 Ca (b) 1,2 1.0 ' " 0.0 . 5 1.0C1c,5 y 0,6 -. . o.............. .." -" "x Figure 4.2. The steady-state behavior of the simple and compound drops when subject to a uniaxial extensional flow. (a) D vs. Ca for the simple (k2=1, Y2=1) and compound drops (k2= 1, k3=1, Y2=1, Y3=1, R3=0.5): * : simple drop (Stone and Leal, 1989a); - - - : simple drop (present results); o : compound drop (Stone and Leal, 1990); -: compound drop (present results). Shapes corresponding to the capillary numbers I-IV are shown. (b) Velocity vectors at steady-state for Ca=0.1. Our results for a simple drop with nondimensional viscosity k2= 1 are in good agreement with those of Stone and Leal (1989a), including the prediction of the critical capillary number, as indicated in Figure 4.2(a). This latter is the value approached asymptotically by the D vs. Ca curves. Beyond the critical capillary number, no steady-state shape is possible and the drop stretches continuously. The critical capillary number of Cacr=0.12 obtained in the present case corresponds to the value obtained by Stone and Leal (1989a). Also shown in Figure 4.2(a) is the result for a compound drop with a core viscosity k3=1 and radius R2=0.5, and a shell viscosity k2=1. The upper solid curve corresponds to the whole compound drop, while the lower solid curve represents the core. It is noted that although the compound drop deforms markedly and becomes unstable at a lower capillary number than the simple drop, the deformation of the core is very limited. Also to be noted is that, in agreement with Stone and Leal (1990), and as deducible from the flow field shown in Figure 4.2(b), the core deforms into an oblate spheroid, while the shell becomes prolate. The D vs. Ca curves for the compound drop model shows behavior in excellent agreement with those of Stone and Leal (1990). The steady-state velocity vectors in the drop are shown in Figure 4.2(b) for a capillary number just below the critical value, Ca=0.1. The recirculating flows inside the compound drop can be seen to cause the oblate and prolate deformation of the steady-state inner and outer interfaces. 4.2.2. Unsteady Response ( Ca > Car ) In this section, we retain the condition that the viscosity of the shell fluid k2= 1, while the value of R2 changes from 0.5 to 0.585 for the reason described in Section 4.1.2, and explore how the presence of a core alters the behavior of the simple drop. In Figure 4.3(a), we show the curves of L vs. time for a simple drop, a compound drop with a core viscosity k3=1 and a compound drop with a highly viscous core k3=100. The deformation of the shell under the imposed extensional flow first proceeds rapidly for all three cases. Thereafter the deformation slows down leading to a linear growth region. After 0 5 10 15 20 25 30 35 40 t 0.12- dL2 dt 0.08 0.06 0.04 0.02 U 0 5 10 15 20 25 30 35 40 t Figure 4.3. Response of the simple and compound drops to the uniaxial extensional flow for Ca=0.14 > Cacr. (a) L vs. time. (b) Rate of deformation dL2/dt. ...... simple drop (X= 1) --- : compound drop (2=1, 3) '/ -.-. compound drop (X2=1, X3=100) / shell core / / . / / / I - -t I. , core : simple drop (X=1) : compound drop (X2=1, X3=1) - compound drop (X2=1, X3=100) shell V1 ./ ï¿½ / "/ / / If - I N - a certain extent of deformation there is a second region of accelerated stretching. When the rate of deformation of the shell (dL/dt) is plotted against time, as shown in Figure 4.3(b), it can clearly be seen that there is a minimum in the deformation rate. The slowdown in the intermediate stage of the drop deformation in each of these cases is due to the formation of large curvatures at the lateral ends of the drop. Thus, during this period the capillary forces opposing deformation are higher, resulting in a lower rate of overall deformation. As the shell continues to deform however, the equator of the shell gets progressively flatter and instability results when the curvature at the equator begins to achieve negative values. At this point, the capillary forces assist the flow field in deforming the drop, at least in the regions close to the equator. Thus, less energy needs to be expended by the flow in deforming the interface and the shell is stretched continuously. This aspect of the dynamics has been demonstrated for simple drops by Stone and Leal (1989a). The only distinction between simple and compound drops in the deformation process is the time scale of deformation. The inclusion of the core enhances the overall stiffness of the compound drop as can be seen from the dL/dt vs. time curves, even though its presence produces a continuous stretching mode of the compound drop at a smaller Cacr than for the simple drop. 4.3. Recovery Analysis We next consider the recovery characteristics of simple and compound drops for two different imposed flow strengths, Ca=0.14 and Ca=0.35. In each of these cases, we first deform the shell to a desired stretched state by imposing an extensional flow at a given capillary number and then switch off the flow to allow recovery. 4.3.1. Effect of the Capillary Number In Figure 4.4, we show the deformed and recovering shapes of the compound drop for the above two capillary numbers. The case of a less viscous shell (k2=1) and highly viscous core (X3= 100) is studied. The deformed shapes under the imposed flows are shown in Figure 4.4(al) for Ca=0.14 and Figure 4(bl) for Ca=0.35. We define that td is the time for deformation to the desired extent and tr is the time for recovery starting from the initial deformed shape. It is noted that the characteristic times for deformation are substantially different for the two cases. Furthermore, although the extent of deformation is the same (L2=3), the shapes are different. For the larger rate of stretching, Ca=0.35, the shell assumes a sharper shape laterally due to stronger imposed shear rate effect. The time scale of deformation of a drop is proportional to the drop viscosity for highly viscous drops, as can be expected from the physical consideration that the larger of the viscosities between the drop and suspending phase should influence the dynamics (Rallison, 1984). This scaling arises naturally by obtaining the integral equation form for the interface velocity as in the boundary integral methods (Rallison and Acrivos, 1978). A natural parameter that emerges as a convenient scaling factor for the time of deformation of a highly viscous simple drop is (1+X) where k is the drop viscosity non-dimensionalized with respect to the suspending phase viscosity. Thus, for the compound drop configuration shown in Figure 4.4, the response of the shell to deformation is much more rapid than that of the core (corresponding to k3/k2= 100). Therefore, at the higher stretch rates, the deformation of the drop is governed by the less viscous shell. It is important to remember that the recovery rate of the drop depends on the history of deformation for both simple and compound drop models. It is recalled that although the flow field itself is quasi-stationary, the time-history of deformation is stored in the drop in terms of its shape. This is obvious from Figures 4.4(al) and (bl). Once the flow is switched off in the recovery phase, the return of the shell to its spherical unstressed condition is controlled by the starting shape at release. In the case of the shell deformed at Ca=0.35, the produced surface curvatures are higher than those of the case at Ca=0.14, hence the capillary forces causing recovery are larger. This accounts for the faster recovery of the shell after high rate of stretching. (bl) td=5.5, tr=0 (b2) tr=l.5 C:*: (b3) tr=6.5 (b4) tr=12.5 (a5) tr=48.1 (b5) tr=28.5 Figure 4.4. Effect of the capillary number or imposed flow strength on the deformation (L2=3) and subsequent recovery of the compound drop (y2=1, y3=1, X2=1, l3=100). (a) Recovery shapes for the compound drop deformed at low strain rate (Ca=0.14). td is the time for deformation to the extent shown. tr is the recovery time starting from the shape in (al). (b) Recovery shapes for the drop deformed at higher strain rate (Ca=0.35). The sequences of figures in (a) and (b) correspond to identical initial extension lengths. (a2) t3.1 (a2) tr=3.1 (a3) tr= 11.1 (a4) tr=19.1 (al) td=35, tr-0 4.3.2. Effect of the Core Viscosity The behavior of a compound drop is now examined for two core viscosity values, X3= 1 and k3= 100. The shell and suspending fluid retain identical viscosities, X2=1= 1. We study the recovery characteristics of compound drops after they have been deformed to the same extent under identical imposed flow (Ca=0.14). Shown in Figure 4.5 are the velocity vectors for the two cases. In Figures 4.5(a 1) and (b 1) we show the velocity vectors just before the flow is switched off and the recovery process begins. It is clearly seen that the less viscous core is only slightly more deformed. In both cases, we have maintained the same surface tension value for the core. Since the radius of the core is 0.585 times the whole drop, the capillary forces on the core are stronger than on the outer shell surface. Thus, the core is harder to deform. Furthermore, during the stretching process, the core lies in a region of low flow strength. These combined factors make the value of the core viscosity plays a relatively minor role in the deformation stage. Consequently, the time taken for the compound drop to deform to the extent shown in Figures 4.5(al) and (bl), namely L2=3, is not very different in the two cases (td=31 for k3=1 and td=35 for k3=100, respectively). It is observed that during recovery, the time of recovery for the more viscous core (tr=48.1) is much longer than that of the less viscous case (tr=24). This is due to the fact that, although in the deformation process the core lies in a region of low flow strength, during recovery this is not the case, as can be seen in Figure 4.5(a2). The capillarity driven flow attempts to deform the core. For the less viscous core, the recovery flow is able to deform it into an oblate shape as shown in Figure 4.5(a3). Notice that although in Figure 4.5(a2) the flow vectors in the core are in the same direction as the flow in the shell, in Figure 4.5(a3) the outer flow sets up a circulation in the opposite direction in the core. The velocity magnitudes in the final stage, i.e., Figures 4.5(a5) and (b5), are small, and the shapes are very close to spherical. (bl) L2=3, td=35, tr=O ------ --------------- ------------------------------- (al) LO=3 td=31 t -0 (a2) t 4 (a3) tr=9 (a4) 13: ,3 , (a 5 ) t 2 4... .. .. . .. .. . (b2) tr=3.1 (b3) tr= 11.1 (b4) t,= 19.1 (b5) tr=48.1 00 05 10 15 20 25 30 Figure 4.5. Dependence of the recovery velocity field on the viscosity of the core. Velocity vectors for the X3=1 and X3=100 deformation are shown in (al 1) and (b 1). The flow fields in the recovery stages are shown in (a2)-(a5) for X3=1 and (b2)-(b5) for X3=100. The coressponding velocity at different recovery times along the x-direction is plotted against x in Figures (a6) and (b6) respectively. At tr=O0, L2=3 and L3=0.618. z - - -- - - - -- - - ---------------------------------------- ------------ 1-1- I'll-, ------------- ------ ------Z Z - -- ---- -------------------- ------------------------------------------------------------------------------------------------ ---- - ------- -------------------- - - - - ------------ - -------. . . . . . . . . . . I - - - - - - - - - - - - - - - - - - - - - - - - - - - -------------------------- ------ - - - - - - - - - - - - - - - - - - -- - - - - - - - - - - - - - - - - - - - - - - - - ::-:--------------- ------- --- ---------- The effect of the highly viscous core is clearly seen by comparing Figures 4.5(a2)-(a4) with Figures 4.5(b2)-(b4). Figures 4.5(a6) and (b6) show the velocity along the x-direction, i.e., the normal velocity of the recovering interfaces along the symmetry axis y=0. It can be observed that the velocities at the enclosed drop are much smaller than at the outer interface. This is particularly the case for the highly viscous inclusion, which acts like a blockage to the flow induced by the capillary forces at the outer shell interface. The U(x,y=0) vs. x curves also give a picture of the flow strength at the different stages of the recovery, indicating for example that the velocity magnitudes shown in Figures 4.5(a5) and (b5) are very small and the recovery is practically complete. The velocity field in the vicinity of the highly viscous core is very weak, causing it to remain undeformed in the recovery process. This enhances the overall viscosity of the compound drop and causes the large difference in the recovery rates mentioned above. 4.3.3. Comparison of Simple and Compound Drops We now compare the recovery results of simple and compound drops for the cases presented above by plotting the L values versus time. In Figure 4.6(a) we show the L vs. tr curves for the capillary number (deformation rate) of 0.14 and deformation extent L2=2. The recovery of a simple drop, a compound with k3=1 and a compound drop with k3= 100 are shown. The difference between the recovery processes of different drops is not noticeable in the initial stage of recovery. After this stage the effect of the core becomes more noticeable in the recovery curves. Thus, the compound does recover much like a simple drop in the initial stage; but as the recovery proceeds, the presence of the core enhances the "apparent viscosity" of the compound drop. Also, as the viscosity of the core increases, the departure from the Newtonian simple drop behavior becomes more pronounced. Figure 4.6(b) shows the recovery characteristics after the drop has been subjected to the larger deformation rate Ca--0.35. In this case, the behavior of a compound drop with an extremely viscous core (k3=1000) is also examined in order to establish a trend. The results 87 (a) 2 ........ : simple drop (X= 1) --: -compound drop (2=1, k3=1) . . : compound drop (2= 1, X3=100) 1.5 % .shell L - core 0.5 -- 0 0 5 10 15 20 25 30 35 40 tr (b) 2 ........ : simple drop (X=l) : compound drop (X2=1, X3=1) . :compound drop (12=1, X3=100) : compound drop (X2=1, X3=1000) 1.5 shell L core 0.5 0 0 5 10 15 20 25 30 tr Figure 4.6. Recovery curves for the simple and compound drops after deformation at different strengths (capillary numbers) of the imposed extensional flow. Shown here are L vs. time curves for (a) Ca=O. 14 and (b) Ca=0.35. L2=2 for both cases before recovery is started. for the two higher viscosities are very close and such highly viscous cores can be taken to approximate a solid inclusion within the compound drop. A noteworthy point in this plot however is the close agreement between the compound drop with the core viscosity of X3= 1 and the simple drop. As can be seen in Figure 4.6(b), the core of the compound drop with k3= 1 is highly deformed by the imposed extensional flow. Thus, in the recovery process the core does not behave like a passive blockage but is itself recovering, a situation no different from the Newtonian simple drop. The only noticeable deviation occurs in the final stages when, as seen in Figure 4.6(b), the core has been driven to an oblate shape and has to recover its spherical shape. In that stage, the dynamics of the shell-core combination deviates from a Newtonian simple drop behavior. Finally, for the two extents of drop deformation studied (L2=2 and L2=3) at Ca=0.14, the recovery curves are very similar. This is true for both simple and compound drops with core size R3 < 0.6. 4.4. Summary From the above study, the following decisions are made regarding the parameters for the 3-layer leukocyte model to be presented in the next Chapter: (1). The capillary number of the flow employed to deform the compound drop does not have a significant impact on the recovery when the core viscosity is sufficiently high. Therefore Ca=0.14 is chosen in what follows. (2). The core viscosity, when high enough, does not change the dynamics markedly, as seen from Figure 4.6(b), where the values of k3= 100 and 1000 show little difference in the recovery curves. (3). The extent of the drop deformation does not have a significant effect on the recovery if the compound drop is sufficiently deformed. As a consequence, an extension value of L2=2.1 will be selected for the leukocyte modeling study in order to match some reported experimental data. CHAPTER 5 LEUKOCYTE MODELING 5.1. Choice of Model Parameters In this Chapter, we examine the effects of a viscous cytoplasm enclosing a highly viscous nucleus. The material parameters for the nucleus, cytoplasm and suspending fluid as well as the surface tensions at the inner and outer interfaces are chosen based on values pertaining to leukocytes. In view of the nebulous picture of the properties of the nucleus and cytoplasm of leukocytes (only apparent values being known based on global properties and there being wide discrepancies in the reported values), we choose to select only representative values for the material properties of the 3-layer model. The "apparent viscosity" of the leukocyte as measured in the recovery experiments by Tran-Son-Tay et al. (1991, 1994b) is around 800 Poise. The suspending phase viscosity is usually in the vicinity of 0.1 Poise. Non-dimensionalizing viscosities with the suspending phase value gives kapp=8000 for the overall cell viscosity. Numerically, dealing with the deformation of a drop with such a large viscosity ratio is extremely taxing, since the deformation time scale is very large (O(Xapp)). In order to speed up the calculations we set the suspending phase viscosity to 1 Poise. Thus the nondimensional viscosity Xapp of the cell (i.e. whole cell) is 800. This is justified from the fact that for viscosity values high enough, as will be shown later, there is no significant change in the dynamics of liquid drops, except that the time scale of deformation is scaled in proportion to the viscosity. The viscosity of the cell is some averaged property, based on the viscosity of its contents. Therefore, with a 21% volume for the nucleus (corresponding to the leukocyte nucleus) and assuming that the nucleus is 100 times more viscous than the cytoplasm, we assign viscosity values X2=40 and X3-=4000. In order to reduce the parameter space, we assume for now that the nondimensional surface tensions, Y2=y3= 1. The flow is applied at a capillary number of 0.14 and shut off when the cell reaches a deformed length of L2=2.1. 5.2. Effect of a Viscous Cytoplasmic Layer In contrast to the cell investigated in Chapter 4.3.2, where the cytoplasm fluid viscosity was equal to that of the suspending fluid (i.e., k2= 1), the cytoplasmic fluid viscosity under study here is much larger than the suspending phase (i.e., X2=40). The recovery curves of the cytoplasm and slightly deformed nucleus are represented in Figure 5.1(a) by the solid lines on the L(t) vs. tr plot. The recovery process is considered completed when the cell length L < 1.05. The recovery curve for a simple Newtonian drop with a viscosity value equals to the cell volume averaged viscosity, namely 832, is also shown for comparison in Figure 5.1(a). The recovery process of a simple drop is much slower than that of a cell with a comparable averaged viscosity value. Thus, the use of the volume averaging procedure is highly inadequate for obtaining the effective or "apparent" viscosity of the cell. In fact, Figure 5.1(a) shows that the recovery behavior of a simple drop with a lower viscosity value, e.g., 240 or 100, does approach better that of the cell under study. Although theoretical estimates exist for the effective viscosity of a dilute suspension of compound drops (Stone and Leal, 1990), the only measure of the effective viscosity of a compound drop is qualitative. The only guidance offered (Johnson and Sadhal, 1985; Stone and Leal, 1990) is that the effective viscosity of the drop is intermediate between that of a solid sphere (which limit is approached as R3- 1), and that of a simple drop (composed only of the cytoplasm, this limit of course being realized as R3-0). Furthermore, the effective or "apparent viscosity" of a compound drop, as measured by the rate at which the drop recovers, depends on the deformation history, in addition to its material properties. 0 200 400 600 800 1000 12 00 1400 1600 tr 0 0.5 1 15 2 2.5 3 ts=tr/( 1 +( [c)2-layer) ....... simple d -- -: simple d -- : simple d -: cell (k2=scaling cell rop (X=100) rop (X=240) rop (X=832) 40, X3=4000, factor Xapp=285) 0 0.5 1 1.5 2 2.5 3 0 0.5 1 1.52 ts=tr/(1 +(c)2-layer) ts=t/( 1 +(tc)2-layer) Figure 5.1. The recovery behavior of the 3-layer leukocyte model compared to several averaged simple drop models. (a) L vs. unscaled recovery time tr. (b) L vs. scaled recovery time ts. (c) Scaled rate of deformation dL/dts vs. scaled recovery time ts. (d) L vs. scaled recovery time ts. All viscosities and scaling factors chosen for 3-layer model are indicated in the legends. ....... simple drop (=-100) - - -: simple drop (X=240) ' -.-.-: simple drop (k=832) --- : cell (k2=40, k3=4000) nucleuslf nucleus 2 ....... simple drop (k= 100) : simple drop (k=240) \ .-. simple drop (=832) - cell (k2=40, k3=4000, 1.5 scaling factor kapp=285) cell L nucleus 0.5 12 (c) 10 8 dL 6 dts 4 2 (d) 2 1.8 ....... : simple drop (X= 100) ---: simple drop (X=240) -.-.-: simple drop (k=832) __--: cell (k2=40, k3=4000, scaling factor kapp-=40) cell A different perspective is obtained by plotting the L(t) curves against a rescaled dimensionless time ts defined as, tr ts - tr 1 + kapp (5.1) where the kapp is either the apparent viscosity of the cell or the viscosity, k, of a simple drop. A similar representation of the recovery results by rescaling the recovery time is also employed by Tran-Son-Tay et al. (1991, 1994a). In that work, experimentally obtained recovery data of leukocytes are compared to the behavior of Newtonian simple drops by fitting the numerically generated curve to the experimental data. Here, for the 3-layer model, the scaling factor is chosen to be kapp=285 based on matching asymptotic behavior of the cell with that of a simple drop. This yields the picture in Figure 5.1(b). The simple drop curves collapse together if the dimensionless recovery time defined in Eq. (5.1) is used. The result is expected because of the similarity in the recovery paths for highly viscous drops. The only driving forces in the recovery process are the capillary forces, which for the same starting shapes, should only yield differences in recovery times. Further insight can be obtained by plotting the scaled rate of recovery dL/dts against ts as in Figure 5.1(c). It is clearly observed that, while forts < around 0.2 the recovery of the cell is extremely rapid, a much lower recovery rate takes over at later times. There is a distinct change in the recovery behavior of the cells. Thus, the rapid recoil ascribed to leukocytes in initial recovery stages can be provided by such a 3-layer model. This rapid phase is of course due to the fact that the nearly undeformed stiff nucleus does not participate in the recovery dynamics, as seen in Chapter 4.3.2. Therefore, the initial recovery stage is dominated by the cytoplasm region. The viscosity of this region is low and hence the time scale of response is fast, based on the scaling in Eq. (5.1). When the deformation of the cytoplasmic layer decreases, the interaction of the cytoplasmic surface with the nuclear surface becomes more important, leading to the dynamics of the cytoplasm-nucleus combination which makes this period one of slow recovery. Therefore, even within the Newtonian fluid consideration there is a possible wide disparity in observed cell apparent viscosity depending on the recovery characteristics. We will tackle this issue in more detail later. If the scaling factor for the cell is chosen to be kapp=40 (the viscosity of the cytoplasm) instead of 285, the earlier phase of the recovery will be matched as seen in Figure 5.1(d). As seen in this figure, the recovery behavior in the initial stage follows the slope of the simple drops, before deviating from those curves at about ts=0.2. Thus, at least in the initial stages, the apparent viscosity of the cell is close to that of the cytoplasm. 5.3. Effects of Nucleus Deformation and Time Scales As mentioned in the previous section, it is possible to capture an initial rapid recoil of a leukocyte with a 3-layer Newtonian model, supporting the suggestion of Hochmuth et al. (1993). This is seen to be possible due to the preferential flow of the less viscous cytoplasmic liquid in the recovery process. However, one important fact has to be considered at this juncture. Under certain conditions, such as long holding times or low deformation rates, the recovery of the leukocytes ensues like a Newtonian simple drop. Based on our findings above, it appears that the existence of a nucleus leads to the distinct fast and slow regimes in the recovery process. In order to examine the conditions under which a cell can yield a Newtonian simple drop-type recovery behavior and to clarify the essential elements of the required 3-layer model, we will study the recovery dynamics of two sets of ratios (y2/2) and (y3/3). 5.3.1. Incompatible Nucleus and Cytoplasm Time Scales We want to investigate the impact of the deformation of the nucleus on the cell recovery. One difficulty encountered in this regard is that, as mentioned before, the extensional flow is weak in the vicinity of the nucleus. Also, it is difficult to obtain sufficient deformation of the nucleus due to its inherent stiffness, on account of the higher viscosity and smaller radius amounting to larger capillary forces. We attempted to increase the strain rate on the cell by raising the Capillary number to 0.35. This failed to achieve the desired deformation of the nucleus, due to the enhanced rapidity of flow of the cytoplasmic fluid. The outer interface deforms rapidly and contacts the nucleus at which time the calculations are terminated. The nucleus itself is slow to deform and remains mildly extended. Therefore, our strategy to test the effect of the nucleus deformation is to start the recovery process from given initial shapes of the cytoplasm and the nucleus, and exploit the quasi-stationary nature of the flow field. In the present computations, the cell has values Y2=1, Y3=1, k2=40, and k3=4000. This corresponds to a highly viscous nucleus. The shape of the drop is taken to be that obtained after deformation to L2=2.1 under Ca=0.14, and the shape of the nucleus is imposed. The velocity field in the recovery period is shown in Figure 5.2(a) for the case where the nucleus is deformed to a very small extent (L3=0.618) which resulted from the application of the extensional flow as in the previous section. Clearly, in the recovery stage, the response of the cell is limited to the outer interface only, while the nucleus itself serves merely as an internal solid-like obstruction to the flow field. It is this characteristic of the flow field which leads to the rapid recoil response of the cell. The flow fields for the larger extent case, L3=1, are shown in Figure 5.2(b). The deformed nucleus is assumed to be a cylinder with a hemisphere cap at both ends. One quadrant of the cell is illustrated in Figure 5.2(bl 1). The dimensions of the cylinder and hemisphere are assigned, for a given extent of deformation, so that the volume of the nucleus is maintained at 21% of the cell. The initial shape of the cell, i.e., just before recovery, is the one given in Figure 5.2(al). It can be seen that due to the higher viscosity of the nucleus, the flow resulting from the recovery of the outer interface is decelerated as it approaches the nucleus, leading to a rapid initial recoil phase. The apparent viscosity values, .app=285 for L3=0.618 and .app=220 for L3= 1, are obtained by using the scaling in Eq. (5.1) when the asymptotic recovery trend of the cell matches with that of the simple drop (Figure 5.3). Figure 5.3 shows the recovery curves of the simple drop and the cell with different initial nucleus deformations. The off-Newtonian trends of the cells are clearly seen when the nucleus is only slightly deformed. In the later |

Full Text |

PAGE 1 COMPUTATIONAL STUDY OF LEUKOCYTE RHEOLOGY BASED ON A MULTILAYER FLUID MODEL By HENG-CHUAN KAN 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 1997 PAGE 2 ACKNOWLEDGMENTS I would like to express my sincere thanks toward my advisor. Dr. Roger Tran-Son-Tay, for his expert guidance, encouragement, and financial support during the course of my graduate studies. He has been a constant source of inspiration for me both professionally and philosophically. The same gratitude is also given to my cochair. Dr. Wei Shyy, for his understanding and help in so many ways. Very special thanks go to Dr. H. S. Udaykumar for his tremendous help and enthusiasm in this work. In spite of spending the countless and frustrated hours in front of a variety of computers with me during the early debugging process, his patience and belief kept this work going. I would like to thank the members of my supervisory committee. Dr. Renwei Mei, Dr. Richard B. Dickinson and Dr. Nicolae D. Cristescu, for their fruitful advices. I wish to thank my colleagues at the Laboratory of Cellular Mechanics and Biorheology for their direct and indirect contributions and entertaining discussions. In particular, I thank Emmanuelle, David, Neal, Jie, Philou, and Philippe for their friendship and assistance. Words are not enough to describe my sincere gratitude to my parents, Pao-Hsiang Wang and Che-Chang Kan, for their persistent support and love. Their emphasis on higher education for their children pave the opportunity for me to complete this work. I also wish to thank my brothers, Po-Wen Kan and Feng-Jung Kan, for their continuous encouragement throughout these years. Finally, I reserve my deepest appreciation to my wife, Min-Hui. I will be forever indebted to her endless sacrifices and love. To her, I dedicate this dissertation. ii PAGE 3 TABLE OF CONTENTS page ACKNOWLEDGMENTS ii ABSTRACT v L INTRODUCTION AND BACKGROUND 1 LL Motivation and Objectives 1 1.2. Importance of Compound Drops Dynamics 2 1.3. Significance of Leukocyte Rheology 2 1 .4. Observed Rheological Behavior of Leukocytes 4 1.5. Current Status of Leukocyte Modeling 6 1.6. Summary 9 2 THEORETICAL ANALYSIS 12 2.1. Background 12 2.2. ELAFINT: A Mixed Eulerian-Lagrangian Solution Algorithm 15 2.2.1. Characteristics of the Numerical Methodology 15 2.2.2. Interface Tracking Procedure 16 2.2.3. Formulation of the Pressure-Based Flow Solver 18 2.2.4. Discretization Procedure of the Flow Field Equations 20 2.3. Momentum Interpolation 23 2.4. Interactions between Interfacial Physics and Flow Field Variables 28 2.4. 1 . Cut-Cell Procedures for Solid-Liquid Boundaries 28 2.4.2. Immersed Boundary Technique for Fluid-Fluid Interfaces 31 2.5. Mass Conservation in Individual Phases 36 2.6. Summary 39 3 ASSESSMENT OF COMPUTATIONAL TECHNIQUE 40 3.1. Flows With Fixed Arbitrary-Shaped Internal Boundaries 40 3.2. Moving Boundary Flow Problems 47 3.2. 1 . Deformation and Recovery of Viscous Drops in Extensional Flows 47 iii PAGE 4 3.2.2. Deformation of Drops in Extensional Flows with Inertia Effects . 58 3.2.3. Deformation of Drops in Constricted Tubes 65 3.3. Summary 71 4 DYNAMICS OF COMPOUND DROPS 72 4.1. Problem Statement and Mathematical Formulation 72 4. 1 . 1 . Computational Setup 76 4. 1 .2. Selection of Model Parameters 77 4.2. Deformation Analysis 77 4.2.1. Steady-State Configuration ( Ca < Cacr ) 77 4.2.2. Unsteady Response ( Ca > Cacr ) 79 4.3. Recovery Analysis 81 4.3.1. Effect of the Capillary Number 81 4.3.2. Effect of the Core Viscosity 84 4.3.3. Comparison of Simple and Compound Drops 86 4.4. Summary 88 5 LEUKOCYTE MODELING 89 5.1. Choice of Model Parameters 89 5.2. Effect of a Viscous Cytoplasmic Layer 90 5.3. Effects of Nucleus Deformation and Time Scales 93 5.3.1. Incompatible Nucleus and Cytoplasm Time Scales 93 5.3.2. Compatible Nucleus and Cytoplasm Time Scales 96 5.4. Discussion And Implications to Leukocyte Rheology 102 6 SUMMARY AND FUTURE RESEARCH 105 6.1. Conclusion 105 6.2. Future Research 106 6.2.1. Experimental Work 107 6.2.2. Numerical Algorithm 107 REFERENCES 109 BIOGRAPHICAL SKETCH 117 iv PAGE 5 Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirement for the Degree of Doctor of Philosophy COMPUTATIONAL STUDY OF LEUKOCYTE RHEOLOGY BASED ON A MULTILAYER FLUID MODEL By Heng-Chuan Kan May 1997 Chairperson: Dr. Roger Tran-Son-Tay Major Department: Aerospace Engineering, Mechanics and Engineering Science Knowledge of the rheological properties of leukocytes (white blood cells) is essential not only for the comprehension of microcirculatory flow dynamics, but also for the understanding of their functions and behavior in health and disease. These properties, together with the leukocyte's structural characteristics, determine the deformability of the cell, especially during large deformations, such as those involved in their release from the bone marrow or extravasation into the interstitium. Existing viscoelastic solid, Newtonian and nonNewtonian liquid drop models used to calculate the cell rheological properties do not adequately explain reported experimental observation. This is because such models do not take into account the morphology of the cell and its multi-layer structure. In this work, the deformation and recovery of a compound liquid drop (three-layer fluid model) are investigated by means of a mixed Eulerian-Lagrangian computational methodology which allows large viscosity and capillarity differences between layers. Analysis of the deformation process confirms published results and clearly indicates that the history of the deformation is stored in the compound drop in terms of the shape of its components. The recovery process is found to provide information on the hydrodynamics of compound drops PAGE 6 that is not supplied by the deformation analysis. It is discovered that a compound drop does not recover like a single phase Newtonian liquid drop except when the core is sufficiently deformed and its material properties are such that the core deformation/recovery time scale is compatible with that of the shell layer. A difference in the core and shell layer time scales causes an initial rapid recoil of the drop during which the shell fluid is the sole participant in the hydrodynamics, followed by a slower relaxation period during which the core and shell layer couple. The large deformation and subsequent recovery processes of compound drops also provide new insights into the complicated behavior of leukocytes. This three-layer fluid model describes the morphology of leukocytes better than existing ones. It consists of an outer interface (cell membrane), containing a shell layer (cytoplasm), and a core (nucleus). Based on the compound drop dynamics, it is concluded that (1) the nucleus plays an essential role in the cell response, and (2) the energy stored in the deformed interfaces can be important. The present investigation indicates that unless the nucleus and its deformation are included in the analyses, neither Newtonian nor non-Newtonian drop models for leukocytes can yield a reliable picture of the hydrodynamics of such cells. vi PAGE 7 CHAPTER 1 INTRODUCTION AND BACKGROUND 1.1. Motivation and Objectives Knowledge of the rheological properties of leukocytes (white blood cells) is essential not only for the comprehension of microcirculatory flow dynamics, but also for the understanding of their functions and behavior in health and disease. These properties, together with the leukocyte's structural characteristics, determine the deformability of the cell, especially during large deformations, such as those involved in their release from the bone marrow or extravasation into the interstitium. This deformability is also essential in the leukocyte's response to disease/infection in humans because it determines the ability of the cell to flow and deform in capillaries and/or migrate in tissue. Therefore, obtaining the correct mechanical description for the leukocytes has been a focus of research in the recent years. Existing viscoelastic solid, Newtonian and non-Newtonian liquid drop models used to calculate the cell rheological properties do not adequately explain reported experimental observation. This is because such models do not take into account the morphology of the cell and its multilayer structure. The three-layer (or compound drop) model has been suggested to be a good model for the leukocytes (Skalak et al., 1990; Dong et al., 1991; Hochmuth et al., 1993). This model gives a better representation of the morphology of the leukocyte and its interior over the existing ones. It consists of a cell membrane (outer interface), containing a cytoplasmic layer (shell), and a nucleus (core and inner interface). However, the three-layer model has not hitherto been approached analytically due to the limitation of available tools. The objectives of this work are: (1) to develop a general computational methodology for investigating the rheological behavior of multilayer drop systems, (2) to test the validity 1 PAGE 8 2 and accuracy of the algorithm, (3) to study the hydrodynamics of compound drops, and (4) to provide a better understanding of leukocyte rheology. 1.2. Importance of Compound Drops Dynamics The deformation and subsequent relaxation of a compound liquid drop is a fundamental fluid mechanics problem of medical importance. The underlying idea that makes the application of compound drops attractive and has spawned several new technologies (Jackson and Lee, 1991; Powell and Mahalingam, 1992; Batich and Vaghefi, 1994; Goeletal., 1994; Martin and Parthasarathy, 1995; Hassan et al., 1996; Thies, 1996) is that a core of active ingredient (whether solid, liquid or gas) can be protected by enclosing it within a shell of a second material. The active ingredient, which must not of course attack the chosen encapsulating material, is released at a determined rate, time and location to achieve certain localized effects. These include controlled release of drugs (Tai and Sun, 1993), flavors (Jackson and Lee, 1991; Gorski, 1994), enzymes (Hassan et al., 1996) or latent heat (Goel et al., 1994; Mulligan et al., 1996). While many aspects of compound drops have been studied in the literature (Johnson and Sadhal, 1985), including in particular deformation in an extensional flow (Stone and Leal, 1990) and stability (Landman, 1985; Oguz and Sadhal, 1987), the range of conditions under which compound drop dynamics has been investigated does not span the parameters of interest in this work. 1.3. Significance of Leukocyte Rheology The persistence and widespread occurrence of diseases associated with defects in the circulation of blood has largely provided the impetus for the study of blood rheology. Human blood is a suspension of erythrocytes (red blood cells), leukocytes (white blood cells) and platelets in plasma. Generally speaking, blood is considered as a Newtonian fluid in large vessels at shear rates greater than 100 s~' . When the shear rate is less than 10 s~' , blood behaves like a Casson fluid (Cokelet et al., 1963). However, for blood vessels less than 500 PAGE 9 3 \im in diameter, blood cannot be considered as a continuous fluid. Blood must be treated as a two phase fluid (core and plasma layers) in vessels with diameters between 30 to 500 ^m. Several effects (e. g., Fahraeus, and Fahraeus-Lindqvist) occur in that diameter range and become stronger as the vessel diameter decreases (Fung, 1993). For blood vessels less than 20 ~ 30 \im in diameter, the presence of individual cells cannot be neglected. The flow is creeping and the rheological properties of individual blood cells have to be taken into account. The rheological behavior of erythrocytes has been intensively studied and found to be adequately modeled by a Newtonian fluid enclosed in a viscoelastic membrane (Evans and Hochmuth, 1976; Evans and Skalak, 1980). On the other hand, the search for a suitable model to describe the rheological behavior of leukocytes continues. Rheological studies of leukocytes have been neglected because they occupy less than 1% of the total volume content of blood. However, since leukocytes have found to have a larger volume and a lower deformability than red blood cells, the number of studies has steadily increased. Leukocytes, because of their relative size and rigidity, play a significant role on the overall circulatory flow. They can affect in particular the apparent viscosity of blood and therefore the pressure drops encountered across vessels and capillaries (Pries et al., 1994). In addition, the ability of the leukocytes to flow and deform in capillaries and/or migrate in tissue is essential in the response to disease/infection in humans. Therefore, knowledge of the rheological properties of the leukocytes is critical for understanding microcirculatory flow dynamics and its effects in health and disease. Given the microscopic size of the blood cells (about 8Â—10 [im in diameter), experimental efforts face obvious hurdles. However, micro-manipulation techniques have been successfully developed not only for deforming but also for evaluating mechanical and chemical effects on blood cells (Schmid-Schonbein et al., 1981; Evans and Kukan, 1984; Needham and Hochmuth, 1990; Usami et al., 1992; Skierczynski et al., 1993; Tran-Son-Tay et al., 1994a). One common way of determining the rheological properties of the leukocyte has PAGE 10 4 been to deform the cell by aspiration into a micropipet (3 ~ 6 |xm in diameter) and then to eject the cell and observe the recovery characteristics. The recovery time and path then provide an estimate of the "apparent viscosity" of the cell. The apparent viscosity is an average rheological property of the whole cell, including all the inclusions. 1 .4. Observed Rheological Behavior of Leukocytes The most populous type of leukocytes in human blood is the neutrophil. The neutrophil contains a relatively small segmented nucleus, occupying about 21% of its volume, and a cytoplasm that is rich in granules. In contrast, in a lymphocyte, another type of leukocyte, the nucleus occupies approximately 44% of the cells (Schmid-Schonbein et al., 1980). A neutrophil is spherical in a passive, non-reactive state. The membrane of the neutrophil is highly ruffled, giving it an excess surface area (estimated at 110%Â— 120% of the undeformed sphere), and permits deformation at constant volume. In spite of the intensive research effort, to date, a variety of behaviors of the leukocyte cannot be explained clearly. For example, when a neutrophil is aspirated into a micropipet it flows into the pipet when the aspiration pressure exceeds a value called the critical suction pressure. Neutrophils rapidly aspirated into a micropipet with a diameter smaller than the cell diameter enter at an initial velocity much faster than the steady-state velocity. Also, as the rate of deformation increases, the cell viscosity is found to decrease, which has been speculated to be due to a shear-thinning effect (Needham and Hochmuth, 1990; Tsai et al., 1993). Needham and Hochmuth (1990) observe that after the cell enters into the pipet, the subsequent deformation occurs at a constant viscosity, independent of the aspiration pressure. However, as the excess aspiration pressure increases, the apparent viscosity of the cell decreases. Furthermore, following small but rapid deformations (in less than 0.3 seconds) of the cell, a rapid initial recovery (around 0.3 seconds) is obtained, followed by a slower phase (Dong et al., 1991). On the other hand, when the neutrophil is aspirated and held inside the pipet for a short time period of 5 seconds or less, and then expelled, the cell recovers faster PAGE 11 5 Authors Experiments Cell Apparent Viscosity Operating conditions ^(Poise) Dong etal. (1988) Aspiration ixeL-uvciy Small deformation 300 Evans and Yeung (1989) Aspiration Large deformation 1000-2000 Needham and Hochmuth (1990) Aspiration High deformation rate Low deformation rate 1000 2000 Dong etal. (1991) Aspiration Recovery High deformation rate Large deformation O AA O AAA zlHJÂ— zUUU Tran-Son-Tay etal. (1991) Dong and Skalak (1992) Recovery Aspiration Large deformation (short holding time) Tyarpe deformation (long holding time) Large deformation 800-1500 1000-2000 50-600 Hochmuth et al. (1993) Tsai et al. (1993) Recovery Aspiration Small deformation Low aspiration pressure High aspiration pressure 600 5000 500 Waugh and Tsai (1994) Aspiration High deformation rate 1000-1500 Table I. A list of the published values of neutrophil cytoplasmic viscosity based on a fluid model under different experimental conditions. than when held for a longer time with an initial rapid recoil (Tran-Son-Tay et al., 1991). Finally, neutrophils that are only slightly deformed in a pipet before ejection, yield much smaller values of apparent viscosity (Hochmuth et al., 1993) than cells that are significantly deformed. This is analog to the finding that the cell appears to be more viscous when it flows into smaller diameter pipets (Evans and Yeung, 1989). It appears that the rate of deformation, the extent of deformation, the excess aspiration pressure and the holding time all affect the cell recovery characteristics, resulting in different apparent viscosities. Table I summarizes the apparent viscosity values for neutrophils re- PAGE 12 6 ported by various groups. These results obtained from the aspiration (flow into a pipet) and recovery experiments. Therefore, these curious properties of the cell have yet to be fully explained by existing theoretical models based on the Newtonian or the non-Newtonian framework. 1 .5. Current Status of Leukocyte Modeling Much effort has been done over the years on obtaining the correct rheological model for leukocytes. In earlier studies, Bagge et al. (1977) show that granulocytes deform as they flow down a tapered pipet. Transit times down the tube were on the order of seconds, for pressure drops of magnitudes found in-vivo. When expelled from the pipet, the granulocytes recover their original spherical shape in times ranging from 20 to 80 seconds. In the work of Schmid-Schonbein et al. (1981), the leukocyte is modeled as a simple viscoelastic solid with a Maxwell element in parallel with an elastic element, as shown in Figure 1.1. Figure 1.1. Illustration of a standard viscoelastic solid model for the leukocyte. The model constants \i and k's represent the viscous and elastic coefficients. The elastic element is thought to impart shape memory, which enables the cell to recoil to the resting shape. However, Evans and Kukan (1984) show that the deformation of a leukocyte into a pipet is a continuous flow process with no approach to a static deformation limit. It is therefore suggested that the leukocyte be modeled as a Newtonian liquid drop with a PAGE 13 7 Figure 1.2. Illustration of a two-layer Newtonian liquid model with constant cortical tension for the leukocyte, y represents the isotropic surface tension coefficient. prestressed cortical tension (Yeung and Evans, 1989), as shown in Figure 1.2, rather than a viscoelastic solid. The shape memory then would be induced in much the same way as in drops, by ascribing a cortical tension (corresponding to surface tension in drops) to the membrane enclosing the cytoplasm. Based on this representation, Yeung and Evans (1989) simulate the flow of a neutrophil into a pipet by imposing a constant "ring" force at the pipet tip and a constant orifice pressure inside the pipet to satisfy the free slip condition. Although the liquid-like behavior of neutrophils is widely accepted, the precise nature of this hquid is not established. In their experiments, Needham and Hochmuth (1990) observed that the cell exhibits a non-Newtonian shear-thinning behavior when flowing into a pipet at high rates of deformation produced by large aspiration pressures. The non-Newtonian nature of the cell is also suggested by its contents, for example actin filaments, microtubules and other suspended particles, even though, their role may increase the cell viscosity only. It is clear from experimental observations, as summarized in the previous section, that the apparent viscosity changes with rate and extent of deformation. Thus, the behavior of neutrophils cannot be described in the same fashion as simple Newtonian liquid drops. These types of behaviors have been thought to indicate a "fading elastic memory" reminiscent of Maxwell liquids. Dong et al. ( 1 988) use a Maxwell liquid model with a constant corti- PAGE 14 Figure 1 .3. Illustration of a two-layer Maxwell liquid model with constant cortical tension for the leukocyte. cal tension, as illustrated in Figure 1 .3, to analyze the behavior of neutrophils flowing continuously into a pipet. The deviatoric stress x is obtained in such a model from: ^u + f^ij = 2Kj (1.1) where fX and k are model constants controlling the viscous and elastic properties of the Maxwell liquid. The strain rate ri is obtained from _ I avj 5Vj ^iJ 2^ ^ (1.2) where v is the fluid velocity and x represents spatial coordinate. However, experimental data could be correlated with their model only if the material properties of the Maxwell liquid were continuously varied as the deformation proceeded. Hence, their model cannot be considered to yield a fundamental description of the leukocyte in terms of a Maxwell fluid. Another attempt to obtain a non-Newtonian description of the leukocyte is made by Tsai et al. (1993) through a power-law representation of the cell apparent viscosity. The expression used by these authors is of the form: where \io is the viscosity at a reference shear rate r\oThe constant b is obtained empirically. PAGE 15 9 The power-law model essentially tries to link the shear rate dependent with the cytoplasmic viscosity to explain the shear-thinning behavior observed in the micropipet experiments. Although the power-law fluid model could be matched, for suitable initial conditions, to a fairly wide region of the deformation path of leukocytes, the rapid recoil region could not be predicted. 1.6. Summary Numerical study of leukocyte modeling has relied on methods based on semi-analytical approaches (Schmid-Schonbein et al., 198 1 ; Tran-Son-Tay et al., 199 1) or finite-elementbased structural analysis (Dong et al., 1988; Dong and Skalak, 1992). It is our opinion that the full scope of the numerical techniques available to the computational fluid dynamics community has not been brought to bear on the study of leukocyte modeling. It is clear from the review of current status of leukocyte modeling in the previous section that the existing viscoelastic solid, Newtonian and non-Newtonian liquid drop models fail to explain the complex rheological behavior of leukocytes. This is because such models do not take into account the internal structure of the cells. As Hochmuth et al. (1993) point out in the conclusion to their review of the state of affairs in leukocyte modeling, it seems that the correct model of the leukocyte can only be obtained when it is true to the morphology of the cell and its contents. However, since the internal structure of the cell is extremely complex, how can a morphologically consistent model be arrived at? One strategy, as suggested in the three-layer representation (Skalak et al., 1990; Dong etal., 1991; Hochmuth etal., 1993) shown in Figure 1.4, is to simplify the situation by factoring in the most important physical components, namely the membrane cortex, the cytoplasm and the nucleus and to temporarily ignore the contributions of the smaller intracellular particles including the microtubules and filaments. This is a good first order approximation since the cytoplasm is mainly composed of water. The cytoplasmic elements are assumed to only produce a higher cytoplasmic viscosity. PAGE 16 10 Middle layer (cytoplasm) Figure 1 .4. schematic of a three-layer or compound drop model for the leukocyte. In the three-layer model or compound drop under study here, the outer layer is the cortical region surrounding the cytoplasm. The cortex is comprised of an actin network. Its thickness is estimated to be of the order of 0. 1 microns (Zhelev and Hochmuth, 1994) compared to the cell diameter of about 8 microns. Thus, the cortical layer can be considered to be a thin shell with constant tension, estimated at 0.024 to 0.038 mN/m (Evans and Yeung, 1989; Needham and Hochmuth, 1990; Zhelev and Hochmuth, 1994). The second layer, namely cytoplasm, is of a more complex nature and under the three-layer viewpoint, some scenarios have been advanced. It could be a two-fluid mixture, with a highly viscous liquid intermixed with a less viscous liquid. Osmotic stresses may exist between these two liquids. When the cell is deformed, the less viscous liquid would flow preferentially. The resulting osmotic stresses would tend to restore the equilibrium state between the two fluids. With regard to the third layer, the nucleus, researchers (Skalak et al, 1990; Dong et al., 1991; Hochmuth et al., 1993) suggest that it can be modeled as a highly viscous, solid-like component which retains its shape when the cell is deformed. Therefore, the leukocyte is modeled as a thin cortical shell with a persistent isotropic surface tension (outer surface) surrounding a thick layer of cytoplasm (shell) with a nucleus inside (core and inner interface). This representation is also known as a compound drop model. PAGE 17 11 The layout of this thesis is as follows. In Chapter 2, the details of the salient features of the computational methodology employed are given. In Chapter 3, several well-defined and challenging flow problems are tested to assess the capabilities of the current numerical technique in terms of accuracy, robustness and flexibilities. The large deformation and recovery behavior of viscous compound drops are given in Chapter 4. This phase of the study elucidates the hydrodynamics of compound drops in relation to simple Newtonian drops and provides a basis for the development of a three-layer leukocyte model. In Chapter 5, the hydrodynamics of the three-layer model and its implications for the rheological behavior of the leukocytes are discussed. Finally, Chapter 6 summarizes the findings of this study and provides a list of suggestions for future research. PAGE 18 CHAPTER 2 THEORETICAL ANALYSIS 2.1. Background The challenges posed by fluid flow problems involving moving boundaries have resulted in a variety of computational techniques directed toward their solution (Shyy et al., 1996). The moving boundaries may represent material discontinuities in the flow field, and their dynamics is coupled with that of the flow. Numerical methods applied to such problems are required to follow the evolution of the often complex shapes assumed by the moving boundaries. To tackle these evolving interfaces, as a natural extension of adaptive grid methods for stationary complex geometries, algorithms based on body-fitted coordinates have been applied to moving boundary problems (Kang and Leal, 1987; Glimm et al., 1988). A boundary-conforming grid arrangement is convenient and accurate; the discontinuity across the interface is always maintained and interfacial conditions are applied at the exact locations without any smearing or redistribution. However, when the interface deformation becomes large, it is difficult to adequately resolve the geometrical complexities while maintaining desired mesh control. Mesh skewness and stretching may impact negatively on accuracy and efficiency of the flow solver. The boundary-fitted grid method presents further difficulties when interfaces merge or fragment. In simulating multi-layer fluid phenomena involving highly deformed interface shapes, purely Eulerian methods (Sethian, 1990; Gunstensen et al., 1991; Kothe and Mjolsness, 1992; Grunau et al., 1993; Sussman et al., 1994) have come to be widely used. In such methods, a fixed Cartesian grid is usually used, so problems associated with grid generation are circumvented. Eulerian methods based on the Volume of Fluid (Hirt and Nichols, 198 1 ), 12 PAGE 19 13 Level-Set (Sethian, 1990) or Lattice Boltzmann Equation (Grunau et al., 1993) have been developed to handle free-surface flows, such as wave-breaking, sloshing and splashing type problems, and multi-phase fluid flows in complex geometries. In purely Eulerian methods, the details of the interface shape are not explicitly needed, since the interface is deduced based on a computed field variable. The interface information can be obtained from an appropriate computed field variable in order to describe the details of the interface shape and curvature. An accurate estimation of the interface curvature is very important when capillarity plays a significant role in the interfacial physics. In tracking highly distorted interfaces with good accuracy, a combination of the strengths of Eulerian and Lagrangian methods has been shown to be useful. This mixed Eulerian-Lagrangian approach has been used in recent years for problems involving fluid-fluid (Unverdi and Tryggvason, 1992; Nobari et al., 1996) as well as solid-liquid interactions (Udaykumar and Shyy, 1996; Juric and Tryggvason, 1996). These methods employ fixed structured grids, and afford sifnplicity and availability of well-tested flow solvers. The only moving component is the interface which is a lower-dimensional surface overlying the fixed grid and is explicitly tracked. One limitation of the methods based on fixed grids is that since they have traditionally employed Cartesian grids, complex flow geometries are difficult to incorporate. A complex geometry facility is required for studying the interaction of the moving boundary with flows in arbitrary domains. The representation of a solid boundary on a fixed Cartesian grid layout calls for special treatment, since the control volumes through which the irregular boundary passes become fragmented into solid and liquid regions. This situation is illustrated in Figure 2.1. Quirk (1992) details some of the procedures involved in dealing with cell fragments or cut cells in the framework of compressible, inviscid flows around stationary obstacles. However, the calculation of fluxes is not presented in detail in that work. Also in connection with compressible, inviscid flows, Pember et al. ( 1 995) use local refinement and a redistribution procedure to enforce conservation in the vicinity of the solid-fluid boundary. Arbitrary-shaped boundaries running through a Cartesian mesh have also been dealt with using PAGE 20 14 (b) Control point (t)(i,j) 2 3 Figure 2. 1 . Illustrations of (a) the situation when a solid object in represented on a Cartesian grid, (b) a typical control volume encountered in the vicinity of the solid-liquid interface (Interfacial cell). finite element methods (Young et al., 1991). In Zeeuw and Powell (1990) and Bayyuk et al. (1993), a Cartesian grid is employed to track the motion of solid objects through an inviscid compressible fluid. Goldstein etal. (1993, 1995) employ a modified immersed boundary technique (Peskin, 1977) to impose no-slip boundary conditions over arbitrary-shaped solid PAGE 21 15 boundaries. This facility is useful for complex geometries, an important consideration for the spectral methods in use therein. Turbulent flows around solid obstacles have been calculated using this approach and the results have been shown to compare favorably with experiment. The method, as used in these works, treats the solid boundary as a series of steps conforming with the Cartesian grid layout, i.e., in piecewise constant fashion. In contrast, a piecewise linear representation of the free surface is used by Miyata (1986) for the simulation of wave breaking; the calculations are performed on a Cartesian grid using a finite-difference technique. In a finite volume setting, Udaykumar and Shyy (1995a) use a cut-cell technique to compute incompressible, viscous flows in arbitrary geometries. In that work, it is demonstrated that by means of a suitably defined consistent mosaic of control volumes, it is possible to obtain accurate solutions for flows on Cartesian meshes. However, the staggered variable layout in Udaykumar and Shyy (1995a) renders the measures adopted to handle the internal boundaries rather tedious. In this work, we present a fixed grid methodology for the simulation of flow through complex geometries with fluid-fluid moving boundaries, considerably simplified by the choice of a collocated variable layout (Lien and Leschziner, 1994). The performance of the flow solver in terms of stability and accuracy is thoroughly tested for various flow conditions. 2.2. ELAFINT: A Mixed Eulerian-Lagrangian Solution Algorithm 2.2. 1 . Characterishcs of the Numerical Methodology There are two distinct aspects to the computational framework, called ELAFINT, which stands for Eulerian-Lagrangian Algorithm for INterface Tracking, reported in Udaykumar et al. (1995a, 1996). These are: (1) A pressure-based flow solver to simulate the incompressible fluid flow over a wide range of Reynolds numbers. PAGE 22 16 (2) A procedure to deal with the interaction of fixed arbitrary geometries and evolving moving boundaries with the flow. The first aspect is straightforward and approached in a conventional manner using the robust pressure-based solution procedure detailed by Patankar (1980) and in the context of more modern developments by Shyy (1994). The scheme to track the moving fronts and describe complex geometries is summarized below, following which the flow solver and its discretization procedure are described. 2.2.2. Interface Tracking Procedure The interfaces are described by means of markers indexed sequentially and separated into objects. There is no restriction on the number of objects or on the end conditions that can be imposed on these curves. As shown in Figure 2.2, the objects can be open or closed. Furthermore, they can be fixed or moving. They can enclose different phases or materials with different properties. Once the interface description is obtained in terms of markers and their positions, the shape characteristics such as normals and curvatures can be extracted. The convention adopted by the algorithm is to consider the normal to the interface as pointing from phase II to phase I, as shown in Figure 2.2, such that when one traverses the interface along the marker sequence, phase II lies to the right. In order to handle multiplej+1 j+2 n PHASE I PHASE II Figure 2.2. Illustration of objects, markers and normal convention in the description of interfaces. PAGE 23 17 valued interfaces, the interface tracking algorithm employs piecewise polynomial fits to compute the necessary shape derivatives at each marker location. Facility exists to construct quadratic or cubic-spline functions Xi(s) andyi(s). The polynomial fits are obtained by parametrizing the curve as a function of arc length s so that, for the ith marker Xi = Xi(s), y,. = yiis) (2.1) Then, the normal and curvature are obtained by taking appropriate derivatives of these piecewise polynomials. Thus, the components of the normal are, ^ '-^ ^ + ^^2)1/2 ' y (x,2-Ky/)l/2 (2.2) The curvature x is, for the two-dimensional case: Â• n = Â— TTT O'.^ + Xs^)y^ (2.3) and the axisymmetric case: Xs |_ ysXsS ~ " y(ys^ + Xs^)y^ + x,2)3/2 (2.4) Interfacial markers are advected to their new positions by a Lagrangian translation. Thus, = + At Uf (2.5) y/+^ = y^l^+AtVf (2.6) where the subscript Â£ indicates the interfacial marker number and the superscript k implies the time level. It is demonstrated (Udaykumar and Shyy, 1995a, b; Shyy et al., 1996) that the interface tracking procedure can handle extreme deformation, including topological changes such as repeated mergers and breakups. It is also established that the interface tracking facility is versatile and accurate. Once the objects are described on the computational domain, the interaction with the flow field needs to be established. This is mainly done by providing the following information: PAGE 24 18 (1) The cell in which each interface marker lies. (2) The markers contained within each interfacial cell (i.e., control volume through which the interface passes). Since a Cartesian grid is employed, establishing these relationships is a straightforward task. Based on the information mentioned above, the flow solver can incorporate either changes in the control volume definitions (as in the cut-cell approach) or appropriate source term contributions (as in the immersed boundary technique). These two approaches enable the effects of the internal stationary as well as moving boundaries to be communicated with the flow field. In turn, based on the flow field developed, the motion of the interface can be determined. In Udaykumar and Shyy (1995b), moving as well as stationary solid-liquid boundaries are handled by the cut-cell technique. In this work, the cut-cell technique is retained for the stationary solid-liquid boundaries, while the immersed boundary technique (Peskin, 1977; Fauci and Peskin, 1988; Unverdi and Tryggvason, 1992) is used to treat moving fluid-fluid boundaries. This, along with the use of a collocated variable arrangement, enhanced the robustness and simplicity of the algorithm. Extensive testing of the algorithm has been performed for a wide variety of flows and the results thereof are presented elsewhere (Shyy et al., 1996; Udaykumar et al., 1996). 2.2.3. Formulation of the Pressure-Based Flow Solver Consider the two-dimensional planar or axisymmetric, incompressible Navier-Stokes flow. The governing equations in conservation forms are: Continuity equation: ^ + ^(PM) + -L^(yÂ» = 0 PAGE 25 19 Momentum equations: d{gu) dt dx ^ dx^ dx' ^ y'" dy^y ^df ^ ^' u (2.8) d{Qv) _^ _d_ dt dx dy ^ ^2 + dx^dx^ ^ y^dy^ ^ dy^ ^ (2.9) Energy equation: dt dx (2.10) where m=0 for planar and 1 for axisymmetric case. In the above, g is the fluid density, u and V are the fluid velocity components, T is the temperature, p is the pressure and fi is the coefficient of viscosity and 5Â„ and 5^ are the source terms in xand y-momentum equations which includes, for example, the buoyancy term under the Boussinesq approximation and the surface tension contributions discussed later. The boundary conditions to be applied on the interface for cases involving no mass exchange between phases are the continuity condition and normal stress balance: where the subscripts 1 and 2 represent the two adjoining fluids and / represents the interface, VÂ„ is the normal velocity component, y is the surface tension and x is the curvature of the interface, with the normal pointing from fluid 2 into fluid 1, as shown in Figure 2.3. (VÂ„), = iVn)2 = (Vn), (2.11) (2.12) PAGE 26 20 normal to interface Underlying fixed Cartesian grid Liquid 2 normal to solid surface Figure 2.3. Ilustration of the interaction of the solid-liquid and liquid-liquid boundaries with Cartesian grid. The sohd-liquid is treated using cut-cell technique, while the liquid-liquid interface is modeled with the immersed boundary technique. 2.2.4. Discretization Procedure of the Flow Field Equations The discretization of the governing equations is carried out using a control volume formulation, which with particular reference to two-dimensional planar flow. Similar treatment applies for the axisymmetric case. Consider the conservation law for the variable 0, defined to be (a) 1 for the continuity equation (b) u and v for the xand ymomentum equations and (c) specific enthalpy (h=CpT) for the energy equation. In the general case, we have The control volumes in the vicinity of a solid boundary passing through the grid are irregularly shaped, and in the general case can assume a five-sided shape as shown in Figure 2.1(b). In order to evaluate the fluxes through the cell faces, we integrate Eq. (2.13) over d(Q(p) dt + V Â• (^M0) = (V Â• rV(^) + s (2.13) PAGE 27 21 the control volume and employ the divergence theorem for two-dimensional problems to obtain dt dA + -0) Â• ndl = ()(rV0) Â• ndl + I SdA I I A (2.14) where the second term on the left hand side is now the line integral of the outward normal convective fluxes and the first term on the right side is the line integral of the outward normal diffusion fluxes through the faces of the control volume. We now proceed to discretize each of the terms in Eq. (2.14) for the control volume shown. Thus 1 k+\^k+\ _ Qkfkk dt (2.15) where the superscript k, k+1 indicate the time levels and dt is the time step size. A^v is the area of the irregular control volume. The time derivative is discretized in the backward Euler form leading to an implicit-in-time solution procedure which relaxes the time step constraints in comparison with explicit schemes. Next, 5 (){qu )-ndi = Y.'<^((^h)'^'di i (2.17)
(x,y) = ax^+by^+cxy+dx+ey+f. The coefficients are obtained by inverting the 6x6 matrix obtained by the known values of (p at six suitably chosen points in the required phase around the end of the probe. The inversion of the matrix is performed by LU decomposition. Thus, the required normal gradient of the field
+ a 4 a + ^^1"Â— otherwise (2.45) where a is the bound region around the interface which is equal to 2h in our computations. Usually this bound is taken to be proportional to the grid size. Based on the Heaviside function, the material properties are assigned as : ^ =/3, +O32-/3,)aG(0)Â« (2.46) where ^ is any material property such as viscosity and density. The subscripts 2 and 1 indicate the two adjoining phases, in line with our convention for the normal. In the RIPPLE approach (Kothe et al., 1992), which is based on a VOF representation of the phase field, 0 is taken as a color function dependent on the fluid fraction field, while in the level-set algorithm the ((> value corresponds to a distance function. In these methods |