You are on page 1of 40



Introduction to Multiphase Fluid Dynamics

Scope of the Book
This book focuses on simulating engineering-level multiphase fluid dynamics. The term multiphase flow was coined by the late Prof. Soo of the University of Illinois in 1965 and includes fluid dynamics motion of various phases. The particle composition can be solid, liquid, or gas, whereas the surrounding fluid can be a liquid or a gas. The flows can also be categorized by concentration. Multiphase flows are considered dispersed if the changes in particles trajectories are dominated by their interaction with the surrounding fluid. In particular, the particle motion is chiefly governed by particle inertia and the fluid dynamic surface forces (e.g., drag force) and body forces (primarily gravitational). This regime is generally associated with the particle volume fractions that are less than 10% and particles sizes which are small compared to that of the overall domain. Volume fractions greater than 10% typically lead to dense multiphase flows, where particle movement is dominated by particle-particle interactions (e.g. by collisions or contact mechanics). On the other hand, volume fractions of 0.1% or even less tend to yield sparse flows, where particles are so dilute that particle-particle interactions are negligible and particles concentrations are too low to affect the surrounding flow. This text will focus on dispersed multiphase flows, which include the subset of sparse regime but excludes the dense regime. Additionally, the surrounding fluid will be generally assumed to have a Newtonian viscosity and be in a continuum with respect to the particle. Furthermore, for sake of brevity the focus will be on fluid dynamics and not on the influence of heat and mass transfer. Dispersed particle dynamics are common in many engineering systems. Because of this, robust and accurate description of multiphase flow is important as knowledge of the characteristics can lead to increases in performance, reduction in cost, and improved safety for engineering systems. In recent years, Computational Fluid Dynamics (CFD) has become an indispensable predictive tool for gathering information to be used for design and optimization for fluid systems. Thus, the combination of CFD and multiphase flow has emerged as an important research area, but one with unique characteristics and issues. For example, multiphase flow approaches include a wide variety of approaches such as variations in reference frame representation (Eulerian or Lagrangian), phase coupling (intra-phase and/or inter-phase coupling), and particle/flow detail (e.g., high resolution around a single particle or bulk description of thousands or millions of particles). These different approaches are associated with large variations in computer usage and predictive fidelity. Thus, simulation of dispersed multiphase flow requires careful consideration of both the flow regimes and the relevant numerical approaches. As such, this text is designed to relate various multiphase fluid physics to the appropriate governing equations and numerical schemes. For additional information on such topics and on topics outside the scope of this book, there are several excellent publications to which the reader is referred. In particular, multiphase fluid physics are discussed in the texts of by Clift et al. (1978), Wallis (1969), Soo (1990), Crowe et al. (1998), Brennen (2005) and Michaelides (2006). Reviews focusing on computational multiphase flow include Elghobashi (1994), Faeth (1987), Shirolkar et al. (1996), Shyy et al. (1997) and Tomiyama (1998), as well as Prosperetti & Tryggvason (2007). 14

Important studies which focus on the mathematical treatments of two-phase flow include Drew & Passman (1998) and Prosperetti (1998) while heat and mass transfer aspects are discussed by Williams (1965), Oran & Boris (1987), Kuo (1986), and Sirignano (1999).



Examples of Multiphase Flow in Engineering Systems

Dispersed particle-laden multiphase flow is important to many engineering applications including aerospace, atmospheric, biological, chemical, civil, mechanical, and nuclear systems. Table 1.1 lists some engineering applications for the four primary conditions of dispersed multiphase flow: solid particles in a gas, liquid droplets in a gas, solid particles in a liquid, and gas bubbles in a liquid. Table 1.2 illustrate different categories of particles as a function of typical size ranges, varying from scales less than a micron to scales greater than a centimeter. The table indicates that multiphase fluid dynamics can span a large range of length scales. To indicate the importance of numerical analysis of the dispersed dynamics, examples from three types of systems which involve multiphase flow are discussed in the following. These three systems include: energy systems (associated with the generation of energy or mechanical work), processing systems (associated with the production and transport of fluids), and environmental systems (associated with natural and biological transport processes). 1.2.1. Multiphase Flow in Energy Systems

Circulating Fluidized Beds Combustion of coal and other solid fuels for power generation is often most efficiently accomplished for power generation via a fluidized bed. In such a system, the solid fuel is fed into the bed after crushing and the injected upward-flowing gas flows causes the particles to be approximately suspended and exposed to high convection velocities and heat transfer rates which, in turn, cause rapid combustion (Basu, 1999). At low flow rates, the fluidized bed is quite compact with particle volume fractions on the order of 20-50%. This condition represents a dense mixture, i.e., particle-particle fluid dynamics and collisions dominate the bed dynamics. However, the emphasis has recently been on circulating fluidized beds (Fig. 1.1) which are more-efficient and generally yield fewer pollutants due to more complete combustion. In these fluidized beds, the upward gas flow rates are much higher such that the particles are not just suspended at the bottom of the chamber, but instead burn as they move upward and then recirculate generally along the side walls of the chamber. To avoid particles exiting the system, vertical U-channels can be employed to capture the particles and force them back down into the combustion region for continued recirculation (Fig. 1.2a). Smaller particles which exit the primary burner can have secondary recirculation in an economizer before they are re-injected into the lower bed. For finer particles which are not re-circulated at all, downstream stages employ other filtering devices based on centrifugal forces and settling velocity, such as the cyclone shown in Fig. 1.2b, whereas later stages use porous materials to remove the finest (smallest) particles. Employing various stages of filtration is intended to allow maximum particle removal while avoiding significant pressure losses. In these circulating fluidized bed systems, the lower-most region is dense but most of the combustion and filtration takes place in regions which can be reasonably described as dispersed. As such, these regions are amenable to the numerical methods discussed in this book. Numerically simulating the filtration systems can be important to allow optimization


with respect to efficient particle removal performance. Such simulations can be quite difficult due to the systems geometric complexity and wide range of scales, as well as the variety of relevant physics including many types of particle-fluid interaction (turbulence, convection, settling, and re-suspension).

Air-Breathing Engines Air breathing engines primarily use internal combustion to produce mechanical power based on chemical reactions of fuels with the oxygen gas. Two types of engines are most prevalent: cyclic compression engines which employ piston motion for intermittent compression; and jet-turbine engines which employ streamwise steady compression by passing the inflow though rotating blades. Piston-based engines are primarily used for ground transport, while turbine engines are primarily used for air transport. However, piston engines have been used for some aerospace operations and turbine engines are increasingly being used as ground-based auxiliary power generation as they are reliable and clean-burning. Air-breathing engines generally involve multiphase flows in three aspects of the process: the presence of particles within the air intake, the injection and distribution of the fuel spray in the combustion chamber, the production and evolution soot particles. With regard to the air intake aspect, particles entering through the air stream can be a problem for aircraft jet engines since mechanical air filter are impractical due to the high flow speeds. To prevent this problem, large particle ingestion causes erosion of the turbo-machinery components during near-ground operations in sandy or dusty regions. This problem can be minimized by flow curvature (Fig. 1.3a) but the performance of such filtering is sensitive to particle concentration, size and trajectory which can be difficult to predict. Also water droplet ingestion at near freezing temperatures can lead to ice formation on compressor blades (prior to any combustion). This can be a concern as the accreted ice can degrade compressor performance and also break-off in pieces which can be hazardous to engine performance (may even lead to flame-out conditions). Therefore, prediction of the particle or drop deposition rates on surfaces is important to predict turbo-machinery erosion and icing problems so that designers can determine effective and efficient countermeasures. The fuel spray is another important multiphase aspect, and results in a complex cloud of droplets created from pressurized ejection (Fig. 1.3b). The size, velocity, and spatial distribution of the spray are important in describing the performance of the system. For example, drops which are too large can impact the chamber side walls or not undergo rapid vaporization and combustion. This can have deleterious effects in terms of combustion efficiency and pollution emission, including an increase in NOx, CO and soot particulates. Modern engines take advantage of advanced spray nozzles, turbulence generation, swirl, and often create special lean zone regions minimize emissions. The third multiphase aspect of airbreathing engines is the production and evolution soot particles resulting from the combustion process which are emitted in the exhaust. Particulates in the exhaust as an emission source is receiving increasing attention because of environmental concerns associated with fine and ultra-fine particles. Soot production is sensitive to the spray dynamics and combustion in the combustion regions, so that these multiphase problems are coupled.


Space Systems Rocket engines flows are a common aerospace application which involves multiphase flow. Rocket engines create high-pressure gasses which are exhausted at supersonic speeds to create thrust, and typically operate at high altitudes insufficient for air-breathing propulsion, e.g. in space. The two most common rocket engines are solid-propellant engines and liquid-propellant engines. Solid-propellant engines allow very high thrust loads useful for take-off or single-use missiles and are based on surface-burning for the combustion process. The solid-propellant fuel often includes a binder with aluminum particles since their high energy release can increase thrust levels. The particles break-off from the internal sidewall surfaces as the binder around them burns away, and they are then transported along the core flow of the combustion chamber (Fig. 1.4). The particles undergo complex changes, since they transition from solid to liquid and then evaporate and combust. However, some conditions can yield particles which do not vaporize and lead to undesired build-up or ablation in the nozzle region. Therefore, simulating the particle trajectories is crucial to model and optimize the rocket performance, especially since the hostile high-temperature environment of the rocket chamber makes experimental access difficult. For liquid-fueled rockets, multiphase flow issues include cavitation and evaporation as the liquid fuel and oxidizer is pumped from tanks and converted into gas. Another aerospace application which involves multiphase flow is space-craft thermal management, needed to handle the severe variations in ambient conditions. In particular, boiling is used in many thermal management systems as it is a highly effective heat transfer mechanism and can be employed in many processes: component cooling, energy-conversion systems, and storage systems for cryogenic materials and life-support fluids. However, the complexity of microgravity conditions complicates the fluid physics since the gravitational pull-off process is eliminated and convection is needed to transport bubbles. Microgravity is also difficult to replicate experimentally for full-scale conditions. As such, computational approaches are powerful assets for design and analysis of the multiphase aspects of such systems.

Nuclear Reactors Nuclear power plants rely heavily on heat removal by water, which will often boil due to the high temperatures. Thus, vapor bubble flow convection (Fig. 1.5) has been widely studied due to its importance in the safety and thermal efficiency of these systems. Simulating these flows can be complex since they involve non-spherical bubble shapes with changing volumes that may be subjected to non-linear lift and wall-interaction forces for pipe flow. Furthermore, the random passage of bubble wakes in the flow causes additional velocity fluctuations beyond that of single-phase turbulence. The phenomenon has been termed pseudo-turbulence or burbulence (since it can occur even in laminar flow conditions and is often most apparent for bubbly systems). These phenomenon are not well understood


individually so that the combination of these aspects has been a great challenge to numerical simulation of these flows. 1.2.2. Multiphase Flow in Processing Systems

Chemical Reactors Multiphase flow is important in many liquid chemical reactors since they often employ bubble columns and air-lifts (Figs. 1.6) to cause high mass and/or chemical transport rates among liquids, gasses, and solids. The reactor processes include fermentation, biooxygenation, production of cell cultures and proteins, waste-stream treatment, etc. The transfer rates are a function of the bubble size distribution, mean and fluctuating velocities and trajectories, bubble coalescence and break-up, as well as particle-bubble interactions. Advanced bubble columns employ impellers and non-linear configurations which further complicate the fluid dynamics. The integrated fluid physics for these phenomena and configurations are not fully understood. As such, CFD plays a critical role in the design. While most systems focus on chemical reactions outside of the bubbles, there is a related process which conducts chemical reactions inside a bubble, called sono-chemistry. This phenomenon is based on acoustically-driven micro-bubble oscillations which momentarily create very high pressures and temperatures during the contraction phase. This is accomplished by driving the bubbles with high-amplitude pressure fluctuations near their natural frequency. As a result, the bubbles will undergo a rapid and violent collapse during the compression portion such that the gas inside the bubble may be driven to temperatures of more than 10,000oK and pressures of more than 10,000 atmospheres. In fact, measurements of inner bubble conditions in a sulpheric acid bath have yielded temperatures as high as 20,000oK (Suslick & Flannigan, 2005). As a result, the gas inside the bubble becomes a plasma and emit light which can be seen with the unaided eye. This emission is termed sono-luminescence and an example is shown in Fig. 1.7a. These intense chemical reactions are valuable for a number of applications including degradation of organic species in contaminated water and the production of protein micro-spheres (Fig. 1.7b). These micro-spheres can be used for drugdelivery, biomedical imaging, and blood substitutes.

Pumping Systems Liquid pumping by propellers is used in many systems to achieve high throughput. However, the associated high blade speeds can cause the local pressure to drop below the local vapor pressure, giving rise to a multiphase flow phenomenon known as cavitation. If their presence is extensive, the emergence of cavitating vapor bubbles can reduce the overall pumping performance. In addition, cavitating bubbles may undergo violent condensation collapse on the blade surfaces which can cause an intense localized pressure wave and cause surface pitting. Therefore, CFD is often used to design advanced propellers systems so that cavitation will be primarily occur in the trailing vortices instead of on the blade surface.


Another manner in which multiphase flow aspects are important to liquid pumping is the use of polymers or bubbles to reduce the mean wall skin friction losses (and thus pumping power required) in a channel flow. In particular, polymer chains added in very small volume fractions (about 10-4) have been found to reduce the mean skin friction by more than 80% by suppressing turbulent transport near the walls (Fig. 1.8). One of the most famous examples of this is the successful use of drag-reducing polymer along the Alaskan oil pipe-line (they have also been used for other flows including fire hose nozzles). Similarly, the injection of microbubbles has been found to reduce the skin friction losses (also shown in Fig. 1.8), though the required volume fractions are typically larger (about 10-2 or more). While not as robust and understood, the bubbly flow drag reduction eliminates some environmental and cost concerns associated with polymers. Therefore, this technique is also being actively explored for implementation into commercial or military vessels as it could be highly beneficial in terms of reduced fuel for transportation. Since the fluid physics of micro-bubble drag reduction is not well understood, multiphase CFD with micro-physics modeling is being employed to determine how best to incorporate phenomenon into full-scale engineering systems. The above drag reduction systems involve the interaction between wall turbulent features and particle motion. This aspect of multiphase flow is also important for analyzing systems which involve transporting particles in pipes. Particle conveying systems are called slurries when they are liquid-driven and pneumatic transport when they are air-driven. In both cases, the dynamics of particles can involve regions where convection effects are important (though, much of the transport involves dense multiphase regions). Predicting the trajectories of particles in pipe flows (including sliding, bouncing and re-suspension dynamics) is important to understanding and minimizing wall deposition, especially for turbulent flows. Even with very small particles and laminar flows, wall deposition can occur due to random Brownian motion resulting from the collisions of the surrounding fluid molecules of the particle surface. This phenomenon can be important in clean rooms where the flow handling (filtering and convection) is designed to avoid particle contamination on relevant surfaces to allow microand nano-scale fabrication.

Sloshing Dynamics Multiphase flow is also important to the transport of liquids in containers where sloshing may occur. Such transport is particularly important for Liquid Natural Gas (LNG) shipping as oceanic vessels have steadily increased in size to respond to mounting international energy demand. The sub-cooled and pressurized LNG tanks can suffer high impact sloshing loads when large unsteady ship motions occur (Fig. 1.9a). Small-scale experiments (Fig. 1.9b) have qualitatively shown that these impact pressures can be substantially cushioned by the presence of dispersed bubbles due to the resulting compressibility of the mixed fluid. However, it is well-known that the experiments can not be properly scaled to quantitatively replicate both wave motion and compressibility. Therefore, attention has been placed on simulations, where challenges lay in successfully predicting the unsteady and wide range of bubbles sizes and locations. In particular, one must pay attention to the processes of formation (by wave cresting and plunging), spreading (by the complex wave motion), and removal (by buoyancy rise out of


the liquid). In the latter two processes, dispersed particle dynamics can be used to simulate the conditions. Cavitation and condensation can also occur due to the rapid changes in pressure, which further complicate the physics and numerics.

Sprayer Systems Sprayer systems include applications such as surface coating/painting, powder production, particle filtering, hazard mitigation, etc. For surface coating, it is desirable to apply a spray with droplets of a particular size range since large droplets tend to cause surface anomalies and small droplets may not deposit on the target zone and instead be carried away with air currents. Thus, numerical tools are used to help design optimum spray conditions and impact rates for a given target distance and target shape. A particularly complex multiphase flow process is that of extreme-durability single-crystal solid metallic coatings. To achieve such coatings, plasma sprays are employed (Fig. 1.10) whereby a metal powder is injected into a supersonic hightemperature gas jet and the particles become molten or semi-molten as they are entrained. The convected particles ideally impinge on the target surface with an appropriate speed and size to liquefy (with little or no splashing) and then crystallize. This results in surface coatings with high toughness, such as used for gas turbine blades to increase their wear resistance. Thus, the plasma spray process has been the subject of substantial CFD to determine and optimize the particle and jet properties for given powder material and plasma temperatures. Spraying systems have also been used to produce powders, i.e. fine solid particles. In an industrial spray dryer (Fig. 1.11a), this is typically accomplished by first emulsifying the material into liquid form and then downward spraying the liquid mixture into a large chamber in which dry-gas is injected. The droplets from the spray then solidify in flight and settle to create a powder at the bottom of the chamber (Fig. 1.11b). Using multiple stages increases the system thermal efficiency and particle size uniformity, and also prevents overheating of powder particles along their trajectories. Achieving a desired powder size and consistency is critical for the production of foods, detergents, pharmaceuticals, ceramics, herbicides, pesticides, etc. However, achieving the proper agglomeration (or avoiding any agglomeration in the synthesis of some micro- and nano-particles) and avoiding sidewall deposition is often difficult and conducted empirically. Because of the complex physics and geometries of these systems, numerical methods have begun to play a larger design role to help improve systemlevel quality and efficiency.


Multiphase Flow in Biological and Environmental Systems

Particle Motion in the Blood Stream and the Body The bloodstream is the most important liquid flow within the human body and it may contain particles by disease or design. One naturally occurring multiphase condition is the transport of cholesterol or calcium particles in the bloodstream. Unfortunately, these particles can undergo wall deposition and thus can accrete on certain vessel surfaces (perhaps due to low


shear-stress and certain bio-chemical conditions). The resulting growth is called plaque deposition and, if severe, can cause a blockage in the blood vessel (called a stenosis) which can lead to poor circulation and serious health problems. Perhaps more worrisome is that large pieces of the plaque deposition can break-off and move into smaller arteries causing a more profound and serious blockage. For example, if a strong blockage occurs in on the brain blood vessels it may result in an intracranial aneurysm, i.e., a stroke. Therefore there is strong interest in predicting and understanding the mechanisms of plaque movement and deposition in bloodstreams. This is a complex problem as the flow is pulsatile with three-dimensional geometries (especially near grafts) and may also be transitional if not turbulent in the large vessels near the heart. Note that the flow in the very small vessels (less than 1 mm) is not easily considered as continuum since the blood cells themselves are discrete bodies (about thick discs about 8 microns in diameter) with concentrations as high as 40% moving in a liquid medium called plasma. Thus, capillary vessels may be thought of as a dense multiphase flow (which is generally beyond the scope of this book) even when there are no solid particles present. Another bio-medical multiphase flow associated with bodily fluids is the artificial introduction of particles for purposes of imaging or drug delivery. For example, encapsulated micro-spheres (which contain a liquid interior inside a shell that can dissolve or be broken with acoustic radiation) can focus drugs in targeted areas. The movement of the particles in the body to areas of importance can be controlled by magnetic force if the particles posses a small charge. The shell can also be chemically engineered to attach to specific molecules for the purpose of drug deliver. Particles can also be injected to serve as a contrast agent, i.e. a diagnostic marker within the body. In particular, micro-bubbles (about three microns in diameter and composed of helium, carbon dioxide or other gasses) can be injected into the body. The small bubble size allows for both very low volume concentrations and very high number concentrations (e.g., a 0.0001% volume fraction corresponds to 20,000 bubbles per cc). The low density of the gas makes them a good contrast-agent for angiography, ultrasound, and MRI imaging. For example, they can be used to can show blood-flow patterns in the heart muscle to determine if there are any regions with poor circulations. The micro-bubbles are often coated with a protein, lipid, or polymer layer termed a shell (Fig. 1.7b). This shell is helpful in that it prevents dissolution of the gas while being thin and elastic enough to oscillate due an applied local sound field (Fig 1.12a). This property allows the bubbles to be good contrast agents for acoustic imaging in non-vascular regions (Fig. 1.12b). It also allows the micro-bubbles to be targeted (e.g., concentrated on a vessel wall) by using acoustic radiation forces and is helpful. The computational approaches for micro-bubbles transport share attributes with those associated with other micro-particles in the blood-stream as discussed above. However, the additional issue of (symmetric and asymmetric) size oscillation is critical to their performance and thus individual bubbles have been examined for their response to pressure oscillations.

Particle Motion in Respiratory Passages The introduction of particles to the gas flow of the human respiratory system is a related and important dispersed multiphase problem which can be addressed by numerical methods.


For instance, oral spray-drug delivery by inhalers has been modeled with particle simulations (Fig. 1.13) to assess the aerosol trajectories. Particular interest is placed on the ability of the aerosol droplets to avoid wall deposition so as to maximize air-borne delivery to the lung cavity. In the case of airborne pollutants, the opposite is desired, i.e., deposition of particles in the mouth and nasal passages is favorable since it avoids particulate transport to the sensitive tissue of the lungs. The probability of particle deposition in the respiratory system varies significantly with particle size. For particles which are greater than one-micron in diameter, inertia impact and curved passages dominate the deposition rates. For smaller particles, say 0.1 to 1 micron in diameter, the terminal velocity caused by gravitational settling becomes dominant. The deposition of ultra-fine particles (0.1 microns in diameter) is dominated by Brownian diffusion, i.e. diffusion by the random molecular collisions of the gas with individual particles. Since this deposition rate is very weak due to the nearly negligible inertia at these sizes, such particles can reach the lungs where they can permanently settle and do significant harm, particularly among children and the elderly. Numerical techniques to simulate depositions for the large range of possible particle must adapt to the controlling physics, e.g. different approaches are generally required for modeling Brownian vs. turbulent diffusion.

Atmospheric Solid Particulates Because of the concerns to the respiratory system of humans (and the potential impact on animals and plant-life), there has been a significant focus placed on modeling the concentration distribution of micron and sub-micron sized particles in the atmosphere. The concern has primarily been associated with man-made emissions, commonly regarded as particulate pollution. As such, air quality is often characterized in terms of the micro-particle concentration (Table 1.3). In 2003, over 25 urban areas were found to have Air Quality Index (AQI) levels above the unhealthy threshold of 100 for more than 10 days. On the other hand, particle concentrations levels far from any pollution sources are typically less than 1 g/m3, which corresponds to AQI levels of less than 10. Another motivation to understand atmospheric transport of particles is to be able to assess the threat of terrorists placing a device which emits hazardous particles. In particular, the consequences of the release of a highly lethal powder (e.g., anthrax) or the use of a dirty bomb (where an explosive device is used to scatter radioactive materials) are now being carefully studied. Numerical techniques are thus being developed to understand the evolution of particle concentration in the atmosphere, both locally and globally. Generally, the particle terminal velocities are small compared to the atmospheric velocities and can be quite reasonably be considered as passive scalars transported with the continuous-phase, i.e. the particle velocity is can be assumed equal to the surrounding gas velocity. Furthermore, the concentrations of the particles are generally so small that their effect on the surrounding air fluid dynamics is quite negligible. Thus once the particle source is identified, the primary numerical issue is the proper modeling of the air flow velocity field. This can be quite complex in that the atmospheric boundary layer includes turbulence which produces highly unsteady threedimensional flowfields. The modeling for these particles can extend to the entire planetary


atmospheric boundary layer as shown in Fig. 1.14a. It can also focus on interaction within urban structures as shown in Fig. 1.14b, where one must consider the effect of the street canyons between buildings that can have increased particulate concentration due to local sources and re-circulating (trapped) air. Furthermore, interaction of such solid particles with water droplets from clouds or rain is a three-phase flow aspect which is often important to understand the evolution of atmospheric particle concentration. There are also natural events which can cause very high particulate concentrations, including particles emitted from a large fire or by gigantic dust storms (Fig. 1.15). Predicting the particle concentration and movement for these types of events is important to environmental modeling and to determine flight space regions suitable for aircraft operation. Furthermore, the deposition of these particles is of interest particularly with respect to the impact on surface water and agricultural areas. When the particles are airborne, the small particle sizes allow use of a passive scalar approximation for the particle velocity, though sometimes a terminal velocity is superimposed for larger particles. In some cases, the particle concentrations can be high enough that they increase the effective density of the particle-laden air and thus cause the mixture to move downward as compare to unladen air flows which are less dense in comparisons. As such, numerical formulations with very small particles in high concentrations can assume a single mixed-fluid for the multiphase flow field whose density is a function of particle concentration. Another natural even involving particle injection into the atmosphere is the case of volcanic eruptions. These can be particularly violent and far-reaching since ash particles can eject upwards at speeds of hundreds of meters per second and disperse on a continental scale as a result. A particularly deadly condition can occur when the cloud being ejected from the volcano reduces in velocity but retains enough particle mass content that it has a significantly increased effective (mixed-fluid) density. In this case, the effective density of particle-laden becomes greater than that of the surrounding unladen air so that the particle cloud is driven along the ground forming what geologists call a pyroclastic flow (Fig. 1.16). Such flows can have temperatures up to 600oC with speeds up to 100 km/hr can propagate up to a distance of 60 kilometers. Thus pyroclastic flows have the capacity to rapidly incinerate and asphyxiate such that they have killed more people than any other aspect of volcanic eruptions (Chester, 1993). The particle sizes are small often enough to be considered fluid tracers but with a terminal velocity, however their high concentration can substantially modifies the particleladen gas flow, e.g. increasing the effective mixed-fluid density by as much as 10%. Because of this, these flows are often referred to as particle-driven or gravity-current flows since the density stratification (compared to neighboring particle-free air) drives the overall flow. Numerical simulation of these flows is complex due to the large range of flow speeds and particle sizes as well as the strong unsteadiness and large scales of the interaction.

Atmospheric Liquid Particles The formation and coalescence of drops in clouds is another area in which multiphase CFD is being applied to understand the controlling processes. The prediction and understanding of rain cloud formation is important to weather prediction and climatology. For example, the


mechanisms which cause clouds to rapidly transition from small droplets (about 1-20 microns) into large rain drops (about 100-1000 m) during the formation of a thunderstorm are underpredicted by conventional collision theory. Numerical simulations are being used to help develop models which may better explain and predict this behavior (Chun et al. 2005). For aircraft, droplets concentrations and size impact the visibility for pilots, but take on additional significance at near-freezing conditions as they can rapidly accrete on wings and other aircraft surfaces (Fig. 1.17a). The accretion, commonly known as aircraft icing, can result in a substantial degradation to the aerodynamic performance, such as a two-fold reduction in lift under some conditions, which can present a substantial safety hazard. The resulting loss of maneuvering control has been responsible for several fatal aviation accidents. The influence of particle inertia, turbulent dispersion and the surface roughness can be significant in terms of the net accretion, and CFD has played a major role in increasing the understanding the associated droplet deposition. In addition, computations are now used to simulate ground-based experiments whereby droplets are injected upstream of a refrigerated wind tunnel test section to study the droplet deposition and icing physics in consistent and well-defined conditions (Fig. 1.17b).

Sediment and Bubble Transport in Bodies of Water Sediment transport in natural bodies of water can play a significant role in the health and maintenance of rivers and coastal areas. For example, silting (accumulation of sediment) can dramatically reduce the capacity of rivers and has been cited as one of the reasons for some of the more devastating floods in Asia and elsewhere. Coastal land erosion is another problem that has become more severe in the last few decades since it can lead to substantial losses of animal habitats and affect human population centers. Much of the important physics of sediment transport for silting and erosion occurs in the dense shear region, which separates the static sediment bed from the water above. In particular, the lifting-up of particles (typically by drag and lift forces) along with the subsequent convection and re-settling are important dispersed multiphase fluid dynamics which control this mobile layer. These aspects are further complicated by the wave-induced and current-induced turbulent boundary layers. Since theoretical analysis is difficult to apply for such conditions, dispersed multiphase numerical techniques are being used to understand these overall processes. Particles in artificial bodies of water are often sought to be filtered out. In particular, the environmental cleaning of waste streams often involves the addition of chemicals to form a precipitate. To remove the unwanted solid precipitate, bubbles are introduced into the sediment stream to achieve air flotation, whereby the bubbles attach to the suspended particles making buoyant agglomerations which are more easily filtered (Fig. 1.18a). The bubbles are generally generated by injecting high-pressure air-saturated water through nozzles, which results in air bubbles coming out of solution in diameters of about 30 to 300 microns. The resulting bubble/particle agglomerates rise to the top of the flotation tank where they are collected in a sludge blanket that is later removed while the particle-free water passes below a baffle (Fig. 1.18b). The sludge blanket is then removed to produce waste-free water. To optimize the system it is desired to use bubbles small enough to filter out the smallest particles,


but large enough so that the agglomerates are sufficiently buoyant to ensure a reasonable flow rate (i.e., processing speed). This multiphase process is quite complex due to the nonlinearities associated with turbulent mixing, the large variety of particle shapes and sizes, and the kinematics associated with collisions that cane lead to break-up or agglomeration. Bubbles are also injected in plumes for reservoirs and wastewater treatments to increase oxygen content and control pH levels. In these large-scale water bodies, the upward rise of the bubbles causes the surrounding liquid to also rise due to turbulent mixing and entrainment. This entrainment is important to overall effectiveness of the plume. For numerical treatment, it is important to take into account this coupling between the dispersed- and continuous-phases flows. This can be complicated by presence of large bubbles which can coalesce, deform due to inertia effects, and exhibit non-linear (zig-zag) trajectories. Interestingly, biological bubble injection has been used by humpback whales and dolphins to catch fish. For example, a whale (or several whales) will dive beneath a group of prey and then slowly spiral upwards while blowing bubbles (Fig. 1.19). This creates a cylindrical wall of bubbles, which traps the prey in a virtual net, allowing a whale to then swim and feed though this region. Leighton (2004) suggested that the reason for the trapping is that the bubbles create an acoustically-insulated environment. Bubbles of specific sizes and concentrations are also released along the hull of military ships for stealth. This process is called masker and is intended to modify or reduce the underwater acoustic detectability of the vessel. To predict this multiphase flow, numerical simulation of the bubble trajectories must account for the presence of a turbulent boundary layer along the ship and a downstream turbulent wake. If the acoustical properties are also to be simulated, the oscillation (dynamics) of the bubble diameters due to the unsteady surrounding pressure field must also be simulated.



Basic Terminology and Assumptions

To describe multiphase flow regimes, a few essential definitions are necessary. A "particle" will be defined herein as a relatively small unattached body immersed in a flow. In terms of shape, often particles are spherical or can be reasonably approximated as such, so that a simple geometric diameter (d) can be used to describe their volume. In terms of size, a macroscopic length scale of the continuous-flow domain (D) can be defined (such as a channel width or a jet diameter) and herein the particle size will be assumed to be small in comparison such that dD (Fig. 1.20a). A particle can be a solid particle (such as dust or soot) or a fluid particle (such as droplet or a bubble), and will generally interact with the surrounding fluid. Since each particle is assumed to be surrounded by a fluid which otherwise fills the domain, the particles in general will be referred to as the dispersed-phase and the surrounding fluid will be referred to as the continuous-phase. Most of the analyses herein will assume a single dispersed-phase and a single continuous-phase, but they are generally extendable to multiple phases. Dispersed flow, the focus of this book, indicates that particle motion is generally dominated by the interaction with the surrounding flow as opposed to the presence of nearly particles. At low particle concentrations, there is one-way coupling whereby the dispersed-phase motion is affected by the continuous-phase but not vice-versa. The dispersed-phase motion is typically controlled by particle inertia, drag and gravitational forces. At higher concentrations, the dispersed-phase can also affect the continuous-phase through by changing the momentum or effective density and this is referred to as two-way coupling. At even higher concentrations, three-way coupling may occur whereby particle-particle fluid dynamic interactions can occur when particles in close proximity, but not touching. In particular, three-way coupling occurs when the flow around a particle will influence the flow around a neighboring particle, e.g. though wakes. Finally, four-way coupling is refers to particle-particle collisions when actual contact is made. As noted in 1.1, dispersed-flow conditions (the focus of this text) arise when particle motion is dominated by the continuous-phase flow interactions. On the other hand, dense flow occurs when the particle motion is dominated by particle-particle interactions, i.e. three-way or four-way coupling effects are stronger than one-way coupling effects. Details and criteria of these different coupling regimes will be discussed in Chapter 2. In the remainder of this chapter, we will consider a single isolated particle whose influence on the surrounding flow will be negligible (consistent with one-way coupling). The continuous-phase surrounding the particle will always be a fluid (gas or liquid) and it will be associated with the subscript "f". It will generally be assumed to behave as a continuum such that a dynamic fluid viscosity (f) and density (f) of this phase can be defined. The subscript "p" will be used when referring to particle properties, e.g., the volume-averaged particle density is p. For a fluid particle (such as a gas bubble in liquid, a liquid drop in a gas, or a liquid drop in another immiscible liquid), the viscosity of the particle (p) becomes finite and may play a role in the particle dynamics. However, if the particle is composed of a solid substance (such as a coal particle), the viscosity of the particle is not defined or can be


considered infinite. A fluid associated with either phase can also be assumed to have either gas or liquid properties (denoted by subscripts g and l). A volume will generally be denoted as , and the overall domain of the continuous-phase will be denoted D while the volume of an individual particle volume will be denoted as p. The length-scale d will refer herein to the particle volumetric diameter, i.e., the diameter of a sphere which has the same net volume as that of the particle (Fig. 1.20b):

6 d 2rp p

1/ 3


The above expression similarly defines the particle radius (rp) which will be equal to the actual radius for a spherical particle but will be a volume-based radius for a non-spherical particle. The volume-averaged particle density is the ratio of the particle mass (mp) to particle volume: p
mp p 1.2

In this text, we will assume the particle density is uniform so that the center of mass is equal to the centroid for spherical particles. The particle density ratio can be defined as * p p f 1.3

We may refer to very high density particles as having * p 1 and very low density particles as having * p 1. Examples of very high density particles include solid particles and droplets in a gas flow (which fall in the direction of gravity). Examples of very low density particles include bubbles and hollow particles in a liquid (which rise opposite the direction of gravity). Examples of particles with density ratios of the order of unity include drops in immiscible liquids. Herein, immiscible indicates that there is no molecular mixing at the interface so that the drop properties (mass, volume, etc.) remain distinct from those of the surrounding liquid. For fluid particles, the ratio of the viscosities can be important and is defined as * p p f 1.4a

The viscosity ratio will tend to exhibit the same tendencies as that for the density ratio, i.e., * liquid drops in a gas correspond to * p 1 and p 1 while gas bubbles in a liquid correspond to
* * Similarly, drops in an immiscible liquid correspond to * * p 1 and p 1. p ~1 and p ~1.

Fluid particles also contain an interface whose curvature causes an interface stress due to surface tension (). For quiescent conditions and convex curvature, this stress is balanced by the jump from the surrounding fluid pressure (p) to particle pressure (pp):

p p = p + (1/ rsurf ,1 + 1/ rsurf ,2 )



In this expression, rsurf,1 and rsurf,1 are the principal radii of the interface curvature in the two orthogonal directions. For a sphere, this result yields:

p p = p + 2 / rp


The density, viscosity, surface tension and other properties for two liquids and two gasses are given in Table A.1. A force or a velocity will be represented as a vector quantity when given in bold-face (e.g., F) and as a tensor quantity when given with an index (e.g., Fi). The scalar magnitude of a vector will be represented by a regular type-face (e.g., F equals |F|). To determine the fluid dynamics forces, the velocities of the multiphase flow with respect to a fixed reference frame should be defined. The particle velocity (v) is defined as the translational velocity of the particle center of mass (xp) as shown in Fig. 1.21a. The continuous-fluid velocity (u) is defined in all areas of the domain unoccupied by particles. To construct the relative particle velocity, the continuous-phase velocity away from the particle can be hypothetically projected to the particle center of mass, e.g. by interpolating the continuous-flow field velocity to xp as if the particle was not present. This theoretical velocity is denoted as u@p and is termed the unhindered velocity. The relative velocity of the particles (w) to the continuous-phase is then based on the unhindered velocity along a particle trajectory
w(t) v(t) - u@p(t)


It is important to note that u@p is un-physical (the continuous-phase does not exist inside the particle) and neglects the presence of the particle, i.e. u@p does not include the deviation of flow paths around the particle itself (though, it can include the fluid dynamic effects of neighboring particles associated with general two-way or three-way coupling). In addition to the above translational velocities, equivalent angular velocities may be defined for both phases. For the continuous-phase flow, vorticity () is conventionally defined as twice the fluid angular rotation: f 2 f 1.6 One may then define the unhindered vorticity, @p, as the fluid vorticity hypothetically projected to xp. Since the particle can rotate with an angular velocity (p) about the particle center of mass (Fig. 1.21a), the relative angular velocity can be defined (similar to Eq. 1.5) as p,rel(t) p(t) - f@p(t) = p - f@p 1.7 Often, the subscript f will be implied for , while the subscript p will be implied for . Forces, velocities, and particle positions can generally vary in space (x) and time (t). Spatial description will be generally represented in Cartesian coordinates for convenience. In particular, the three Cartesian components of the particle velocity and position are based on the unit vectors in the x, y, and z directions as 1.8a v(t) vx(t) ix + vy(t) iy + vz(t) iz
xp(t) xp(t) ix + yp(t) iy + zp(t) iz



However, sometimes the velocities, positions and forces will be given in spherical coordinates fixed to the particle to examine flow inside or in the immediate vicinity of a particle. Temporal variation can be considered in terms of various reference frames including a fixed Eulerian reference frame, a particle-Lagrangian reference frame, and a fluid-Lagrangian reference frame (Fig. 1.21b). For example, the Eulerian derivative of an arbitrary parameter q is based on temporal changes at a fixed point in the domain, x (though, the domain itself may be moving at some velocity) and thus is given as
q q(x, t + t) q(x, t) t t t 0


In comparison, the particle-path and fluid-path Lagrangian time derivatives are defined along a path specified by v and u respectively, i.e.
q dq q(x + vt , t + t) q(x, t) = + v q dt t t t 0

1.10a 1.10b

Dq q(x + ut , t + t) q(x, t) q = + u q Dt t t t 0

This assumes that the variable q can be considered as a continuum throughout the domain (this may not be possible with certain properties) so that these two derivatives can be related for a scalar field as:
dq Dq = + w q dt D t For a vector field q, the Lagrangian derivatives are related as: dq q Dq = + ( v ) q = + ( w ) q dt t Dt 1.12 1.11

For the declaration of several other variables and symbols used in the text, the reader is further referred to the nomenclature. However, new variables will generally be defined as they are introduced.



Description of Uncoupled Continuous-Phase Flow

In this section, we consider the condition where the continuous-flow is uncoupled, i.e., there is a negligible impact of the particles on the surrounding flow evolution given by u@p. This neglects disturbances in the immediate vicinity of the particle surface as well as macroscopic changes due to the collective action of the particles. This uncoupled condition is consistent with a small number of particles in the domain, such that the continuous-phase behaves as in a single-phase flow and whose conservation equations are discussed in Appendix A. This condition also serves as a baseline to which we can add the influence of particles on the surrounding flow, i.e., the coupled conservation equations to be discussed in Chapters 5, 7 & 8. In both cases, a continuum flow assumption is often made as discussed below.

Continuum Approximation

To establish the conventional derivative-based conservation equations for the continuousphase which will be introduced in the next sub-section, it is important that the surrounding fluid can be considered as a continuum. The continuum approximation is reasonable if the relevant length-scales for the flow-field (e.g. D) are much larger than the average molecular spacing (lm-m). In this case, a large number of molecular interactions occur over the relevant length-scale, such that only the averaged effect is important and the individual discrete (random) interactions are not significant. The breakdown of continuum occurs quite differently in a liquid than in a gas. For liquids, the molecules have an inter-particle spacing which is on the order of the molecular length-scale (less than a nanometer), whereas gasses typically have an inter-particle spacing which may be on the order of a thousand times the molecular diameter. Therefore, the continuum approximation is more likely to be violated with gas flows than with liquid flows for a given particle size. To determine the continuum criteria for a gas, the average mean-free path must be obtained. This is defined as the average distance traveled by gas molecules between collisions with each other. The mean free-path is inversely proportional to the product of the number of molecules per unit volume and their effective cross-sectional area (which is based on the molecular diameter). However, the molecular diameter is not a directly measurable quantity, so a more practical means of assessing the mean-free path is to relate it to the average rms (root-meansquare) speed of the molecules (Vm-m) and the measured gas viscosity and density:
l m-m = 2 g g Vm-m


For a perfect gas, the average speed of the molecules is related to the gas temperature and gas constant (defined in Appendix A.1), such that
Vm-m =

R g Tg


and is therefore roughly one-third greater than the gas speed of sound (ag). Based on this, the molecular mean-free path for air at STP (standard temperature and pressure conditions of 20 oC


and 1 atm.) is about 0.07 m and thus the continuum assumption is generally reasonable except at very small scales. Since gasses in general have similar dynamic viscosities which are only a moderate function of absolute temperature (approximately square root dependence), substantial increases in the mean-free path primarily occur at low gas densities such as those found in vacuum conditions (where the path increase is inversely proportional to the pressure decrease) or very high altitudes (where lm-m in the atmosphere is about one meter at an altitude of 130 km). To characterize when effects of molecular collisions can be neglected and a continuum can be assumed, it is helpful to parameterize the competing effects. This is accomplished with the Knudsen number (Kn), which is defined as the ratio of the mean-free path to the relevant flow length-scale. In multiphase flow, one can define both a macroscopic (continuous-phase) Knudsen number and a particle Knudsen number as follows
l m-m D l Kn p m-m d Kn D

1.15a 1.15b

In general, very small values of the Knudsen number (e.g., less than 10-2) are consistent with the continuum approximation and a no-slip boundary condition on a solid surface. This is because a very large number of molecular collisions occur over the relevant length-scale, such that the effects of their discrete interactions are negligible. Boundary conditions for small but finite Knudsen number (about 10-1) as illustrated in Fig. 1.22 can be incorporated by assuming a finite slip velocity at solid surfaces. For Knudsen numbers of order unity, the flow can not be considered as a continuum with respect to the relevant length-scales and non-continuum aspects should be incorporated directly. As the Knudsen number becomes much greater than unity, a free-molecular flow regime is established, which is dominated by individual molecular collisions. Based on this definition, continuum approximation criteria for the macroscopic fluid dynamics and the local fluid dynamics over the particle can be defined as KnD 0 Knp 0
macroscopic continuous-phase acts as a continuum flow around the particle acts as a continuum

1.16a 1.16b

For air at STP, the Knp criterion suggests that particles which are on the order of 5-10 microns or larger can be reasonably considered to obey the continuum approximation. While this assumption will generally be used in this text, the non-continuum effects for micron-sized or sub-micron particles will be discussed in Chapters 6 and 8. The assumption of a continuum for the macroscopic continuous-phase flow will be used throughout the text.

General transport equations

Assuming a continuum fluid with respect to the uncoupled flow, the transport equations can be written by assuming no mass, momentum, or energy transfer between the fluid system


and the particles or between the fluid system and the surroundings. The general (threedimensional, unsteady, compressible, viscous) conservation equations of mass, momentum, and energy are given by Eqs. A.1, A.2 and A.4. If the flow is compressible, these partial differential equations (PDEs) are supported by an equation of state to relate pressure changes to density changes discussed in Section () A.1. The Mach number is often used to assess the flow compressibility, whereby higher Mach numbers indicate larger fluid dynamic changes in pressure and density. Various inviscid formulations are discussed in A.2, while viscous formulations are discussed in A.3 and A.4 in terms of the three primary viscous regimes: laminar, transitional, and turbulent flow.



Equation of Motion for an Isolated Spherical Particle

This section will discuss general and specific equations of motion for an isolated, spherical particle subjected to drag and gravitational forces and simple flows. The resulting basic equations introduced in this section will be useful to illustrate the most common aspects of multiphase fluid dynamics discussed in Chapter 2-4. 1.5.1. General Particle Equation of Motion

The overall particle translational equation of motion simply specifies that the rate of change of the particles linear momentum is equal to the net sum of the forces acting on the particle. This yields a general Lagrangian equation of momentum given by


dv = Fbody + Fsurf + Fcoll dt


The right hand side (RHS) includes the forces associated with these temporal changes. Body forces (Fbody) are those proportional to the particle mass, surface forces (Fsurf) are those proportional to the particle surface area and related to the surrounding fluid stress, and collision forces (Fcoll) includes the effects of other particles or walls which may come in contact with the particle. These forces control the particle velocity, which is needed for the companion Lagrangian equation of particle position d xp dt


Integration of this equation yields the particle trajectory as a function of time. In the following, the RHS of Eq. 1.17 will be constructed in terms of particle and continuous-phase characteristics under specific assumptions. In this text, the body forces are assumed to be represented solely by the gravitational force (Fg) which acts in the direction of the gravity acceleration vector (g), such that Fbody Fg = mp g 1.19 This assumes that other body forces (such as electromagnetic forces) are negligible. In this chapter, collision forces will be neglected since we have assumed an isolated particle. The surface forces are also idealized based on several key assumptions. First, the surface force is assumed to be the linear sum of various fluid dynamic forces as follows Fsurf = FD + FL + F + FH + FS + FBr + FT 1.20

These individual components include forces due to: drag (FD) which resists the relative velocity, lift (FL) which arise due to particle spin or fluid shear, virtual-mass (F) which is related to the surrounding fluid that accelerates with the particle, history (FH) which takes into account unsteady stress over the particle,

fluid-stress (FS) which stems from the fluid dynamic stresses in the absence of the particle, Brownian motion (FBr) random motion from discrete molecular interactions, and thermophoresis (FT) force due to molecular interactions along a temperature gradient. Second, the fluid surrounding the particle is assumed to: have constant density so that compressibility effects are negligible, to have weak temperature gradients ( Tg  Tg / d ) so

that thermophoresis forces are negligible, and to be in continuum (Knp0.) so that Brownian motion is negligible. Third, we assume that the surrounding flow-field velocity gradients and particle spin are small so that particle lift can be neglected (FL=0). Fourth, the particle is assumed to be spherical with constant diameter and density. Using the above assumptions, a simplified surface force expression can be obtained as discussed in the following two sections. Note that Chapter 6 will consider more general conditions so that these assumptions can be relaxed. In particular, it will address issues associated with various shapes and particle surface conditions, flow compressibility, non-continuum effects, lift forces, particle-particle interaction, particle-wall interactions, etc. Note also that we have not included two effects associated with continuous-phase flow rotation: centrifugal and Coriolis forces. The centrifugal force is associated with particles being thrown out of rotating flow regions. It is not a surface force but instead represents the inertia associated with a particle which tends to translational motion and thus resists rotational motion. Similarly, the Coriolis force is not a surface force but is instead a relative acceleration to a continuous-phase system which is uniformly rotating. As such, both effects are already implicitly incorporated into Eq. 1.17 so long as a Cartesian coordinate system is employed, but would need to be added if one chooses a rotating coordinate system.
1.5.2. Quasi-Steady Drag Force

The quasi-steady drag force arises from pressure and viscous stresses applied to the particle surface. The resulting force resists the relative velocity so that it is defined to act in the direction opposite of the particle relative velocity (w). In this section, In this section, the force will be considered for steady continuum conditions. Furthermore, the particle and the far-field velocities are both assumed to be steady, i.e. w is constant, in order to obtain the quasi-steady drag force. Furthermore, the continuous-phase flow in the absence of the particles influence is assumed to be spatially uniform, i.e. .u@p=0, such that w=v-u. The dependence of the drag force on the magnitude of the relative velocity is primarily dictated by the particle Reynolds number (Rep), defined as Rep f w d f 1.21

The particle Reynolds number is the non-dimensional ratio of fluid inertial forces to viscous forces with respect to the fluid dynamics in the vicinity of the particle. The Reynolds number can be used to categorize the flowfield regime seen by the particle (Fig. 1.23). The fluid


dynamics around the particle, in turn, determine the quantitative relationship between the relative velocity and the drag force exerted on the particle.

Creeping Flow: Stokes Solution

When the inertial effects are negligible, the viscous and pressure fields dominate the flowfield solution resulting in a fully attached laminar flow over the particle (e.g., Fig. 1.23a). The dominance of viscosity also means that the shear stresses are felt far from the particle surface. As a result, this low Reynolds number condition is often called creeping flow. For this condition, Stokes (1851) derived the drag force in the limit of negligible convection terms and by assuming incompressible continuum flow of constant density and viscosity around the particle. Simplifying the incompressible Navier-Stokes momentum equation based on these assumptions and additionally neglecting gravitational and unsteady terms, yields a balance between the pressure gradient and the viscous stresses (Eq. A.56a). Considering this equation with respect to the flow around a particle yields the Stokes equation: p = f 2u
for Rep0


The neglect of gravity eliminates the hydrostatic pressure gradient associated with buoyancy, but this effect will be re-introduced in the next section. For a particle, application of this PDE and the proper surface boundary conditions can describe the local velocity field of the continuous-phase fluid around the particle. As such, the resulting velocity field is affected by the particle displacement and thus is different from the unhindered velocity (i.e. uu@p). The PDE of Eq. 1.22 can be expressed in terms of axisymmetric flow in spherical coordinates (see A.3) with a control volume is consistent with the particle center as shown in Fig. 1.24, and moves at the speed of the particle. In this case, the radial and tangential velocity components (ur and u) are defined throughout the surrounding flow (rrp) and the swirl velocity is assumed to be negligible (u=0). From this, Eq. 1.22 can be decomposed into radial and tangential momentum equations as discussed in A.3 1 p 2 u r 2 u r 2u r 1 2 u r cot u r 2 u 2u cot = 2 + 2 + 2 + 2 1.23a f r r r r r r 2 r r 2 r2 u 1 p 2 u 2 u 1 2 u cot u 2 u r 1.23b = 2 + 2 2 + 2 + 2 + r 2 rf r r r r sin r 2 r By introducing the Stokes spherical stream function which satisfies continuity, the radial and tangential velocity components of the continuous-phase fluid can be expressed as ur = 1 r sin

1.24a 1.24b

u = -

1 r sin r

Substituting these velocity components into Eq. 1.23 yields the stream-function based momentum equations


p = 2 f 2 ) ( r r sin

1.25a 1.25b

1 p f = ( 2 ) r r sin r

For convenience, these expressions use a specialized spherical differential operator 2 defined as
2 q sin 1 q q = 2 + 2 r r sin


This special operator is similar to, but different than, the true Laplacian operator of Eq. A.51d. The pressure term can be eliminated by cross-differentiating the two equations of Eq. 1.25 with respect to and r, yielding a single fourth-order differential equation for the stream function
2 (2 ) = 0


The boundary conditions are now needed to allow solution of this PDE. The boundary conditions for a solid non-rotating sphere are given by a no-slip condition on the particle surface. In a coordinate system which is moving with the particle centroid, this boundary condition corresponds to ur=u=0 at r=rp. Based on Eq. 1.24, these boundary conditions can be written in terms of the stream function
= =0 r
as r = rp


Far from the surface, the boundary conditions correspond to uniform flow in the horizontal direction (x-direction) with a velocity of -w (which is based on u@p). Thus, the radial and tangential components can be given as as r ur = w cos 1.29a 1.29b as r u = -w sin Note that this coordinate system translates the continuous-phase velocities to ones relative to the velocity of the particle centroid. Based on Eq. 1.24, this yields a boundary condition on the stream function given by
= w r2 sin2 + const.
as r


The constant in the above expression is arbitrary and so can be set to zero. Introducing a separation of variables =f(r)sin2(), which satisfies the governing equation (Eq. 1.27) and the above boundary conditions (Eq. 1.30), the stream function can be obtained as (White, 1991):
3 wr 2 3rp rp 2 1 sin + 3 2 2r 2r


The first term of this expression corresponds to the far-field condition (Eq. 1.30), the second term is called the Stokeslet (a viscous correction), and the third term is called the doublet (which is an inviscid correction). The viscous correction decays much more slowly with


increasing distance from the particle surface than the inviscid correction. From Eq. 1.24, the radial and tangential velocities can be obtained (with free-stream, doublet and Stokeslet contributions) as 3rp rp3 u r = ( w cos ) 1 + 1.32a 2r 2r 3 3rp rp3 u = ( w sin ) 1 3 1.32b 4r 4r The resulting streamlines are symmetric both front-to-back and top-to-bottom, which is consistent with the flow seen in Fig 1.23a and the condition of reversibility associated with Eq. A.56a. This indicates that the particle influence is felt as far upstream as it is felt downstream.
To obtain the drag coefficient, both the pressure and shear stresses must be obtained on the particle surface and integrated over its area. To obtain the pressure field, the pressure changes can be expressed in terms of radial and tangential pressure gradients dp = p p dr + d r 1.33

The derivatives of the pressure on the RHS of this equation are given by Eq. 1.25 in terms of the spherical differential operator which, in turn, can be obtained from Eqs. 1.26 and 1.31 as 2 = 3wrp sin 2 2 r 1.34

The pressure can then be obtained by substituting this expression into Eqs. 1.25 and 1.33 and integrating to yield p = p 3wrpf 2r 2 cos 1.35

This can be evaluated on the sphere surface by setting r=rp. The result indicates that the pressure is symmetric top-to-bottom but is anti-symmetric in the streamwise direction. Since w<0 (Fig. 1.24), the pressure is higher on the front of the particle (=0) than on the back (=), giving rise to a drag in the positive x-direction. Once the velocity and pressure distributions are known, the drag force can be obtained by integrating the stresses on the particle surface. Integrating the streamwise component of the pressure (pcos) with a differential area (R2sind) over the upper surface (=0 to ) and assuming symmetry for the lower surface (by introducing a factor of 2), the resulting force is

Fpressure = 2 rp2

p cos sin d = 2rp f w


This is known as the pressure drag or sometimes called form drag. The shear stress is given by Eq. A.51b as


u 1 u K r f r + r r r r


Given the fluid velocity from Eq. 1.32, this yields the local tangential shear stress on the particle surface. Using surface integration (in a manner similar to Eq. 1.36) of this stress yields the friction drag:

Fshear = -4 rp f w


This is equal to twice the pressure drag, affirming the importance of the viscous stresses. The total drag force is simply the sum of the form and friction drags:

FD = - 6 rp f w


This is often referred to as the Stokes drag, owing to his derivation and is in good agreement with the experimental results for small but finite Reynolds numbers (e.g. within 2% accuracy for Rep<0.1). The Stokes solution can be expressed in terms of the drag coefficient (CD) which is defined in terms of the total drag force normalized by the Bernoulli dynamic pressure (Eq. A.36) and the cross-sectional area: F 24 for Rep0 CD 2 D = 1.39b 1 2 Re p d 2 f w 4


Therefore, the creeping flow conditions yields a drag force linearly proportional to relative velocity and a drag coefficient inversely proportional to Reynolds number. For linearized convection, the theoretical Oseen approximation can be employed up to Rep of about 1 (6.2.1), but unfortunately, there is no closed-form analytical solution for the viscous flowfield once the particle Reynolds number is greater than unity.

High Reynolds Number Flow: Newton-Based Drag

At the other extreme of very high Reynolds numbers (Rep1), some theoretical analysis is also possible if we assume the viscous effects are confined to a small regions near the surface so that much of the external flow can be analyzed with inviscid potential flow theory. For incompressible flow, the equations of motion for mass and momentum in spherical coordinates are given by Eqs. A.50 and A.52. Considering only the flow outside of the boundary layer, an inviscid boundary condition can be approximately applied at the particle surface, i.e. ur=0, while the far-field boundary conditions are the same as Eq. 1.29. The stream function solution turns out to be a linear combination of a uniform stream and a point-doublet (White, 1991), which can be given in terms of radial and tangential velocity components by Eq. 1.24:


r 2 rp 2 2 wr sin =1 1.40a p 2 rp2 r rp3 u r = w cos 1 3 1.40b r rp3 u = w sin 1 + 3 1.40c 2r For the velocities, the first term in the parentheses is the far-field solution while the second term is the inviscid doublet. In contrast to the viscous solution, the tangential surface velocity is finite and given by u = 3 2 w sin . Using this velocity and the Bernoulli relation (Eq. A.36), the surface pressure is given as: 9 p = p + 1 f w 2 1 sin 2 1.41 2 4 This distribution is shown in Fig. 1.25 by defining a pressure coefficient as the difference between the static and the far-field pressures normalized by the dynamic pressure (fw2): p p Cp 1.42 1 w2 2 f The inviscid pressure is thus symmetric about the fore and aft surface, with the stagnation pressure recovered at both =0 and (where Cp=1). Integrating this symmetric inviscid pressure distribution over the particle surface leads to the solution FD=0. This is known as dAlemberts paradox, i.e., that the net force acting on a body moving though an inviscid fluid is zero. This is, of course, inconsistent with the to drag observed on bodies at high Reynolds numbers because viscous effects, while generally negligible in the majority of the flowfield for high Rep values, become important in the thin region near the particle surface defined as the boundary layer. The boundary layer is important since it incorporates the finite shear stress that acts on the body owing to the no-slip condition and since it determines if and where the flow will separate from the particle surface. At high particle Reynolds numbers, the boundary layer always separates and yields a turbulent wake. Once the flow separates, the concept of a thin boundary layer is lost and the outer potential flow solution is also unreasonable. Prior to separation, the attached boundary layer can generally have two states: laminar or turbulent.

The condition where the boundary layer is thin and laminar before separation is termed the sub-critical range, and corresponds approximately to 2000<Rep<300,000. In this range, laminar separation occurs just before the mid-plane (~80o) and produces a fully turbulent wake downstream (Fig. 1.23c). This creates a pressure distribution similar to the inviscid solution in the upstream attached region (given by Eq. 1.41) but yields a nearly constant pressure zone in the separated region, which is lower than that of the front stagnation region (Fig. 1.25). This results in a net pressure drag opposing the relative velocity. The portion of drag associated with friction is related to the shear stress acting on the surface. Before the separation point, the velocity field just above the boundary layer closely follows the inviscid solution (Eq. 1.40). The no-slip condition within the attached boundary layer results in a shear stress component of drag opposing the relative velocity. After the separation point, the wall shear stress is nearly negligible due to lower velocities and lower gradients in the detached


flow zone and drag in this region primarily stems from the surface pressure. Thus, friction and form drag are both important in this Reynolds number range and are related to viscous effects close to the particle surface. Since both the pressure and shear stresses are related to the local dynamic pressure and since the separation point is approximately constant for 2,000<Rep<300,000, it can be expected that the resulting drag in this range will be simply proportional to the free-stream dynamic pressure and the particle frontal area. This is consistent with measurements that have shown that drag coefficient is approximately constant in this Rep range (Fig. 1.26) where its value has been called the critical drag coefficient with a value approximately given by: for 2,000<Rep<300,000 CD,crit 0.40 0.45 1.43 The occurrence of a nearly constant drag coefficient for large particles was first recognized by Sir Isaac Newton when he analyzed drag resistance for projectiles in 1710. As such, the subcritical Rep range is often called Newtons regime. The maximum Reynolds number for the above condition is termed the critical Reynolds number, Rep,crit (about 300,000), which is generally defined when CD drops below 0.3. At Rep>Rep,crit, the attached boundary layer becomes turbulent and causes the flow separation point to be shifted to a further downstream position (~120o). The reduced separation region associated with a wire tripping is shown in Fig. 1.23d, which is significantly smaller than that associated with the sub-critical flowfield shown in Fig. 1.23c. This allows the pressure to follow the potential flow solution over a larger portion of the surface and thus substantially increases the pressure in the aft region as compared to the sub-critical case (Fig. 1.25). As a result, the drag coefficient in this super-critical regime drops dramatically, which creates a phenomenon termed the drag crisis (Fig. 1.26). This reduction in drag was first observed by G. Eiffel, who was both an avid fluid dynamicist and the famous French architect. In the super-critical regime, the drag becomes is again sensitive to the Reynolds number since the turbulent boundary layer separation point varies with Rep. Because most multiphase flows are well characterized by Rep<Rep,crit, the super-critical region is usually not applicable unless the particles are very large and moving at very high speeds. For example, a 75 mm (3) diameter baseball moving at 36 m/s (80 mph) in air at STP has an Rep of about 1.8x105 consistent with sub-critical conditions for a smooth sphere. However, the stitching on the surface prematurely trips the attached boundary layer to become turbulent so that the drag crisis occurs with a reduced resistance and thus faster speeds. A similar drag reduction occurs with the addition of dimples on a golf ball, thus allowing longer trajectories.

Drag at Intermediate Reynolds numbers

For intermediate particle Reynolds numbers (about 0.1<Rep<2,000), the flowfield in the rear of the particle undergoes several transitions: from an attached laminar wake (Fig. 1.23a), to a separated but still laminar wake (Fig. 1.23b), to an unsteady transitional wake. As discussed in the previous two sub-sections, creeping flow conditions are recovered below this


range of Reynolds numbers while turbulent wakes (e.g. sub-critical and super-critical flow) are recovered above this Reynolds number range. Thus, the intermediate range is characterized by wakes which are range from steady laminar conditions to transitional conditions depending on the Reynolds number. For intermediate flow with Rep<20, there is no flow separation but viscous effects are strong so that the Cp distribution tends lie in between the inviscid solution and the creeping flow solution of 6 cos / Re p , given by Eq. 1.35. As the Reynolds number increases further, about 25<Rep<2,000, a separation point occurs at the rear stagnation point (=180o) and moves upstream as the Reynolds number increases, reaching a value of about 108o Rep=100. The Cp distribution at Rep=100 (Fig. 1.25) is similar to the potential flow distribution region in the fore portion of the sphere but does not reach the inviscid minimum due to viscous effects. After the flow separation, the pressure only changes weakly indicating a lack of recovery which contributes to the overall drag. Based on measurements of the drag coefficient for intermediate Reynolds numbers (Fig. 1.26), the drag coefficient transitions from the creeping flow value given by Eq. 1.39b to the nearly constant critical drag coefficient of Eq. 1.43. While there is no general analytical solution for intermediate Reynolds numbers, empirical drag coefficient expressions can be used to describe the overall variation, e.g. the expression by White (1991): CD = 24 6 + + 0.4 Re p 1 + Re p
for Rep<2x105


The first term on the RHS is the Stokes limit, the last term on the RHS is the Newton limit, and the intermediate term is simply an empirical correction. The agreement with experimental data for this expression is good for a large range of conditions as shown in Fig. 1.26. Also shown in this figure is the Putnam (1961) fit, which is another empirical expression which is suitable for low and intermediate values of Reynolds numbers:
CD = 24 4 + 1/3 Re p Re p

for Rep<1000


A similar fit (which would be difficult to distinguish from the Putnam fit in Fig. 1.26) for small to moderate Reynolds numbers is the Schiller-Naumann (1933) expression, which is CD = 24 1 + 0.15Re0.687 p Rep

for Rep < 800


This is the most commonly used drag coefficient in multiphase flows since it is quite accurate and since many particles are constrained to Rep values in this range. A more detailed comparison of these drag coefficient expressions is given in 6.2. Unlike the analytical creeping flow drag, such empirical drag coefficients are non-unique and only valid for conditions for which they have been calibrated. There are, unfortunately, many other empirical expressions which can arise in computational multiphase fluid dynamics. One may quantify the departure of the drag from the Stokes solution by normalizing the drag force by the creeping flow solution for a sphere in an incompressible uniform continuum,


FD CD 3d f w (24 / Re p )


This ratio is called the Stokes correction factor and allows the drag force to be generally written as

FD = - 3 d f f w


If only Reynolds number effects are considered (i.e. a solid sphere in incompressible flow), this correction factor can be denoted as fRe, which can then be obtained from Whites fit as: FD (Re p ) Re p /4 Re f Re = 1+ + p 1.49 FD (Re p 0) 60 1 + Re p To account for conditions that may include continuous-fluid compressibility, non-spherical particle shapes, etc., a set of detailed models for CD and f will be discussed in 6.2. However, the departures from the above expressions in many cases are often not very large, i.e., of order unity. For example, White (1991) notes that the drag coefficients for non-spherical particles tend to be in the ball park of those for a sphere with the same volume. This is also generally true with respect to differences between solid and fluid particles, where the latter is next discussed.

Drag for Spherical Fluid Particles

Compared to the solid sphere discussed above, a fluid sphere (drop or bubble) that allows internal recirculation will generally have a reduced drag. This condition for Rep0 permits an analytical solution if the interface between the phases can be considered clean (i.e. fully mobile) and that the surface tension is uniform over the particle (reasonable if the effect of contaminants is negligible). The free-stream boundary conditions of Eq. 1.25 and the stream function PDE of Eq. 1.27 still apply (for both and p) as was the case for a rigid sphere in creeping flow. However, the mobile interface leads to modified boundary conditions and a modified stream function solution outside of the particle () as well as introduces a separate stream function defined inside the fluid particle (p). In particular, the surface conditions at r=rp for steady flow are related to the external and internal stream functions:
= p p = r r 1 1 p f 2 = p 2 r r r r r r p 2 f 1 1 p 2 + 2 = pp p 2 r r sin r r sin rp 1.50a 1.50b 1.50c 1.50d

The first expression (Eq. 1.50a) simply indicates that the surface of the sphere is a streamline (no mass transfer across the particle surface). The second and third expressions (Eqs. 1.50b and 1.50c) indicates that the tangential velocity and the tangential stress are continuous across


the interface, where the latter is obtained by setting Kr=Kr,p via Eq. 1.37 and the stream function definition of Eq. 1.24b. This is sometimes called the stress-free conditions since there is assumed to be no contamination which interferes with the stress balance given by the individual viscosities. The fourth equation gives the normal stress balance across the interface arising from surface tension and pressure stresses (Eq. 1.4c) as well as viscous stresses (A.51a) which are expressed in terms of stream function using Eq. 1.24a. The analytical solution for the continuous-phase and particle-phase streamlines for the above boundary conditions is generally called the Hadamard-Rybczynski solution as it was independently derived by both in 1911. Details of the stream function derivation are given by Batchelor et al. (1967) and with the solution:
* 2 3rp3* wr 2 3rp (2 + 3 p ) p = + sin 1 * 3 * 2 2r(3 + 3 p ) 2r (3 + 3 p )


This expression reverts to the solid sphere expression of Eq. 1.31 as * p . The internal motion associated with p is called the Hills vortex and the theoretical result and observed streamlines are shown in Fig. 1.27 for a drop falling in a different liquid. The pressure distribution external to the particle of fixed volume can be obtained using the relationships of Eqs. 1.25 and 1.51 to yield:
p = p 3wrp f (2 + 3* p) 2r 2 (3 + 3* p) cos


The pressure within the fluid particle (pp) is associated with a discontinuity at the interface due to the surface tension and shear stress by Eq. 1.50d so that the net change from the far-field is: p p = p + 5wrp p 2 2 cos rp 2r (1 + * p) 1.53

In the limit of no relative velocity, the internal pressure is simply specified by the surface tension jump, i.e. Eq. 1.4c. Thus, a spherical shape is reasonable for creeping flow conditions consistent with the analysis of Saito (1913). The pressure drag and shear stress (both the tangential and the deviatoric normal components) can be integrated over the particle surface to determine the drag in steady flow conditions yielding

FD = - 6 rp f w

2 + 3* p 3 + 3* p


If we define the Stokes correction for internal circulation (f*) as the ratio of the drag force acting on fully mobile fluid sphere to the drag force on a solid sphere at Rep0, we have
f * FD (* p , Re p 0) FD (* p , Re p 0) = 2 + 3* p 3 + 3* p



Therefore, the drag of a fluid sphere at this condition with * p 0 (e.g. a gas bubble in a liquid) is 2/3 of that for a solid sphere. However, contaminants in the system often cause the interface to become immobile (or partially immobile) such that this drag reduction is not realized (or only partially realized). Furthermore, spherical liquid droplets falling in a gas will tend to be characterized by * p so that the correction of Eq. 1.55 will be negligible. The impact of contamination, deformation and impact of finite density and viscosity ratios on drag at finite Reynolds numbers will be discussed in 6.2.3 and 6.2.4. Furthermore, low density ratio fluid particles tend to be especially influenced by the acceleration-based forces in addition to the velocity-based drag force.


Surface Forces due to Accelerations

Assuming again a solid spherical particle in a continuum incompressible flow, the surface forces associated with accelerations of the far-field continuous-phase velocity and the particle velocity are considered in this section. This will give rise to three types of acceleration forces: the fluid stress force, the added mass force, and the history force. The first acceleration force is the fluid-stress force (FS) which arises due to the continuousphase pressure and shear stresses that will exist in the absence of the particle. For example, consider the deformable fluid element in the steady diverging nozzle flow of Fig. 1.28. As the element moves downstream, it will decelerate according to the area increase. This deceleration along the fluid element Lagrangian path (xf) is simply controlled by the fluid stresses (pressure and shear) integrated over the fluid element surface. The surface-area integration (over Af) can be converted to a volume integration (over f) by applying the divergence theorem. For example, the pressure stress which acts opposite the outward unit normal vector n of the fluid element surface yields a fluid stress force based on pressure divergence within the volume: ( pn ) dA f = ( p ) df = f p 1.56 The RHS expression assumes that the pressure gradients are approximately uniform in the vicinity f (reasonable for small length-scales compared to changes in the fluid properties). A similar conversion can be applied to the viscous stresses if the differential volume is assumed small, so that the net fluid stress over the volume of a fluid element is given by: FS,f = f ( p + K ij ) 1.57 Therefore, the unhindered stresses (which neglect the influence of the particle) give rise to a fluid-stress force for either a fluid element or a solid spherical particle which inhabits the same space (Maxey and Riley, 1993). Within the fluid volume, one may also relate the pressure and viscous stresses by the single-phase momentum equation (Eq. A.5b): Du g ( p + f K ) f = f f 1.58 Dt These unhindered fluid stress forces (which neglect the influence of the particle within f) will act on a fluid element, e.g. to decelerate it along its path as shown in Fig. 1.28.


Now consider a spherical particle which occupies this volume (p=f). It will also be subjected to these forces (6.5) so that the fluid stress forces on a particle can be similarly expressed in terms of the Lagrangian fluid acceleration and hydrostatic force (Maxey, 1993): D u@p - g FS = f p 1.59 Dt The fluid-stress force is thus related to the convective and gravitational accelerations of the continuous-phase field, but is independent of the volume shape (i.e. particle shape). If the continuous-phase fluid acceleration is negligible, then the fluid stress force reduces to the buoyancy force Fbuoy - f p g 1.60 In other words, the buoyancy force is equal to the hydro-static pressure gradient integrated over the particle surface and acts in the opposite direction of the gravitational (body) force. It can also be recognized from Eq. 1.59, that the fluid stress force is proportional to the displaced continuous-phase mass m f @ p f p 1.61 Since the fluid stress is proportional to the continuous-phase density, it will tend to be more important for particles with low density ratios, i.e. * p 1. The second acceleration force is the added-mass or virtual mass force (F) which arises when the particle acceleration is different from that of the continuous-phase field, i.e. when there is a non-zero relative acceleration (dw/dt ) . This relative acceleration gives rise to an additional force because a portion of the surrounding fluid (in the vicinity of the particle) is carried along with the particle with the same acceleration. This is due to the particle boundary conditions which prevents mass flux through the surface. The surrounding fluid mass which moves the particle is referred to as the added-mass or the virtual mass (m). It is schematically described in Fig. 1.29 as an integrated effect that is strongest near the particle surface and weakens further away. For a solid sphere, the added-mass can be related to the displaced fluid mass using by defining c as the virtual mass coefficient: m = c m f@p 1.62

Analytically, the force require to accelerate the added-mass associated with a sphere at the relative particle acceleration is defined as the added-mass force, i.e.
F = - c f p

dw dt


The added-mass coefficient, c, is typically taken as since this corresponds to both the value for both the Stokes regime and inviscid flow (Maxey, 1993), and moreover it has been found to be reasonable for a wide range of Reynolds numbers for spheres (see 6.4). The third acceleration force arises when the viscous stresses on the particle surface are unsteady and is called the history force (FH) since it will depend on the temporal development of this flow. This force can be significant if the accelerations are large and the viscous effects are strong but is zero if inviscid flow is considered (unlike the added mass force which is


present for both viscous and inviscid flow conditions). The history force arises due to the unsteady development of the viscous region due to the particle no-slip condition. It is generally treated separately from the quasi-steady drag force FD, but is, in fact, the unsteady portion of the drag force so that the two forces together are the total drag. To aid in the understanding of this force, Crowe et al. (1998) discuss the equivalent history force on a flat plate impulsively accelerated to a fixed velocity at some initial time (t=0) which was beforehand quiescent (u=0). The creeping flow analytical solution to the time development of the velocity filed is often referred to as Stokes first problem and is derived from the one-dimensional unsteady flow equation f u x 2u x p = + f t x x 2 1.64

The flowfield nearest the plate will be accelerated first leading to high initial shear rates, while further regions will take some time to accelerate according to viscous diffusion of the shear. The resulting boundary layer on the plate will grow with time. This growth leads to a constantly reducing surface velocity gradient such that the frictional force on the plate continually decays (even though the relative velocity of the plate to the far-field is steady). For a more general motion, the plate friction can be expressed as an integral of the relative 1/ 2 accelerations weighted by lapse , where lapse is the time since an acceleration occurred, so that their influence is continually reduced as time progresses. Based on these same principles, the derivation for the history force on a spherical particle with general acceleration was first formulated by Basset (1888) for creeping flow conditions (Rep0). Assuming the particle is initially at a quasi-steady state and the continuous-flow is spatially uniform away from the particle, the theoretical formulation is given by t dw 1 w 3 FH = d 2 f f d + t =0 1.65 0 2 t d t Thus the influence of previous accelerations decays as t-1/2 consistent with that of the flat plate. Unfortunately there is no analytical solution for the history forces at higher Reynolds numbers, but fortunately it tends to be weak (and thus negligible) at high Rep and has been reasonably modeled for intermediate Rep (6.6). Furthermore, the history force is generally insignificant if the particle accelerations are mild (compared to gravitational acceleration) or the particle density ratio is high ( * p 1).


Simplified Equations of Motion for an Isolated Particle

Translational Equations of Motion

Depending on the assumptions used, a variety of equations of motion can be put forth. Based on the above surface and body force expressions at creeping flow (Rep0), the particle


equation of motion for the translational acceleration of an isolated, rigid, non-rotating, spherical particle in a continuous-phase velocity field is: (p + c f)p

dv = -3 d f w + p (p-f) g + dt t dw d u @p D u @p 3 2 f p c + d f f 0 dt Dt 2 d

w t =0 1 d + t t


This is result is the classic Basset-Bousinesq-Oseen (BBO) equation. It is worth noting that this expression is non-empirical, i.e., it can be obtained directly from the incompressible continuum flow equations (Eq. A.48) applied to the particles exterior surface. If one wishes to employ a similar equation of motion at larger Reynolds numbers, the creeping drag force is generally replaced by the empirical drag force (Eq. 1.48). The added mass and fluid stress expressions are unchanged since they are also theoretically correct in the potential flow limit, and have been found to be generally reasonable for a wide range of Rep conditions. However, the history and lift forces are generally neglected for non-creeping conditions since they are somewhat controversial and often weaker for higher Reynolds numbers (although significant advances have been made for some conditions as discussed in 6.3 and 6.6). Thus, a corresponding non-creeping equation of motion can be given as (p + c f)p
d u@p D u @p dv = -3 d f fw + p(p-f)g + fp c + dt dt Dt


This equation of motion is used often by many multiphase researchers in computations for general Rep conditions with the above assumptions. For a more general equation of motions, Chapter 6 will discuss incorporation of non-spherical and compressible flow effects, lift due to vorticity and/or particle rotation, history force, thermophoresis, etc. If one explicitly assumes a spatially-uniform surrounding-fluid velocity (consistent with the drag force derivation and the neglect of the lift force), the temporal derivatives of the fluid velocity become equal based on Eq. 1.11 so that the above equation of motion becomes: (p + c f)p u@p dv = -3 d f fw + p(p-f)g + fp (1 + c ) t dt 1.68

The particle acceleration terms have been placed on the LHS (left-hand-side) which suggests that an effective mass of the particle can be defined as the sum of the particle mass and the virtual mass: meff mp + m = (p + f c) p 1.69 One can consider two limits of this effective mass expression. If the particle density is much greater than that of the surrounding fluid ( * p 1), as is typically the case of a liquid or solid particle in a gas, then the effective mass is approximately equal to mp. However, if the particle density is much less than that of the surrounding fluid ( * p 1), as is typically the case of a gas bubble in a liquid, then the effective mass is approximately equal to m.


As with the effective mass, the effects related to gravity were similarly grouped together in Eqs. 1.66-1.68. These terms can be combined to form an effective gravitational force which is proportional to the difference between the particle mass and the displaced mass as noted by the Archimedes principle: Fg,eff p (p- f) g = (mp - mf@p) g 1.70 * If the particle density is much greater than that of the surrounding fluid ( p 1), then the effective gravitational force can be approximated as the body force of the particle (Fg,effFg). However, if the particle density is much less than that of the surrounding fluid ( * p 1), then the effective gravitational force can be approximated as the buoyancy force acting on the particle (Fg,effFbuoy). With the above notations, the isolated spherical particle equation of motion in a spatiallyhomogeneous surrounding flow-field becomes meff u@p dv = FD + Fg,eff + p f (1 + c) t dt 1.71

The particle acceleration can then be written in terms of the density ratio ( * p ) and the acceleration parameter (R) defined below: u@p 18 fw dv = - 2 *f + (1 R ) g + R dt t d p + c

1.72a 1.72b

+ f c 1 + c = f * p + c p + f c

The acceleration parameter relates the importance of the fluid acceleration to the particle acceleration. Assuming an added-mass coefficient of , R approaches 3 for very buoyant particles ( * However, R p 1) indicating the fluid acceleration effect can be substantial. approaches 0 for very dense particles ( * p 1) indicating the fluid acceleration effect will be small. If the surrounding flowfield is steady (or R is negligible) and a particle is initially in equilibrium with the flow (v=u), then it will accelerate due to the effective gravitational forces. For very heavy particles ( * p 1) the initial particle acceleration will be g, which is the well known ballistic result. For very buoyant particles ( * p 1) the initial particle acceleration will be -2g, due to the added mass effect. Eventually such a particle will reach a relative equilibrium speed which is called the terminal velocity (wterm) and is independent of the initial conditions and is equal to the steady-state particle velocity in quiescent conditions. The terminal velocity can be obtained based on the balance between drag and effective gravitational forces of Eq. 1.72a and then can be expressed in terms of either the Stokes correction factor or the drag coefficient:


w term

( =

f ) gd 2

18f f term

1.73a 1.73b

w term =

4 p f g d 3 f CD,term

Note that the direction for wterm will be in the direction of g for * p >1 or in the direction of g for * p <1. Based on Eq. 1.73, the magnitude of the terminal velocity increases when either the particle diameter or density difference increases, and the degree of sensitivity depends on whether creeping, intermediate, or Newtonian flow conditions apply since f and CD can be a function of Reynolds number for different flow regimes. Let us consider some typical terminal velocities for a 500 m spherical particles at STP (Table A.1). If it is a solid particle of density of 1000 kg/m3 falling in air, it will have a terminal velocity of 1.9 m/s. However, if it is a non-deformable gas bubble with a no-slip surface rising in water, it will have a much lower terminal velocity of 0.055 m/s. In both cases, intermediate drag expressions are appropriate. On the other hand, a Stokes drag assumption requires a small particle Reynolds number. Example conditions for Rep,term=0.1 include a 40 micron coal particle falling in air or a 75 micron air bubble in water. The focus for the next two chapters will primarily be on the translational particle momentum. In particular, Chapter 2 will primarily examine particles which approximately move at terminal velocity conditions (where Eq. 1.73 is applicable). This will allow some understanding of particle characteristics in relatively simple conditions. Chapter 3 will then consider the dynamics of particles in the context of the surrounding fluid, e.g. turbulence (where Eq. 1.66 and 1.67 may be applicable). However, some brief comments on the heat and mass transfer of a particle are offered in the next section.


Heat and Mass Transfer for an Isolated Spherical Particle

While heat and mass transfer are not the focus, simplified equations for heat and mass transfer are introduced herein based on Sirignano (1999) to provide the reader with the basic relationships. For more detailed discussion on heat and mass transfer, the reader is referred to texts of Clift et al. (1978), Soo (1990) Crowe et al. (1998), Sirignano (1999), Brennen (2005) and Michaelides (2006). These equations presented herein are idealized for the simple case of a isolated spherical particle with uniform internal temperature and density. The continuousphase surrounding the particle is assumed to have uniform temperature and radiation and chemical reactions are neglected. Furthermore, the kinetic energy associated with the mass efflux velocity on the particle surface associated with mass transfer is assumed to be negligible, i.e. small compared to changes in the internal energy (temperature) of the particle (Crowe et al. 1998). As such, the relations are only reasonable for small non-combusting particles of uniform composition.


For heat transfer, the Lagrangian ODE related to the conservation of energy states that the increase in particle internal energy minus the energy required for any phase change (associated with mass transfer) is equal to the heat transferred from the continuous-phase to a single p : particle, Q

p 1.74  p =Q h phase m dt In this expression, cp,p is the specific heat (at constant pressure) of the particle and hphase is the enthalpy required for the particle matter to change to a phase consistent with the surrounding fluid phase (e.g. equal to the heat of vaporization for a droplet in a gas). The heat transfer rate is related to the thermal conductivity of the continuous phase (kf) and the temperature difference between the far-field fluid and the particle surface (Tf-Tp). Based on the above assumptions one can then define the non-dimensional heat transfer using these parameters and the particle diameter as the Nusselt number: p Q Nu 1.75 dk f ( Tf Tp ) mp
A more formal and general definition of the Nusselt number is given by Sirignano (1999) which allows internal temperature variations, whereas this form assumes that the particle internal temperature is uniform and equal to Tp. Further assuming the specific heat of the particle is constant, the temperature evolution given by Eq. 1.74 can be written in terms of the temperature difference and the change in particle diameter: dTp 6Nuk f 3h dd = Tf Tp ) + phase ( 1.76 2 dt p d c p,p dc p,p dt For creeping flow conditions with no mass transfer, the temperature field around the particle is governed only by thermal diffusion (no heat convection) which allows an exact solution, similar to the result found for Stokes drag (Clift et al. 1978). This limiting heat transfer solution corresponds to Nu=2.

d (c p,p Tp )

If the surrounding fluid temperature is approximately constant, then the ODE without mass transfer can be integrated analytically (Crowe et al. 1998). The result indicates that the relative temperature difference will decay exponentially based on a thermal response time (t): 12k f t Tp (t) Tf t = exp = exp 1.77a 2 dc Tpo Tf T p p,p
T p d 2c p,p 12k f


This indicates, that the particle mass can be considered constant for short-times (tt) and can be considered to be approximately in equilibrium with the surrounding fluid for long-times (tt). The thermal response time is an indication of the time it takes for a particle to reach equilibrium if there is an initial difference with the surrounding continuous-phase temperature. For the particle mass transfer, the change in particle diameter can be expressed in terms of  p (Sirignano, 1999): the mass transfer rate from the continuous-phase to a single particle m

1.78 dt d 2 The rate at which mass is transferred to the particle from the surrounding fluid will be governed by two effects: 1) how close the surrounding fluid is to saturation for the particle species and 2) how fast the particle mass species can move (diffuse) through the surrounding fluid. The first effect can be quantified by the Spalding transfer number, B, which is defined in terms of the difference between the particle species mass fraction just outside of the particle surface (assumed to be saturated) and the particle species mass fraction in the far-field. If the surrounding fluid is a gas, these mass fractions are proportional to the vapor pressure normalized by the local pressure, which is assumed to be constant. Since the particle is assumed to have uniform temperature, this is simply equal to the saturated vapor pressure based on the droplet temperature (Tp), which is denoted as pvapor,surf. The corresponding farfield vapor fraction is denoted as pvapor,f, and will be zero at a minimum (if there is negligible amounts of the particle species far from the particle) or will be the saturated vapor pressure based on the gas temperature (Tf) at a maximum. The resulting definition of the Spalding transfer number also takes into account the difference between the particle species vapor molecular weight (MWvapor) and the far-field gas molecular weight (MWf): ( p vapor,surf / p ) ( p vapor,f / p ) Bvapor = 1.79 ( MWf / MWvapor ) ( pvapor,surf / p ) As such, the Spalding number will be zero if the surrounding gas has the same vapor pressure as the saturated value found on the particle surface and there will be no mass transfer. If there is a vapor pressure difference, the Spalding number sign will be opposite to that of the mass transfer rate. For example, an ambient vapor pressure is higher than the particle surface vapor pressure (which can occur if Tp<Tf and both are saturated) there will be condensation and mass will move from the surrounding gas to the particle (positive mass transfer) which is consistent with a negative Spalding number. On the other hand, an ambient vapor pressure lower than the saturated conditions (which will occur if Tp>Tf) will yield evaporation and mass will move from the particle to the surrounding (negative mass transfer) which is consistent with a positive Spalding number. The second mass transfer effect is the diffusivity of the particle species while in the surrounding fluid, which is denoted by f,p, For a droplet evaporating or condensing in a gas, f,p is given by the molecular diffusion of the droplet vapor within the surrounding fluid, which indicates its net diffusivity (Eq. A.42). For a bubble condensing in a liquid or a particle melting in a liquid, f,p becomes the molecular diffusion of the particles species in a liquid based on conditions just outside of the particle surface. Higher diffusivity allows the species to be move faster toward the particle (positive mass transfer) or faster away from the particle surface (negative mass transfer). One may then define a normalized mass transfer as the Sherwood number (based on the surrounding fluid density and the particle diameter): p m Sh 1.80 df Bf ,p

d ( pd )

p m


A more formal and general definition of the Sherwood number is given by Sirignano (1999). Based on the Sherwood number definition, the mass transfer expression of Eq. 1.78 can be written as: p Bf ,pSh m dd = f = 1.81 pd p d 2 dt For creeping flow (Rep0) and saturated conditions (B0), an exact solution for the mass transfer rate can be obtained using the Stokes solution of velocity field coupled with solution of the concentration diffusion equation for negligible convection. The resulting analytical solution corresponds to Sh=2 (Clift et al. 1978). If the Spalding number and the diffusivity are also constant (consistent with weak changes in temperature) and the particle density is constant, this yields a particle surface area which changes linearly in time, a result which is often referred to the d2 law: 4 f Bf ,pSh t d2 = d2 t = d2 1 1.82a t =o t o =  p m
m  pd 2 t =o 4 f Bf ,pSh


Thus, the mass-transfer time-scale ( m  ) indicates the time it will take for the particle mass to disappear (for negative mass transfer corresponding to B>0) or the time it will take for the particle surface area to double (for positive mass transfer corresponding to B<0). For finite Reynolds and Spalding numbers, these heat and mass transfer rates can be enhanced (beyond their creeping values of 2) due to convection effects. The most common approach to take this into account is to employ the empirical Ranz-Marshall (1952) relationships which assume that the mass transfer Spalding number based on is equal to the heat transfer Spalding number: 1/3 ln (1 + B ) Sh = ( 2 + 0.6 Re1/2 1.83a p Sc f ) B 1/3 ln (1 + B ) 1.83b Nu = ( 2 + 0.6 Re1/2 p Prf ) B For these expressions, the Schmidt number of the surrounding fluid is defined by Eq. A.45 and the Prandtl number by Eq. A.46. While these expressions are quite common and useful, there are more sophisticated expressions for Sh and Nu which take into account variable Spalding numbers, conditions of high mass transfer, effects of particle spacing, etc. (Sirignano, 1999). For bubbles, evaporation and condensation can be related to vapor pressure while the temperature is assumed to be approximately constant for both phases, owing to high thermal inertia of the surrounding liquid (Owis & Nayfeh, 2003).