Document Sample

APPENDICES: FLUID DYNAMICS AND COMPUTATION FOR BUBBLES, DROPS AND PARTICLES E. Loth Draft for Cambridge University Press April 20, 2009 1 Appendix A: Single-Phase Flow Equations and Regimes Some of the common assumptions and formulations for the surrounding-fluid conservation equations will be reviewed in this appendix, including the PDEs for various flow regimes. Derivation details of these equations are available in Liepmann & Roshko (1957), White (1991) and Schlichting (1979), while the numerical discretization (in space and time) of these equations will be discussed in Appendix B. In the first section of this appendix (A.1), the continuous-phase equations are discussed along with relevant equations of state and flow properties. Afterwards, the equations will be considered for idealized conditions of inviscid compressible flow (A.2) and viscous incompressible flow (A.3 and A.4). The appendix concludes with a description of statistical properties and averaging techniques for turbulent flow (A.5). Since no particles are considered in either Appendix, u=U and p=P for the continuous-phase. For convenience, the appendices use u and p to be consistent with most of the transport equations in the main text. However, one may substitute U and P to describe local flow around a particle and V and Pp to consider flow inside a particle. A.1. Continuous-Flow Governing Equations General Conservation Equations The Partial Differential Equations (PDEs) which govern the uncoupled continuous-phase flow are based on the conservation of three quantities: mass, momentum, and energy. These equations can be obtained by equating the time rate of change of a quantity to the flux across the surface of that quantity and any source or sink of that quantity within the volume. If q is an arbitrary quantity, then the Reynolds Transport Theorem states that the time-rate of change of q within an arbitrary volume is given by the inward flux through the surface of the volume and the internal generation rate of the quantity. This concept is shown in Fig. A.1 and can be expressed as: q t d q u n dA q d A A.1 In this equation, u is the velocity of the quantity transport, n is the outward normal to the surface, and q is the source strength within the volume. If q and u are finite and continuously differentiable in space, then the first term on the RHS of this equation is the inward flux to the control volume. The opposite of this term is the outward flux. This can be transformed into a volume integral by the Gauss Divergence theorem for a vector q: A q n dA q d A.2 Rewriting the vector as (qu) yields a relationship between outward flux and internal gradients: q u n dA (qu) d A A.3 Therefore, the PDE for the transport of q in an infinitesimally small fluid volume is 2 q (qu) q A.4 t This is often referred to as the conservation equation if q =0. The conservation of mass with no sources or sinks can be written in terms of the density by noting that it is linearly proportional to the fluid mass if the differential fluid volume is held fixed. Based on the conservation equation and the definition of the substantial derivative, the conservation of mass in differential form can be expressed in either the Eulerian derivative or the Lagrangian derivative of the density f f u 0 A.5a t D f f u 0 A.5b Dt The first equation is particularly convenient for Eulerian formulations of the continuous-phase flow. The second equation makes use of the relationship defined in Eq. 1.13b. The conservation of momentum of the fluid in a differential volume is dictated by equating the change in momentum with the forces applied to the mass in the control volume (Newton’s second law). The applied forces include body forces and surface forces. The most common body force is the gravitational force mf g = f g , but others include electrostatic and magnetic forces. The surface forces are based on the fluid pressure (p) and viscous stresses (Kij) acting as area-integrated quantities. These forces can be converted into volume-integrated quantities by again using the divergence theorem. Based on this approach, the conservation equation (assuming no source or sinks for momentum in the control volume) can be written in the following differential forms as f u f uu f g - p K i j A.6a t Du f f g - p K i j A.6b Dt Note that the conversion of Eq. A.6a to A.6b is based on the continuity equation (Eq. A.5b) and the expansion of the spatial velocity gradient term f uu f u u u f u A.7 Note that the continuity and momentum equations given by Eqs. A.5a & A.6a are in the so- called “conservative” form since all derivative terms are not multiplied by dependent variables, whereas Eqs. A.5b & A.6b are in the “total derivative” or “non-conservative” form. To complete these equations it is important to relate the viscous stresses to the velocity field (the pressure will be later related to density by an equation of state). For a Newtonian fluid, these stresses are linearly proportional to the deformation. The viscous stresses can be expressed in Cartesian coordinates in terms of the shearing and compressive velocity gradients as well as the shear and bulk viscosities: u u j u k K i j f i x j x i bulk,f x k A.8 Often the last term is negligible compared to the shear stresses especially for low-speed flows, although this compressive viscous stress is needed to resolve the features within a shock-wave 3 or other highly compressive phenomena. To simplify, we may invoke the Stokes’ hypothesis to relate the bulk viscosity to the shear-based dynamic viscosity: bulk,f 2f / 3 . Further employing the divergence definition, the viscous stresses can be expressed as u u j 2 u k u i u j 2 K i j f i ij f ij u A.9 x j x i 3 x j x i 3 x k In the equation, ij is the Kronecker delta defined as 1 if i=j i j A.10 0 if i j This tensor is symmetric such that K ij K ji . White (1991) notes that the Stokes’ hypothesis may not be correct but since generally u is typically small (e.g. zero for incompressible flow), the exact nature of the bulk viscosity is not critical, although it can be important for diffusion effects related to sound attenuation and shock waves. One may also combine the viscous and the pressure stresses with a surface normal (nj) to define a fluid stress vector as: Tj pn j K ijn i A.11 This vector represents the force per unit area acting on the surface defined by nj. For “inviscid” flow, the vector notation form is T pn . To develop a transport equation for fluid energy, one may define the total energy per unit mass (etot) as the internal energy due to random molecular motion (e) plus the kinetic energy due to the mean molecular motion: e tot e + u 2 / 2 A.12 The conservation of energy can then be developed based on the first law of thermodynamics which states that the energy change of a system is equal to the heat added to the systems and the work done by the system. The heat flux across the control volume border surface can be related to the gradient of the temperature and thermal conductivity ( k f T ), where k f is the thermal conductivity coefficient of the fluid. The work is related to the power needed to overcome the fluid dynamic stresses (while shaft work is neglected). Using the Reynolds Transport Theorem, the energy conservation can be written as: f e tot f e tot u f g u + K i j u - pu k f T A.13 t This equation is in conservative form and indicates that the transport of total energy (LHS) is related to the work and energy applied to the fluid system including (RHS): including body force work (1st term), viscous stress work (2nd term), pressure stress work (3rd term) and heat transfer due to the thermal conductivity of the fluid (4th term). This transport equation can also be written in terms of the substantial derivative of the internal energy as u = K ij pij i + k f T De f A.14 Dt x j In contrast to Eq. A.13, the pressure and body force work do not appear as identifiable terms. However, it can also be written in conservative form by employing Eq. A.5a: f e u f eu = Kij pij i + k f T A.15 t x j 4 As will be discussed in §B.3.4, conservative forms of the transport equations are preferred for compressible flows. Fluid Properties At NTP, several of the fluid properties used in the above PDEs (viscosity, density, conductivity, etc.) are given in Table A.1 for air, methane, water and ethanol. Often the viscosity may be assumed constant, but for significant temperature variations (especially experienced in compressible flows), it is a function of temperature. For 250 K to 750 K, air viscosity can be estimated as f fo T / To 1.82x105 kg /(m-s) T / 293K 2/ 3 2/ 3 A.16 In this equation, the subscript “o” indicates NTP and absolute temperatures are used. This form is generally reasonable for other gasses up to about 400 K. White (1991) provides detailed viscosity relationships for several gasses and liquids. Equation of State and Compressibility Relationships In many flows, density changes of the continuous-phase occur due to fluid dynamics. In these “compressible” flows, it is important to determine the relationship between pressure and density, and this is generally referred to as the equation of state for a fluid. This relationship can also be used to determine the speed of sound in the fluid since acoustic signals are transmitted via compressible pressure waves. In addition, the speed of sound will help determine the character of the governing PDEs. In the following we will consider the equation of state for both a gas and for a liquid. In order to relate the variations in pressure, temperature, and energy it is useful to define enthalpy (h), which is proportional to the internal energy per unit mass p h e A.17 f In general, enthalpy and specific energy increase with temperature and the proportionalities are the specific heat at constant pressure (cp) and specific heat at constant volume (c): h cp A.18a T p e c A.18b T If these specific heats are independent of temperature for a particular fluid, then the fluid is called “calorically perfect” and the internal energy and enthalpy become linear functions of temperature. The ratio between these two values is the specific heat ratio cp / c A.19 For liquids at moderate pressures and temperatures, this ratio is approximately equal to unity ( c p c ). For a diatomic calorically perfect gas, the specific heat ratio is constant and equal to 5 7/5, which is a reasonable approximation for air at NTP. At higher temperatures, the gas is generally no longer calorically perfect and is reduced. The specific heat ratio can also be used to define the speed of sound as the rate of propagation of weak (isentropic) pressure pulses p a2 A.20 s The relationship between pressure and density for a gas is straightforward for the “perfect gas” assumption, whereby the molecules are assumed to be far enough apart that intermolecular forces are negligible and only molecular collisions determine their motion. In this case, pressure can be related to the volumetric concentration of molecules and their average kinetic energy. The latter can be related to the absolute gas temperature (measured in Kelvin) based on the statistical mechanics and the kinetic theory of gasses. The result is a linear relationship between the pressure and the product of temperature and density for a perfect gas. The proportionality constant is the specific gas constant, R g, and Dalton showed that this is equal to the universal gas constant (R univ) divided by the molecular weight of the gas (MWg). The resulting equation of state is known as the perfect gas law: R pg R g univ A.21 MWg g Tg Based on this relationship and the definition of enthalpy, the gas constant is equal to the difference between the specific heats: cp c R g A.22 For a calorically perfect gas, the ratio of these specific heats (Eq. A.19) is constant such that the local speed of sound of Eq. A.20 for a perfect gas can be expressed as a g R g T A.23 Most gasses at room temperature will have a speed of sound of a few hundred m/s (Table A.1). In the case of a liquid, the equation of state relating pressure to density is often given in terms of the empirical Tait equation which includes a liquid compressibility pressure constant (Bl) and a compressibility exponent (Kl) as well as the pressure and density at NTP: K p+Bl po l l A.24 (1+ Bl )po lo The speed of sound for liquid (al) from the Tait equation based on isentropic conditions becomes K l (p Bl p o ) al2 A.25 l Note that this equation reverts to the calorically perfect gas conditions for B=0 and K=For a liquid these parameters are non-zero. For example, the compressibility of water at normal temperature is reasonably describe with Bl=3000 and Kl=7.15. Typical liquid sounds speeds 6 are on the order of a 1000 m/s (Table A.1) and Thompson (1972) gives B and K for other fluids, all of which tend to much higher than that of typical gasses. The qualitative differences between the pressure and density relationship for a gas and a liquid can be seen in Fig. A.2. For gasses or liquids, the effect of compressibility on fluid motion is typically characterized by the Mach number, which is defined as the local speed of the fluid normalized by the speed of sound: Mu/a A.26 This parameter indicates the relative speeds of convection as compared to acoustic signals. If one defines the total or stagnation condition as that associated with an isentropic deceleration to zero velocity, the changes in flow properties can be related to the local Mach number as 1 1 1 T p 1 2 1 M for a perfect gas A.27 Ttot tot gas p tot gas 2 In these relationships, Ttot, tot and ptot are the total temperature, pressure and density. This equation indicates that changes in flow properties become more pronounced as the Mach number become significant. This is qualitatively similar to the influence of Mach number of liquids. Therefore, a gas or liquid flow with M0 throughout will tend to have negligible density changes due to velocity variations, while flows with M~1 or greater can be expected to have large changes in flow properties. The Mach number is useful in identifying flow regimes and the mathematical character of the PDEs as will be discussed in the next section. A.2. Inviscid Flow Equations and Mach Number Regimes In many flow conditions, the quantitative effects of friction over surfaces and dissipation of vorticity are generally not of primary importance. In such cases, one may neglect viscosity in the governing continuous-flow equations, yielding “inviscid” flow. When neglecting viscosity for the continuous-flow description, the wall boundary condition can be described as a “slip” condition. As such, only the velocity perpendicular to the wall is prescribed (and set equal to the wall speed if the wall moves relative to reference coordinate system). Let us first consider the very low Mach number regime where the density will not be altered by fluid dynamics. Incompressible Flow The M0 range corresponds to “incompressible flow” indicating changes in the flow velocity will not significantly affect the fluid density. Since liquid flows tend to have smaller flow speeds and higher acoustic speeds than gas flows, the assumption of liquid incompressibility is quite common. If the flow is adiabatic and homogeneous, then the temperature everywhere is constant (“isothermal flow”) and it no longer influences the mass and momentum equations. Thus, a continuous-phase energy equation (as in Eq. A.12) is often not needed for incompressible flow if heat transfer is negligible. If the mechanisms which cause gradients in density (density stratification, heat transfer, or mechanical compression) are 7 weak, then the continuous-fluid density (and viscosity) can be considered everywhere constant along streamlines, i.e. f const. A.28 This allows the conservation equations to be correspondingly simplified. The conservation of mass indicates that the velocity divergence is zero: u 0 A.29 This is often referred to as the continuity equation and in Cartesian coordinates becomes: u x u y u z 0 A.30 x y z This result is also valid when the flow is stratified but non-mixing since incompressibility will enforce that density is constant along streamlines Df /Dt=0 which can be combined with Eq. A.5b to yield Eq. A.29. Three effective descriptors of fluid fields are: 1) the stream function and streamlines (which serves to indicate fluid paths), 2) the vorticity field (which serves to indicate rotational and shear regions), and 3) the velocity potential (which can be related to the pressure field). These three descriptors are defined in the following but are all related. The first is the stream function () which is only defined for two-dimensional flow. For Cartesian coordinates, it is: ux A.31a y -u y A.31b x This definition automatically satisfies the continuity equation for incompressible flow (Eq. A.30). Stream-lines are instantaneous curves along which the stream function in two- dimensional flow is constant and are quite useful for visualizing the flow angles since the streamlines are always parallel to the flow direction. In addition, the mass flow between a pair of streamlines is constant so that a reduction in gap between them indicates a local increase in velocity and a decrease in pressure. Pathlines are the actual (history-based) path of a fluid element and will be equal to the streamlines if the flow is steady. The second descriptor is the continuous-phase fluid vorticity which is a vector defined as twice the angular velocity of a fluid particle (Eq. 1.7): u x i x y i y z i z A.32 The RHS expression is the vorticity in Cartesian coordinates and its magnitude (always taken as positive) can be defined with the Cartesian form or in tensor form as: u u y u x u z u y u x 22 2 z 2 2 2 y z z x x y x y z A.33a 1 u i u j u j u i 2 x j x i x i x j A.33b Two fundamental vortical flow conditions with constant vorticity are: “linear shear” flow and the “simple vortex” flow. For “linear shear” flow, the flow lines are parallel and the velocity 8 field has a uniform spatial gradient. An example is shown in Fig. A.3a with u y u z 0 where the vorticity magnitude given by the streamwise velocity gradient: u shear x const. for linear shear flow A.34 y For “simple vortex” flow, the flow lines are circular with a tangential velocity that increases linearly with radius, i.e. solid-body rotation. An example of this is shown in Fig. A.3b with u r u z 0 where the vorticity is given by: 2u vortex const. for simple vortex flow A.35 r A flow with no vorticity flow throughout (=0) is called “irrotational”. The third descriptor is the velocity potential defined as u A.36 Since q 0 for any scalar, any flow for which a velocity potential can be defined (i.e. “potential flow”) must also be irrotational. The velocity potential is always perpendicular to surfaces and streamlines and thus is helpful to visualize flux surfaces in a flowfield. Furthermore, unlike the stream function, a velocity potential can be defined for three- dimensional flow. For steady flows, this equation can be combined with the flow continuity (Eq. A.30) to yield: 2 0 A.37 Thus, the governing PDE of the steady velocity potential is given by a linear Laplace equation. If we consider a two-dimensional flow which is irrotational (=0), Eqs. A.31 and A.33 can be used to show the resulting stream function will also satisfy Laplace’s equation: 2 0 A.38 This property allows linear super-position of solutions, i.e. if 1 and 2 satisfy Eq. A.37, so then does 1+2 , and a similar statement can be made with respect to two stream function solutions. This feature is helpful for analytical solutions of flows in terms of or . If vorticity is important, the stream function equation for two-dimensional flow in the x-y plane yields a Poisson equation: 2 z A.39 This relationship may be combined with a vorticity transport equation to describe two- dimensional incompressible flow, independent of whether viscosity is included. As a result approaches have been used in both analysis and CFD to describe a flowfield with given boundary conditions, since the velocity and pressure fields can be obtained thereafter. However, most current CFD approaches use the equations for pressure and velocity directly (i.e.; “primitive” variable approach) to more generally allow three-dimensional flow fields. The primitive velocity formulation is generally given by Eq. A.6b coupled with Eq. A.29. If viscous effects are neglected and steady flow is assumed, the momentum conservation (Eq. A.6b) can be written as: 1 2 f u2 p f g p f gi z A.40 9 This equation can be integrated to yield a relationship for the dynamic pressure (pdyn) as: pdyn 12 f u2 p f gz const. A.41 This is called the Bernoulli equation and is valid for steady inviscid incompressible flows where fgz accounts for the hydrostatic pressure effect. If this term is small compared to the dynamic pressure for a given flowfield (e.g. gD u 2 ), the relationship becomes: D p 12 f u 2 p tot const. for weak hydrostatic gradients A.42 In this case, the total pressure is uniform throughout the domain and equals the pressure at the stagnation points. To apply a similar concept when hydrostatic effects should instead be retained, one may define a fluid dynamic pressure (p′) based on the RHS of Eq. A.41 so that its gradient can be expressed in terms of the velocity field by Eq. A.40: p p f gz p tot 12 f u 2 A.43a p 1 2 f u 2 A.43b Even when unsteady viscous effects are important, the fluid dynamic pressure is convenient since substitution into the momentum equation removes the gravitational term, e.g. Eq. A.6b can be rewritten as Du f p K i j A.43X Dt This also eliminates the need to include hydrostatic effects in the pressure boundary conditions. Once this fluid dynamic pressure distribution is known, then Eq. A.43a can be used to determine the local pressure. Compressible Flow If the condition M0 is not satisfied (say, M>0.1), then the governing equations should be written in compressible form since local changes in velocity will significantly affect the density. This compressibility condition also mandates that the energy equation be considered along with the momentum and mass conservation equations. The resulting set of compressible flow equations for inviscid adiabatic conditions is usually arranged in conservative form. Rearranging Eqs. A.5a, A.6a and A.13 and neglecting hydrostatic pressure gradients (Eq. A.43b) in the momentum and energy equations (since gD u 2 for gas most flows), these D equations can be written as: f f u j 0 A.44a t x j f u i f u i u j pij 0 A.44b t x j f e tot f u je tot pu j 0 A.44c t x j 10 This formulation is convenient for numerical approaches as it prescribes linear derivatives of the conservative variables (f, fui, and fetot), to which finite volume approaches can be applied in a straightforward fashion. If one assumes that the flow is irrotational, then a velocity potential can be again defined for compressible flow by Eq. A.36. If we further assume steady isentropic conditions (no shocks or viscous effects) and a calorically perfect gas, then the equations can be linearized based on small velocity perturbations from a free-stream velocity in the x-direction (ux∞). The equations of continuity, momentum and energy can be employed to yield the linearized steady compressible velocity potential equation: 2 2 2 2 (1 M ) 2 2 2 0 A.45 x y z In this equation, the free-stream Mach number is defined as M∞=ux∞/a∞. This can then be evaluated for the different Mach number regimes in terms of flow compressibility: M«1, M<1, M~1, and M>1. In addition, we can identify the mathematic character of the governing PDE (Eq. A.45) by considering the slope of the characteristic directions (Anderson et al. 1997): dy dz 1 A.46 dx dx M 1 2 Characteristics for the four Mach number regimes are shown in Fig. A.4 and discussed below. For M«1, we note that the pressure waves travel at the speed of sound since the convection speed is approximately negligible in comparison. However, the pressure waves are generally so weak that they have a negligible quantitative effect on the velocity field. Therefore, the incompressibility assumption is generally invoked under these conditions (Eq. A.45 reverts to the incompressible form given in Eq. A.37) unless acoustic signals need to be captured. The M<1 range corresponds to “subsonic flow” indicating that the convection speed is less than the acoustic speed, but is generally not small in comparison. As such, the pressure signals which move according to the combination of the local sound speed and convection speed (u+a) can travel in all directions. This is consistent with the fact that Eq. A.45 has imaginary characteristic slopes for M∞<1 so that it is an elliptic PDE such that the flow solution is coupled in all directions. As such, all points within the domain can affect all other points in the domain, i.e. both the region of influence and domain of dependence for a given point include the entire domain. We shall see later in Appendix B that the numerical solution should mimic this behavior. The M~1 condition corresponds to “transonic flow” indicating that the convection speed is approximately equal to the acoustic speed. If M∞=1, the pressure signals can only move downstream and Eq. A.45 becomes a parabolic PDE in the flow direction. For flow primarily in the x-direction, this is consistent with characteristic curves in the x-y plane having a slope (dy/dx) equal to ∞. Under this condition, the region of influence for a given point is confined to the downstream domain. The M>1 condition corresponds to “supersonic flow” whereby pressure signals are further limited to only move in the downstream within a cone of influence. The cone angle is consistent with the slope of the characteristic curves in the x-y plane as M 1 1/ 2 2 . Such signal propagation indicates a hyperbolic PDE at this condition. 11 A.3. Incompressible Viscous Flow Equations Diffusion Parameters Viscous flows require the inclusion of shear stresses in the PDE formulation. Viscosity is the mechanism which characterizes the rate of momentum transport due to molecular interactions. For Newtonian flows, the magnitude of the viscosity is given by the linear constant of proportionality between the velocity gradients and the shear stress (Eq. A.8). Thus, the dynamic viscosity of a fluid (f) is the rate of momentum transport in the fluid per unit volume. Similarly, the kinematic viscosity corresponds to the rate of momentum transport per unit mass and is related to the dynamics viscosity by: f f A.47 f The presence of viscosity causes flow structures with velocity gradients to diffuse and dissipate. Viscosity also causes a no-slip wall condition if the fluid can be considered as a continuum, i.e. the relative flow along the wall must be zero. For example, if the wall is static then the flow must have zero velocity at the wall. There are two other important diffusion properties of a fluid: mass diffusion and thermal diffusion. The equations which govern these effects are discussed below. The mass diffusion is analogous to the momentum diffusion concept except it employs the concentration gradient instead of the velocity gradient. If one applies Fick’s law, the mass diffusion of a species per unit area is linearly proportional to the concentration gradient: mi f i A.48 A The LHS is the mass flux per unit area of species in the direction of decreasing concentration, i is the mass of species i per unit volume of the mixture, and f is the mass diffusion coefficient based on the rate of molecular exchange between the species and the overall mixture (White, 1991). The negative sign indicates that diffusion occurs in the direction of the lower concentration region yielding a tendency towards uniform concentration. If we consider the diffusion of only one species, and define its local volume fraction to be , then the component density can be related to the overall fluid density by k k f A.49 The transport equation of the species volume fraction can then be obtained for both variable and constant density conditions as (Hinze, 1975): f k f k u i f ,k f k A.50a t x i x i k k u i f ,k k for constant density A.50b t x i x i 12 This mass diffusion assumes that the species are miscible within the overall mixture. This is a good assumption for gasses since their molecules can readily form mixtures and this is also true for many liquids. However, some liquids (such as oil and water) are not miscible. If the tracer species is not miscible, then f=0. Note that the mass diffusion of a species based on f is analogous to the momentum diffusion based on f (Eq. A.47). Because of this, it is natural to consider the ratio of the momentum to mass diffusion, which is defined as the (laminar) Schmidt number f Sc A.51 f For air and many gasses, the Schmidt number is of order unity indicating that the mass and momentum diffusion rates are similar. However, water and many liquids have a Schmidt number on the order of 103, indicating such fluids have a mass diffusion which is much slower than their momentum diffusion. Similar to the above, thermal diffusion in fluids is assumed to be linearly proportional to the local temperature gradient. This is referred to as Fourier’s Law and the constant of proportionality for thermal diffusion is kf, which relates temperature gradients into heat diffusion (Eqs. A.13-A.15). The ratio of momentum diffusion to thermal diffusion is given by the (laminar) Prandtl number f f c p,f Prf k f / f c p,f A.52 kf In this expression, cp is the heat capacity at constant pressure. For gases, the Prandtl number is of order unity (about 0.7 for air at NTP) indicating that thermal and momentum diffusion rates are generally similar. For liquids, the Prandtl number varies considerably with both the fluid type and the temperature: Pr ranges from 1.7-13.7 for water, 0.004-0.03 for liquid metals, and 50-100,000 for oils. Governing Equations Viscous flows can be classified as either compressible or incompressible. However, to focus on the influence of viscosity and simplify matters, only incompressible flow will be considered. For isothermal homogenous flows, the viscosity and density can be considered constant throughout. As noted above, if an incompressible flow is homogenous, adiabatic and non-reacting, the temperature is approximately constant. For such an isothermal flow, the influence of the energy equation on the mass and momentum equations can be neglected. For incompressible flow, the mass conservation equation is simply given by Eq. A.29. For the momentum conservation, the assumption of incompressibility allows the bulk stress to be neglected in Eq. A.8. Further assuming constant viscosity, the gradient of the viscous stress tensor of Eq. A.9 in tensor and vector notation is given by: u i u j K i j f f u u 2 x j x j x i A.53 13 Since the velocity divergence is zero for incompressible flow (Eq. A.29), the momentum equation of Eq. A.6b in vector form becomes: Du f p f g f 2u A.54 Dt The Cartesian tensor forms of the continuity and the momentum equations are thus: u i 0 A.55a x i u i u p 2ui f f u j i f g i - f A.55b t x j x i x 2j In the second equation, gi is the component of gravitational acceleration in the xi direction. In the limit of quiescent flow (ui=0), this yields p / x i f g i A.56 This is the hydrostatic equation It is often convenient to consider incompressible flow in spherical coordinates. Defining r as the radial coordinal, as the polar angle, and as the azimuthal angle about the axis =0, the incompressible continuity equation in spherical coordinates becomes 1 r ur 2 1 u sin u 0 A.57 r 2 r r sin The continuity equation in cylindrical and general curvilinear orthogonal coordinates is also given by Batchelor (1967) as are the momentum equations in spherical coordinates: u r u r u u r u u r u u 2 2 ur gr t r r r sin r r A.58a 1 p 2 2u r 2 u sin u f u r 2 2 f r r r sin u u u u u u u r u u cot 2 ur g t r r r sin r r 2 1 p 2 u r 1 u A.58b f u 2 2 2 u 2 cos f r r r sin u u u u u u u r u u u cot ur g t r r r sin r r A.58c 1 p 2 2 u r u u u 2 f cot f r sin r sin 2sin The viscous terms make use of the spherical Laplacian operator 1 q 1 q 1 2q 2q 2 r 2 2 sin 2 2 A.59 r r r r sin r sin 2 For reference, the viscous stress tensor in spherical coordinates is also given as: 14 u r K rr 2f A.60a r u 1 u r K r f r A.60b r r r u 1 u K 2f r A.60c r r 1 u u r u cot K 2 f A.60d r sin r r 1 u r r u K r 2 f A.60e 2r sin 2 r r sin u 1 u K 2 f A.60f 2r sin 2r sin It can be important to identify such stress on a spherical interface. In either set of coordinates, the flow variables can be non-dimensionalized by defining D and uD as respectively the reference macroscopic continuous-phase length and speed (e.g. diameter and mean speed in a pipe flow). Thus, Eq. A.55 can be rewritten in vector notation with non-dimensional values u*=u/uD, x*=x/D, t*=tuD/D and p*=p/f(uD)2 as: u* 0 A.61a Du* 1 g 1 * p* 2 u* Dt FrD g Re D A.61b In the RHS of Eq, A.61b, g/g is the unit vector in the direction of gravity. The latter equation includes the non-dimensional macroscopic Froude and Reynolds numbers: u2 FrD D A.62a gD Du Re D f D A.62b f This Froude number definition of Eq. A.62a is consistent with that of Oguz & Prosperetti (1990) and is sometimes called the “second” Froude number, whereby the “first” Froude number is simply the square root of this value. Based on Eq. A.41, it is proportional to the ratio of dynamic pressure to hydrostatic pressure and thus can reflect the importance of convection to gravitational forces. If this number is large (FrD»1), then gravitational effects can generally be ignored with respect to overall flow development. The Reynolds number of Eq. A.62b represents the overall ratio of convective stresses to viscous stresses for the macroscopic continuous-phase flow. The uncoupled viscous flow physics will thus depend on the Reynolds number, the boundary conditions, and the domain geometry. However, as will be discussed in the next section, the viscous effects can generally not be ignored even when ReD»1 because local viscous regions can affect local skin friction. This, in turn, can control flow separation characteristics and affect the entire flow field. In 15 local regions where viscous effects are strong, the shear stresses in the momentum transport of Eq. A.61b render the PDE to be elliptic and consistent with molecular diffusion occurring in all directions. The relevant transport equations are discussed in the next section with respect to the different Reynolds number regimes. A.4. Reynolds Number Regimes for the Continuous-Phase Overview of Regimes Viscous flow behavior is generally classified into three different regimes according to viscous flow physics: laminar, transitional, or turbulent flow. This division is based on the appearance of undamped flow instabilities as the Reynolds number increases. When the instabilities are fully damped by viscosity, the flow is said to be “laminar”. When these instabilities are partially damped, perturbations can give rise to simple (e.g. two-dimensional or axisymmetric) unsteady flow features and the flow is called “transitional”. The onset of transitional flow occurs at the critical Reynolds number (ReD,crit), but this is often quickly followed by turbulent flow (at ReD,turb) whereby strong flow instabilities are observed. This results in a three-dimensional, unsteady flow with non-linear complex coupling between the various instability modes. The values of ReD,crit and ReD,turb are generally obtained from experiments for a given geometry, but for some simple conditions (e.g. flat plate boundary layer flow) it is possible to at least qualitatively predict the critical Reynolds number. The instabilities can occur in wall-shear flows due to the Tollmien-Schlichting instability and in density gradient flows due to the Rayleigh-Taylor instability. Another example is Kelvin-Helmholtz instability in free-shear flows, which result from velocity gradients between parallel flowing streams. Perturbations can be damped by viscosity at low Reynolds numbers, but a higher values will tend to grow and eventually lead to eddy and braid features (Fig. A.5). For transitional flows, the structures will exhibit detailed features and become increasingly non-linear, though still remain approximately two-dimensional (Fig. A.5b). As the Reynolds number increases further, significant three-dimensionality and increased flow complexity will occur. As such, transitional flows can sometimes be intermittent with some laminar-like flow regions (organized with simple features) but also intermingled with bursts of turbulent-like flow (quasi-random with complex features). For ReD>ReD,turb (the minimum Reynolds number for fully-developed turbulence), the flow instabilities become sufficiently pronounced that non–linear chaotic behavior is sustainable. The result is a flow-field with vortices and other features which occur and interact over a wide range of length and time scales. As an example, a fully turbulent free-shear layer is shown in Fig. A.5c, where it can be seen that the large scale structures of Fig. A.5b are present but are combined in a complex way with a spectrum of other structures. These structures are three- dimensional, unsteady, and effectively stochastic at the smaller scales. At the smallest scales, molecular diffusion effects govern the velocity, concentration, and thermal gradients. The Reynolds numbers associated with transition and turbulence are a function of the flow-field geometry. For example, ReD,crit for flat plate boundary layer is on the order of 200,000- 400,000 based on plate length, whereas ReD,crit is about 2,000-3,000 for a round pipe flow and even less for a free shear layer. Fully-developed turbulence generally occurs soon after transition as shown in Fig. A.6 for a pipe flow with ReD,turb≈4,000-6,000. 16 The influence of the various Reynolds number regimes is illustrated in Fig. A.7 for the case of species mass diffusion and convection in a pipe flow with Schmidt number of unity (so that mass and momentum diffusion are of the same order) over a range of Reynolds numbers. Here it can be seen that increasing the Reynolds can eventually lead from steady flow to turbulent flow. In this sense, these figures are similar to the sketches of dye dispersion for the celebrated pipe flow stability experiments conducted by Reynolds in 1883. For very low Reynolds number (ReD«1), mass diffusion will dominate the convection and thus a mass source will tend to propagate upstream nearly as fast as it does downstream (Fig. A.7a). For ReD~1, the convective and viscous diffusion rates will be similar order so that upstream diffusion is significant but limited (Fig. A.7b). As the Reynolds number increases further but the flow remains laminar, the convection effects will dominate such that the transverse diffusion is limited and the mass is mostly carried in the flow direction (Fig. A.7c). Further increases in Reynolds number will cause the flow to be unstable and the resulting unsteady structures of transitional flow distribute the injection concentration with significant lateral perturbations (Fig. A.7d). A further increase in Reynolds number leads to turbulent flow, where the injected mass mixing becomes non-linear and distributed over a wide range of length and time scales (Fig. A.7e). The transverse fluctuations in this regime result in significant lateral mixing, which can be seen to be much greater than that observed for high Reynolds number laminar flow (Fig. A.7c). Laminar Flow Momentum Equations The momentum conservation given by Eq. A.55b assumed constant viscosity and constant density and under those assumptions is valid for all Reynolds numbers for both laminar and turbulent flow. Two more simplifying assumptions are made for the following discussion of laminar flow: the flow is assumed to be steady and gravitational forces are neglected. The first assumption is often reasonable since flow instabilities tend to be damped by viscous forces and so the absence of any time-dependent external forcing or time-dependent boundary conditions will yield a steady flow. Using the definitions of Eq. A.62, the second assumption is reasonable for large ReD and large FrD so that the pressure distribution is dominated by hydrodynamic variations (Eq. A.42). This second assumption is also reasonable if pressure distribution is dominated by viscous forces, i.e. small ReD and small ReD /FrD . Laminar flow may also be classified according to the magnitude of the Reynolds number including ReD0, ReD«1, ReD~1, and 1«ReD<ReD,crit. Each of these regimes may have a unique set of momentum equations. Based on Eq. A.61b and ReD0, the LHS convection term becomes negligible compared to the viscous stress terms. However, special consideration should be given to the pressure term. In particular, regions of the flow where viscous effects dominate suggest that this term be non-dimensionalized as p*=pD/fuD. Employing this alternate definition of p* into Eq. A.61b shows that the pressure term should be retained in the limit of Re D 0 . The resulting steady-state momentum PDE (neglecting hydrostatic effects or by using the fluid dynamic pressure of Eq. A.43) becomes: 17 p 2ui f for ReD 0 A.63 x i x 2 j This often termed the “creeping” flow equation and is typically used for ReD<0.01. For such flows, the pressure term balances the shear stress and neither, but neither are important when the velocity field only varies linearly in space. Also, creeping flows are “reversible” in the sense that if a flowfield u satisfies this equation, then so does –u. For linearized flow with ReD«1 (e.g., 0.01<ReD<0.1), the convective terms are small but finite. As such, the convection velocity may be approximated by the free-stream velocity while the spatial gradients are based on the local fluid velocity. The resulting PDE becomes: u i p 2ui f u j f for ReD 1 A.64 x j x i x 2 j This linearization is especially convenient for analytical solutions. In either creeping or linearized flow, viscous effects on a body will directly influence the velocity field far away and in all directions since the momentum diffusion is slow or negligible. For the intermediate condition of ReD~1, the convective terms are generally on the order of the viscous terms and the pressure terms so that Eq. A.55a yields: u i p 2ui f u j f for ReD~1 A.65 x j x i x 2 j If hydrostatic forces are important, the fluid dynamic pressure (Eq. A.43a) can be used in this equation so that the gravitational terms are implicitly included. In the other limit of ReD»1 with laminar flow, the convection terms will generally dominate the viscous terms in the overall flow field. However, viscous stresses can be important in certain parts of the flow, such as along a wall where a no-slip wall condition applies. Since the Reynolds number is large, the viscous effects will only be felt a small distance away from the wall, which leads to an important characteristic of many high Reynolds number laminar flows: the boundary layer approximation. As postulated by Prandtl in 1904, the boundary layer is a thin region where the viscous effects are significant. A consequence of this flow feature is that the velocity gradients in the streamwise direction are often much smaller compared to those in the transverse direction. As an example, the momentum equations for a two-dimensional attached boundary layer flow along a flat wall in the x-direction become: p 0 A.66a y for boundary layers along x u u p 2u x f u x x f u y x f A.66b x y x y 2 This simplification allows the PDEs to be parabolic in the direction of flow, which substantially reduces the numerical solution of the flow and thus is often employed. However, there are several conditions where Eq. A.66 will not apply for ReD»1. Obviously, three-dimensional geometries or boundary conditions will require a three-dimensional description of the flow. In addition, complex surface features (such as corners or other rapid curvature changes) and adverse pressure gradients can cause flow separation. In this case, the 18 local velocity gradients will be significant in all directions. Also, while laminar flows generally damp out initial perturbations, their ability to do so is reduced as the Reynolds number increases since this corresponds to decreasing viscosity. If this leads to flow separation, then the boundary layer assumptions needed for Eq. A.66 are no longer appropriate. If the flow is also unsteady, then the full Navier-stokes equations (Eq. A.55) must be used. Transitional and Turbulent Momentum Equations While the governing equations are known for transitional and turbulent flow (e.g. Eq. A.55) the flow is too non-linear, unsteady and three-dimensional to allow an analytical solution. Therefore, the equations must be solved numerically in their unsteady form to obtain an exact solution. Since this is often not practical, another option is to consider only the time-averaged flow equations since the mean flowfields for many configurations (e.g. flat wall boundary layers or pipe flows) are repeatable. The mean flow distribution is related to the statistical properties of the various turbulent scales (including velocity, length and time scales) and so these are also important to the overall flow development. The mean properties and flow equations for turbulence are discussed in the next section. A.5. Average Properties and Equations for Turbulent Flow To consider the time-scales associated with turbulent flow and other non-linear eddy- containing flows, it is helpful to review some of their basic characteristics. Turbulence is a complex flow phenomenon that occurs at high Reynolds numbers (as discussed in §A.4). In general, turbulent flows have a wide range of time- and length-scales. It is also three- dimensional, unsteady and inherently unstable. Finally, turbulence is effectively stochastic since the state of the flow at a given time can not be used to deterministically predict the detailed state at a significantly later time. However, as with many processes with random-like features, turbulent flows and other non-linear eddy-laden flows can be characterized in terms of their mean statistics, which are generally reproducible for the same geometry and flow boundary conditions. In the following section, several key statistical parameters will be defined for turbulence including: mean and rms of the unsteady velocity components, turbulent kinetic energy and dissipation, mean correlations as a function of time and space shifts, integral time and length scales, and micro time- and length-scales. These turbulence scales are then used to develop transport equations to numerically solve for the mean velocity and other statistical properties. These transport equations are sub-divided into two primary categories: those which only consider the time-averaged fields (also called Reynolds-averaged) and those which try to fully- or partially-resolve the unsteady turbulent features. A.5.1. Averages of Fluctuating Velocities The most common statistical treatment is a simple time-average at a fixed (Eulerian) point in space for a fluid property q: 19 T 1 T q(x) q(x, t)dt as T A.67 0 In this equation, the overbar indicates a time-average over a sufficiently long period of time (T) such that the average converges (ergodically stationary). This averaging is often called Reynolds-averaging and indicates that the instantaneous value of a property can be decomposed into time-averaged and fluctuating components q(x, t) q(x) q(x, t) A.68 Specifically, the instantaneous velocity (as well as density, temperature, pressure, etc.) can be decomposed as: u i (x, t) u i (x) u (x, t) i A.69 The differences between the mean and instantaneous velocity fields are shown in Fig. A.8 for a turbulent boundary layer. Here it can be seen that all the turbulent structures are removed for the mean velocity field. An illustration of the time-history of this decomposition at a single point is shown in Fig. A.9a. There are several ways to characterize the turbulent velocity fluctuations. The Eulerian mean of the fluctuating components is identically zero ( u =0 for i=1,2,3). However, its i variance about the mean results in a finite root-mean-square (rms) value 1/2 1 u rms (x) u u ( x) u ( x, t) dt 2 as A.70 0 i, i i i The rms is often used to denote the strength of the turbulence. An example of the time-history of a component of the instantaneous fluctuations is shown in Fig. A.9b. The continuous-phase turbulent kinetic energy (k) has units of energy per unit mass and is defined as one-half of the combined strength of the turbulence in all three Cartesian directions k(x) 1 uu (x) 1 u u (x) uy uy (x) u u (x) 2 i i 2 x x z z A.71 One may define a velocity fluctuation strength as: u u rms 2 3 k A.72 This can be used to non-dimensionalize the velocity fluctuations seen in Fig. A.9b, whereas a turbulent integral time-scale (to be defined in §A.5.3) can be used to non-dimensionalize time. Another property associated with the turbulence fluctuations is the turbulent energy dissipation () which is related to the conversion of velocity fluctuations into thermal energy. The viscous stress work of Eq. A.13 indicates that this conversion will be related to the local and instantaneous velocity gradients. If we consider only the energy conversion associated with the fluctuations, turbulent dissipation can be expressed based on the magnitude of the turbulent velocity gradients (White, 1991): 2 u u u f i i 15f 1 A.73 x j x j x1 This definition reflects the conversion of turbulent fluctuations into heat at the smallest scales, and is the conventional expression used in many texts. The approximation assumes isotropy of the turbulence, a concept which will be discussed in §A.5.2. Note that Schlichting (1979) and 20 others include additional terms and thus define the overall dissipation somewhat differently. Methods to estimate turbulent dissipation from velocity measurements are discussed by Elsner & Elsner (1996). For reference, the dissipation and kinetic energy can be related to the macroscopic length and velocity scales for fully turbulent pipe flows as k 0.004u 2 and 0.005u 3 / D though there are also effects of ReD (Hinze, 1975). D D There will also be additional energy dissipation into heat associated with the mean velocity gradients but this is generally orders of magnitude smaller for turbulent flows except in near- wall regions. This is because the magnitude of the instantaneous spatial gradients will far exceed that of the mean velocity gradients. This can be especially seen in the turbulent boundary layer of Fig. A.8 which shows both the instantaneous streamwise velocity field, u x (x, y, z, t) , and its time-averaged value, u x (x, y) . This figure also shows that the mean velocity field effectively removes all the complex eddy features. A.5.2. Homogeneous Isotropic Stationary Turbulence Based on the above statistics, one can determine if it is reasonable to assume homogeneous isotropic stationary turbulence (HIST) based on the mean properties. The term “homogeneous” indicates that the kinetic energy is independent of position (no gradients of the velocity rms values in any direction). The term “isotropic” indicates that the average velocity fluctuations in all three directions are equal at a given position. A homogeneous isotropic turbulent flow thus has the same rms of the velocity fluctuations at all locations and in all directions. It also has no cross-correlations of velocity fluctuations and these properties may be summarized as: k u 0 rms for homogenous turbulence A.74a urms u u u u u u 2 k 2 x x y y z z 3 for isotropic turbulence A.74b uu i j i j 0 for homogenous isotropic turbulence A.74c The term “stationary” indicates that the mean statistics do not vary in time. In summary, the rms of the velocity fluctuations in HIST can be represented by a single value which is steady and independent of direction and spatial location, i.e., u rms ≠F(x,t,i). i, The concept of HIST is “historical” and is a useful foundation for much of the turbulent flow analysis since it allows discussion of basic turbulence concepts along with key analytical simplifications (Pope, 2000). However, it is important to note that HIST is an idealization and that it cannot be realized in actual flow conditions. This is because turbulence is generally created by solid surfaces or instabilities associated with flow gradients and turbulence will decay due to viscous dissipation unless continually supplied with new energy. Despite this, portions of some flows may be reasonably approximated as containing HIST. For example, the boundary layer statistics shown in Fig. A.10. Near the wall, the boundary layer has large (y) gradients in the turbulence and the magnitudes of the turbulence vary depending on the direction, i.e. represents a NAsT flow. Far from the wall (y/≈1), the turbulence strength is nearly the same in all directions and does not have significant spatial variations, i.e. represents an approximately HIST flow. Another example of a nearly HIST flow is the wake far downstream of a planar grid. Snyder & Lumley (1971) made careful measurements for such a 21 flow and found that the rms of the three fluctuations components of velocity were approximately equal indicating isotropic turbulence. Although there was decay of turbulent energy in the streamwise direction, it was limited to small fractions of the overall level over a typical turbulent eddy length scale (to be defined below) indicating nearly homogeneous turbulence In contrast, non-homogeneous anisotropic turbulence (NAsT) contains non-zero gradients of the turbulent kinetic energy ( k 0 ) as well as variations in magnitude among the three velocity components ( u x,rms u y,rms u ). For example, the mean transverse profiles of the z,rms three fluctuating velocity components near the wall (y/≈0.05) in Fig. A.10 indicate that the spatial variations in the velocity fluctuations are quite large and that the stream-wise turbulence intensity is significantly greater than the transverse intensity. This flow is “nasty” in that it greatly complicates the analysis of the fluid physics and multiphase interactions. Other examples of NAsT flows include pipe flows (where the gradients and anisotropy are most profound along the wall) and free-shear flows such as shear layers and jets (where the gradients and anisotropy are distributed throughout the flow but tend to be more gradual). A.5.3. Integral Scales in Turbulent Flows For either HIST or NAsT, the characterization of the turbulence can also be assessed in terms of its correlation over time and space. This is most often quantified via the velocity correlation tensor (also called the temporal correlation function) between velocity component fluctuations defined for a time-shift ( and a spatial vector shift (x) as: u (x, t)u j ( x x , t ) ij (x, x , τ) i NSI A.75 u (x)uj,rms (x x ) i,rms Note that Eq. A.75 is written with no summation of the indices (NSI). This correlation function considers the mean relationship of velocity perturbations with temporal and spatial shifts. With no temporal and spatial shifts, the correlation function becomes unity, i.e. ij (0, 0) 1 . Finite spatial or temporal shifts will result in correlation values less than unity due to changes in the fluctuations over time or space. This decay occurs due to the break-up, merging, distortion, dissipation, etc. of the turbulent structures. For long times or large distances, the correlation will completely decay, i.e. ij (, 0) ij (0, ) 0 . For homogeneous isotropic turbulence (Eq. A.47), the temporal correlations are only a function of the shifts and are independent of location and time. Furthermore, their diagonal components are uniform so that: τ 11 τ 22 τ 33 τ f (x, t) A.76 As such, the correlation can be expressed as a single quantity (no longer a tensor) as a function of the temporal and spatial shifts. For the rest of this section, HIST will generally be assumed so that the temporal correlations can be given as a single isotropic quantity. However, as will be shown later in this section, that the spatial correlations are generally anisotropic even for HIST flows. The three most common temporal shifts are Lagrangian, Eulerian and “moving-Eulerian” and are illustrated in Fig. A.11. The Lagrangian shift follows the instantaneous fluid path so that the associated spatial shift, xL. is based on integration of u over the temporal shift . The 22 Eulerian shift is simply associated with a fixed point in space, i.e. x E=0. The moving-Eulerian temporal shift is associated with the mean fluid path so that the spatial shift is based on the mean velocity vector, i.e. xmE = u . This resulting correlation, mE , is a more convenient to measure than the Lagrangian correlation since xmE is fixed and known in advance for a given time shift. The correlation functions for these three shifts in HIST are thus summarized as: L τ u (x, t)u (x uDt, t ) u u t i i i,rms i,rms A.77a t E τ u (x, t)u (x, t ) u u i i i,rms i,rms for i=1,2, or 3 A.77b mE τ u (x, t)u (x u, t ) u u i i i,rms i,rms A.77c A typical Lagrangian correlation function is shown in Fig. A.12, where it can be seen that the correlation decays from unity and tends to zero for long times, as expected. For intermediate times, the correlation is finite and can be related to the normalizing time-scale, L. This time- scale is therefore useful to define the temporal decay of turbulence and is called the Lagrangian “integral time-scale”. Integral time-scales in general are defined in the following. If above correlation functions are integrated over time, the result is simply the “integral” time-scale of the turbulence and thus can be defined for the Lagrangian (L), moving Eulerian (mE) and Eulerian (E) reference frames as follows: L L d t uDt, d A.78a t 0 0 mE mE d u, d A.78b 0 0 E E d 0, d A.78c 0 0 For the Lagrangian correlation, the integral time is also called the “eddy turnover time” since it is associated with the average time it takes for an eddy structure to change as it convects. A reasonable assumption for HIST at short-times and weak fluctuations ( u rms u ) is that the moving Eulerian and the Lagrangian correlations are approximately equal: L E A.79 This approximation is supported by data from Elghobashi & Truesdell (1984), Squires & Eaton (1990), Loth & Stedl (1999), and Coppen et al. (2001) that showed that the moving-Eulerian time-scale is only somewhat longer than the Lagrangian time-scale (about 10-20%). However, non-homogeneous anisotropic turbulence will generally yield integral time-scales that are a function of position and velocity component as shown in Fig. A.13a. If one estimates the correlation function to have an exponential decay, the relationships with the integral time-scales become simply L (τ) exp( / L ) A.80a E (τ) exp( / E ) A.80b mE (τ) exp( / mE ) A.80c 23 Thus, can be viewed as the time for a fluid velocity perturbation to decay (on average) to about 1/e (0.37) of the original value. This can be seen to be a reasonable description as shown by the exponential curve-fit and the experimental data in Fig. A.12. However, experimental data often tends to show a slight negative-loop, e.g. at /L of 4-6 in Fig. A.12. This is not predicted by the exponential decay form of Eq. A.80 so other correlation functions are sometimes employed to incorporate this negative-loop, such as the exponential-cosine function also shown in Fig. A.12. However, such functions do represent a clear improvement in terms of overall experimental agreement. Furthermore, the integral time-scale is generally more important than the shape of the correlation function. As such, the simpler exponential function is the most common function. In a manner similar to defining the integral time-scale, we may also define an integral length-scale as the characteristic distance for spatial decorrelation with a spatial shift parallel to mean velocity. For HIST, this length-scale is a vector defined by the streamwise spatial shift ( x ) but with zero temporal shift (=0): u (x, 0)u (x x , 0) i ii (x , 0)dx i i dx NSI A.81 0 0 u u i,rms i,rms As such, the length-scale depends on the direction of the turbulent fluctuations. The integral time-scale based on streamwise velocity fluctuations is defined as: u (x, 0)u (x x , 0) (x , 0)dx dx A.82 0 0 u,rms u,rms The corresponding correlation function is typically approximated by exponential decay (Fig. A.14a) since this is consistent with most experimental measurements: (x ) exp x / A.83 An example flow for which the integral-scales can be readily observed is the free-shear flow (also called the “mixing-layer”) of Fig. A.5c. Most of the energy is contained in the coherent large-scale features since they are responsible for the largest velocity fluctuations. However, the integral length-scale is averaged overall energy levels and so tends to be somewhat smaller, e.g. is roughly 1/3 of the shear-layer thickness (Wygnanski & Fiedler, 1970). The correlation function and integral length-scale based on lateral fluctuations (normal to the mean velocity) are still defined with a streamwise spatial shift and are given by: (x ) u (x, t)u (x x , t) u ,rms u ,rms A.84a (x , 0)dx 0 A.84b Unlike the streamwise correlations, the lateral correlation does not have an exponential decay because of an additional constraint associated with mass conservation. To demonstrate this constraint, let us denote the two lateral directions x 1 and x 2 . Next consider, the fluctuations associated with velocity in the first of these directions ( u 1 ) over the surface given by spatial 24 shifts in the other two directions as shown in Fig. A.14b. Since this plane is parallel to the flow, the spatially-averaged mass flux through this plane should be zero: u 1 (x, t) dx dx 2 0 A.85 This mass continuity criterion means that a positive u 1 on the plane should be counter-acted somewhere else with a negative u 1 . If this integral is weighted by u 1 at a streamwise shifted point on this plane and then averaged, the correlation integral over the plane will also be zero, u 1 (x, t)u 1 (x x , t) dx dx 2 0 A.86 For HIST, Csanady (1963) and Tennekes & Lumley (1972) refer to this as the “continuity” effect, and note that this property dictates a relationship between the streamwise and lateral spatial correlations: x (x ,0) (x ,0) (x ,0) A.87 2 x Employing the streamwise correlation given by Eq. A.83, the lateral spatial correlation for zero-time shift must then have the form (x , 0) 1 x exp( x ) A.88 This results in a substantial negative-loop, i.e. this correlation becomes negative for significant streamwise spatial shifts (Fig. A.14b). More importantly, the corresponding integral-scale in the lateral direction will be half that in the longitudinal direction, i.e. 1 2 A.89 Experiments in a variety of free shear flows generally support this ratio. For example, measurements in the center of mixing-layers correspond to 0.5 (Wygnanski & Fiedler, 1970) while those in circular jets correspond to 0.45 (Wygnanski & Fiedler, 1969). This ratio is reasonable even in NAsT flows such as that shown in Fig. A.13b. However, typical anisotropic free-shear and boundary-layer flows tend to have a length-scale ratio that is somewhat larger, about 2-3 (Hinze, 1975). The anisotropic aspect for the integral length-scale will be shown to be important for particle diffusion when the particle paths are no longer reasonably approximated by the fluid-tracer path. Based on arguments of dimensionality, the lateral (or streamwise) length-scale for HIST should be proportional to the Lagrangian time-scale and the turbulence intensity. One may define a constant of proportionality as c so that: c u L A.90 The parameter c has been termed the “turbulence structure parameter” since it relates all three integral scales. Free-shear measurements have shown that c is typically in the range of 0.5-3 (Hinze, 1975), though non-uniform value of c will appear near the edges of the shear- layers or regions where the turbulence is significantly non-homogeneous and anisotropic. One may also expect that the integral time scales reflects the balance between turbulence production and dissipation (Tennekes & Lumley, 1972) and this relationship can be quantified 25 by defining an integral time-scale coefficient which can be related to the lateral length-scale by Eqs. A.72 and A.90: L 3 c A.91 k 2 c k 3/ 2 The value of c is of order unity for most free-shear flows. For example, c0.2-0.3 is reasonable for the grid-generated wake (Snyder & Lumley, 1971) and for pipe flows away from the near-wall region (Oesterle & Zaichik, 2004). Similar values are reasonable for the outer portion (y/>0.15) of a turbulence boundary layer for which the integral time-scale is approximately constant (Hinze, 1975) as is the integral length-scale (Fig. A.13b). However, one can expect some variation with different flow geometries as well as some sensitivity to Reynolds numbers. In near-wall regions, the non-homogenous nature of the flow can lead to substantial spatial variations of c and c. For example, the inner portion (y/<0.15) of a turbulence boundary layer is typically characterized by changes in turbulence intensity (Fig. A.10), integral time- scale (Fig. A.13a) and length-scale (Fig. A.13b) whose combination indicates that c of Eq. A.90 must increase rapidly for locations closer to the wall. Similarly, the integral time-scale shown in Fig. A.13a is not linearly proportional to k/so that c not be considered a constant in this region either. Thus, care must be taken when using Eq. A.90 and A.91 for wall- bounded flows. A.5.4. Micro-Scales and Energy Cascade in Turbulent Flows Free-Shear Micro-scales The Kolmogorov scales are associated with the smallest flow structures in a turbulent flow whereby the velocity gradients are dissipated into heat through viscosity. The corresponding Kolmogorov length-scale () and velocity scale (u) are thus defined to have a Reynolds number of unity: Re=u/f=1. Using these characteristics, one can relate the Kolmogorov scales to viscosity and dissipation by 1/ 4 3 / f A.92a 1/ 4 u f A.92b Thus, the Kolmogorov length and velocity micro-scales act as a filter beyond which the turbulent energy spectrum dissipates due to viscosity. These scales in turn can be used to specify the Kolmogorov time-scale 1/ 2 f A.93 u Since these scales are not influenced by flow domain or large-scale structures, it is reasonable to expect that this micro-turbulence is approximately isotropic and homogenous even for NAsT conditions. This is Kolmogorov’s first similarity hypothesis and is applicable if the ratio of an 26 integral length scale to Kolmogorov length scale is large. Based on Eqs. A.91 and A.92a, the ratio of the length-scales is proportional to an integral scale Reynolds number (Re) u 4/ 3 ~ Re A.94 f Thus, the disparity between these two scales increases with increasing Reynolds number, which itself is much greater than unity for turbulent flow. Therefore, high Reynolds numbers will yield Kolmogorov velocity and length scales that are much smaller than the corresponding integral scales, which are typically much smaller than those of the overall flow evolution, i.e., u«u «uD and ««D. These length scales are shown schematically in Fig. A.15 for a particular flow where it should be noted that choice of the macroscopic domain scale is subjective but refers to a large-scale which controls the overall flowfield changes. It should be kept in mind that the integral and Kolmogorov scales will vary, i.e. they are generally not uniform throughout the domain. At high Reynolds numbers, there is a wide range of scales associated with turbulence as shown in Fig. A.5c and A.7e. This spectrum is often analyzed in terms of the wave-number (n) and the turbulent energy per wave number ( k , “specific turbulent energy”). These variables can be related to the total kinetic energy and a general wave-length as: n 2/l A.95a k k n dn A.95b 0 The conceptual relationship between these two parameters is shown in Fig. A.16. The largest scales (smallest wave-numbers) are determined by the macroscopic length-scale (D), while the most energetic structures (highest k ) typically occur at the integral length-scale (). The energy cascades down to the Kolmogorov length-scales, beyond which viscous dissipation effectively filters the turbulence. If the overall Reynolds number is high enough, the cascade of energy through successively smaller scales gives rise to an “inertial sub-range” which exists between the integral-scales and the micro-scales. This stems from Kolmogorov’s second similarity hypothesis whereby he assumed a range of wavelengths in the inertial sub-range of length-scales l such that «l«. Kolmogorov proposed that the inertial scales will be independent of the viscosity since «l, and that the specific turbulent energy will be dependent on the local wave-number and the overall turbulence dissipation if the cascade of energy is in equilibrium. Dimensional analysis can thus be used to relate the proportionality of these variables for both the kinetic energy per wave number and the kinetic energy at a given wavelength: k ~ 2/3n 5/3 ~ 2/3l 5/3 A.96a kl k / l ~ l 2/3 A.96b Thus, the turbulent energy content decreases as the wave-number becomes higher (wave-length becomes shorter), and the inertial sub-range gives rise to the classic -5/3 slope when considered on a log-log plot (Fig. A.16). While a formal inertial sub-range may require Reynolds number on the order of ReD~108 or more (Tennekes & Lumley, 1972), even moderate ReD experiments 27 are found to have a measured spectrum for which a portion has nearly a -5/3 slope, indicating an inertial sub-range. One may also define two particular turbulent length-scales in between the integral length- scale and the Kolmogorov length-scale. The specific kinetic energy weighted by the wave- number and normalized by the total kinetic energy yields an intermediate length-scale, linter (Spelt & Biesheuvel, 1997): 1 nk n dn/ k n dn A.97 l in ter 0 0 This may be thought of as the length-scale which has the mean (vs. peak) kinetic energy. Another length-scale of similar order is the Taylor micro-scale, lTaylor, which is defined by the mean turbulent velocity gradients of Eq. A.73, which can thus be related to Kolmogorov time- scale of Eq. A.93: u1 / x1 u 15f / 15 u 2 l Taylor u / A.98 The Taylor micro-scale can also be defined based on the curvature of the auto-correlation for small temporal shifts (Tennekes & Lumley, 1972). In HIST with moderate Reynolds numbers, the Taylor length-scale tends to be somewhat smaller (Poorte & Biesheuvel, 2002): l Taylor l int er A.99 Since these two length-scales involve turbulent intensity (associated with integral scales) and fluid viscosity (associated with the micro-scales), they will be bounded by the micro- and integral-scales for high Reynolds number HIST, e.g. «lTaylor « Unlike the large disparity in velocity and length scales noted between the macroscopic and microscopic conditions, the turbulent time-scales are often on the order of the macroscopic features and can even be larger. For example, turbulent pipe flows have integral time-scales that can be approximated as L 1.2D where D is based on pipe diameter and mean speed (Westerweel et al. 1996; Schneider & Merzkirch, 2001)Also, for the channel flow of Kulick et al. (1994), the integral time-scale is only about three times the Kolmogorov time-scale. Therefore, a time-scale comparison can be put forth as < < D. Wall-Bounded Micro-scales In wall-shear turbulence, the smallest scales are often characterized by the mean wall stress in the sub-layer region, which is related to the mean shear rate and kinematic viscosity. If y is the distance normal to the wall and the boundary surface is along the x-direction, the friction velocity, length and time scales can be given as: K xy,wall u f f u fr f , fr , fr A.100 f y wall u fr 2 u fr The non-dimensional velocity, length, and time friction scales are therefore given as 28 u y yu fr t tu 2 u , y = , t fr A.101 u fr fr f fr f The wall-shear rate is generally linear in the laminar sub-layer which extends to a y+ of about 5, so that in this region, the mean velocity simply equals the distance from the wall when both are normalized, i.e. u =y+. These “inner” scales are most appropriate close to the wall (e.g. y+<20), while the boundary layer thickness and other “outer” scales are most appropriate in the middle and outer edges of the boundary layer (e.g. y>0.15). The range between the inner and outer scales can be significant, e.g., + (≡yfr) can be in the hundreds or even thousands. A.5.5. Time-Averaged Navier-Stokes Equations and Techniques To understand and model the behavior of turbulent flow, it is often helpful to consider the time-averaged transport equations. If these equations are written in terms of flow variables which are individually time-averaged (Eq. A.67), the resulting equations are called the Reynolds-Averaged Navier-Stokes (RANS) equations. This approach is especially common for flows with constant density. For variable density flows (e.g. stratified flows or compressible flows), a modified time-averaged approach is used whereby most of the variables are density-weighting before they are time-averaged. This is referred to as Favre-averaging and reduces the complexity of the averaged equations. However, both the Reynolds-averaged and the Favre-averaged transport equations require “closure” models for solution. The objective of the closure models is to relate the correlations based on turbulent fluctuations (e.g., the kinetic energy defined in Eq. A.71) to gradients of only mean flow properties to allow a computationally efficient approach to turbulence. In the following, the Reynolds- and Favre- averaging concepts are discussed and respectively applied to the transport equations for constant-density flow and variable-density flow. Reynolds-Averaged Mass and Species Transport For the RANS approach, the unsteady velocity field is separated into steady and fluctuating components as noted in Eq. A.68, i.e. in tensor notation q i = qi q i A.102 The time-averaging can be applied to the various continuous-phase PDEs including conservation of mass, momentum, and energy. If we assume that the averaging is based on an infinitesimally small control volume which is fixed in time and space we note that the spatial and time derivatives of any property q can be written as q q A.103a x x q 0 A.103b t 29 For some flows, it is sometimes important to incorporate macroscopic unsteadiness of time- scale D while averaging over the higher-frequency turbulence of time-scales than and smaller. For example, an airfoil which has a time-dependent angle change which is much slower than that of the turbulent time-scales may satisfy the criterion: D . In order to separate these frequencies numerically, one may redefine the Reynolds-average as the time- average for a flow-field which has previously been subjected to a low-pass frequency filter. This effectively allows for macroscopic unsteadiness of the mean velocity and nullifies Eq. A.103b. However, the remainder of §A.5.5 will neglect macroscopic unsteadiness so that Eq. A.103b is applicable. The Reynolds-averaged transport equations of mass, momentum and turbulence quantities are discussed in the following. For these PDEs, the flow is assumed to have a uniform constant density. This assumption is consistent with §A.3 and §A.4 and is most common and simplest form of RANS. In this context, application of RANS to the instantaneous mass continuity equation (Eq. A.29) yields zero-divergence for the time-averaged velocity, and also (by Eq. A.102) zero-divergence for the fluctuation velocities: u u i 0 x i A.104a u u i 0 A.104b x i Because the resulting continuity equation for the mean flow is effectively the same as that for incompressible steady flow, they are often solved with similar numerical techniques, e.g. the pressure-based techniques discussed in §B.3.3. Next consider the time-average of the species transport based on Eq. A.50b, where the mean convective term can be broken into two terms: u i u i u f A.105 x i i t x i t x i x i x i Note that an unsteady term for the mean concentration (first term on the LHS) has been retained even though it can be set as zero by Eq. A.103b. Retaining this term is thus not needed for physical reasons, but is instead included to allow the PDE to be conveniently solved with a time-marching scheme. In particular, one may employ a “guessed” initial concentration field along with steady boundary conditions and integrate this PDE for long periods. In this limit, the unsteady effects will be swept out of the computational domain resulting in a converged steady-state solution. Further details of the “pseudo-unsteady” formulation will be discussed in §B.3.2. One may note that Eq. A.105 is similar to Eq. A.50b, which is its laminar flow counterpart equation. A key difference is that A.105 includes a term with gradients of u which i represents turbulent diffusion. Such diffusion is illustrated in Fig. A.17a, where eddies cause fluids to spatially “mix” in a vertically-averaged sense (but remain separate in a local sense since there is molecular diffusion). If the black fluid is defined by =1 and the white fluid by =0, the vertical spatial average of the distribution in Fig. 16a can be considered analogous to a time-average, i.e. =½ corresponds to 50% probability of having white or back fluid at a given vertical point or given point in time. 30 If there is no mean convection ( u x 0 ) and there is a net gradient in along the x-axis, black fluid that moves into the pure white fluid region is associated with a movement to the left, u 0 , and a positive perturbation to the concentration 0 . Similarly, white fluid that x moves into the pure black fluid region is associated with u 0 and 0 . As such both x types of mixing lead to u 0 , i.e. turbulent mixing is in the opposite direction of the mean x concentration gradient so that the gradient is reduced over time. The strength of this mixing can be expected to be proportional to the strength of the turbulence mixing and the magnitude of mean concentration gradient, i.e. is strongest where this gradient is largest (Fig. A.17b). A such, one may expect that turbulent mass diffusion can be reasonably described by Fickian diffusion (analogous to Eq. A.50b) yielding: u turb A.106 x i i This is a heuristic model where the turbulent transport is related to the mean concentration gradient and a turbulent diffusivity (turb). Similar to the relation of Eq. A.51 for molecular diffusivity, this diffusivity can be related to a turbulent viscosity (turb) and a turbulent Schmidt number (Scturb) as: turb turb / Scturb A.107 The turbulent Schmidt number is generally taken as 0.7 based on experiments (Faeth, 1987) while the turbulent viscosity can be expected to be proportional to turbulence intensity. Combining Eqs. A.105-A.107 yields the “pseudo-unsteady” turbulent mass transport: u i f turb A.108 t x i x i Sc Sc turb x i This expression shows that the turbulent and molecular diffusion are both proportional to concentration gradient and have a similar form. For most regions in turbulent flows, the turbulent diffusion typically dominates laminar diffusion, so that f can sometimes be ignored with respect to distribution. However, it should be noted turbulence that can not cause any local mixing of the instantaneous , which is instead controlled by molecular diffusion and the instantaneous gradients in (based on Eq. A.50b). As an example of the difference between diffusion in laminar flows and turbulent flows, consider the diffusion of miscible cream on top of coffee. Assume that both liquids have the same density and the initial condition is a horizontal interface, e.g. with the cream on the top and the coffee below (Fig. A.18). As shown on the LHS of Fig. A.18, if the flow in the cup is stirred at a low speed such that the flow remains laminar with no instabilities or eddies. In this case, molecular diffusion (associated with Θf) acts along the horizontal mixture region only which will require a long before the two liquids are homogeneously mixed. However, higher stirring speeds can cause the flow to be turbulent as shown on RHS of Fig. A.18. This creates a cascade of vortices through which the interface between the two liquids is rapidly stretched, convoluted and distributed throughout the domain. This mixing associated with Θturb causes the regions with only pure back or pure white regions to be broken up into many small such regions. As a result, molecular mixing only has to occur over very short distances to mix these regions and can operate simultaneously throughout the domain. Therefore, mixing in turbulent flow can be considered to be a two-stage process: 31 1) eddy-mixing whereby a wide range of vortex structures causes pockets of high and low concentration to be transported stretched, convoluted and distorted at integral scales, so as to yield a large amount of interfacial area with high concentration gradients, 2) molecular diffusion whereby the high concentration gradients and high interfacial areas lead to rapid local molecular mixing (associated with f) at small scales. As a result, the turbulent fluid will be fully mixed must faster than that for laminar flow where the molecular diffusion alone is limited to operate in a horizontal region initially with typically weaker concentration gradients. Favre-Averaged Mass and Species Transport As mentioned above, Favre-averaging employs a density-weighting to simplify the overall PDE equations when the density can no longer be considered a constant. The density-weighted time-average (i.e., the Favre-average) for a property q is defined as q q f A.109 f Using this notation, an arbitrary property may be decomposed in terms of mean and fluctuating components as q q q A.110 Favre averaging has the same properties given in Eq. A.103 and also leads to the identities: q f q / 0 A.111a qq A.111b The relationship between density-weighted average and the time-averaged can be obtained by combining the above Favre-averaging properties with Eq. A.102: q q q q q / f q q f A.112 Therefore, the difference between time-averaged quantities and Favre-averaged quantities is related to the density fluctuations. Applying the Reynolds- and Favre- averaging to the variable-density mass transport PDE (Eq. A.5a and A.44a), yields the following two PDEs (again retaining the unsteady terms for the pseudo-unsteady formulation as discussed with Eq. A.105: f f u j uj f t x x A.113a j j f f u j A.113b 0 t x j While both equations are simple, the Reynolds-average formulation for variable density includes a fluctuating correlation on the RHS which requires turbulence modeling. In comparison, the Favre-formulation for mean mass conservation does not include this term and thus is more straightforward to solve. Application of Favre-averaging to the concentration transport Eq. of A.50a yields: 32 f k t x i f k u i f ,k f k x i A.114 This is again simpler than its Reynolds-averaged counterpart, which would contain additional correlations based on density fluctuations. Such simplifications are the primary reason for the Favre-averaging popularity for stratified flows, where the fluid is considered incompressible but has variable density. Examples of such low-speed flows (with small Mach number) include a gas jet exhausting into a lower density gas or mixing of oceanic salinity currents. Additional details of the time-averaged equations for variable density flows with effects of compressibility are discussed by Rodi (1980) and Speziale et al. (1991). Time-Averaged Momentum Transport and the Assumption of Boussinesq For the momentum transport of Eqs. A.6 and A.8, the “pseudo-unsteady” formulations for constant-density Reynolds averaging and for variable-density Favre-averaging are given by: u f i f u i u j f g i p f i u i u j f uuj A.115a t x j x i x j x j x i x j f u i f u i u j u u j 2 p f g i f i 3 u k f uu j A.115b t x j x j x i x i x j i The RHS of these equations includes effects of gravitational forces, pressure gradients, mean viscous stress and a final term which is based on a correlation of the velocity fluctuations multiplied by the density, e.g. f uuj . This tensor is often referred to as the “turbulent stress” i or the “Reynolds stress” and has units of force per unit area. The diagonal velocity components of this tensor (i=j), contain the turbulent kinetic energy as shown by Eq. A.74b. The off-diagonal velocity components (i≠j) are responsible for turbulent transport much in the same way as u (of Eqs. A.105 & A.106) was responsible for turbulent mass transport. i The importance of the off-diagonal components of the Reynolds stress for transporting momentum can be illustrated by an incompressible 2-D shear flow layer with mean flow in the x-direction, as shown in Fig. A.19. In this schematic, one may note that an upward turbulent fluctuation (associated with u y >0) will tend to move low-speed fluid into the high-speed (where it will be associated with u x <0). Similarly, a downward fluctuation velocity will tend to entrain high-speed fluid from above and so will be associated with u y <0 and u x >0. Thus both scenarios transfer momentum via u u <0 and require a gradient in the mean velocity. x y The transfer is the reason for the downstream thickening of a boundary layer (Fig. A.8b), a shear-layer (Fig. A.5c), a jet, or a wake. Note that flows with no mean velocity gradient will not have this correlation, u u =0, and will also not have a mean turbulent transport of x y momentum There are two principle approaches to model the velocity correlation in terms of mean flowfield properties: turbulent viscosity models and Reynolds stress models. The first approach is the most common approach and employs the Boussinesq (1877) assumption 33 whereby turbulent momentum transport is related to turbulent viscosity and mean velocity gradients so that the time-averaged and Favre-averaged versions of the turbulent stress tensor are given by: u u j 2 f uuj K turb i k x j x i 3 f ij A.116a i ij u u j 2 u k 2 f uu K turb i k i j ij x j x i 3 ij x k 3 f ij A.116b This relationship employs the Kronecker delta of Eq. A.10 and is analogous to the diffusion of momentum in laminar flow (Eq. A.9). It is also similar to Eqs. A.106 and A.108 whereby turbulent mass transport is related to turbulent viscosity and mean concentration gradients. For algebraic turbulence models, the kinetic energy is not computed explicitly so that the last term on the RHS of these two formulations is often neglected. In this case, it can be seen that mean velocity gradients are needed to produce a turbulent Reynolds stress. Note that the more commonly used two-equation turbulence models include a transport PDE of kinetic energy so that the last RHS term of these two equations is retained. Note that the kinetic energy in the Favre-formulation is defined as: k 1 uu 2 i i A.117 This is the counterpart of the Reynolds-averaged kinetic energy given by Eq. A.71. Substituting the turbulent viscosity model into the incompressible RANS momentum equation of A.115a yields u u i u j p u u j 2 f i f f g i f turb i f kij x j x i 3 A.118 t x j x i x j Therefore, the Boussinesq assumption effectively assumes that the laminar and turbulent momentum transport have the same form. It is instructive to consider the limit of Eq. A.118 for isotropic turbulence with no mean velocity gradients so that the empirical Boussinesq assumption can be removed. While such a condition is not common (since mean velocity gradients are generally needed to produce turbulence), it highlights the role of the turbulent stresses. If one neglects the hydrostatic pressure gradient this yields p 2 k f A.119 x i 3 x i This result is similar to the balance of pressure and viscous stresses of Eq. A.63 for creeping flow, and indicates gradients in the kinetic energy (non-homogenous turbulence) will yield a mean pressure gradient. This equation also indicates that turbulence does not affect the fluid momentum for the idealized condition of HIST. Another interesting limit is that of very high Reynolds numbers. Since turbulence is often a result of essentially inviscid instabilities and since the turbulent viscosity will generally dominate the laminar viscosity, there is a temptation to remove f from Eq. A.118. However, viscous effects become on the order of convective effects in small localized regions, and these in turn influence the overall fluid dynamics. For example, wall-bounded flows generally yield high velocity gradients near the surface which, in turn, dictates the wall shear-stress magnitude and the primary production of vorticity. The local wall-stress also determines whether and 34 when boundary layer separation occurs over a body, which in turn can modify the overall pressure distributions and the surrounding flow-field. In a pipe flow, the Reynolds number determines the mean wall-shear skin friction (Fig. A.6) and thus axial pressure drop. Even if there is no wall-surface in the flow domain (e.g. a jet, a mixing layer, or a wake flow), the dissipation of vorticity and kinetic energy is directly controlled by viscosity. To employ the above equations, a model for the turbulent viscosity is needed to “close” the equations as discussed in the following section. However, it is important to note that Eqs. A.106 and A.116 are empirical approximations since the Boussinesq assumption does not theoretically arise from the Navier-Stokes equations. In fact, Tennekes & Lumley (1972) referred to the concept of eddy viscosity as the “gradient-transport fallacy” because turbulence generally involves complex development from one region to another. As such, models of turbulent viscosity are ultimately a “necessary evil” to provide closure to the time-averaged momentum equations. A.5.6. Techniques for Closing the Turbulence Equations For closure of the mean momentum equations, it can be seen that models are needed for the Reynolds stress given in Eq. A.116. There are two main types of models: 1) “turbulent viscosity” whereby models are employed for the terms on the RHS of Eq. A.116a and A.116b and 2) Reynolds stress models whereby the velocity correlation tensor is modeled directly. Turbulent viscosity models are more common, but both of these approaches require consideration of turbulence characteristics as discussed below. Turbulent viscosity models seek to relate the turbulent mixing characteristics to mean flowfield properties. This can be accomplished in several ways, the simplest of which is the mixing-length models. These are essentially algebraic relationships based on empirical coefficients and mean velocity gradients and include Prandtl, Cebeci-Smith, and Baldwin- Lomax formulations (Chung, 2002). However, these algebraic formulations assume that the turbulent shear stress is only a function of the local mean velocity gradients. This can be a severe limitation in many flows. For example, consider turbulence generated by a screen or grid in a channel of constant cross section. Away from the wall boundary layers, the mean velocity will be approximately uniform and the turbulence slowly decays as the flow moves downstream. However, the lack of mean velocity gradients is a problem for algebraic turbulence models since Eq. A.116 does not allow prediction of turbulence nor of its decay because such models do not include a transport PDE for the turbulent kinetic energy. Models which do have such turbulent transport PDEs (and so can handle grid-generated turbulent decay) include the so-called one-equation and two-equation turbulence models discussed in this section. A well-known one-equation transport formulation is the Spalart-Allmaras (1994) method. However, two-equation techniques are the most commonly used. This is particularly true for two-phase flows because they provide both a length and time-scale of the local turbulence, which are needed to model turbulent diffusion of particles as discussed in §3.5.3. The three most common two-equation approaches are the k- (which often gives good performance for free-shear flows), the k- formulations (which often gives good performance for boundary-layer flows), and the Shear-Stress Transport (SST) model which employs the k- approach close to the wall and the k- far from the wall, and uses a blending function in 35 between (Menter, 1993). All three approaches are similar and can provide information needed for turbulent length- and time-scales. Herein, we will discuss the k- approach (Launder & Spalding, 1972) as it is the simplest of the two-equation models and perhaps the most common for multiphase flows. For constant density flows, this approach approximates the turbulent viscosity in terms of the local kinetic energy (Eq. A.71) and dissipation (Eq. A.73) as follows: t f turb f c k 2 / A.120 This form introduces an empirical eddy-viscosity constant (c) generally taken as 0.09. Using Eqs. A.106 and A.107, the concentration–velocity correlation can be related to the turbulence properties and the mean concentration gradient as: c k 2 u turb A.121 x i Scturb x i i For variable-density flows, the turbulent viscosity is approximated similar to Eq. A.120: turb f turb cf k 2 / A.122 The variable-density version of Eq. A.116b similarly employs the Favre-averaged kinetic energy defined by Eq. A.117 and Favre-averaged dissipation defined by: u u f f f i i A.123 x j x j This is the mean dissipation per unit volume and is the counterpart of Eq. A.73. To complete the k- formulation, transport equations are needed for the kinetic energy and dissipation. To obtain a PDE for the kinetic energy for incompressible flow, one may take the dot product of ui and Eq. A.54 and time-average the result. This provides a Reynolds-averaged transport equation for the turbulent kinetic energy: k k u j k 1 1 ui uuj f 2 uujuj up t x i x i x i x i i i i IV A.124 I II III This equation represents the balance of four aspects of the turbulence. The LHS of this equation (I) represents the change in turbulent kinetic energy along the average continuous- phase path (i.e. mean convection) and includes the unsteady terms to allow a “pseudo- unsteady” formulation. The first term on the RHS (II) is the production term, whereby it can be seen that the turbulence generation occurs in regions of mean shear. The second term (III) is the diffusion of the turbulence (and includes sub-terms associated with molecular diffusion, turbulent diffusion, and pressure diffusion). The last term (IV) is the turbulent viscous dissipation as defined in Eq. A.73. The production term (associated with mean shear) is dominant at the lower frequencies while the dissipation term (associated with instantaneous shear) is dominant at the higher frequencies. A balance occurs at intermediate wave-numbers based on the inertial sub-range (Eq. A.96). The resulting cascade of turbulent energy is shown in Fig. A.16. Note that the degree of separation for D, and is a function of the domain geometry and domain Reynolds number. A similar transport equation can be obtained for dissipation using Eq. A.73. 36 Note that the introduction of a transport equation for the kinetic energy has introduced further turbulence correlations which in turn need to be closed. As such, further empirical relations are needed. A typical set of incompressible k- transport equations obtained with some terms modeled and some terms neglected (Shyy et al. 1997) is given by: k k K u j k ui f turb ij A.125a t x i f x i x i x i 1.44 u turb 2 ui j f 1.92 x i x i 1.3 x i K ij A.125b t x i f k k For the RHS terms of these two equations, the first term represents production and employs Eq. A.116a, the second term represents diffusion and the third term represents dissipation. Note that the pressure fluctuation term is often neglected since it is difficult and controversial to model. However, some formulations for this term are discussed by Chung (2002). The variable-density kinetic energy and dissipation PDEs can be similarly obtained starting with Eq. A.6a and assuming weak correlations between the viscosity and mean velocity fluctuations and neglect the pressure fluctuation term. Both of these assumptions are reasonable unless the flow is supersonic. Using the notation of Eq. A.116b, the PDEs become: u k K u f k p k u f f turb f i f j t x i x i x i x i x i ij i A.126a f f u i u j p 1.44 K u t x i x i x i ij i k 2 A.126b f f f turb 1.92f x i 1.3 x i k These equations include a new production term based on pressure gradient (second term on each RHS) which is related to pressure work. The pressure gradient is multiplied by a time- average of the Favre-perturbation which can be related to the density correlations using Eq. A.112 and then modeled by the density gradients using A.121: c k 2 f u ui u i u / f turb f f A.127 x i Scturb x i i i f Note that this term is identically zero for constant density flows. The remainder of each RHS includes diffusion terms (positive) and a dissipation term (negative). An important benefit of the k- formulation is that it can be used to specify turbulent length and time-scales which are useful in multiphase turbulent dispersion. For example, Gosman & Ioannides (1981) define the turbulent mixing length as: mix c3/4k3/21 turbc1/4k 1/2 A.128 This mixing length is proportional to the integral-length scale by a constant (c) which, in turn, can be approximated by the other integral-scale constants of Eqs. A.90 and A.91: mix c 5cc A.129 37 This constant should be on the order unity from dimensional arguments, and a typical value of c1.3 is obtained for the wake flow of Snyder & Lumley (1971) if one relates kinetic energy and dissipation to integral length-scales. Within a wall-bounded boundary layer, Eqs. A.126 and A.127 require modification and the mixing length in this region is related to the perpendicular distance to the wall (White, 1991): mix 0.41y wall 1 exp y / 25 wall A.129X In this expression 0.41 is the von Karman constant and 25 is a constant for flat-plate flow. Based on the above, the turbulent stress tensors of Eq. A.116 can now be written entirely in terms of the mean velocity, pressure, kinetic energy, and dissipation (as well as mean density and viscosity for stratified flows). As such, this represents closure of the continuity and momentum conservation. However, this closure of the equations is based on empirical coefficients (c, Scturb, 1.44, 1.92) that have been obtained by calibration with experimental data and are not universal. In general, they tend to allow a reasonable description of the mean velocity for free-shear flows and attached boundary layers but only qualitative descriptions of the k and , and are particularly troublesome in flow separation regions (Shyy et al. 1997). Furthermore, the k- formulation cannot predict the anisotropic aspects of turbulence, which are known to occur in many flows (e.g. near the wall Fig. A.10). To capture such behavior, Reynolds-stress models must be used instead. Reynolds stress models avoid the problems associated with defining a turbulent viscosity, i.e. the RHS models of Eq. A.116 are not needed. Instead transport equations are employed for all the components of the turbulent stress tensor. Since the velocity correlation is a symmetric tensor, this leads to six independent transport equations for the different components in a 3-D flow. The time-averaged PDE for these components can be given in terms of the mean transport (LHS) and production, diffusion and dissipation terms (RHS) in a “pseudo-unsteady” formulation: uuj uuj u i u j 1 u uj uk uju uu p i i i t x k x k x k x j x i k i k A.130 uuj u uj f x k i x k 1 uuju pu jk pujik 2 f i x k x k i k i Additional transport equations can be written for the tensor dissipation terms, and other terms such as the triple velocity correlation, the pressure terms, etc. However, this only creates additional correlations which further complicate the closure problem since empirical modeling is eventually needed for all the remaining turbulent correlations. Even without additional transport equations, Eq. A.130 requires more terms to be modeled (thus more empirical constants to be introduced and calibrated or estimated) than that for isotropic k- transport (e.g. Eq. A.125). The increased empiricism of the Reynolds-stress transport approach tends to hinder its robustness so that the predicted anisotropy may only be qualitative. As such, many researchers feel that the increased complexity of the Reynolds-stress model is not worth the extra computational expense of solving more transport equations. If computational expense is not a critical issue and a more accurate approach is desired, one can employ resolution of some or all of the turbulent features as discussed in the next section. 38 A.5.7. Fully-Resolved and Partially-Resolved Techniques To incorporate turbulent dispersion physics, improve velocity field predictions, and reduce empiricism associated with RANS modeling, one may choose to instead resolve all of the turbulent features and scales. This simply entails use of the un-averaged Navier-Stokes equations and is generally termed Direct Numerical Simulation (DNS), for which the resulting flow can be three-dimensional and unsteady (e.g. Fig. A.8a). The DNS approach can be used for laminar, transitional and turbulent flows. However, its computational cost can become excessive even at modest Reynolds number (§4.3.1 and Fig. 4.21). Therefore, intermediate approaches which a significant portion of the turbulent scales have become popular for engineering simulations. Such “partially-resolved” approaches aim to simulate the dynamics of the larger eddies which contain the majority of the turbulent energy, whereas the smaller fine scales are modeled is an empirical manner. The most common of these is the Large Eddy Simulation (LES) approach. A schematic comparison of DNS, LES and RANS is shown in Fig. A.16, where it can be seen that DNS resolves the full spectrum of turbulent scales while LES resolves the portion of the spectrum that has the majority of the kinetic energy. Finally, RANS does not resolve any of the spectrum, but simply represents the turbulence in terms of a single integral-scale mixing length. The LES technique and other partially-resolved techniques such as hybrid RANS/LES and Proper Orthogonal Decomposition treatments, are discussed in the next sub-sections. Large-Eddy Simulations Simulations which resolve some of the individual spatio-temporal features of the turbulent eddy structures will be considered herein as “partially-resolved”. This category primarily includes LES where the turbulence is only resolved up to some cut-off wave number of unfiltered solution (1/G), beyond which a sub-grid-scale turbulence model is employed. Note that in all the partially-averaged techniques, the emphasis is on resolving some form of the larger-scale eddy structures while those at the Kolmogorov scale are always filtered out. For LES, the governing equations are obtained by a low-pass spatial filtering of the Navier- Stokes equations such that all the velocity components are separated into their resolved (unfiltered) and unresolved (filtered) components. The spatial-averaged value of any flow variable (q) can be defined as a transfer from unfiltered space ( x ) to the filtered space (x): q x G x, x q x dx A.131 Here G is the filter function whose integral at any point is unity and can have Gaussian, Fourier, or Box functions. The latter is defined as a sharp cut-off in terms of the filter width (G): 1/ G if |x x| G / 2 G(x) A.132 0 if |x x>G / 2 A common prescription for the filter-width is to set it equal to the discretization length-scale (Garnier et al. 2002): G A.133 39 Ideally the spatial filtering is applied at sufficiently small-scales with high-order spatial discretization so that it occurs within the inertial sub-range so that the modeled turbulence is nearly homogeneous and isotropic. However, there may be practical limitations to increasing the resolution since solving the unsteady and three-dimensional LES flow is much more computationally intensive than solving a steady two-dimensional RANS flow. Use of the box filter is convenient as it yields u 0 so that the spatial-averaging is similar i to the time-averaging, i.e. ui u j ui u j uuj i A.134 If one interprets uuj as the unresolved (sub-grid) Reynolds stress, the LES governing i equations will be the same as that for the “pseudo-unsteady” formulation of RANS. For example, Eqs. A.113, A.115b, A.116b and A.126 can be considered as LES variable density transport PDEs. For many other LES filters, u 0 and u u , giving rise to additional terms i in the LES transport equations. However, these terms can be expressed in terms of the resolved flowfield so that they do not require empirical relations (Chung, 2002). The last terms on the RHS of Eq. A.115b become the sub-grid turbulent stress tensor for an LES formulation for incompressible or compressible flow. As with the similar term which arises in RANS, it can be modeled in a variety of ways. In free-shear flows, it is sometimes ignored altogether (e.g. MILES and VLES techniques) since the resolved-scales dominate and the effect of non-linear spatial schemes tend to mask its effect anyway (Yan et al. 2002). However, this term must be modeled for wall-bounded LES flows, and this is done with a Smagorinsky sub-grid stress model based on a sub-grid turbulent viscosity (): u u j ij uuj i u u A.135a x i x i 3 i k k 1 u i u j u i u j c G 2 2 A.135b 2 x j x i x j x i In these equations, G is the filter length-scale and c is the Smagorinsky sub-grid constant, with typical values ranging 0.2 (Piomelli, 1997) to 0.065 (Urbin & Knight (2001) for free-shear flows. This is effectively an algebraic model for the sub-grid stresses. To improve accuracy, the sub-grid kinetic energy (k) can be modeled with a one-equation transport like Eq. A.125a. In this case, the sub-grid turbulent viscosity can be obtained by assuming that the sub-grid mixing length is proportional to the filter-scale by a constant (Yoshizawa & Horiuti, 1985): ck k G 0.05 k G A.136 The turbulent dissipation can be similarly modeled as k 3/ 2 / G A.137 This indicates that all the turbulent dissipation for LES takes place at the sub-grid-scales. However, the total turbulent kinetic energy is the sum of the resolved and sub-grid quantities: k k k res A.138 The resolved kinetic energy can be obtained by taking a time-average of the unsteady resolved velocity fluctuations. 40 Near the wall-region, a RANS-type wall-damping function can be employed as a function of y+ (defined in Eq. A.101) so that the turbulent viscosity is zero in the sub-layer. This can be implemented by modifying the sub-grid stress constant, about 2 0.4y wall 3 c min 0.2, 1 exp y / 25 A.139 G This equation includes the von K rm n constant (0.4) and a damping constant (25) and reverts to the free-shear value far away from the wall (Balaras et al. 1996). More sophisticated sub- grid stress models are discussed by Piomelli (1997) and Chung (2002) including Favre- averaged formulations (to include density fluctuations) and dynamic sub-grid models (whereby c is a function of the resolved velocity gradients). However, Urbin & Knight (2001) showed that the effect of the sub-grid model is often weak and may even be neglected altogether when non-linear schemes are used to stabilize the temporal integration. This is because laminar viscosity is important in near-wall regions while outer regions are primarily controlled by inviscid instabilities. Hybrid RANS/LES Approaches Several hybrid RANS/LES approaches have emerged which intend to treat the separated flow regions with an LES approach while treating the attached flow regions with a RANS approach. Perhaps the most common of these is the Detached Eddy Simulation methodology (DES) developed by Spalart et al. (1997). The basic concept was to allow for a one-equation RANS treatment in the attached boundary layer regions (where the approach is reasonably robust and where LES would be prohibitively expensive computation-wise) and LES treatment in the separated and free-shear flow regions (where LES is reasonable and generally quite accurate). This blending is achieved by modifying the wall-distance to spatially separate the RANS and LES regions: ywall,DES min cDES, ywall A.140 Typical values for cDES are 0.65 or less. In near-wall regions, ywall,DES becomes the physical wall distance so will be governed by RANS-based turbulence model. In regions far from the wall, ywall,DES approaches a constant so that the turbulent equations yield a Smagorinsky-type viscosity. While the LES and DES equations share similarities, a key difference is that the LES grid resolution in the boundary layer region must be very fine to resolve boundary layer structures (e.g. Fig. A.8a), while DES grid resolution in the attached boundary layer must be intentionally coarser so that region will be described with a RANS approach (e.g. Fig. A.8b). An example of a DES computation is shown in Fig. A.20 for flow over a cylinder. Turbulent structures associated with the separated wake region (ywall»Δ) are three-dimensional and unsteady while the attached boundary layer region is RANS-like (steady and two-dimensional). One problem with the one-equation DES method is that it does not provides independent transport variables for turbulence intensity and turbulence dissipation. This can be a problem for multi-phase flows which require both properties (i.e. both a length-scale and a time-scale of the turbulence) and can be solved by using a two-equation hybrid model. Such a hybrid model has been developed by Nichols & Nelson (2003) based on the SST scheme, which was noted to be a popular and robust RANS technique in the previous section. Ideally, a hybrid turbulence 41 model produces both resolved and sub-grid (unresolved) turbulence which together mimic the actual turbulence, i.e. k = k res + k hy A.141 For RANS regions, the flow is steady (kres=0) and the hybrid model should represent the RANS turbulence (khy=kRANS). For LES regions, the flow is three-dimensional and unsteady so that most of the kinetic energy is resolved (kres»k) and the hybrid model represents the sub- grid turbulence (khy=k). For the Nichols-Nelson model, the conventional SST equations are used with the mean/resolved flow field to determine the RANS-scale turbulent properties of kRANS, RANS, RANS, and RANS. The large-scale length scale is based on the mean/resolved vorticity as: RANS max 6.0 RANS /, k 3/2 /ε RANS RANS A.142 This length-scale is used along with the sub-grid length-scale (, defined by Eq. A.139) to compute the hybrid turbulent kinetic energy: k hy 1 RANS 4/3 2 4/3 = 1 tanh 4/3 A.143 RANS 2 4/3 k RANS 2 The hybrid eddy viscosity can then be calculated using hy RANS k hy k RANS 1 k hy k RANS A.144a min c hy k hy , RANS A.144b The quantity chy determines the transition between RANS and LES behavior is similar to ck in Eq. A.136 and is taken as 0.0854 by Nichols & Nelson (2003). This hybrid turbulent viscosity (hy) is then applied to model the Reynolds stress for the mean/resolved mass and momentum flow equations. However, there is only one version of the turbulent dissipation (RANS=hy) since dissipation is assumed to take place at the smallest scales (res=0). For the RANS portion of the spectrum, this hybrid model behaves like a standard SST model and subsequently transitions to a non-linear k-equation model in the LES portion. Sample turbulence predictions for this hybrid approach are shown in Fig. A.21, which reasonably reproduces Eq. A.141. Proper Orthogonal Decomposition Another partially-averaged technique which is designed predict the primary of the spatio- temporal features is Proper Orthogonal Decomposition (POD). In this approach, the turbulence is decomposed into a set of fundamental modes which can evolve. These simulations employ a low-order construction of the turbulent flow-field typically using spectral or pseudo-spectral functions (Joia et al. 1997) which are tracked in time as 3-D dynamical features. As such, they tend to employ a more modest number of degrees of freedom (compared LES) while simulating the large-scale non-linear flow physics. Unfortunately, POD models typically require input from a detailed realization (experiment or simulation) in order to solve for the best fit of their lower-order dynamical system. Thus, they are empirical in the sense that they cannot quantitatively self-determine the continuous-flow vortex structures for a generic set of boundary and initial conditions. Furthermore, their prediction accuracy degrades 42 as the flow evolves for long times. However, once a POD is constructed for a particular flow it can be reasonably rendered many times to test the transport of a variety of particle conditions for moderate periods. 43 Appendix B: Single-Phase Flow Discretization Techniques Herein we consider computational methodologies for single-phase flow. These continuous- phase methodologies can be used directly for particle flows with one-way coupling but should be modified for two- and three-way coupling with appropriate sink and source terms as discussed in Chapters 7 and 8. We will limit our discussion of single-phase numerical methods to the most common and/or simplest approaches for the sake of brevity. The reader is referred to the texts focused on computational fluid dynamics for a more in depth discussion, e.g. Oran & Boris (1987), Hirsch (1988), Anderson (1995), Anderson et al. (1997), and Chung (2002). Throughout this appendix, there will be comments on numerical techniques appropriate for certain flow regimes (incompressible flow, compressible flow, turbulent flow, etc.) for which PDEs and equations of state were outlined in Appendix A. With respect to the reference frame of the continuous-phase, the primary choices include Eulerian, Lagrangian and the Arbitrary Lagrangian-Eulerian (ALE) approach. With respect to spatial discretization, the primary choices are local or global discretization. For local discretization, the domain is broken up into many smaller discrete volumes; for global discretization, the domain is decomposed into a range of wavelengths. The Eulerian local discretization approach is the most commonly used due to its simplicity and efficiency, and herein will be the primary focus. The associated grids and variable distributions for this approach are discussed in the next section. B.1. Computational Grids and Local Shape-Functions Grid Types and Shapes The organization of the computational nodes and volumes in an Eulerian reference frame for multi-dimensional problems can generally be categorized into three types of grids: structured grids, unstructured grids and hybrid grids (some regions structured and some regions unstructured as shown in Fig. B.1). Structured grids in their simplest form have edges that fall along physical or computational coordinate lines. For example, Cartesian meshes can have rectangular elements in 2-D (Fig. B.1) and hexahedron (“brick”) cells in 3-D (Fig. B.2a). The grid lines can also be curved to accommodate non-orthogonal computational domain boundaries. Structured grids allow direct indicial means of numerical communication along gridlines, i.e. predictable i,j,k indices along physical boundaries and principal directions. This allows significant numerical efficiency for many numerical methods in terms of CPU cost per cell. The simplified grid connectivity and geometry also allows the more straightforward Finite-Difference approach to be employed. Structured grids discretized along coordinate systems are associated with specific spatial resolutions. For 3-D Cartesian grids these are given by side lengths of a brick (x, y, and z) which are equal for an isotropic grid. Increasing grid resolution is often desirable but will increase the computational cost. For example, reducing grid scales in each direction by a factor of two will increase the number of elements in a 3-D domain by about eight-fold. Furthermore, the time-step size tends to scale with grid length-scale (to be discussed in §B.3.2) so that the number of time-steps for a given time will double indicating the total number of operations will increase by a factor of sixteen. As such, there is an important balance between 44 controlling computational resources for practical applications and ensuring independence on spatial and temporal resolution. This issue becomes paramount when selecting a numerical technique since each scheme has different grid resolution requirements for satisfying this independence. If the domain has regions which require much higher physical resolution than others (e.g. shock waves) or is geometrically complex, then structured grids can become inefficient because it is difficult to localize such resolution increases and there can be difficulty of grid generation for complex boundaries. Conversely, unstructured grids do not use grid lines and instead the domain is discretized one element at a time, typically starting from boundaries where grid resolution is specified. They can allow local resolution of features as they are not constrained by grid line connectivity. They are often composed of triangular elements in 2-D can be composed of tetrahedral elements in 3-D (or the other shapes in Fig. B.2) each of which can readily adapt to complex domains. Because of these attributes, unstructured grids are becoming increasingly important in computational fluid dynamics. However, typical unstructured grid elements are not as accurate in regions where grid aspect ratios are very large, e.g. boundary layers. In such cases, hybrid grids (Fig. B.1) which utilize regions of both structured and unstructured cells allow for the benefits of both types of grids. However, a downside of both hybrid and unstructured grids is that the lack of grid lines typically increases the computational time per node. Thus, they may be inappropriate for simple geometries for which there is no significant saving in the total number of elements compared to structured grids. Additionally, unstructured and hybrid grids do not allow implementation of third-order or higher spatial discretization schemes which can be beneficial for flows with a wide range of structures, e.g. direct simulation of turbulence. Therefore, there are many conditions for which a simple structured grid is appropriate. Shape Functions If a computational domain is discretized into a finite number of nodes, the flow variables at these nodes can be used to interpolate the flow variables at all other parts of the domain. Discrete computational methods use shape functions to describe the variations within the computational domain. Often, the shape functions can be defined “locally” (within a single element) as a summation: N q x = q i i xi , x B.1 i 1 In this equation, q represents the field variable, qi are the discrete values of this variable (unknowns), i are the shape (or trial) functions, and N is the number of nodes for a single element where the discrete values are stored. For example, N=2 for a linear one-dimensional element. In some cases, the discrete values are stored at cell centers so that N is the number of cells which border each other. If there is only one field property to be solved, e.g. stream function, then q is a single scalar field. However, in general the unknowns can represent multiple fields, i.e. density, three components of velocity and energy. To overview shape functions, it is helpful to first consider their general characteristics. For discrete volumes, shape functions should have the property of conservation and nodal identity 45 N B.2a i 1 i 1 for cell conservation i (x j ) ij for nodal identity B.2b In these equations, N is the number of shape functions (and nodes) for an individual element, xj are the node coordinates, and ij is the Kronecker delta defined in Eq. A.10. Many computational approaches assume linear shape functions as they are reasonably straightforward to implement. Figure B.3a shows a 1-D domain with two nodes per element. The grid spacing and a dimensionless coordinate distance (x*) may be defined from the local coordinates as: x* x 2 x / x 2 x1 x 2 x / x B.3 Applying, Eqs. B.2a and B.2b to this element yields two linear shape functions which happen to represent the non-dimensional sub-tended length defined by x and the opposing node: 1 x * B.4a 2 1 x * B.4b This result combined with Eq. B.1 recovers the nodal values at x*=0 and x*=1 and in between yields a linear interpolation, e.g. the pressure half-way between two nodes is given by N p= pi i p11 p 2 2 p1 p2 1 2 B.5 i 1 x*=1/2 To allow second-order polynomial variations within a cell, one may use quadratic elements. Such a 1-D element requires three nodes per element and is shown in Fig. B.3b. The quadratic shape functions are determined based on Eq. B.2 and the three local coordinates as: x* x 3 x / x 3 x1 B.6a 1 2x * 1 x * 1 , B.6b 2 4 1 x * x * B.6c 3 2x * 1 x * B.6d Additional nodes can be added for higher-order elements. An alternative method to increase the order of the element, which does not require adding nodes, is to store both the values and their derivatives at the nodes. For two-dimensional elements, the linear shape function for each node can be obtained by defining a subtended area based on the opposing node as shown in Fig. B.4a-b. The shape functions are then simply these areas normalized by the total element area: i =A i /A A * i B.7 Note that the element area of a triangle can be given in terms of the local nodal locations as 46 1 x1 y1 1 1 A 1 x 2 y2 x 2 y3 x 3 y 2 x1y3 x 3 y1 x1y 2 x 2 y1 B.8 2 2 1 x3 y3 The opposing triangular area for a given node can then be simply obtained by replacing the associated nodal coordinate with the interior nodal coordinate. For example, A2 12 xy3 x3 y x1y3 x3 y1 x1y xy1 B.9 From this the associated shape function id given by A2/A. For a quadrilateral element, A can be computed by dividing the element into two triangles based on two of the opposing nodes. The four opposing areas for Ai can be formed by using two interior lines which connect the mid-points of the opposing faces. An alternative approach which gives the same result is to map a quadrilateral as a square by employing the local coordinate system (x*,y*) by defining the local nodes as x*=±1 and y*=±1. For example, if Fig. B.4b was generalized to have non- orthogonal sides, the normalization length would vary linearly from y1-y4 on the left face to y2- y3 on the right face. Based on such a local coordinate system, the four shape functions are: i = 1 x*1 y* B.10 Shape functions which can be defined by the local coordinate system are sometimes called a “natural” or “local” element, and triangular elements can be similarly posed (Eq. B.4). Three-dimensional linear shape functions are similarly equal to the subtended volume associated with the opposing nodes (Fig. B.4c-d) so that i =i / * B.11 Similar to a triangle, the tetrahedral volume (and a subtended volume) can be given by the local nodal locations: 1 x1 y1 z1 11 x2 y2 z2 B.12 61 x3 y3 z3 1 x4 y4 z4 Brick elements, like quadrilateral elements, are generally mapped into a local coordinate system with x*=±1, y*=±1, z*=±1, so that i = 1 x*1 y*1 z* . If the coordinates are orthogonal, the computational volume is simply x y z B.13 Higher-order elements can also be formed with such multi-dimensional elements by using polynomial expressions which satisfy Eq. B.2 and additional nodes per element. Given the nodal values and the shape functions, any variable can thus be described throughout an element by Eq. B.1. Based on Eqs. B.3, B.7, B.12, or B.13, one may define length scales associated with various computational grids as 47 x for 1-D grids B.14a A 1/ 2 for 2-D grids B.14b 1/ 3 for 3-D grids B.14c This length-scale can be used to gauge the computational spatial resolution (which helps determine accuracy) and allowable time-step values (which are sensitive to grid resolution) as will be discussed in §B.3.2. Note that all the above shape functions can be considered “local” since they are valid over a small discrete volume of the overall computational domain () and each element contains N nodes and shape functions. To allow the local shape functions to be used throughout the computational domain, one may define a global shape function Dj associated with each global node xj, where j=1,Nf and Nf is the total number of nodes in the domain. A global shape function will equal to the local shape function (Dj=i) when x is within a “corresponding element”, i.e. an element which includes a local node xi used to define the local shape function. Outside of the corresponding element, the global shape function is zero. For linear shape functions in 1-D, this can be expressed as: 1 x i x /x for x i x x Di x i , x B.15 0 for x i x x This is illustrated in Fig. B.5 where it can be seen that the D2 is non-zero within element 1 and element 2 but zero within element 3. As such, the global shape functions have a tent-like shape which are unity at the associated node then decrease linearly to zero at the neighboring nodes beyond which they are zero (and this is also true for 2-D and 3-D linear global shape functions). Based on this definition, Eq. B.1 can be expanded to encompass the entire computational domain as: Nf q x = q i Di xi , x B.16 i 1 Since it is equally valid to consider a shape function in its local or global form, this choice is typically based on numerical convenience. There are other types of shape functions which are intrinsically “global”, i.e. are only defined over the entire computational domain (D). The next section discusses these and how global shape functions are used for different numerical methods. B.2. Weighted Residual Methods The central idea of the weighted residual method is to minimize a residual that relates the exact solution of the flow variables (qexact) to the discrete solution (q). The exact solution satisfies the governing Partial Differential Equations (PDEs) throughout the computational domain and all the Boundary Conditions (BCs) along the surface of the domain. The residual L is defined as the PDE operator which is zero for the exact solution at all points. For example, the incompressible PDE of Eq. A.54 can be written in terms of the flow variables and space as 48 Du L x, t f p f g f 2u B.17 Dt If one substitutes the exact solution, then L(qexact)=0. However, if one substitutes the numerical solution, then this residual will typically be non-zero, L(q)≠0, owing to (hopefully small) numerical errors. The weighted residual method examines the integral of the numerical residual throughout the domain and seeks to determine the unknowns which minimize this integral. This is accomplished by integrating the product of the residual and a set of “weighting functions”, wj, over the computational domain: w L q d j D 0 B.18 Note that wj are sometimes also called the “test functions” and generally j=1,Nf. The various weighted-residual techniques differ primarily in their choice of the weighting and shape functions as described in Table B.1. As noted in the previous section, a key classification for these functions is local vs. global. Techniques which use local shape functions (discussed in §B.1) will be discussed first, followed by those which use global shape functions. The three most commonly used local weighting functions are: collocation, finite volume, and Galerkin. The collocation method is the simplest since the weighting function is simply unity at the nodes and elsewhere zero (Table B.1). As such, the integral of Eq. B.18 becomes Nf L q qi Di , xi L q, x j 0 B.19 i 1 xi =x j This effectively forces the PDE to be specified as zero only at the nodal locations (j=1,Nf). At each nodal location, a spatial discretization of the PDE operator is typically carried out by considering the finite difference of the variables between neighboring nodes (as discussed in §B.3). As such, the collocation method is more commonly refereed to as the “Finite Difference” (FD) method. The spatial discretization results in Nf equations to be solved for convergence or for each time-step (since Nf is the number of nodes in the domain). Note that the FD method does not require knowledge of the shape function spatial distributions between nodes for solution of the qi values. Instead, these spatial distributions are only needed to interpolate values of q at non-nodal locations afterwards. For the Finite Volume (FV) technique, the weighting function is unity over each local computational cell volume (,k) and zero elsewhere, where there are N computational cells with index k. Thus, Eq. B.18 becomes an integral over each of these cells: Nf N L q qi Di , xi w jdD L qi i xi d,k 0 i1 B.20 i 1 This approach thus requires integration over each computational element for which the flow variables spatial distributions are governed by the local shape function distributions (to be discussed in §B.4). The shape function can be specified as constant, linear, quadratic, etc. depending on the order of spatial representation desired. The Galerkin method is similar but uses the shape functions themselves as the weighting function, so that the system of equations for all elements becomes N N jL qii d,k 0 j1 i1 B.21 49 This integrand includes products of shape functions as will be discussed in §B.5. Other weighting methods which use a local weighting include the least-squares method, where the weighting function is equal to the PDE operator. For intrinsically global shape functions, the computational domain is often broken into many wave-lengths based on Fourier series (sine waves), Chebyshev polynomials, Legendre polynomials, etc (Chung, 2002). This is typically called the Spectral Element method. The two most common weighting functions to use with this method are the Collocation and Galerkin versions (Table B.1). The shape function and weighting function influences the characteristics of the method and dictates its typical usage (Table B.2). The node-based finite difference approaches are the simplest to formulate, but may suffer from conservation errors if there are strong gradients (such as shock waves). This problem can be eliminated by using volume-based methods such as the FV or FE techniques. Spectral Element approaches allow high spatial accuracy for a large range of scales with high computational efficiency and, thus, are quite popular for simulating fully-resolved turbulent flows with DNS. However, spectral approaches can be complex to implement for complex geometries and can yield substantial errors when the flow field is discontinuous (i.e. due to shock waves, reaction fronts, field discontinuities, etc.) so that they are not as widely used as FD, FV and FE approaches. B.3. Finite Difference Methods B.3.1. Spatial Discretization The Finite Difference method sets the residual to zero at discrete locations L(xi)=0. To enforce this, the PDE derivatives must be written in terms of discrete quantities of the dependent and independent variables using finite discretization. In particular, the derivatives are expressed in terms of algebraic equations by the use of Taylor series expansions. For a 1-D example, consider a property q (representing velocity, pressure, etc.) defined at both points x and x+x, where x is a discrete distance which can be related to the derivatives expanded about point x, i.e. q x 2 2q x 3 3q x 4 4q q(x x) q(x) x ... B.22 x x 2! x 2 x 3! x 3 x 4! x 4 x This can be rearranged to represent the first derivative at point x as q q(x x) q(x) x 2q x 2 3q ... B.23 x x x 2! x 2 x 3! x 3 x If we combine the higher-order derivatives on the RHS as a discretization error of order x and use the notation q(x)=qi and q(x+x)=qi+1, the first derivative with the leading error term can be expressed as: q qi 1 qi x 2q q q ... i 1 i + O (x) B.24 x i x 2 x i 2 x 50 This finite difference equation is a “forward-difference” since it uses a position forward of x to obtain the gradient and is “first-order” accurate since the discretization error is linearly proportional to the grid length. The terms of Eq. B.23 which are neglected in the RHS of Eq. B.24 are considered the “truncation terms”, and indicate the scheme’s accuracy. Similarly, a first-order accurate backward difference-expression which instead uses q(x-x)=qi-1 is given by q q q i 1 i + O (x) B.25 x i x In both cases, the truncation error is linearly proportional to the discretization, i.e. first-order accurate. It should be kept in mind that there is another type of numerical error called “round-off” or “machine” error which is related to the finite amount of significant figures that a computer can use to represent a single number. For many computers, single-precision floating-point representation corresponds to approximately 7 digits while double-precision corresponds to approximately 16 digits. Such errors can become significant when subtracting terms which are nearly equal so that their difference approaches the precision of the numbers. However, this error is generally over-shadowed by the truncation term accuracy so that it does not manifest itself unless extremely small grid resolution is used so that the differences between two values as in Eq. B.25 also becomes small and subject to round-off error. In general, finite difference equations may be generated for any order derivative by using additional Taylor series and number of property locations. For example, the forward and backward Taylor series (Eq. B.23 and its counterpart) can be added to give a second-order accurate central difference scheme for the first derivative: q q i+1 q i 1 = O (x 2 ) B.26 x i 2x The differences between the computational stencils of Eq. B.24 and B.26 are shown in Fig. B.6a, where it can be seen that the stencil width has increased to 2x. In some case, it is helpful to have second-order accurate schemes which are not centered. In particular, one-sided schemes are helpful to impose gradient boundary conditions. For example, consider a boundary condition at node 1 given by (q/x)1=const. (e.g. a symmetry boundary condition is given when the constant is zero). In this case, the derivative should be related to nodal values at i=1, i=2, i=3, etc. For such conditions, one may derive one-sided schemes by combining a Taylor series expansions about node i for nodes i+1 and another for node i+2 to create a second-order accurate forward-difference scheme (and similarly for nodes i-1 and i-2 for a backward difference scheme): u 3u i 4u i+1 u i 2 = O (x 2 ) B.27a x i 2x u 3u i 4u i-1 u i 2 O (x 2 ) B.27b x i 2x The forward-difference schemes can be used to represent a gradient boundary condition at i=1, while the backward difference scheme can be used for a gradient boundary condition at i=N. An alternative second-order accurate approach is to define a “ghost node” just outside the 51 domain so that a centered scheme can be used, e.g. specify an additional node at N+1 and then discretize (q/x)N based on Eq. B.26. Note that fixed boundary conditions which simply specify a constant, e.g. q1=const., can be applied directly without any discretization error. A second-order accurate central-difference can be obtained for the second derivative by subtracting the forward and backward Taylor series (Eq. B.23 and its counterpart). Doing so, eliminates the first derivative terms and yields an expression involving nodes at i+1, i and i-1: 2 q q i+1 2q i q i 1 1 4q q 2qi qi 1 2 = x 2 ... i+1 O x 2 x i x 2 12 x i 4 x 2 B.28 Finite difference schemes for numerical implantation are obtained by replacing all the derivatives of the PDE with such expressions and ignoring the discretization errors (also called “truncation errors”). This, along with suitable boundary conditions, yields a set of simultaneous algebraic equations for the q nodal values throughout the computational domain. For both the first and second-derivatives, the order of accuracy can be increased further (e.g. to third-order or fourth order) by employing additional Taylor series terms and expansions about other points (e.g. qi+2, qi-2, etc.). A summary of many such finite differences is given by Anderson et al. (1997). Such higher-order schemes can allow substantial reductions in error or computational requirements. To demonstrate this, consider a one-dimensional domain of length D with uniform discretization so that the number of nodes is given by N f 1 D / x B.29 Therefore, the error of the scheme given in Eq. B.24 is proportional to N f1 i.e. doubling the number of grid points will tend to half the error. In contrast, the error for the scheme of Eq. B.26 is proportional to N f2 so that doubling the grid will lead to a four-fold reduction in error. One may also compare the number of nodes required for either scheme if they are specified to have equal error. In this case, the number of grid points for the first-order scheme will equal the square of the number of grid points for the second-order scheme. Thus, the second-order accurate scheme will require less grid points and this corresponds to a reduction in the amount of required computational resources both in terms of both memory and number of operations. While such benefits of higher-order schemes can be important, there are some conditions for which they are not helpful. For example, involving more nodal values to increase the accuracy also increases the computational stencil width which can increase the complexity and cost of the computational operation. In addition, regions with discontinuities (e.g. shock waves) can give non-physical oscillations with higher-order schemes. As such, many schemes are second-order which yields a balance between accuracy and complexity. Another way to increase accuracy is to have non-uniform mesh spacing when the gradients are non-uniform. The motivation for this is that the truncation error is proportional to grid size and spatial gradients (e.g. Eqs. B.24 and B.28). As such, error distribution can be balanced by using a smaller x where the gradients are high and a larger x where the gradients are low. In this way, the overall error throughout the domain is minimized for a fixed number of nodal values by increasing the grid resolution in the areas where significant changes are occurring. This non-uniform distribution can be imposed with a linear grid stretching factor (): x i 1 x i 1 x i (1 ò )x i (1 ò )(x i x i 1 ) B.30 52 With such stretching the second derivative finite-difference expression of Eq. B.28 will have a modified scheme and will include an additional truncation term: 2q 2q 2 2 ò q i 2 1 ò q i 1 x i2 4 q x i 3q i 1 ò ... x 2 i x i2 2 3ò ò 2 12 x 4 i 3 x 3 i B.31 This additional truncation term is negligible as →but can yield a first-order scheme if its absolute value is large. In practice, the stretching rate is limited to small values, e.g. || ≤0.15, to retain approximately second-order accuracy. For discretization in multiple dimensions, the above finite difference expressions can be linearly combined, e.g. the 2-D Laplacian at node i (in the x-direction) and j (in the y-direction) can be written as q i+1,j 2q i, j q i 1, j q i,j+1 2q i, j q i, j1 q 2 i,j = (x) 2 (y) 2 O (x 2 , y 2 ) B.32 This can be similarly extended to 3-D. Expressions can also be obtained for PDEs in cylindrical or spherical coordinates, and complex computational domains can employ grid discretization along transformed coordinate systems which conform to the domain geometry. In such cases, the discretization and solution occurs in the transformed coordinates (which need not be orthogonal, but should not be substantially skewed either), and then the results are converted into the physical (e.g. Cartesian coordinate system) using the transformation Jacobian. For very complex domains, it may be best to discretize with unstructured triangular or tetrahedral grids in which case the finite volume or finite element approaches are needed. B.3.2. Temporal Discretization Temporal discretization is critical for unsteady computational fluid dynamics for which the temporal integration should be accomplished with sufficient accuracy so that the solution is independent of the time-step employed. However, time-stepping is also used for many steady flows which use a “pseudo-unsteady” formulation by retaining an artificial unsteady terms. For example, the time-averaged turbulence PDEs of §A.5.5 include an unsteady term on the LHS. Such formulations start with an initial “guessed” solution and employ time-stepping with steady boundary conditions until a converged steady-state flow solution is reached, once the time-derivative terms becomes negligible. A similar approach can be used for the solution of other steady flows, e.g. a steady laminar flow can be similarly solved numerically by employing the unsteady PDEs of §A.1 (e.g., Eqs. A.5, A.6, etc.) with steady boundary conditions. For such “pseudo-unsteady” formulations, it is not critical that the temporal integration be accurate since any unsteady features will not be relevant to the converged solution. Instead it is more important that pseudo-time-marching allows rapid progression to convergence while avoiding numerical instabilities. This technique is useful for complex three-dimensional flows for which iteration to the exact solution can be impractical. In addition, pseudo-time-marching can make use of computational frame developed for truly unsteady flows. Time-marching schemes for both steady and unsteady flows can generally described as either “explicit” or “implicit”. The difference between these two types of schemes is based on 53 how a new time-step unknown at a given point is determined. Explicit schemes provide a solution for this new unknown independent of the solution of any of the other new unknowns at other grid points. In contrast, an implicit scheme relates the solution for the new unknown to other new unknown values at other grid points. In 1-D, the functional dependence of these two schemes can be expressed as: uin 1 F uin1 , uin , uin1 ,etc. for explicit schemes B.33a F uin1 , uin 1 , uin1 ,etc. F u in1, u in , u in1,etc. 1 1 for implicit schemes B.33b In the following, some simple explicit and implicit schemes are discussed. Explicit Schemes To consider various finite-difference expressions, let us consider Eq. A.55b under the assumption of one-dimensional flow with linearized inertia and with constant density and pressure so that the non-linear convection, gravitational and pressure terms can be neglected. The result can be expressed in terms of two simple PDEs if one alternatively neglects viscous or convective effects: u u 2u u f 2 B.34 t x x For any ODE or PDE discretization, it is important to evaluate all terms at a consistent time- level and a consistent spatial location. For the above equations, the simplest temporal integration scheme is to evaluate all terms at the current time-step and a nodal location: n u u 2u n n u f 2 B.35 t i x i x i The LHS time-derivatives require discretization which can be obtained with a forward- differencing in time (similar to that as used for the spatial discretization in Eq. B.23): n u u in 1 u in t 2 u u in 1 u in n 2 + O (t ) + O (t) 2 B.36 t i t 2 t i t In this equation, un+1 represents the solution at t+t and un at time t and the resulting scheme is first-order accurate in time. The RHS spatial derivatives of Eq. B.35 can be evaluated with second-order spatial (central) difference schemes about location xi: u t(u in1 u in1 ) f t(u in1 2u in u in1 ) u in 1 =u in O ( x 2 , t 2 ) B.37 2x (x) 2 The resulting one-step finite difference equation is “explicit”, i.e. un+1 can be obtained directly as a function of the previous unknowns, which is numerically convenient. Note that this scheme is first-order accurate in time with respect to acceleration (Eq. B.36) second-order in time with respect to velocity (Eq. B.37) because of the multiplication by t. 54 While this explicit scheme seems reasonable to use for solution of Eq. B.34, it has major problems related to numerical stability. The diffusion-based discretization of Eq. B.37 (associated with f) is only conditionally stable and thus the time-steps must be carefully chosen (to be discussed in the next sub-section). A more severe problem occurs with the convection-based portion (associated with u∞) in that it is unconditionally unstable and thus never appropriate! In the following, explicit convection schemes are considered which are conditionally stable and thus appropriate. The unconditional instability of the convection-based portion may be traced to a fatal flaw in its mathematical character. To see this, consider u∞>0. Only upstream information should influence u in 1 so that use of u in1 and u in is appropriate but use of the downstream value u in1 suggests that information is moving in the direction of -u∞, which is unphysical and this leads to instabilities. This problem may be conditionally solved by “upwinding” the spatial discretization. In particular, one may use Eq. B.24 or Eq. B.25 depending on the flow direction, such that only physically relevant upstream quantities are employed (for f=0): u in 1 =u in u (u in1 u in )t / x O (x, t 2 ) for u in 0 B.38a u in 1 =u in u (u in u in1 )t / x O (x, t 2 ) for u in 0 B.38b This first-order accurate upwinding scheme is “non-linear”, in that the scheme itself depends on the local flow solution. The spatial schemes can also use the higher-order discretization of Eqs. B.26 and B.27a while maintaining upwind behavior: u in 1 =u in 1 u (3u in 4u in1 u in 2 ) t / x O ( x 2 , t 2 ) 2 for u in 0 B.39a u in 1 =u in 1 u (3u in 4u in1 u in 2 )t / x O (x 2 , t 2 ) 2 for u in 0 B.39b This yields an explicit scheme which is second-order accurate in space and first-order accurate in time. One drawback of this method is that it employs an expanded computational stencil width which can make the scheme awkward to use near boundaries. A second-order accurate scheme which avoids the above increased stencil width is the Lax- Wendroff method, which is second-order accurate in time. This scheme differentiates the inviscid PDE of Eq. B.34 with respect to time to relate the second derivative in time to the second derivative in space u n n u 2 u u 2 u n n n 2 2 u u u 2 t t i t i t B.40 x i x t i x i The last term on the RHS comes about from re-substituting the unsteady term of the PDE for the convective term. This result can be used to replace the term which is second derivative in Eq. B.36 with a spatial derivative. Retaining the resulting term, shift the truncation error to a term associated with the derivative in time so that the scheme is second-order accurate in time: n u u n 1 u in t 2 2 u n i u 2 + O (t 2 ) B.41 t i t 2 x i 55 Replacing the new spatial derivative with a second-order accurate central difference and using central differencing for the convective spatial derivative (in Eq. B.37) yields the one-step inviscid Lax-Wendroff method: u t u t 2 n 1 u =u (u n 1 n i 1 u ) 1 (u in1 2u in u in1 ) O (x 2 , t 3 ) n i 1 B.42 x x i i 2 2 The increased time-accuracy associated with this Lax-Wendroff method (there is also a two- step version) importantly results in a scheme which is conditionally stable, and one that can be extended to include the viscous term discretization of Eq. B.37. Time-Step Limits The convection-based stability constraints for the upwind schemes (Eqs. B.38a-B.38b or Eqs. B.39a-B.39b) and the Lax-Wendroff scheme (Eq. B.42) can be obtained by examining the conditions for which a numerical error in the solution will decay (or grow) over a single time- step. For incompressible flow, linear instability analysis of the inviscid version of Eq. B.35 (f=0) yields the Courant–Friedrichs-Lewy (CFL) criterion: t conv x / u B.43 Such a constraint is equivalent to limiting the travel of information within a single time-step to the physical domain for which this information is available (consistent with one grid spacing). This constrain indicates that reduced spatial resolution or increased convection speeds will tend to reduce the numerical time-step allowed. If the convection becomes non-linear (u∞ is replaced by u), then this time-step may need to be reduced further in regions of large gradients. This is typically accomplished by multiplying the RHS of Eq. B.43 by a coefficient (cCFL) which can be adjusted between 0 and 1 and is typically set at the maximum value which ensures numerical stability. This constant is sometimes called the Courant or CFL number. If the inviscid PDE of Eq. B.34 is extended to include finite compressibility, information may travel due to both convection and acoustic signals. In this case, it can be shown that the Eigenvalues of this system yield wave speeds of u∞, u∞+a, and u∞-a, where a is the speed of sound. The CFL condition must be modified to account for the fastest possible wave-speed t conv x / ( u a) B.44 The inclusion of the acoustic speed indicates the importance of compressibility but leads to a more restrictive time-step constraint. This restriction becomes impractical for flow at very low Mach numbers ( u a ) since the flow will develop based on convection scales while the time-step will be restricted on the much finer acoustic scales. Therefore, compressibility is ignored in such cases (e.g. acoustic waves are neglected) and specialized incompressible schemes are employed (e.g. §B.3.3). The diffusion time-step limit for the non-convective portion of Eq. B.37 can also be obtained from stability analysis which yields: t visc ( x) 2 / 2 f B.45 This indicates that the viscous time-step is even more sensitive to grid spacing. 56 The techniques for the diffusion terms and convection terms of Eq. B.35 can be combined by satisfying the combined stability limits. For example, one may combine central differencing for the diffusive terms with upwinding for the convective terms to get a one-step explicit scheme which is second-order accurate in space and first-order accurate in time. The resulting scheme has a stability constraint which can be approximated as t min t visc , t conv B.46 However, this criterion may be overly conservative and also needs to be considered for both “local time-stepping” for steady-flows and “global time-stepping” for unsteady flows. For flow fields which are steady, “local time-stepping” is a common solution technique whereby an unsteady time-marching scheme is coupled with an initial guess for the flow conditions and steady boundary conditions is used to solve a “pseudo-unsteady” formulation noted in §A.5.5. Temporal integration over a long-time, barring numerical instability, will yield a flow field which converges to the steady-state solution since the q/t terms will eventually vanish. In this case, temporal accuracy is neither important nor desired. Instead, one wishes to arrive to the converged steady-state solution as soon as possible. To do so, each node in the domain can use its own time-step which is maximized for the stability. Anderson et al. (1997) proposed a local time-step constraint for 3-D flows which have convection, compressibility and viscous effects as: 1 Re min u x x, u y y, u z z B.47a 1 u uy u 1 1 1 t conv x z a x y z x y z 2 2 2 B.47b cCFL t conv i, j,k t i, j,k local time - stepping for steady flows 1 2 / Re i, j,k B.47c The last equation includes a Courant number (cCFL), which is nominally unity but can be reduced for non-linear flows to improve stability. This time constraint is local in the sense that it is applied at each node and at each time in order to accelerate the solution to convergence. Unsteady problems require that the time-step be universally synchronized for the entire domain in order to ensure temporal accuracy of the fluid dynamics. In this case, we must use a single “global” time-step for all nodes which satisfies stability constraints throughout the domain: t min t i, j,k global time-stepping for unsteady flows B.48 This global time-stepping can lead to quite small convection time-steps throughout to accommodate regions with small grid resolution (Eq. B.47a). In regions dominated by viscous effects, where the flow velocities and grid resolutions are small, the Re can become very small, which will reduce the overall time-step even further. This can lead to impractically small time-step, but this is typically avoided by using an implicit scheme for the diffusion term, as discussed in the next sub-section. 57 Implicit Schemes As mentioned above, time-step constraints on the diffusion terms with explicit schemes can be prohibitive. Therefore a common approach is to treat such terms with an implicit approach so that the scheme will be unconditionally stable. In this case, the time-step is limited only by accuracy and not by stability, and thus generally allows greater time-steps to be used. A popular implicit time-integration technique which is 2nd-order accurate and is unconditionally stable is the Crank-Nicolson method. For this method, the time derivative is considered at tn+1/2 and is discretized using a central-difference scheme: n 1/2 u i u in 1 u in O (t 2 ) B.49 t t Any spatial derivatives are then approximated as an average at the initial and final time-step values base on a Taylor series expansion in time (similar to that in space for Eq. B.22): 1 u u n 1/2 n 1 u n O ( x 2 ) B.50 x i 2 x i x i To apply the Crank-Nicolson method to non-convective version of Eq. B.35, the LHS is replaced with Eq. B.49 and the RHS is discretized with Eq. B.50, after which the spatial derivatives are approximated with central differencing in space (Eq. B.281) yielding: u in 1 u in f u in11 2u in 1 u in1 u in1 2u in u in1 1 O (x , t ) 2 2 B.51 t 2 x 2 x 2 This second order-accurate scheme and can be written in terms of the new velocity values as f t u in11 2u in 1 u in1 u in 1 h i u i 1 , u i , u i 1 n n n n 1 B.52 2 x 2 Unlike Eq. B.37, this scheme is implicit due to the multiple n+1 terms on the LHS. In general, this may be solved using an iterative approach whereby the unknown variables are continually updated until they are converged for a given timestep. However, the 1-D discretization of Eq. B.52 happens to yield a tri-diagonal system which can be solved directly using the form: a i u in11 bi u in 1 ci u in1 h in 1 for i=1,N B.53 The RHS is a known function of the values associated with the previous time-step as in Eq. B.33b. This linear system of equations can be solved efficiently with the Thomas algorithm which modifies the coefficients with two down-sweeps and one up-sweep. If there are no ghost nodes, a1=cN=0, this scheme can be expressed as 58 bimod bi a i ci 1 / bi 1 for i=2,N B.54a h imod h in a i d i 1 / bi 1 for i=2,N B.54b uN h mod i /b mod N B.54c u i h imod ci u i 1 / bimod for i=N-1,1 B.54d This form is convenient in that it avoids any iteration, though it is only directly applicable to 1- D discretizations. For 2-D and 3-D flows, iterative techniques are common because the system is no longer simply tri-diagonal since neighboring j (and k) nodes are coupled to the neighboring i nodes. However, the direct tri-diagonal scheme can be modified to handle multi-dimensional systems using the Alternating-Difference Implicit (ADI) approach. In the ADI approach, the differencing is first applied implicitly in one-direction (and thus solvable with the Thomas algorithm) and explicitly in all other directions. This procedure is then alternated among the other directions so that each of the directions is treated implicitly for one time-step and then the whole process is repeated. For example, consider the 2-D PDE given by u 2u 2u f 2 2 B.55 t x y The 2-D ADI approach used two alternating time-steps (the first is implicit in the x-direction and the second is implicit in the y-direction): n 1/2 u i u in 1 u in u in1, j 2u i,1 u in1, j u i, j1 2u i, j u i, j1 1 n 1 n n n f j B.56a t t x 2 y 2 n 3/2 u i u in 2 u in 1 u in1, j 2u i,1 u in1, j u i,21 2u i, 2 u i,21 1 n 1 n n n f j j j j B.56b t t x 2 y 2 The ADI technique approximates the Crank-Nicolson scheme in accuracy and is widely used for linear 2nd order derivative terms in space. However, iterative techniques can be applied to improve temporal accuracy and handle non-linear terms (e.g. if viscosity is not constant). Perhaps the simplest iterative technique is the Jacobi method which uses fixed-point iteration based on the values from the previous iteration. Application of this technique to Eq. B.55 evaluates the newest (k+1) iteration of an unknown based upon the surrounding values at the previous (k) iteration: f t n 1,k t n 1,k x 2 u i 1, j u in1,k f 2 u i,1 u i,1 h n 1, j y j n 1,k j u i,1,k 1 n B.57 j 2 1 f t / x 2 f t / y 2 In this equation, hn refers to an explicit RHS operator. This scheme can be iterated to a convergence tolerance (maximum net change of u between iterations) or an iteration count (kmax). The convergence can be accelerated by the Gauss-Seidel method which uses the most recently updated values (be they k or k+1) as one sweeps from i=1,Nx and j=1,Ny (where Nx 59 and Ny are the number of nodes in the x and y directions). This method can be coupled with successive over-relaxation (SOR) with a tunable parameter (cSOR) as: f t n 1,k t n 1,k 2x 2 u i 1, j u in1,k f 2 u i,1 u i,1 h n 1, j 2y j n 1,k j n 1,k 1 u i, j 1 cSOR u i,1,k n B.58 1 f t / x f t / y / cSOR 2 2 j If cSOR=1, this result is the standard Gauss-Seidel while cSOR=2 corresponds to an instability limit for linear systems. In the range 1<cSOR<2, changes will be accelerated compared to standard Gauss-Seidel and optimum values for convergence acceleration for linear systems are on the order of 1.5-1.8. Non-linear systems often require a lower value due to stability issues. Convergence of any of the Gauss-Seidel methods can be improved further by sweeping forwards and then backwards. For example, a sweep of Eq. B.58 with i=1,Nx and j=1,Ny can be alternated by a sweep for i=Nx,1 and j=Ny,1. The resulting scheme is called the Symmetric SOR Gauss-Seidel method and is quite popular for handling the viscous terms in many CFD approaches. Other commonly used iterative methods include the Conjugate Gradient and the Generalized Minimum Residual methods (Chung, 2002). Such methods are robust (nearly always converge) but are also more complex and require more memory. At the other extreme, the Gauss-Seidel SOR method is relatively simple and efficient and uses the minimum memory. The ADI method is intermediate to these two extremes, i.e. moderately robust and moderately expensive to employ. Other implicit schemes have also been developed for the terms associated with convection and pressure for incompressible flows, as will be discussed in §B.3.3. However, compressible flows typically use explicit schemes for the convection and pressure terms, as is discussed next. B.3.3. Incompressible Flow Implicit Schemes While explicit methods are generally favored for convection-dominated flows where the flow Mach numbers are significant, low-speed flows often use a different class of time- marching called “implicit methods”. In this case, the numerical schemes for a given space- time stencil are formulated as a function of multiple unknowns at the n+1 time-step. This is consistent with the elliptic character of the pressure field since disturbances at any location are felt all other locations instantaneously once a→∞. While this yields a more costly simultaneous evaluation of all the unknowns at n+1 time-level, it has the substantial advantage of increased numerical stability so that time-step constraints such as Eqs. B.43 and B.45 are no longer applicable. Most of the implicit schemes are unconditionally stable for linearized PDEs so that a time-step constraint is based more on convergence rates for steady-state problems. For unsteady incompressible flow problems, time-steps are typically constrained by Eq. B.47c in order to properly capture convection physics. However, this is still much larger than a time- step based on the acoustic speed or pseudo-compressibility. There are a great many implicit schemes which have been used including those of Laasonen, Crank-Nicolson, Beam-Warning, Briley-McDonald, and MacCormack (Chung, 2002). These schemes are generally second- order accurate in space and time but sometimes add fourth-order diffusivity to avoid dispersion errors for non-linear problems. Multiple dimensions can be addressed efficiently by using “Alternating Difference Implicit” (Eq. B.56) or the more general “Approximate Factorization” schemes. 60 A key issue for viscous incompressible flows is the proper treatment of the continuity though the solution of the pressure field. Some of common implicit schemes include Artificial Compressibility, Marker-and-Cell (MAC), Semi-Implicit Method for Pressure-Linked Equation (SIMPLE), SIMPLER (SIMPLE Revised), and Pressure Implicit with Splitting Operators (PISO). A common aspect of such incompressible flow techniques is the use of a staggered grid, whereby the fluid velocities are defined on a stencil that is shifted by a x/2 and y/2 from the pressure stencil. This combination avoids the unphysical checker-board type oscillations which can occur if a single unified grid is used for all variables. The MAC scheme is relatively straightforward to implement and is generally robust. The PISO scheme (Issa et al. 1986) is particularly popular since it allows large time-steps and avoids iterations so that it is generally efficient. Both of these schemes are overviewed in the following. Like most of these techniques, the PISO scheme uses a momentum equation which substitutes A.55a into Eq. A.55b, as well as a Poisson equation for pressure. The latter is obtained by taking the spatial gradient of the momentum equation, reapplying the continuity equation (A.55a), and assuming the temporal and spatial derivatives are independent. In vector notation, the PISO governing equations become: u i f u i u j p f f 2 u i - f g i B.59a t x j x i 2p u f u i u j 2p f i f 2 u i B.59b x i 2 t x i x i x j If hydrostatic forces are important, the fluid dynamic pressure (Eq. A.43) can be used in this equation so that the gravitational terms are implicitly included. The PISO scheme uses the following five steps to update the momentum and pressure fields, where only the temporal discretization is presented: uiu j pred1 f pred1 n p n t u i - u i f x j f u 2 pred1 i - x i f g i B.60a uiu j pred1 f u in 2 pred1 p 2 pred1 f f u i t x i x i x j B.60b uiu j pred1 f pred2 n p pred1 t u i - u i f x j f 2 u ipred1 - x i f g i B.60c uiu j pred 2 n 1 f u in 2 pred2 p 2 p2 pred 2 f f u i t x i x i x j B.60d f u i u j pred 2 f n+1 n p pred 2 t ui - ui x j f 2 u ipred2 - x i f g i B.60e Note that u associated with continuity equation appears as the first term on the RHS of Eqs. B.60b and B.60d. For increased accuracy and stability, the terms in the square brackets can be 61 split into diagonal and non-diagonal terms (Chung, 2002). The PISO method is well-suited for unsteady flow but also is quite reasonable for steady-flow, which is a reason for its popularity. However, the SIMPLE method (whose first three steps are similar) tends to be faster for steady-state problems for which the momentum coupling to the pressure variations is strong. Other techniques for incompressible flows include methods which successively employ a scheme for the convection and viscous terms followed by an implicit iterative scheme to include pressure corrections. This broad class of solvers includes the projection method and the Marker-and-Cell (MAC) method (Chung, 2002). For these methods, the momentum equation (e.g. Eq. B.59a) can be written in terms of convection and diffusion acceleration functions as: u 1 B.61a Fconv Fdiff g p t f Fconv uu B.61b Fdiff f 2u B.61c To discretize this equation, the projection method uses central-differencing in time on Eq. B.61a and then breaks up the change in velocity by using a predicted velocity (upred) and a corrected velocity (un+1): n 1/2 u u n 1 u n u n 1 u pred u pred u n t t O t 2 t t O t 2 B.62a u u pred n t n 1/2 n Fconv Fdiff1/2 g O t 2 B.62b u n 1 u pred t 1 f p n 1 O t 2 B.62c Thus, upred is the velocity at time tn that neglects the pressure gradient while un+1 includes this effect. Note that un and pn are assumed to satisfy both the above momentum equation and the incompressibility condition (Eq. A.29) at time tn and a similar statement can be made for un+1 and pn+1. For steady-flows where the time-derivative is non-physical and just used for convergence, the predictor step for the velocity can be accomplished with fast 1st-order explicit schemes. For unsteady problems, higher-order schemes are instead preferred since temporal accuracy is critical. A common approach for the non-linear convective portion is the 2nd-order accurate explicit Adams-Bashforth method. This method is based on a Taylor series expansion and a backward-difference for the temporal derivative about tn: t Fconv n F n F n 1 n 1/2 Fconv = Fconv n 2 t O t 2 Fconv conv conv O t 2 n 2 B.63 Since the viscous terms are linear, they can often be treated implicitly, e.g. with a Crank- Nicolson scheme (Eq. B.52) based on two Taylor series expansions about tn+1/2: n n Fdiff1 Fdiff t 2 2Fdiff1/2 n F n 1 Fdiff n n Fdiff1/2 = 2 2 t 2 ... diff 2 O t 2 B.64 62 Combining Eqs. B.62b, B.63 and B.64 yields a prediction of the uncorrected velocity: 2 n upred un 1 t 3Fconv Fconv Fdiff Fdiff tg n 1 pred n B.65 Note that the viscous portion is implicit (since Fdiff is based on upred) but can be solved with conventional iterative methods, such as the ADI scheme of Eq. B.56. For the correction step of B.62c, the pressure field is first determined to be consistent with condition of incompressibility, i.e. un 1 0 . To do this, one may take the divergence of Eq, B.62c to obtain a semi-discrete equation for pressure: f 2 pn 1 upred B.68 t This is an implicit equation for the pressure which requires a Poisson solver. Compared to the viscous terms (which include local diffusive effects), the pressure field requires more care for convergence as its influence is global (owing to the infinite acoustic speed for incompressible flow). Therefore, advanced iterative solvers are often employed for Eq. B.68 including pre- conditioning and multi-grid methods (Wesseleing & Oosterlee, 2001). Once the pressure field is obtained to within a suitable convergence level, the predicted velocity is corrected for the (previously missing) pressure gradient term using Eq. B.62c, i.e.: t n 1 u n 1 u pred p B.69 f The MAC method thus explicitly satisfied the momentum equation with a pressure that indirectly satisfies the continuity equation. A variant of the above technique is to use the old pressure gradient for the predictor step and then use just the change in pressure (from n to n+1) for the correction. For example, Eq. B.65 can be transformed as upred un 1 t 3Fconv Fconv Fdiff Fdiff tg pn / f 2 n n 1 pred n B.70 This scheme is second-order accurate for the convection and viscous terms but only first-order accurate for the pressure gradient. To achieve second-order accuracy for the pressure gradient and to satisfy continuity, a pressure increment can be defined through a Poisson equation based on Eq. B.68: p p n 1 p n B.71a 2 p f / t u pred B.71b u n 1 u pred t / f p B.71c Similar to the conventional MAC approach, the predicted velocity (Eq. B.70) is used to solve a Poisson equation (Eq. B.71b) with an implicit iterative scheme from which the corrected velocity can be computed (Eq. B.71c). This variant is similar to the SIMPLE method described earlier. The timestep for integration is generally given by Eqs. B.43, B.45, and B.46. Other incompressible flow techniques include the artificial compressibility method and vortex methods, though these are not as commonly used. One problem with all Finite Difference (both compressible and incompressible) methods for complex domain geometries (e.g. Fig. 1.13) is that it is difficult to discretize the domain for 63 the spatial derivative terms in Eq. B.88 with an adaptive 3-D coordinate system. This problem can be solved to some degree by overlapping grids (Chimera schemes) but these are limited to rather simple geometries. A more robust solution is to use unstructured grids, e.g. triangles in 2-D (Fig. B.1) and tetrahedron in 3-D (Fig. B.2); these require FV or FE approaches which are discussed in the next two sections. B.3.4. Compressible Flow Flux-Based Schemes Many of the schemes used today for compressible flow are designed to ensure conservation by employing flux-based methods. For examples, consider the 1-D convection equation for a variable q given by q uq B.72 t x For transport equations in conservative form, q can represent , u, or etot. A general scheme can be constructed using central-differencing in space and time as t uq i+1/2 uq i-1/2 O (x 2 , t 3 ) n+1/2 n+1/2 qin+1 qin B.73 x The terms in the parentheses are the “interface fluxes” between two cells. A common method to obtain these for uniform mesh spacing is to assume a linear interpolation between the cell- centered values: 2 qi+1/2 1 qi+1 qin+1/2 n+1/2 n+1/2 B.74a ui+1/2 1 n+1/2 2 u n+1/2 i+1 uin+1/2 B.74b In these equations, the values at the half time-step (n+1/2) can be obtained by first using a forward-differencing from the initial time-step values with t/2. Substituting the RHS expressions into Eq. B.73 yields a central-difference scheme, similar to Eq. B.26, but with conservation of the fluxes enforced directly. This approach can be similarly extended to multi-dimensional system of equations with the pressure terms, e.g. for the 3-D Euler equations (Eq. A.44) can be written as q Q x Q y Q z S B.75 t x y z The vectors of conserved quantities and fluxes for inviscid compressible flow are then: u x u y u z u u 2 p u u u u x x x y x z q u y , Q x u x u y , Q y u y p , Q z u y u z B.76 2 u z u x u z u y u z u z p 2 e u e p u e p tot x tot u y e tot p z tot 64 Here, Qx, Qy and Qz represent the flux vectors and S represents the source vector which can be due to viscous effects, etc. Extending Eq. B.73 to a two-dimensional scheme using central- differencing in time and space is thus given by t t Qx i+1/2, j Q x i-1/2, j Q y i, j+1/2 Q y i, j-1/2 n+1/2 n+1/2 n+1/2 n+1/2 qi, j qi, j n+1 n B.77 x y This can be extended to three-dimensional conditions and to include source terms. The fluxes at the half time-steps can be computed by using a variety of explicit and implicit schemes, as discussed by Chung (2002). Use of these fluxes in 2-D or 3-D requires definition of the cell volume so as to obtain the flux areas. For grids which are not simply rectangles in 2-D or “brick” elements in 3-D (Fig. B.2a), the finite-difference model is generally inadequate since careful integration along the cell surfaces is needed. In this case, finite volume or finite element discretizations (§B.4) are typically used. To ensure stability, one may not simply compute the fluxes with central differencing (recall the unconditionally unstable scheme realized with the inviscid version of Eq. B.37). One solution is to use the upwinding (similar to Eqs. B.39a and B.39b) for the convected quantities q. For example, the x-flux terms for the continuity equations can be computed with a 1st order upwind scheme as 1 i u x,i 1 u x,i for ui 0 B.78a u x i1/ 2 = i1/ 2 u x i1/ 2 = 12 2 i 1 u x,i 1 u x,i for ui 0 B.78b This is referred to as the donor-cell scheme and q denotes a convected quantity (Prosperetti et al. 2007). A more accurate solution is to compute these flux terms using inviscid solutions to the discrete local flow (Gudunov, 1959). The first-order accurate version of this technique assumes that the flow properties are piecewise constant within a cell so that there is a discontinuity at the interface given by the left and right states: qi+1/2,left qi+1 B.79a qi+1/2,right qi B.79b The flow will develop in time in accordance with an exact, but local solution to the inviscid Riemann problem. Figure B.6b shows an example where the adjoining cells have a pressure difference (but there is no velocity difference) which produces a local 1-D shock wave and contact discontinuity into the low-pressure cell and a local 1-D expansion wave into the high- pressure cell. A second-order accurate description can include construction based on piecewise linear variations of the conservative variables for which the discontinuities in pressure at the cell interfaces are generally smaller. In either case, the exact Riemann problem based on the interface differences can be obtained analytically for the general case where the velocities of the left and right states are not equal. The solutions for the four primary cases (shock-shock, shock-expansion, expansion-shock, and expansion-expansion) are outlined by Knight (2006) as a function of the “contact pressure”, which is an implicit function of the left and right states and must be obtained iteratively (e.g. by the Newton method). From this solution of the interface state condition after the interaction, the Riemann-based fluxes given by 65 and uq i-1/2 , can then be used to obtain the cell values at i, i+1, etc. at the next time uq i+1/2 n n (as in Eq. B.73). At the next time-step, a whole new set of local Riemann problems is solved for each set of adjoining cells, and so forth. This has the advantage of yielding the exact solution if there is a 1-D shock wave. However, there are some downsides: reconstruction can be expensive and use of three local 1-D Riemann reconstructions for a 3-D problem or a viscous problem introduces some errors. To address the first downside, there are two well- known groups of schemes: a) Roe’s scheme, which uses an approximate Riemann solution to avoid having to iterate to solution at each time-step and thus is a commonly employed variant of Gudunov’s scheme; and b) flux-vector splitting methods which decompose the fluxes on each side into two components, each associated with left and right running waves (Knight, 2006). While donor-cell schemes and Gudonov schemes can be stable, the above techniques will require special treatment if there are high gradients in the flow variables coupled with weak or negligible diffusion terms, e.g. shock waves in an inviscid flow. In such conditions, low-order upwind schemes (such as Eq. B.38) can be used for the flux differentiation of Eq. B.73. However, such schemes can be highly diffusive causing discontinuities to be “smeared” over an increasing number of cells. This diffusion is avoided with higher-order schemes (e.g., 2nd order Gudonov schemes or 2nd order upwind schemes based on Eq. B.77). However, these can lead to non-physical oscillations (termed “ringing” and illustrated in Fig. B.7) due to the absence or weakness of the stabilizing diffusion terms. Therefore, a practical solution for conditions with sharp gradients (especially when viscous effects or the effects of sources and sinks on the transport become important) is to combine low-order and high-order schemes in a non-linear manner, e.g. include numerical dissipation only in regions of sharp gradients. Such schemes can be represented with a slope limiter coefficient (clim) for selection of the convected variables: q i 1 clim,i 1/2 q i+1 q i 2 for u i+1 0 q i+1/2 B.80 1 q i+1 2 clim,i 1/2 q i+1 q i for u i+1 0 As such, clim=1 advances the scheme to a second-order central difference scheme with no limiting (Eq. B.77), while eliminating the mid-point terms (clim=0) reverts the scheme to a first- order upwind difference scheme (Eq. B.24 or B.25). To determine clim , one may use many different “limiting” schemes. Some of the more popular ones include: Essentially Non-Oscillatory (ENO), min-mod limiter, Monotone Upward-Centered Scheme for Conservation Laws (MUSCL), Total Variation Diminishing (TVD), and Flux Corrected Transport (FCT). These methods all employ low-order non- oscillatory schemes in high gradient regions while allowing high spatial resolution in other regions (Knight, 2006). For these schemes, clim is related to a smoothness coefficient: csmooth,i 1/ 2 qi qi-1 / qi+1 qi for u i+1 0 B.81 csmooth,i 1/ 2 q i+1 q / q i+2 i qi+1 for u i+1 0 This parameter is of order unity for smooth regions, e.g. regions where the gradient is approximately constant. Some of the limiter schemes are shown in Table B.3, and it can be seen that they yield cmin1 as csmooth1 for “smooth” conditions and cmin0 as csmooth0 for 66 regions with high variations in the gradients. For example, limiters for 2nd order schemes will be Total Variation Diminishing if 0 clim min(2,csmooth ) . Note that four of the limiters can have cmin2, which is effectively allowing anti-diffusion. This helps to create sharp discontinuities but can be related to instabilities in some cases. Application of B.81 and B.78a to Eq. B.76, thus yields t n+1/2 n+1/2 qin+1 qin x u i+1/2 qi+1/2 u n+1/2 i-1/2 qi-1/2 n+1/2 B.82 The Flux Corrected Transport scheme is a variant of these approaches and employs a linear combination of a Low-Order (LO) scheme (such as a 0th-order or 1st-order accurate scheme) with a High-Order (HO) scheme (such as a 2nd order accurate scheme) to compute the new value of q. The weighting of the contributions is defined to be the maximum HO contribution such that no new undershoots or overshoots are introduced in the variable q compared to the values at the previous time-step. This method can be summarized as t t n+1/2 n+1/2 qin+1 qi,LO n+1 clim uq HO uq LO clim uq HO uq LO B.83 x i 1/ 2 x i 1/ 2 In this approach, 0 clim 1 , and its value is determined by comparing the LO and HO solutions to q with the previous time-step values of q to identify the maximum value for mono- tonicity (Chung, 2002). Typically such methods allow physical discontinuities to be smoothly and monotonically resolved over 3-4 grid cells, whereas low-order schemes would spread the change over many more grid cells and would diffuse further in time (Fig. B.7). A detailed discussion and comparison of these methods is given by Knight (2006). For unsteady problems, one also prefers to employ higher-order temporal accuracy. Since the typical one-step schemes above are first-order (except for Lax-Wendroff which is limited to linear inviscid PDEs), most time-dependent numerical solutions use multi-step schemes. Perhaps the two most popular of these schemes are the MacCormack scheme and the Runge- Kutta scheme. The MacCormack (1969) scheme is a two-step scheme (predictor-corrector) which uses forward-differencing in space for the predicted variables and backward differencing for the corrected variables. It is particularly well-suited to non-linear problems written in conservative form, e.g. the compressible Navier-Stokes form in Eqs. A.5a and A.6a. For example, the 1-D ODE for an unknown conservative variable q and flux Q can be written q / t Q/x B.84 This form can be generalized to vectors q and Q corresponding to conservation equations for mass, momentum and energy. The corresponding MacCormack predictor-corrector scheme (with coordinate node index i) is then q ipred =q in (Qin1 Qin )t / x B.85a q in 1 = q in q ipred (Qipred Qipred )t / x / 2 1 B.85b This combination of averaging turns out to be second-order accurate in time and space and any non-linear derivatives in Q can be decomposed into derivatives of each variable, i.e. if Q≡q1q2+q3 then the forward differencing operation becomes 67 n Qin1 Qin =q1,i q n 1 q n q 2,i q1,i 1 q1,i q3,i 1 q3,i 2,i 2,i n n n n n B.86 A similar equation can be written for the backward differencing of the predictor method. For multi-dimensional problems, forward-differencing is used for the first-step and backward- differencing is used for the spatial derivatives of the corrector step. For even higher-order accuracy, the four-step Runge-Kutta scheme can be used as: q ipred1 =q in Qin t / 2x B.87a q ipred2 =q in Qipred1 t / 2x B.87b B.87c q ipred3 =q in Qipred2 t / x B.87d q in 1 =q in Qin 2Qipred1 2Qipred2 Qipred3 t / 6x B.87e Qi Qi 1/ 2 Qi 1/ 2 Qi 1 Qi 1 / 2 B.87f This scheme uses three intermediate predictions of the variable q (noted by superscripts pred1, pred2 and pred3) before finally advancing to the next time-step (n+1). The first two steps can be seen as forward-difference operators to obtain improving mid-point estimates, the third step obtains an estimate of the final value, and the fourth step combines these expressions together for a fourth-order accurate prediction. The approximation used for the Runge-Kutta spatial flux difference on the RHS of Eq. B.87f assumes a simple linear variation. A similar assumption is made for the MacCormack and Lax-Wendroff schemes. However, this scheme requires approximately two-fold the storage of the MacCormack and other two-step schemes. B.4. Finite Volume Method While spatial discretization of the differential form is called the Finite Difference (FD) method, the discretization of the integral form is called the Finite Volume (FV) method. As noted in Table B.2, FV (or control volume) methods are formulated from the inner product of the governing PDE with a unit function integrated over a discrete volume. This process results in spatial integration of the governing equations. Using discrete control volumes within the domain allows the use of unstructured grids with triangular and tetrahedral elements, which are especially useful for grid adaptivity and complex domain geometries. Because of this, CFD codes have shifted their focus from FD to FV methods over the last decade or so. Another important attribute of FV methods is that conservation (of mass, momentum and energy) can be ensured within each computational cell and the computational domain. This is especially important for flows with contact discontinuities (e.g. shocks, slip layers, and scalar or reaction fronts), for which some finite difference schemes can have significant conservation errors. This is because FD methods only satisfy the differential form at a point instead of the integral for a volume. To illustrate the FV method, consider the conservative form of the compressible Euler equations (Eq. A.44) as a vector of PDEs: q Q x Q y Q z = B.88 t x y z 68 The conservative unknown variables (q) and the fluxes (Qx, Qy, Qz) are given by u x u y u z u p u 2 u u u u x x x y x z q = u y , Q x = u y u x , Q y = p u y 2 , Q z = u y u z B.89 u z u z u x u z u y p u z 2 e e u pu e u pu tot tot x x e tot u y pu y tot z z This conservative form can be maintained even with inclusion of viscous terms (Eqs. A.5a, A.6a & A.15). Writing the PDEs in finite-volume form for a discrete volume (Δ) yields: q Q Q y Q t xx y zz d 0 B.90 If we define a volume-average of the conservative variables as q and employ second-order central-differencing in time, the PDE can be expressed as: n 1/2 t Q x Q y Q z x n 1 n 1/2 q q q n d B.91 y z This scheme can be applied to either cell-centered or node-centered (also called vertex- centered) schemes depending on the choice of the discrete volume (Fig. B.8). Cell-centered schemes are more common since their volumes are more straightforward to define and their unknown variables are stored at either the cell centers (Fig. B.8b) or at the nodes themselves. To illustrate implementation, consider a cell-centered scheme which stores the unknowns at the nodal locations 1,2,3 for two-dimensional triangular cells (Fig. B.4a). The RHS integral of Eq. B.91 becomes a volume per unit width integral over the cell area (A) which can be converted into a line integral over the cell edge lengths using Green’s theorem Q x Q y dA Q x dy Q y dx B.92 x y Assuming a linear variation of the fluxes along the cell edges and using area of Eq. B.8, t Q x 2 Q x1 Q Q Q Q q y2 y1 x3 x 2 y3 y2 x1 x3 y1 y3 2 2 2 B.93 t Q y2 Q y1 Q Q Q Q - x 2 x1 y3 y2 x 3 x 2 y1 y3 x1 x 3 2 2 2 Each cell-averaged difference for a cell (k) is then extrapolated to the cell’s nodal locations (i) using a discrete version of Eq. B.88 (Chung, 2002). Then the contributions from the K cells surrounding a vertex are averaged to find the net nodal correction and the new value 1 K qin 1 = qin qi,k K k 1 B.94 To improve stability and accuracy, a second (corrector) step is generally added, e.g. a four-step Runge-Kutta temporal integration. Another popular approach is to use temporal integration to predict qn+½ as a first step, and use the resulting Qn+½ in the RHS of Eq. B.91 to get qn+1 with 69 second-order accuracy (this is effectively a two-step Runge-Kutta method). For supersonic (highly compressible) flows, upwinding of the convective terms is often important for stability. Implicit FD methods (§B.3.2) can also be used for FV methods to further increase stability, especially for incompressible flows. B.5. Finite Element Method For the Galerkin formulation of the finite element method (commonly used in CFD), the shape functions are used to weight the residual (Eq. B.21). However, some care must be taken with respect to temporal integration as is discussed in the following. For comparison with the FV method, the FE method is applied to the PDEs of Eq. B.88. If the shape function is a function of space only (and not time), the resulting FE form for a single computational element and forward-difference time marching yields: Q x Q y n q q jdA t x y jdA n 1 n B.95 For computational efficiency, it is convenient to assume the same shape functions also apply to the fluxes, i.e. Q qi i iQ qi iQi B.96 Applying this expression and substituting the discrete form of the variables (Eq. B.1), yields i n dA q i j n 1 i qin t i x jdA Q x,i Q y,i y n B.97 For a single element, i and j range from 1 to N (the number of nodes per element), thus yielding a set of N for each element and for each variable in the vector q. The first term in parentheses on the LHS of Eq. B.97 is the called the mass matrix, which does not appear in the FD or FV methods (see Eqs. B.37 & B.91). The mass matrix causes the RHS to be coupled, but this is sometimes eliminated by using a “lumped” mass matrix whereby the area integral is simply replaced by A. This decouples the qn+1 variables so that the result is fully explicit and thus faster computationally. While lumping reduces the temporal accuracy, it is reasonable for local-times for which only a steady-state solution is required. For unsteady problems, an iterative procedure can be used for Eq. B.97 using successive improvements on qn+1 with a (low-cost) lumped mass matrix for increments. Typically three to four iterations are suitable for convergence with second-order accurate temporal schemes. The term in brackets on the RHS of Eq. B.97 requires differentiation of the shape functions and a non-unique solution requires them to be linear or higher-order. For a linear 1-D element in a local coordinate system with nodes 1 and 2 (N=2), Eq. B.4 yields 1/x*=1 and 2/x*=-1 within the element. Note that adjoining elements which share a node will also share the q at the node, which must satisfy the FE equations for both elements. To demonstrate this coupling with Eq. B.97, consider a three element (four-node) discretization of a 1-D domain with the linearized advection equation: u/t=u∞u/x. The resulting set of equations (one for each unknown/node in the domain) is obtained by global assembly as 70 2 1 0 0 u1 1 u1 n n 1 1 0 0 u1 n n 1 n 0 u n x 1 2 1 0 u 2 u 2 u t 1 0 1 2 B.98 6 0 1 2 1 u 3 1 u 3 n n 2 0 1 0 1 u3 n n 1 n 0 0 1 2 u 4 u 4 0 0 1 1 u n 4 With mass lumping of the LHS matrix (terms in a row are summed to replace the diagonal and other terms are reset to zero), a general interior node equation reduces to u in 1 u in (u t / 2x)(u in1 u in1 ) B.99 Comparing this with the inviscid version of Eq. B.37, it can be seen that this one-step forward- difference, central-space discretization yields the same scheme as when coupled with the finite- difference method (and the FV method), and thus also the same problem of numerical stability. To avoid these, unconditional stability (at least for linear systems) may be obtained with implicit time-stepping, e.g. the first-order accurate backward differencing scheme or the second-order iterative Newton-Raphson scheme. Otherwise, the same type of explicit time- stepping solutions used for FD and FV are also appropriate for FE, e.g. use of a Lax-Wendroff scheme (termed “Taylor-Galerkin” when used with an FE approach) or upwinding when viscous stress terms are also important. The inclusion of viscous stress terms in the FE approach with linear elements requires that the discrete form of the equations be modified. This is necessary because direct application of Eq. B.97 would result in terms such as 2/x2 which become zero. This problem can be solved with higher-order elements, but is more generally addressed by employing integration by parts. For example, the RHS terms with the second derivative of velocity are revised as: 2u u j u x 2 jd x x d x jdA B.100 i j i u i d jdA x x x Thus, any second derivative in the control volume is reduced to a first derivative so that a linear shape function can be employed. Furthermore, the boundary surface integral can be used to explicitly impose boundary conditions. Because of this, integration by parts is sometimes used on all the RHS terms. A similar result can be obtained through the variational principle approach (Chung, 2002). In either case, the resulting system of equations can be then be solved explicitly or implicitly. B.6. Spectral Element Method Spectral Element methods differ from the above FD, FV and FE methods in that they use smooth basis functions which are global, i.e. finite over the entire domain instead of just over a small discretized region. Spectral approximations differ primarily based on the choice of the basis functions and the locations over which the PDE is evaluated. 71 The most common basis functions are the Fourier and Chebyshev expansions, for which fast Fourier transforms and real cosine transforms exist. Fourier basis functions, also called spectral series, are simply the sum of linear sinusoidal functions (Fig. B.9) Nf q= ci (t) sin ix * B.101 i 1 The spatial resolution is increased by increasing the number of terms in the series (Nf). Note that this form differs from other weighted residual methods defined by Eq. B.1 in that spectral methods solve for basis function coefficients rather than for the values at nodal locations. For complex geometries or domains where harmonic series can not be readily applied (e.g. a flat wall extending to infinity), Chebyshev polynomials are preferred (Orszag, 1980). The terms are based on orthogonal cosine terms in Fourier space, so that the first terms are given by q=c1 c2 x* c3 2(x* )2 1 c4 4(x* )3 3x* c5 8(x* ) 4 3(x* ) 2 1 ... B.102 For an unsteady problem the polynomial coefficients are solved as a function of time. Because of their simple and global description for the basis functions, spectral methods are best for problems with either square or circular boundaries. When the boundaries have other shapes, a conformal mapping procedure can be used. The above representations are generally employed in two forms: fully-spectral and pseudo- spectral. The fully spectral method is related to the Galerkin form whereby the PDE is weighted by the basis functions themselves over the entire domain Nf sin jx* L ci sin ix * d 0 B.103 i 1 This is similar to Eq. B.21 for discrete methods. With either Fourier or Chebyshev basis functions, this technique can take advantage of Fast Fourier Transforms to solve for the system of coefficients. The pseudo-spectral method uses collocation at particular points, i.e. Nf L ci sin ix* 0 B.104 i 1 x =x j This is similar to Eq. B.19, and the locations employed are equally spaced if Fourier series are employed (or are the roots of the polynomials if Chebyshev polynomials are employed). In comparison to the spectral method, the pseudo-spectral method can be very efficient (twice as fast as the fully-spectral method) since direct integration is not required. However, it may suffer from aliasing or conservation errors if not independently filtered. B.7. Lattice Boltzman Method The Lattice Boltzmann (LB) method is a powerful technique for the computational modeling of a wide variety of complex fluid flow problems. Like FD, FV and FE methods, it uses a discrete domain. However, discretization is not based on macroscopic continuum equations but is instead based on microscopic models. In particular, it considers a collection of “fluid particles” each with a velocity distribution in specific directions. The fluid particles move in discrete time and can collide with each other under prescribed interaction forces. The 72 rules governing the collisions are designed such that the time-average motion of the fluid particles is consistent with the Navier-Stokes equation. The collision principle allows boundary conditions for interfacial dynamics and complex boundaries to be implemented in a straightforward approach. The LB space is discretized consistent with a kinetic equation for the velocity of each fluid particle, i.e. the coordinates of the neighboring points around x are x+ui, where ui are the discrete velocities. As shown in Fig. B.10, a 1-D lattice for a given volume has only three directions (+1, 0, and -1), but a 2-D lattice can have 5 or 9 directions, while 27 directions are possible for 3-D meshes. Each direction is associated with a certain fluid particle mass/volume (i) which can be summed to give the macroscopic density, momentum and kinetic energy for a given mesh volume as N = i B.105a i=1 N u = iui B.105b i=1 3R g T N (ui u) 2 = i 2 i=1 2 B.105c In these equations, N is the number of directions of particle velocity directions at each node, e.g. N=3 for 1-D. To obtain the above macroscopic properties, the LB method uses a convection/collision equation for the fluid particles at each node and in each direction, which can be written as: i (x ui t, t t) i ( x, t) F,i ( x, t) for i=1,N B.106 Therefore a set of transport equations are solved along discrete paths as compared to conventional finite difference equation for which there is only one transport equation for density. Thus the LB technique has more differential equations, but each is limited to specific communication paths with other nodes, allowing computations to be done independently for each path. For the RHS, F,i is the collision operator which represents the rate of change of each fluid particle density resulting from collisions. It may be expressed as F,i = i ieq / B.107 In this expression, ieq is local equilibrium distribution function and is the equilibrium relation time-scale to be consistent with pressure and shear stresses. The collision operator is required to satisfy conservation of total mass and total momentum at each lattice: N N F,i 0, i=1 F i=1 ,i ui 0. B.108 In 3-D flows, the collision process combined with the high number of components per cell volume make this method resource intensive (memory and CPU) compared to other discrete techniques. However, since it only requires nearest neighbor information it avoids the Poisson equation for pressure (generally needed for incompressible flows) making it highly 73 parallelizable. More work is needed to make it robust for highly compressible flows or flows with large density ratios. B.8. Direct Simulation Monte-Carlo Method For flow conditions where the Knudsen number (Eq. 1.19) becomes finite, the governing equations must adapt as shown in Fig. B.11. The above techniques coupled with the equations in Appendix A are appropriate for continuum flow (Kn0), and can even be used for small but finite Kn values so long as accommodation coefficients are included (Eq. 6.98a). However for Knudsen number of order of unity or greater, it is important to simulate the finite molecular spacing and collision effects. The most common numerical approach for this is the Direct- Simulation Monte-Carlo Method (DSMC) which utilizes many random “representative” molecular trajectories to represent the much larger number of actual molecules present in a system (Bird, 1976). The DSMC method follows such molecules as they interact with each other and any surface or flow boundaries. To focus on nearby interactions, the domain is divided into cells and evaluated over discrete time intervals. A representative molecule is considered within a cell when it is closest to the point which specifies that cell. As such, cells need not be structured and may be simply a distribution of points in space, which is convenient for complex domains. Furthermore, the molecular properties need not be stored with respect to a cell to allow for efficient data storage. The cell population is chosen to be statistically significant, e.g. on the order of 20 representative molecules. As with continuum CFD, the DSMC cell sizes should be such that the flow properties do not substantially vary across a single cell. Similarly, the time-step is chosen to be small compared to the mean collision time per molecule (t«m-m) so that the molecular movement and collisions can be decoupled and treated sequentially within the time-step period. For unsteady flows, initial conditions are implemented by specifying random molecular velocities that are consistent with the mean temperature and mean flow speed. At each time- step, the molecules are first moved according to their individual velocities or introduced at an inflow boundary. Appropriate action is taken if a molecule crosses a domain boundary (e.g. leaves for outflow, returns with the same tangential velocity for symmetry, etc.). Then a number of molecule-molecule collisions are simulated which are appropriate to the time-step, mean collision rate, and the number of representative molecules in a cell. For each specified collision to occur within a cell, two pairs of molecules are randomly selected to see whether they should collide (Bird, 1976). This is determined by using an acceptance-rejection method which is based on a random number generator and a Probability Distribution Function determined from the relative velocities (neglecting relative positions) between the two molecules. If a pair is rejected to collide, then another pair is selected and examined until a collision is “accepted”. At this point, the pre-collision velocities are replaced by the post- collision velocities using an inter-molecular force model, e.g. the inverse power law approximation or a more sophisticated model which includes the attraction forces. To compute mean properties, a sampling time is used which is much greater than the time-step thus ensuring a large number of representative collisions have occurred. As with the (deterministic) LB method, this technique is computationally intensive but is highly parallelizable as it depends only on nearest neighbor information. 74 References Abrahamson, J. (1975) “Collision rates of small particles in a vigorously turbulent fluid” Chem. Engng. Sci. Vol. 30, pp. 1371-1379. Achenbach, E. (1972) “Experiments on the flow past spheres at very high Reynolds numbers,” Journal of Fluid Mechanics, Vol.54, Part.3, pp. 565-575. Achenbach, E. (1974) “The effects of surface roughness and tunnel blockage on the flow past spheres,” Journal of Fluid Mechanics, Vol. 65, Part 1, pp. 113-125. Acrivos, A. Mauri, R. and Fan, X. (1993) “Shear-induced resuspension in a Couette device,” Int. J. Multiphase Flow, Vol. 19, pp. 797–802. Adams, P.J. and Seinfeld, J.H. (2003) “Disproportionate impact of particulate emissions on global cloud condensation nuclei concentrations,” Geophysical Research Letters, Vol. 30, No. 5, pp. 43-1 - 43-4 Aerosty, J. (1962) “Sphere drag on a low density flow”, Report HE-150-192, Univ. of California Berkeley, January. Ahmadi, G. & Ma, D. (1990) “A thermodynamical formulation for dispersed multiphase turbulent flows: I basic theory” Int. J. of Multiphase Flow, Vol. 16, pp. 323-340. Ahmadi,G. (2005), Clarkson University, personal communication. Ahuja, V. Cavallo P.A. and Hosangadi, A. (2000) “Multi-phase flow modeling on adaptive unstructured meshes,” AIAA Fluids 2000 and Exhibit, Paper No. AIAA- 2000-2662, June 19-22, Denver, Colorado. Al Taweel, A.M. and Landau, J. (1977) "Turbulence modulation in two-phase jets," Intl. J. Multiphase Flow, Vol. 3, pp. 341-351. Albertson, L. (1952) “Effect of shape on the fall velocity of gravel particles,” Proc. 5th Hydr. Conf. Studies in Engineering, University of Iowa, Iowa City, Iowa. Alexander, Z. (2004) “High order computation of the history term in the equation of motion for a spherical particle in a fluid,” Journal of Scientific Computing, Vol. 21, pp.129-143. Alger, G.R. & Simons (1968) J. Hyd. Eng. Div. ASCE 94, p. 721. Aliseda, A., Cartellier, A., Hainaux, F. and Lasheras, J.C. (2002) “Effect of preferential concentration on the settling velocity of heavy particles in homogeneous isotropic turbulence,” Journal of Fluid Mechanics, Vol. 468, pp. 77-105 Ambari, A. Gauthier-Manuel, B. & Guyon, E. (1984) “Wall effects on a Sphere Translating at Constant Velocity” Journal of Fluid Mechanics, Vol. 149, pp. 235-353. Amsden, A.A. (1993) “KIVA-3: a KIVA program with block-structured mesh for complex geometries,” (Tech Report LA-12503-MS), Los Alamos National Laboratory. Amsden, A. A., O'Rourke, P. J. and Butler, T. D. (1989), "KIVA-II: a computer program for chemically reactive flows with sprays," Los Alamos National Laboratory report LA-11560-MS. Anderson, D.A., Tannehill, J.C. & Pletcher, R.H. (1997) Computational Fluid Mechanics and Heat Transfer, Taylor & Francis, Philadelphia. Anderson, J.D. (1995) Computational Fluid Dynamics, McGraw-Hill, New York. Apte, S.V., Mahesh, K. and Lundgren, T. (2008) “Accounting for finite-size effect in disperse two-phase flow,” International Journal of Multiphase Flow, Vol. 34, pp. 260-271. 75 Armenio, A. & Fiorotto, V. (2001) “The importance of the forces acting on particles in turbulent flows,” Phys. Fluids, Vol. 13, pp. 2437–2440. Aroesty, J. (1962) “Sphere drag in a low density supersonic flow,” Univ. Calif. Tech. Rept. HE-150-192, Berkeley, 1962. Ashgriz, N. and Poo, J.Y. (1989) “Coalescence and separation in binary collisions of liquid drops” Journal of Fluid Mechanics, Vol. 221, pp. 183-204. Asmolov, E.S. (1999) "The inertial lift on a spherical particle in a plane Poiseuille flow at large channel Reynolds number," J. Fluid Mech. Vol. 381, p. 63. Asmolov, E.S. & McLaughlin, J.B. (1999) "The inertial lift of an oscillating sphere in a linear shear flow field," International Journal of Multiphase Flow, Vol. 25, pp. 739- 751. Auton, T.R., (1987) "The lift force on a spherical body in a rotational flow," J. Fluid Mech. Vol. 183, pp. p. 199. Auton, T.R., Hunt, J.C.R. and Prud'Homme, M. (1988) "The force exerted on a body in inviscid unsteady non-uniform rotational flow," J. Fluid Mech. Vol. 197, pp. 241-257. Aybers, N.M. and Tapucu, A. (1969) “The motion of gas bubbles rising through stagnant liquid,” Wärme Stoffübertrag. Vol. 2, pp.118-128. Azuma, H. and Yoshihara, S. (1999) "Three-dimensional large-amplitude drop oscillations: experiments and theoretical analysis," J. Fluid Mech. Vol. 393, pp. 309- 332. Baba, J. and Komar, P.D. (1981) “Measurements and analysis of settling velocities of natural quartz sand grains,” J. Sed. Petrol. Vol. 51, pp. 631–640. Babcock & Wilcox (2005) personal communication. Bagchi, P. and Balachandar, S. (2002a) “Effect of free rotation on motion of a solid sphere,” Physics of Fluids, Vol. 14, No. 8, pp. 2719-2737. Bagchi, P. and Balachandar, S. (2002b) “Shear versus vortex-induced lift on a rigid sphere at moderate Re,” J. of Fluid Mech. Vol. 473, pp. 379-388. Bagchi, P. and Balachandar, S. (2002c) “Steady planar straining flow past a rigid sphere at moderate Reynolds number,” J. of Fluid Mech. Vol. 466, pp. 365. Bagchi, P. and Balachandar, S. (2003) “Effect of turbulence on the drag and lift of a particle,” Physics of Fluids, Vol. 15, No. 11, pp. 3496-3513. Bailey, A.B. & Hiatt, J. (1972) “Sphere drag coefficients for a broad range of Mach and Reynolds numbers,” AIAA J., Vol.10, pp 1436-1440. Balachandar, S. and Ha, M.Y. (2001) “Unsteady heat transfer from a sphere in a uniform cross-flow,” Physics of Fluids Vol. 13, No. 12, pp. 3714-3728. Balachandar, S. and Maxey, M. (1989) “Methods for evaluating fluid velocities for spectral simulations of turbulence,” J. of Computational Physics, Vol. 83, pp. 96-125. Balaras, E., Benocci, C. & Piomelli, U. (1996) AIAA J. Vol.34, pp. 1111-1119. Barkla, H.M. & Auchterlonie, L.J. (1971) “The Magnus or Robbins effect on rotating spheres,” J. Fluid Mech. Vol. 47, pp. 437-448. Barth, T. & Ohlberger, M. (2004) “Finite volume methods: foundation and analysis,” Encyclopedia of Computational Mechanics. Eds. E. Stein, R. de Borst & T.J.R. Hughes, John Wiley & Sons. Barton, I.E. (1996) "Exponential-Lagrangian tracking schemes applied to Stokes law," AMSE Journal of Fluids Engineering, Vol. 118, pp. 85-89. 76 Basset, A.B. (1888) "On the motion of a sphere in a viscous liquid," Phil .Trans. R. Soc. London, Vol. 179A, pp. 43-63. Basu, P. (1999) “Combustion of coal in circulating fluidized-bed boilers: a review,” Chemical Engineering Science, Vol. 54, pp. 5547-5557. Bataille, J., Lance, M., and Marie, J.L. (1991) "Some aspects of the modeling of bubbly flows," Phase-Interface Phenomena in Multiphase Flow, Hewitt, G.F., Mayinger, F. and Riznic, J.R. eds., Hemisphere Publishing Corporation. Batchelor, G.K. (1967) An Introduction to Fluid Dynamics, Cambridge University Press, Cambridge, England. Batchelor, G.K. (1972) "Sedimentation in a dilute dispersion of spheres," J. Fluid Mech. Vol. 52, pp. 245-268. Beard, K.V. (1976) “Terminal velocity and shape of clouds and precipitation drops aloft,” Journal of Atmospheric Science, Vol. 33, pp. 851-864. Beard, K.V. & Prupacher, H.R. (1969) "A determination of the terminal velocity and drag of small water drops by means of a wind tunnel," J. of Atmos. Sci. Vol. 26, pp. 1066- 1071. Benjamin, T.B. (1987) "Hamiltonian theory for motions of bubbles in an infinite liquid," J. of Fluid Mechanics, Vol. 181, pp. 349-379. Berlemont, A., Achim, P. and Chang, Z. (2001) “Lagrangian approaches for particle collisions: the colliding particle velocity correlation in the multiple particles tracking method and in the stochastic approach,” Physics of Fluids, Vol. 13, pp. 2946-2956. Biesheuvel, A. & Spoelstra, S. (1989) “The added mass coefficient of a dispersion of spherical gas bubbles in liquid,” International Journal of Multiphase Flow, Vol. 15, pp. 911-924. Biesheuvel, A. & van Wijngaarden, L. (1984) “Two-phase flow equation for a dilute dispersion of gas bubbles in liquid,” Journal of Fluid Mechanics, Vol. 148, pp. 301- 318. Beard, K.V. (1976) “Terminal velocity and shape of cloud and precipitation drops aloft,” J. Atmos. Sci. Vol. 33, pp. 851-864. Bel Fdhila, R. & Duineveld, P.C. (1996) “The effect of surfactant on the rise of a spherical bubble at high Reynolds and Peclet numbers,” Phys. Fluids, Vol. 8, pp. 310- 321. Benjamin, T.B. (1987) "Hamiltonian theory for motions of bubbles in an infinite liquid," J. of Fluid Mechanics, Vol. 181, pp. 349-379. Benson, C.M., Levin, D.A., Zhing, J., Gimelshein, S.F. & Montaser, A. (2004) " Kinetic model for simulation of aerosol droplets in high-temperature environments," J. of Thermophysics and Heat Transfer, Vol. 18, pp. 122-134. Berelemont, A., Desjonqueres, P. and Gouesbet, G. (1990) "Particle Lagrangian simulation in turbulent flows," International Journal of Multiphase Flow, Vol. 16. Besnard, D.C., Kataoka, I. & Serizawa, A. (1991) "Turbulence modeling and multiphase turbulence transport modeling," ASME Fluids Engineering Meeting, FED-Vol. 110, pp. 51-57. Bhaga, D. & Weber, M.E. (1981) “Bubbles in viscous liquids: shapes, wakes and velocities,” J. of Fluid Mechanics, Vol. 105, pp. 61-85. 77 Bhargava, C., Loth, E. & M. Potapczuk (2005) "Numerical Simulation of Icing Clouds in the NASA Glenn Icing Research Tunnel", AIAA Journal of Aircraft Vol. 42, pp. 671- 684, May-June 2005. Biesheuvel, A. & Spoelstra, S. (1989) “The added mass coefficient of a dispersion of spherical gas bubbles in liquid,” Intl J. Multiphase Flow, Vol. 15, pp. 301-318. Bilger, R.W. (1976) “Turbulent diffusion flames,” Progress in Energy and Combustion Science, Vol. 1, pp. 87–109. Bird, G.A. (1976) Molecular Gas Dynamics, Clarendon Press, Oxford. Bird, G.A. (1994) Molecular Gas Dynamics and the Direct Simulation of Gas Flows, Oxford Science Publications. Blanco, A. & Magnaudet, J. (1995) “The structure of the axisymmetric high-Reynolds number flow around an ellipsoidal bubble of fixed shape,” Phys. Fluids, Vol. 7, pp. 1265–1274. Bocksell, T. (2004) "Numerical simulation of turbulent particle diffusion," Ph.D. Thesis, University of Illinois at Urbana-Champaign, Aeronautical & Astronautical Engineering. Bocksell, T. and Loth, E. (1999) "A CRW model for diffusion of bubbles and particles," ASME Summer Fluids Engineering Meeting, San Francisco, June, FEDSM99-7366. Bocksell, T. and Loth, E. (2001) "Discontinuous and continuous random walk models for particle diffusion in free-shear flows," AIAA Journal, Vol. 39, pp. 1086-1096, June. Bocksell, T. and Loth, E. (2006) “Stochastic modeling of particle diffusion in a turbulent boundary layer,” Intl. J. of Multiphase Flow, Vol. 32, pp. 1234-1253. Bonometti, T. & Magnaudet, J. (2006) “Transition from spherical cap to toroidal bubbles,” Physics of Fluids, Vol. 18, pp. 052102. Bosse, T., Kleiser, L. & Meiburg, E. (2006) “Small particles in homogeneous turbulence. Settling velocity enhancement by two-way coupling,” Physics of Fluids, Vol. 18. Boussinesq, J. (1877), "Théorie de l’Écoulement Tourbillant", Mem. Présentés par Divers Savants Acad. Sci. Inst. Fr., Vol. 23, pp. 46-50. Boussinesq, J. (1903) Theorie Analytique de la Chaleur (L'Ecole Polytechnique, Paris), Vol. 2, pp. 224. Box, G.E.P & Muller, M.E. (1958) “A note on the generation of random normal deviates,” Annals Math. Stat, Vol. 29, pp. 610-611 Boysan, F., Ayers, W.H., & Swithenbank, J. (1982) “A fundamental mathematical modeling approach to cyclone design,” Trans. Inst. Chem. Eng., Vol. 60, pp. 222. Brady, J.F. (1993) “The rheological behavior of concentrated colloidal dispersions,” Journal of Chemical Physics, Vol. 99, pp.567. Bragg, M.B. (1982) "A similarity analysis of the droplet trajectory equation," AIAA Journal, pp. 1681-1686. Bragg, M.B. (2005), University of Illinois at Urbana-Champaign, personal communication Braslow, A.L., Hicks, R.M., Harris, R.V. (1966) “Use of grit-type boundary layer- transition trips,” NASA SP-124. Brazier-Smith, P. R., Jennings, S. G. and Latham, J. (1972) “The interaction of falling water drops: coalescence,” Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences Vol. 326, No. 1566, pp. 393-408 78 Brenn, G., Valkovska, D. and Danov, K.D. (2001) “The formation of satellite droplets by unstable binary drop,” Physics of Fluids, Vol. 13, pp. 2463-2477. Brennen, C.E. (1995) Cavitation and Bubble Dynamics, Oxford University Press. Brennen, C.E. (2005) Fundamentals of Multiphase Flow, Cambridge University Press. Brennen, C.E. and H. Winet (1977) “Fluid mechanics of propulsion by cilia and flagella,” Annual Rev. of Fluid Mechanics, Vol. 9, pp. 339-398. Brenner, H. (1964) “The Stokes resistance of an arbitrary particle—IV arbitrary fields of flow,” Chemical Engineering Science, Vol. 19, pp. 703-727. Brenner, H. & Cox, R.G. (1963) “The resistance to a particle of arbitrary shape in translational motion at mall Reynolds numbers,” J. Fluid Mech. Vol. 17, pp. 561-595. Briggs, L.J. (1959) “Effects of spin and speed on the lateral deflection (curve) of a baseball and the Magnus effect for smooth spheres,” Am. J. Phys. Vol. 27, pp. 589-596. Brock, J.R. (1962) “On the theory of thermal forces acting on aerosol particles,” J. Colloid Sci. Vol. 17, pp. 768-780. Broday, D., Fichman, M., Shapiro, M. and Gutfinger, C. (1998) “Motion of spheroidal particles in vertical shear flows,” Physics of Fluids, Vol. 10, pp. 86-100. Bröder, D., Sornmerfeld, M. and Tisler, S. (1998) “Analysis of the hydrodynamics in a bubble column by phase-doppler anemometry,” Proc. 3rd lnt. Conf. Multiphase Flow, Lyon (France), June 1998. Brown, F.N.M. (1959) Personal communication, courtesy of the University of Notre Dame. Brucker, C. (1999) “Structure and dynamics of the wake of bubbles and its relevance for bubble injection,” Physics of Fluids, Vol. 11, pp. 1781-1796. Builtjes, P.J.H. (1975) "Determination of the Eulerain longitudinal integral length scale in a turbulent boundary layer," Applied Sci. Res. Vol. 31, pp. 397-399. Burger, J. M. (1942) “On the influence of the concentration of a suspension upon the sedimentation velocity,” Proc. Kon. Nederl. Akad. Wet. Vol. 45, pp. 126. Burton, T.M. & Eaton, J.K. (2005) “Fully resolved simulations of particle-turbulence interaction,” Journal of Fluid Mechanics, Vol. 545, pp. 67-111. Bush, R.H. and Mani, M. (2001) “A two-equation large eddy stress model for high sub- grid shear,” 15th AIAA Computational Fluid Dynamics Conference, June, Anaheim, CA, AIAA 2001-2561. Carlson, D.J. & Hoglund, R.F. (1964) “Particle drag and heat transfer in rocket nozzles,” AIAA J. Vol. 2, pp. 1980-1984. Chahed, J., Roig, V. & Masbernat, L. (2003) "Eulerian-Eulerian two-fluid models for turbulent gas-liquid bubbly flows," International Journal of Multiphase Flow, Vol. 29, pp. 23-49. Chahine, G.L., Delepoule, E. and Hauwert, P. (1993) Cavitation and Multiphase Flow Forum, Washington D.C., ASME FED- Vol. 153, pp. 39-45. Chang, E.J. and Maxey, M.R. (1994) "Unsteady flow about a sphere at low to moderate Reynolds number. Part I. Oscillatory motion," J. Fluid Mech. Vol. 277, pp. 347-379. Chang, H.T., Hourng, L.W. & Chien, L.C. (1996) “Application of flux-vector splitting scheme to a dilute gas-particle JPL nozzle flow,” Intl. J. for Num. Methods in Fluids. Vol. 22, pp. 921-935. 79 Chang, Y.C., Hou, T.Y., Merriman, B. & Osher, S. (1992) “A level-set formulation of Euleran interface capturing methods for incompressible fluid flows,” Journal of Computational Physics, Vol. 124, pp. 449-464. Chapman, B.K. & Leighton, D.T. (1991) “Dynamic viscous resuspension,” Int. J. Multiphase Flow, Vol. 17, pp.469-483. Chen, C.P & Wood, P.E. (1985) "Turbulence closure models for dilute gas-particle flows," Canadian J. of Chem. Engr. Vol. 63, pp. 349-360. Chen, J.-H., Wu, J.-S. & Faeth, G.M. (2000) "Turbulence generation in particle-laden flows" AIAA J, Vol. 38, pp. 636-642. Chen, R.C. and Wu, T.L. (1999) “The flow characteristics between two interactive spheres,” ASME Fluids Engineering Conference, San Francisco, FEDSM99-7776, June. Chen, X.-Q. (1997) "Efficient particle tracking algorithm for two-phase flows in geometries using curvilinear coordinates," Numerical Heat Transfer. Chen, X.-Q. & Pereira, J.C.F. (1994) "Comparison of the various particle dispersion models in plane mixing shear layers," Numerical Methods in Multiphase Flows, FED- Vol. 185, pp. 217-225. Chen, X.-Q. & Pereira, J.C.F. (1995) "Eulerian-Eulerian predictions of dilute turbulent gas-particle flows," Numerical Methods in Multiphase Flows, FED-Vol. 228, pp. 265-272. Chen, Y. and Heister, S.D. (1994) "A numerical treatment for attached cavitation," Journal of Fluids Engineering, Transactions of the ASME, Vol. 116, pp. 613-618. Cheng, N.S. (1997) "Simplified settling velocity formula for sediment particle," Journal of Hydraulic Engineering, Vol. 123, pp. 149-152. Cherukat, P., McLaughlin, J.B. and Graham, A.L. (1994) "The inertial lift on a rigid sphere translating in a linear shear flow field," International Journal of Multiphase Flow, Vol. 20, No. 2, pp. 339-353. Chester, D. (1993) Volcanoes and Society, London, Edward Arnold. Chhabra, R.P. (1995) Powder Technology. Vol. 85, pp. 83-90. Chhabra, R.P., Agarwal, L. & Sinha, N.K. (1999) "Drag on non-spherical particles: an evaluation of available methods," Powder Technology. Vol. 101, pp. 288-295. Christiansen, E.B. & Barker, D.H. (1965) “The effect of shape and density on the free settling of particles at high Reynolds number,” AIChE J. Vol. 11, pp. 145-151. Chun, J., Koch, D.L., Rani, S., Ahluwalia, A. and Collins, L. (2005) “Clustering of aerosol particles in isotropic turbulence,” J. Fluid Mech. (submitted). Chung, J., Choo, Y., Reehorst, A., Potapczuk, M. and Slater, J. (1999) “Navier-Stokes analysis of the flowfield characteristics of an ice contaminated aircraft wing,” NASA/TM—1999-208897. Chung, T.J. (2002) Computational Fluid Dynamics, Cambridge University Press, Cambridge, U.K. Clamen, W.H. & Gauvin, A. (1969) “Effect of turbulence on the drag coefficients of spheres in a supercritical flow regime,” AIChE J. Vol. 15, p. 184. Clark, N.N., Gabriele, P., Shuker, S. & Turton, R. (1989) Powder Technology, Vol. 59, p. 69. Cliffe, K.A. and Lever, D.A. (1985) “Isothermal flow past a blowing sphere,” International Journal of Numerical Methods in Fluids, Vol. 5, pp. 785-802. 80 Clift, R. and Gauvin, W.H. (1970) “The motion of particles in turbulent gas streams,” Proc. CHEMECA '70, Butterworth, Melbourne, Vol. 1, pp. 14-28. Clift, R., Grace, J.R., and Weber, M.E. (1978) Bubbles, Drops and Particles, Academic Press, New York. Climent, E. & Maxey, M.R. (2003) “Numerical simulation of random suspensions at finite Reynolds numbers,” International Journal of Multiphase Flow, Vol. 29, pp. 579-601. Cohen, E.G.D., Verberg, R. & de Schepper, I.M. (1997) "Newtonian viscosity and visco- elastic behavior of concentrated neutral hard-sphere colloidal suspensions," Int. J. of Multiphase Flow, Vol. 23, pp. 797-807. Coimbra, C.F.M. & Kobayashi, M.H. (2002) "On the viscous motion of a small particle in a rotating cylinder," J. Fluid Mech. Vol. 469, pp. 257-286. Collins, J.P., Ferguson, R.E., Chien, K., Kuhl, A.L., Krispen, J. and Glaz, H.M. (1994) "Simulation of shock-induced dusty gas flows using various models," AIAA Fluid Dynamics Conference, Colorado Springs, CO, AIAA-94-2309. Colucci, P.J., Jaberi, F.A., Givi, P. & Pope, S.B. (1998) “Filtered density function for large eddy simulation of turbulent reacting flow,” Physics of Fluids, Vol. 10, pp. 499- 515. Comer, J.K. and Kleinstreuer, C. (1995) “A numerical investigation of laminar flow past non-spherical solids and droplets,” ASME J. Fluids Engineering, Vol. 117, pp. 170- 175. Cooper, C.A. and Young, R.J. (1999) “Investigation of tructure/property relationships in particulate composites through the use of Raman spectroscopy,” Journal of Raman Spectroscopy, Vol. 30, pp. 929-938. Coppen, S., Manno, V. and Rogers, C.B. (2001) "Turbulence characteristics along the path of a heavy particle," Computers and Fluids, Vol. 30, pp. 257-270. Coppola, G., Sherwin, S. J., & Peiro, J. (2001) “Nonlinear particle tracking for high-order elements,” Journal of Computational Physics, Vol. 172, No. 1, Sept., pp. 356-386. Couder, J. (2005) personal website, http://beaufix.ipgp.jussieu.fr/~couder/flowEnglish.html Cox , R.G. (1969) “The deformation of a drop in a general time-dependent fluid flow,” J. Fluid Mech. Vol. 37, pp. 601–623 Crowe, C.T. (1967) “Drag coefficient of particles in rocket nozzle,” AIAA J. Vol. 5, pp. 1021-1022. Crowe, C.T. (1982) “REVIEW - numerical models for dilute gas-particle flows,” J. Fluid Eng. Vol. 104, pp. 297-303, Sept. Crowe, C.T., Babcock, W.R. & Willoughby, P.G. (1972) “Drag coefficient for particles in rarefied low-Mach number flows,” Progress Heat Mass Transfer, Vol. 6, pp. 419- 431. Crowe, C.T., Chung, J.N. and Troutt, T.R. (1988) “Particle mixing in free shear flows,” Progress in Energy and Combustion Science, Vol. 14, No. 3, pp. 171-194. Crowe, C.T., Sommerfeld, M. and Tsuji, Y. (1998) Multiphase Flows with Droplets and Particles, CRC Press, Boca Raton, FL. Crowe, C.T., Troutt, T.R. and Chung, J.N. (1996) “Numerical models for two-phase turbulent flow,” Annu. Rev. Fluid. Mech. Vol. 28, pp. 11-43. 81 Csanady, G. T. (1963) “Turbulent diffusion for heavy particles in the atmosphere,” J. Atmospheric Sciences, Vol. 20, May, pp. 201-208. Cuenot, B., Magnaudet, J. and Spennato, B. (1997) “ Effects of slightly soluble surfactants on the flow around a spherical bubble,” Journal of Fluid Mechanics, Vol. 339, pp. 25-53. Curtis, J. (2000) Personal communication, Purdue University, Chemical Engineering. Czys, R.R. and Ochs, H.T. (1998) "The influence of charge on the coalescence of water drops in free fall," Journal of the Atmospheric Sciences, November, pp 3161-3168. D’Annibale, F. (2005) Private Communication, Institute of Thermal Fluid Dynamics, Italy. Dallavalle, J.M. (1948) Micrometrics, Pitman Publishing, New York. Dandy, H.A. and Dwyer, D.S. (1990) “Some influences of particle shape on drag and heat transfer,” Physics of Fluids A, Vol. 2, No. 12, pp. 2110-2118. Dao-Thi, M.-H., Maes, D., Lisgarten, J., Zegers, I. & Wyns, L. (2005) “Protein crystallization in microgravity conditions,” Personal Communication, Free University of Brussels. Darton R.C. and Harrison, D. (1974) “The rise of single gas bubbles in liquid fluidized beds,” Transactions of Institution Chemical Engineers, Vol. 52, pp. 301–306. Davies, J.M. (1949) “The aerodynamic of golf balls,” J. Applied Physics, Vol. 20, pp. 821-828. Davies, R. & Taylor, G.I. (1950) “The mechanics of large bubbles rising through extended liquids and though liquids in tubes,” Proc. R. Soc. London, Ser. A, Vol. 200, p. 375. Davis, R.H (1991) “Sedimentation of axisymmetric particles in a shear flows,” Physics of Fluids A, Vol. 3, No. 9, pp. 2051-2060. Davis, R.H. and Acrivos, A. (1985) "Sedimentation of non-colloidal particles at low Reynolds numbers," Annual Reviews of. Fluid Mech. Vol. 172, pp. 91-118. Davis, R.H, Rager, D.A. and Good, B.T. (2002) “Elastohydrodynamic rebound of spheres from coated surfaces,” Journal of Fluid Mechanics, Vol. 468, pp. 107-119. Davis, R.H, Serayssol, J.M. & Hinch, E.J. (1986) “The elastohydrodynamic collision of two spheres,” Journal of Fluid Mechanics, Vol. 163, p. 479. DeAngelis, B., Loth, E., Lankford, D. and Bartlett, C.S. (1997) "Computations of turbulent droplet dispersion for wind tunnel icing tests," AIAA Journal of Aircraft, Vol. 34, No.2, pp. 213-219, March-April. DeCroix, D.S. (2003) “Visualizing chemical dispersion in populated regions,“ E- newsletter for Tecplot Users, Issue 18 DeGaspari, J. (2003) “Penetrating vision,” Mechanical Engineering, May 2003. Desjardins, O., Fox, R.O. & Villedieu, P. (2007) “A quadrature-based moment method for dilute gas-particle flows,” 6th International Conference on Multiphase Flow, Leipzig, Germany. Dennis, S.C.R., Singh, S.N. and Ingham, D.B. (1980) "The steady flow due to a rotating sphere at low and moderate Reynolds numbers," J. Fluid Mech. Vol. 101, pp. 257- 279. Deshpande, M., Feng, J. and Merkle, C.L. (1993) "Navier-Stokes analysis of 2-D cavity flows," Cavitation and Multiphase Flow Forum, Washington D.C., ASME FED-Vol. 153, pp. 149-155. 82 di Felice (1994) "The voidage function for fluid-particle interaction systems," Int. J. Multiphase Flow, Vol. 20, pp. 153-159. Dorgan, A. & Loth, E. (2004) "Simulation of particles released near the wall in a turbulent boundary layer, " Intl. Journal of Multiphase Flow, Vol. 30, p. 649-673. Dorgan, A. & Loth, E. (2007) "Efficient calculation of the history force at finite Reynolds numbers," Intl. Journal of Multiphase Flow (to appear). Dorgan, A. & Loth, E. (2009) "Dispersion of bubbles in a turbulent boundary layer," ASME Fluids Engineering Division Summer Conference (submitted). Dorgan, A.J., Loth, E., Bocksell, T.L. & Yeung, P.K. (2005) "Boundary layer dispersion of near-wall injected particles of various inertias," AIAA Journal, Vol. 43, pp. 1537- 1548. Drazin, P.G. & Reid, W.H. (1981) Hydrodynamic Stability. Cambridge University Press, Cambridge, U.K. Dressel, M. (1985) “Dynamics shape factors for particle shape characterization,” Part. Charct. Vol. 2, pp. 62-66. Drew, D. A. (1983) “Mathematical modeling of two-phase flow,” Annual Review of Fluid Mechanics, Vol. 15, pp. 261-291. Drew D.A & Lahey, R.T. (1987) “The virtual mass and lift force on a sphere in rotating and straining inviscid flow,” International Journal of Multiphase Flow, Vol. 13, pp. 113–121. Drew, D.A. and Lahey, R.T. (1993) "Analytical modeling of multiphase flow," Particulate Two-Phase Flow, Ed. by M.C. Roco, Butterworth-Heinemann, Boston, MA. Drew, D.A. and Passman, S.L. (1998) “Theory of multi-component fluids,” Applied Mathematical Sciences, Vol. 135, Springer, New York. Drew, D.A. & Passman, S.L. (1999) Theory of multicomponents fluids, Springer, NY. Druzhinin, O.A. (1995) "On the two-way interaction in two-dimensional particle-laden flows: the accumulation of particles and flow modification," Journal of Fluid Mechanics, Vol. 297, pp. 49-76. Druzhinin, O.A. and Elghobashi, S. (1998) "Direct numerical simulations of bubble-laden turbulent flows using the two-fluid formulation," Physics of Fluids, Vol. 10, No. 3, pp. 685-697. Druzhinin, O.A. and Elghobashi, S. (1999) "On the decay rate of isotropic turbulence laden with micro-particles," Physics of Fluids, Vol. 11, No. 3, pp. 602-610. Duineveld, P.C. (1995) “Rise velocity and shape of bubbles in pure water at high Reynolds number,” Journal of Fluid Mechanics, Vol. 292, pp. 325-332. Dukowicz, J.K. (1980) “A particle-fluid numerical model for liquid sprays,” Journal of Computational Physics, Vol. 35, pp. 229-253. Dukowicz, J.K. (1982) “An exact solution for the drag of a sphere in low Reynolds number flow with strong uniform blowing or suction,” Physics of Fluids, Vol. 25, pp. 1117-1118. Dukowicz, J.K. (1984) “Drag of evaporating or condensing droplets in low Reynolds number flows,” Physics of Fluids, Vol. 27, pp. 1351-1358. Durst, F., Milojevic, D., and Schonung, B. (1984) "Eulerian and Lagrangian predictions of particulate two-phase flows: a numerical study,” Appl. Math. Modeling. Vol. 8, April, pp. 101-112. 83 Eaton, J.K (1999) “Local distortion of turbulence by dispersed particles,” AIAA Fluid Dynamics Conference, Norfolk, VA, AIAA-99-3643, June. Eaton, J.K. and Fessler, J.R. (1994) "Preferential concentration of particles by turbulence," Intl. J. of Multiphase Flow, Vol. 20, pp. 169-209. Edge, R. M. & Grant, C.D. (1971) “Terminal velocity and frequency of oscillation of drops in pure systems,” Chemical Engineering Science, Vol. 26, No. 7, pp. 1001- 1012. Edge, R.M. & Grant, C.D. (1972) “Motion of drops in water contaminated with a surface- active agent,” Chemical Engineering Science, Vol. 27, No. 9, pp. 1709-1721. Eisenklam, P., Arunchalam, S.A. and Weston, J.A. (1967) “Evaporation rates and drag resistance of burning drops,” Eleventh Symposium on Combustion, Combustion Institute, Pittsburgh, pp. 715-728. Elamworawutthikul and Gould (1999) “The flow structure around a suspended sphere at low Reynolds number in a turbulent freestream,” 3rd ASME/JSME Joint Fluids Engineering Conference, July, FEDSM99-6996 Elghobashi, S. (1994) "On predicting particle-laden turbulent flows," Applied Scientific Research, Vol. 52, pp. 309-329. Elghobashi, S. (2006) Multiphase Flow Handbook, CRC Press, C. Crowe edit., Boca Raton, FL, USA Elghobashi, S. & Abou-Arab, T.W. (1983) "A two-equation turbulence model for two- phase flows," Physics of Fluids, Vol. 26, pp. 931-938. Elghobashi, S. and Truesdell, G.C. (1984) "Direct simulation of particle dispersion in a decaying isotropic turbulence," Journal of Fluid Mechanics, Vol. 242, pp. 655-700. Ellingsen, K. and Risso, F. (2001) “On the rise of an ellipsoidal bubble in water: oscillatory paths and liquid-induced velocity,” Journal of Fluid Mechanics, Vol. 440, pp. 235-268. Elsner, J.W. & Elsner, W. (1996) "On the measurement of turbulence energy dissipation," Meas. Sci. Technol. Vol. 7, pp. 1334-1348. Eppinger, K. (1995) "Study of the movement of bubbles in homogeneous isotropic turbulence," Ph.D. Thesis, Instiut National Polytehnique de Toulouse, France. Epstein, N. and Pruden, B.B. (1999) “Liquid fluidization of binary particle mixtures – II stratification by size and related topics,” Chemical Engineering Science, Vol. 54, pp. 401-415. Epstein, P.S. (1929) "Zur theorie des radiometers," Z. Physik Vol. 54, pp. 537-563. Erbe, E. (2005) Personal Communication, USDA Beltsville Agricultural. Ervin, E.A. (1993) "Computations of bubbles and drops in shear flow," PhD Dissertation, Dept. of Mechanical Engineering, University of Michigan, Ann Arbor, MI. Estrade, J.-P., Carentz, Hervé, Lavergne, G. and Biscos, Y. (1999) “Experimental investigation of dynamic binary collision of ethanol droplets – a model for droplet coalescence and bouncing,” International Journal of Heat and Fluid Flow, Vol. 20, pp. 486-491. Fadlun, E.A., Verzicco, R., Orlandi, P. and Mohd-Yusof, J. (2000) “Combined immersed-boundary finite-difference methods for three-dimensional complex flows simulations,” Journal of Computational Physics, Vol. 161, pp. 35-60. 84 Faeth, G. M., (1987) “Mixing, transport and combustion in sprays,” Progress in Energy and Combustion Science, Vol. 13, pp. 293-345. Fan, L.-S. and Tsuchiya, K. (1990), Bubble Wake Dynamics in Liquids and Liquid-Solid Suspension, Butterworth-Heinemann, Boston. Fan, L.-S. and Tsuchiya, K. (1993) “Chapter 23: bubble flow in liquid-solid suspension,” Particulate Two-Phase Flow, edited by M.C. Roco, Butterworth- Heinemann, Boston. Fan, R., Marchisio, D.L. & Fox, R.O. (2004) “Application of the direct quadrature method of moments to polydisperse gas-solid fluidized beds,” International J. of Multiphase Flow, Vol. 139, pp. 7-20. Fang, H., Wan, R. and Lin, Z. (2002) “Lattice Boltzmann model with nearly constant density,” Phys. Rev. E, Vol. 66, pp. 66-69. Faxen, H. (1922) “Resistance to the movement of a rigid sphere in a viscous fluid bounded by two parallel flat walls,” Ann. Phys. Vol. 68, pp. 89–119. Felton, K. and Loth, E. (2001) "Spherical bubble motion in a turbulent boundary layer," Physics of Fluids, Vol. 13, No. 9, 2564-2577, September. Feng, Z. C. and Leal, L.G. (1997) "Nonlinear bubble dynamics," Annual Review of Fluid Mechanics, Vol. 29, December, pp. 201-243. Feng, Z.-G. and Michealidis, E.E. (2001) "Drag coefficients of viscous spheres at intermediate and high Reynolds numbers," Journal of Fluids Engineering, Vol. 123, December, pp. 841-849. Ferrante, A. & Elghobashi, S. (2003) "On the physical mechanism of two-way coupling in particle-laden isotropic turbulence," Physics of Fluids, Vol. 15, pp. 315-329. Ferrante, A. & Elghobashi, S. (2004) "On the physical mechanisms of drag reduction in a spatially developing turbulent boundary layer laden with microbubbles," J. Fluid Mechanics, Vol. 503, pp. 345-355. Ferrante, A. & Elghobashi, S. (2005) "On Reynolds number effect on drag reduction in a microbubble-laden spatially developing turbulent boundary layer," J. Fluid Mechanics, Vol. 543, pp. 93-106. Ferrante, A. & Elghobashi, S. (2007) “On the effects of finite-size particles on decaying isotropic turbulence,” International Conference on Multiphase Flow, Leipzig, Germany.CP03. "On the Ferry, J. and Balachandar, S. (2001) "A fast Eulerian method for disperse two-phase flow," Intl. J. of Multiphase Flow, Vol. 27, pp. 1199-1226. Ferry, J. and Balachandar, S. (2002) “Equilibrium expansion for the Eulerian velocity of small particles,” Powder Technology, Vol. 125, No. 2-3, pp. 131-139 Fessler, J.R., Kulick, J.D. & Eaton, J. K. (1994) “Preferential concentration of heavy particles in a turbulent channel flow,” Phys. Fluids, Vol. 6, pp. 3742-3749. Field, S.B., Klaus, M., Moore, M.G. and Nori, F. (1997) "Chaotic dynamics of falling disks," Nature, Vol. 388, pp. 252-254, July. Foerster, S.M., Louge, M.Y., Chang, H. & Khedidja, A. (1994) “Measurement of the collision properties of small spheres,” Physics of Fluids, Vol. 6, pp. 1108-1115. Foote, G.B. (1974) “The water drop rebound problem: dynamics of collision,” Journal of the Atmospheric Sciences, Vol. 32, pp. 390-402. Ford, B. and Loth, E. (1998) "Forces on ellipsoidal bubbles in a turbulent free shear layer," Physics of Fluids, Vol. 10, No. 1, pp. 178-188, Jan. 85 Forsberg, F., Goldberg, B.B., Liu, J.-B., Merton, D.A., Rawool, N.M. and Shi, W.T. (1999) “Tissue-specific US contrast agent for evaluation of hepatic and splenic parenchyma,” Radiology, Vol. 210, pp. 125-132. Fortes, A.F., Joseph, D.D. and Lundgren, T.S. (1987) “Nonlinear mechanics of fluidization of beds of spherical particles,” J. Fluid Mechanics, Vol. 177, pp. 467-483. Fox, R.O. (2003) Computational Models for Turbulent Reacting Flows, Cambridge University Press. Fraser, A. (2005) Personal Communication about Bad Rain, Pennsylvania State University. Freidman, P.D. & Katz, J. (2002) "Mean rise rate of droplets in isotropic turbulence," Physics of Fluids, Vol. 14, pp. 3059-3073. Galdi, G., Vaidya, A., Pokorny, M., Joseph, D.D. and Feng, J. (2002) “Orientation of symmetric bodies falling in a second-order liquid at non-zero Reynolds number,” Mathematical Models and Methods in Applied Sciences, Vol. 12, No. 11, pp. 1653- 1690. Gañán-Calvo, A.M. and Lasheras, J.C. (1991) "The dynamics and mixing of small spherical particles in plane, free shear layer," Physics of Fluids, Vol. 3, No. 5, pp. 1207-1217. Ganser, G.H. (1993) “A rational approach to drag prediction of spherical and non- spherical particles,” Powder Technology, Vol. 77, pp. 143. Garg, K. and Nayar, S.K. (2004) “Photometric model of a raindrop,” Proceedings of IEEE Conference on Computer Vision and Pattern Recognition, Washington DC, June. Garner, F.H. & Hammerton, D. (1954) “Circulation inside gas bubbles,” Chem. Eng. Sci. Vol. 3, No. 1, pp. 1–11. Garner, F.H. and Lihou, D.A. (1965) DECHEMA-Monograph, Vol. 55 pp. 155-178. Garnier, E., Sagaut, P. and Deville, M. (2002) “Large eddy simulation of shock/boundary-layer interaction,” AIAA Journal, Vol. 40, No. 10, pp. 1935–1944. Gauthier, D., Zerguerra, S. and Flamant, G. (1999) “Influence of the particle size distribution of powders on the velocities of minimum and complete fluidization,” Chemical Engineering Journal, Vol. 74, pp. 181-196. Gauthier, G., Gondret, P. and Rabaud, M. (1998) “Motions of anisotropic particles: application to visualization of three-dimensional flows,” Physics of Fluids, Vol. 10, No. 9, September, pp. 2147-2154. Gavze, E. and Shapiro, M. (1997) "Particles in a shear flow near a solid wall: effect of nonsphericity on forces and velocities," International Journal of Multiphase Flow, Vol. 23, No. 1, pp. 155-182. Geisler, R. (2005) “Sonolumineszenz,”,personal webpage: http://www.physik3.gwdg.de/~rgeisle/nld/sl.html. Gelfand, B.E. (1996) “Droplet breakup phenomena in flows with velocity lag,” Progress in Energy and Combustion Science, Vol. 22, pp. 201-266. Gibilaro, L.G. (2001) Fluidization-Dynamics, Oxford Press, Oxford, U.K. Gidaspow, D. (1994) Multiphase Flow and Fluidization, Academic Press, Boston. Glowinski, R., Pan, T.W., Hesla, T.I., Joseph, D.D. and Périaux, J. (2001) “A fictitious domain approach to the direct numerical simulation of incompressible viscous flow 86 past moving rigid bodies: application to particulate flow,” Journal of Computational Physics, Vol. 169, pp. 363–426. Gogus, M., Ipecki, O.N. and Kokpinar, M.A. (2001) "Effect of particle shape on fall velocity of angular particles," Journal of Hydraulic Engineering, Vol. 127, No. 10 October, pp. 860. Golman, A.J., Cox, R.G. & Brenner, H. (1967) "Slow viscous motion of a sphere parallel to a plane wall," Chem Engr. Sci. Vol. 22, pp. 637-651. Gondret, P., Lance, M. & Petit, L. (2002) “Bouncing motion of spherical particles in fluids,” Phys Fluids, Vol. 14, No. 2, pp. 643–652. Gore, R.A. and Crowe, C.T. (1989) "Effect of particle size on modulating turbulent intensity," International Journal of Multiphase Flow, Vol. 15, pp. 279. Gorial, B.Y. & O’Callaghan, J.R. (1990) "Aerodynamic properties of grain/straw materials," J. Agric. Engrg. Res. Vol. 46, pp. 275-290. Gosman, A.D. and Ioannides, E. (1981) "Aspects of computer simulation of liquid-fueled combustors," Aerospace Sciences Meeting, Reno NV, AIAA-81-0323, Jan. Govier, G.W. and Aziz, K. (1972) The Flow of Complex Mixtures in Pipes, Van Nostrand-Reinhold, New York. Grace, J.R. (1973) “Shapes and velocities of bubbles rising in infinite liquids,” Trans. Inst. Chemical Engineering, Vol. 51, pp. 116-120. Grace, J.R. and Sun, G. (1991) “Influence of the particle size distribution of powders on the performance of fluidized bed reactors,” The Canadian Journal of Chemical Engineering, Vol. 69, pp. 1126-1134. Grace, J.R., Wairegi, T. & Nguyen, T.H. (1976) “Shapes and velocities of single drops and bubbles moving freely through immiscible liquids,” Trans. Inst. Chemical Engineering, Vol. 54, pp. 167-173. Graham, D.I. (1998) "Improved eddy interaction models with random length and time scales," Intl. Journal of Multiphase Flow, Vol. 24, No. 2, pp. 335-345. Graham, D.I. and James, P.W. (1996) "Turbulent dispersion of particles using eddy interaction models," International Journal of Multiphase Flow, Vol. 22, No. 1, pp. 157-175. Greenspan, H.P. and Nadim, A. (1993) “On sonoluminescence of an oscillating gas bubble,” Phys. Fluids A, Vol. 5, 1065–1067. Griffith, R.M. (1962) “The effect of surfactants on terminal velocity of drops and bubbles,” Chemical Engineering Science, Vol. 17, pp. 1057-1070. Groszmann, D.E., Fallon, T.M. & Rogers, C.B. (1999) “Decoupling the roles of inertia and gravity on the preferential concentration of particles,” 3rd ASME/JSME Joint Fluids Engineering Conference, FEDSM-99, San Francisco, CA. ASME, New York. Haberman, W.L. and Morton, R.K. (1953) “Experimental investigation of drag and shape of air bubbles,” David W. Taylor Model Basin Report, Vol. 802. Haberman, W.L. and Morton, R.K. (1954) “Experimental study of bubbles moving in liquids,” American Society of Civil Engineers -- Proceedings Separates, Vol. 80, Separate No. 387, 25p. Haider, A. and Levenspiel, O. (1989) "Drag coefficient and terminal velocity of spherical and non-spherical particles," Powder Tech. Vol. 58, pp. 63-70. 87 Hallouin, E., Gondret, P., Lance, M., and Petit, L. (1998) "On the motion of a sphere toward a plane surface: from lubrication to bouncing regime," International Conference on Multiphase Flow, Lyon, France, June. Hamed, A. (1988) "Effect of particle characteristics on trajectories and blade impact patterns," AMSE Journal of Fluids Engineering, Vol. 110, pp. 33-37. Hamed, A. and Jun, Y.D. & Yeuan, J.J. (1995) "Particle dynamics simulations in inlet separator with an experimentally based bounce model," J. of Propulsion and Power Congress & Expo. Vol. 11, pp. 230-235. Hamed, A. and Kuhn, T.P. (1993) "Effects of variational particle restitution characteristics on turbo-machinery erosions," ASME Intl. Gas Turbine and Aeroengine Congress & Expo. Cincinnati, OH, May 24-27, 93-GT-124. Hanratty, T. J. and Bandukwala, A. (1957) “Fluidization and sedimentation of spherical particles,” AIChE Journal, Vol. 3, pp. 293-296. Happel, J. & Brenner, H. (1973) Low Reynolds Number Hydrodynamics, Noordhoff. Harding, S.E. (1995) “On the hydrodynamic analysis of macromolecular conformation,” Biphysical Chemistry, Vol. 55, pp. 69-93. Harmathy, T.Z. (1960) “Velocity of large drops and bubbles in media of infinite or restricted extent,” AIChE J, Vol. 6, pp. 281-288. Hartunian, R.A. and Sears, W.R. (1957) “On instability of small gas bubbles moving uniformly in various liquids,” Journal of Fluid Mechanics, Vol. 3, Part 1, pp. 27-47. Haworth, D.C. & Pope, S.B. (1987) “A generalized Langevin model for turbulent flows,” Physics of Fluids, Vol. 30, p. 1026. Healy, D.P. & Young, J.B. (2005) "Full Lagrangian methods for calculating particle concentration fields in a dilute gas-particle mixture," Proc. Royal. Soc. A, Vol. 461, pp. 2197-2225. Helenbrook, B.T & Edwards, C.F. (2002) "Quasi-steady deformation and drag of uncontaminated liquid drops," Intl. Journal of Multiphase Flow, Vol. 28, pp. 1631- 1657. Henderson, C.B. (1976) “Drag coefficients of spheres in continuum and rarefied flows,” AIAA J. Vol. 14, pp. 707-708. Hermsen, R.W. (1979) “Review of particle drag models,” JANAF Performance Standardization Subcommittee 12th Meeting Minutes, CPIA, p. 113. Heron, I., Davis, S. & Bretherton, F. (1975) “On the sedimentation of a sphere in a centrifuge,” Journal of Fluid Mechanics, Vol. 68, pp. 209-234. Hersch, A.S., Friichtenicht, J.F. & Slattery, J.C. (1969) “Drag of microscopic spheres in free-molecule flow,” Sixth International Symposium on Rarefied Gas Dynamics, Vol. 1, Academic Press, New York, pp. 757-766. Herzhaft, B. & Guazzelli, E. (1999) “Experimental study of the sedimentation of dilute and semi-dilute suspensions of fibres,” J. Fluid Mech. Vol. 384, pp. 133-158. Hetsroni, G. (1989) “Particles-turbulence interaction,” Int. J. Multiphase Flow, Vol. 5, pp. 735-746. Hill, R. and Power, G. (1956) "Extremeum principles for slow viscous flow and the approximate calculation of drag," Quart. J. Mech. Appl. Math. Vol. 9, pp. 313-319. Hinze, J.O. (1975), Turbulence, McGraw-Hill Book Company, New York. Hirsch, C. (1988) Numerical Computational of Internal and External Flows, John Wiley & Sons, New York. 88 Hobson, E. W. (1931) The Theory of Spherical and Ellipsoidal Harmonics, Cambridge Univ. Press, republished 1955 by Chelsea Publishing Comp., New York. Hoerner, S.F. (1965) Fluid-Dynamic Drag, published by author, Midland Park, N.J. Hosangadi, A., York, B.J., Sinha, N. and Dash, S.M. (1993) “Progress in transient interior ballistic flowfield simulations using multi-dimensional upwind/implicit numerics,” AIAA Paper 93-1915. Hsiao, C.-T. and Chahine, G.L. (2001) “Numerical simulation of bubble dynamics in a vortex flow using Navier-Stokes computations and moving chimera grid scheme,” CAV2001. Hu, S. & Kintner, R.C. (1955) “The fall of single drops through water,”, AIChE J, Vol. 1, pp. 42. Iliopoulos, I. & Hanratty, T.J. (1999) “Turbulent dispersion in a non-homogeneous field,” J. Fluid Mech. Vol. 392, pp. 45–71. Ishii, M. (1975) "Thermo-fluid dynamic theory of two-phase flow," Eyrolles, Paris, France. Ishii, M. and Chawla, T.C. (1979) "Local drag laws in dispersed two-phase flow," ANL- 79-105. Ishii, M. and Zuber, N. (1979) “Drag coefficient and relative velocity in bubbly, droplet, or particulate flows,” AIChE Journal, Vol. 25, pp. 843-855. Issa, R.I., Gossman, A.D. & Watkins, A. P. (1986) “The computation of compressible and incompressible recirculating flows by a NON-iterative implicit scheme,” J. Comp. Physics, Vol. 62, pp. 66-82. Jacobson, S. and J.R. Brock (1965) “The thermal force on spherical sodium chloride aerosols,” J. Colloid Sci. Vol. 20, pp. 544-554. Jahn, T.L. & Votta, J.J. (1972) “Locomotion of protozoa,” Annual Rev. of Fluid Mechanics, Vol. 4, pp. 93-116. Jayaweera, K.O.L.F. & Mason, B.J. (1965) “The behavior of freely falling cylinders and cones in viscous fluid,” J. Fluid Mech. Vol. 22, pp. 709– 720. Jeffrey, G.B. (1922) Proc. Roy. Soc. Ser A, Vol. 102, pp. 161-179. Jeffrey, R.C. & Pearson, J.R.A. (1965) “Particle motion in a laminar vertical tube flow,” J. Fluid Mechanics, Vol. 22, pp. 721-735. Jimenez, J.A. & Madsen, O.S. (2003) “A simple formula to estimate settling velocity of natural sediments,” J. Waterway, Port, Coastal & Ocean Engineering, pp. 70-78. Johnson, K.L. (1985) Contact Mechanics, Cambridge University Press, UK. Joia, I.A., Ushijima, T. and Perkins, R.J. (1997) "Numerical study of bubble and particle motion in turbulent boundary layer using proper orthogonal decomposition," Applied Sci. Res. Vol. 57, pp. 263-277. Joseph, D.D. (1994) "Interrogation of numerical simulation for modeling of flow induced microstructure," ASME/FED Liquid-Solid Flows, Vol. 189, pp. 31-40. Joseph, D.D. (2006) "Rise velocity of a spherical cap bubble," J. Fluid Mechanics, Vol. 488, pp. 213-223. Joseph, D.D., Belanger, J. and Beavers, B.S. (1999) "Break-up of a liquid drop suddenly exposed to high-speed airstream," Intl. J. of Multiphase Flow, Vol. 25, pp. 1263-1303. Joseph, G.G. & M.L. Hunt, (2004) “Oblique Particle-Wall Collisions in a Liquid”, J. Fluid Mechanics, Vol. 510, pp. 71-93. 89 Joseph, G.G., Zenit, R., Hunt, M.L. & Rosenwinkel, A.M. (2001) “Particle–wall collisions in a viscous fluid,” J. Fluid Mech. Vol. 433, pp. 329–346. Kaftori, D., Hetsroni, G. and Banerjee, S. (1995), “Particle behavior in the turbulent boundary layer, part I. motion, deposition and entrainment,” Phys. Fluids A. Vol. 7, pp. 1095-1106. Kallio, G.A. and Reeks, M.W. (1989) “A numerical simulation of particle deposition in turbulent boundary layers," International Journal of Multiphase Flow, Vol. 15, No. 3, pp. 433-446. Kameda, M. & Matsumoto, Y. (1996) “Shock waves in a liquid containing small gas bubbles," Physics of Fluids, Vol. 8, pp. 322-335. Kameda, M., Shimaura, N., Higashino, F. & Matsumoto, Y. (1998) “Shock waves in a uniform bubbly flow," Physics of Fluids, Vol. 10, pp. 2661-2668. Kane, R.S. & Pfeffer, R. (1973) “Heat transfer coefficients of dilute flowing gas-solids suspensions," NASA-CR-2266. Kariyasaki, A. (1987) "Behavior of a single gas bubble in a liquid flow with a linear velocity profile," Proc. of the 1987 ASME-JSME Thermal Engineering Joint Conference, ASME, New York, NY, Vol. 5, pp. 261-267. Kashiwa, B.A. and Gore, R.A. (1987) "A four equation model for multiphase turbulent flow," Proc. of the 1987 ASME-JSME Thermal Engineering Joint Conference, ASME, New York, NY, Vol. 5, pp. 261-267. Kassinos, S.C. and Reynolds, W.C. (1997) "Particle and one-point structure-based modeling of slow deformations of homogeneous turbulence," Proc. of the 11th Symp. on Turbulent Shear Flows, Institut National Polytechnique, University of Joseph Fourier. Keller, J.B. & Kolodner, I. I. (1956) “Damping of underwater explosion bubble oscillations," J. Appl. Physics, Vol. 27, pp. 1152-1161. Kendoush, A.A. (2003) “The virtual mass of a spherical-cap bubble," Physics of Fluids, Vol. 15, No. 9, pp. 2782-2785. Kenning, V.M. and Crowe, C.T. (1997) “On the effect of particles on carrier phase turbulence in gas-particle flows," International Journal of Multiphase Flow, Vol. 23, No. 2, pp. 403-408. Kersey, J., Loth, E. & Lankford, D. (2008) “Effect of evaporating droplets on shock wave attenuation", AIAA Aerospace Sciences Meeting, Reno, NV, AIAA-2008-0793. Kim, J. & Choi, H. (2002) "Laminar flow past a sphere rotating in the streamwise direction," J. Fluid Mech. Vol. 465, pp. 354-386. Kim, I., Elghobashi, S.E. and Sirignano, W. (1993) "Three- dimensional flow over two spheres placed side by side," J. Fluid Mechanics, Vol. 246, pp. 465-488. Kim, I., Elghobashi. S., and Sirignano, W.B. (1998) "On the equation for spherical- particle motion: effect of Reynolds and acceleration numbers," J. Fluid Mech. Vol. 367, pp. 221-253. King, C.J., Hsueh, L. and Mao, K.W. (1965) “Liquid phase diffusion of nonelectrolytes at high dilution,” Journal of Chemical and Engineering Data, Vol. 10, No. 4, pp. 348-350. Khishmatullin, D. (2005) Personal Communication, Duke University. 90 Khismatullin, D.B., Renardy, Y. and Cristini, V. (2003) "Inertia-induced breakup of highly viscous drops subjected to simple shear," Phys. Fluids, Vol. 15, No. 5, pp. 1351 - 1354. Klaseboer, E., Chevaillier, J.-P., Mate, A., Masbernato, O. & Gourdon, C. (2001) “Model and experiments of a drop impinging on an immersed wall,” Phys. Fluids, Vol. 13, No. 1, pp. 45-57. Klebanoff, P.S. (1955) “Characteristics of turbulence in a boundary layer with zero pressure gradient,” NACA Report 1247. Knight, D.D. (2006) Elements of Numerical Methods for Compressible Flows, Cambridge Univ. Press. Koch, D.L. & Hill, R.J. (2001) “Inertial effects in suspension and porous-media flows,” Annual Review of. Fluid Mech. Vol. 33, pp. 619-647. Kochevsky, A. (2005) “Capabilities of numerical simulation of multiphase flows in centrifugal pumps using modern CFD software,” ArXiv physics/0509193, Cornell University, Ithaca, N.Y. Koebe, M., Bothe, D. & Warnacke, H.-J. (2003) “Direct numerical simulation of air bubbles in water glycerol mixtures: shapes and velocity fields,” Proceeding of the 4th ASME-JSME Joint Fluids Engineering Conference, FEDSM2003-45154. Kojima, E., Akehata, T. & Shirai, T. (1968) J. Chem. Engr. Japan, Vol. 1, p. 45. Kok, J.B.W. (1993) “Dynamics of a pair of gas bubbles moving through liquid. Part I. theory,” Eur. J. Mech. B: Fluids, Vol. 12, pp. 515-540. Koren, Lewis, van Brummelen, and van Leer (2002) “Riemann-problem and Level-Set approaches for homentropic two-fluid computations,” Journal of Computational Physics, Vol. 181, pp. 654-674. Kruis, F.E. & Kusters, K.A. (1997), “The collision rate of particle in turbulent flow,” Chemical Engineering Communications, Vol. 158, pp. 201-230. Kubota, M., Akehata, T. and Shirai, T. (1967) “The behavior of single air bubbles in liquids of small viscosity,” Kagaku Kogaku, Vol. 31, p. 1074. Kulick, J.D., Fessler, J.R., and Eaton, J.K. (1994) “Particle response and turbulence modification in fully developed channel flow,” J. Fluid Mech. Vol. 277, pp. 109-134. Kuo, K.K. (1986) Principles of Combustion, John Wiley & Sons, New York. Kurose, R. and Komori, S. (1999) “Drag and lift forces on a rotating sphere on a linear shear flow,” J. Fluid Mech. Vol. 384, pp. 183-206. Kurose, R., Misumi, R. & Komoni, S. (2001) “Drag and lift forces on a spherical bubble in a linear shear flow,” International Journal of Multiphase Flow, Vol. 27, pp. 1247- 1258. Kussin, J. and Sommerfeld, M. (2001) “Investigation of particle behavior and turbulence modification in particle laden channel flow,” International Congress for Particle Technology, Session 12-046 Nuremberg, Germany, 27-29 March. Labous, L., Rosato, A.D. & Dave, R.N. (1997) “Measurements of collisional properties of spheres using high-speed video analysis,” Physical Review E. Vol. 56, pp. 5717- 5725. Lackme, C. (1974) “Two regimes of a spray column in countercurrent flow,” AIChE Symposium Series, Heat Transfer R and D, Vol. 70, pp. 57-62. Lai, R.Y.S. & Mockros, L.F. (1972) “The Stokes flow drag on prolate and oblate spheroids during axial translatory accelerations,” J. Fluid Mech. Vol. 52, pp. 1-15. 91 Lamb (1945) Hydrodynamics, Dover, New York. Lance, M. and Bataille, J. (1991a) “Turbulence in the liquid phase of a uniform bubbly air-water flow,” J. Fluid Mech. Vol. 222, pp. 95-118. Lance, M. and Bataille, J. (1991b) “Homogeneous turbulence in bubbly flows,” J. Fluid Engineering, Vol. 113, pp. 295-300. Langmuir, I. and Blodgett, K.B. (1946) “A mathematical investigation of water droplet trajectories,” Army Air Force Tech. Report 5418, Contract W-33-038-ac-9151. Lasheras, J, private communication, (1998) University of California, San Diego. Lasso, I.A. & Weidman, P.D. (1986) “Stokes drag on hollow cylinders and conglomerates,” Phys. Fluids, Vol. 29, pp. 3921–3934. Launay, K., Huilier, D. and Burnage, H. (1998) "An improved Lagrangian method for predicting the long-time turbulent dispersion in gas-particle flows," ASME Summer Fluids Engineering Meeting, FEDSM98-5012, Washington, D.C., June. Launder, B.E. & Spalding, B. (1972) Mathematical Models of Turbulence, Academic Press, New York. Lawrence, C.J. & Mei, R. (1995) “Long time behavior of the drag on a body in impulsive motion,” J. of Fluids Mechanics, Vol. 283, pp. 307-327. Leal, L.G. (1980) “Particle motions in a viscous fluid,” Annual Review of Fluid Mechanics, Vol. 12, pp. 435-476. Lee, V. & Loth, E. (2007) “Local adaptive timestepping for Lagrangian particle tracking,” AIAA Aerospace Sciences Meeting, Reno NV, January, AIAA-2007-0335. Legendre, D., Daniel, C. & Guiraud, P. (2005) “Experimental study of a drop bouncing on a wall in a liquid,” Physics of Fluids, Vol. 17. Legendre, D. and Magnaudet, J. (1997) “A note on the lift force on a bubble or drop in a low Reynolds number shear flow,” Physics of Fluids, Vol. 9, pp.3572-3574. Legendre, D. and Magnaudet, J. (1998) “Lift force on a bubble in a viscous linear shear flow,” J. of Fluids Mechanics, Vol. 368, pp. 81-126. Legg, B.J. (1983) “Movement of plant pathogens in the crop canopy,” Phil. Tans. Roy. Sci. Long , Vol. B302, pp. 559-574. Legg, B.J. & Raupach, M.R. (1982) “Markov chain simulation of particle dispersion in inhomogeneous flow,” Boundary Layer Meteorology , Vol. 24, pp. 3-13. Leighton, D. & Acrivos, A. (1987) “The shear-induced migration of particles in concentrated suspensions,” Journal of Fluid Mechanics, Vol. 181, pp. 415-439. Leighton, T G. (2004) “From sea to surgeries, from babbling brooks to baby scans: bubble acoustics at ISVR,” Proceedings of the Institute of Acoustics, Vol. 26, Part 1, 357-381. Leith, D. (1987) “Drag on non-spherical objects,” Aerosol. Sci. Tech. Vol. 6, pp. 153-161. L’Esperance, D., Trolinger, J.D., Coimbra, C.F. & Rangel, R.H. (2006) “Particle response to low-Reynolds-number oscillation of a fluid in micro-gravity,” AIAA J. Vol. 44, pp. 1060-1064. Letan, R. and Kehat, E. (1967) “The mechanics of a spray column,” AIChE Journal, Vol. 13, pp. 443-449. Levich, V.G. (1949) ”Motion of a bubble at lagre Reynolds numbers,” Zhur. Eksp. i. Teoret. Fiz, Vol. 19, p. 18. Levich, V.G. (1962) Physico Chemical Hydrodynamic., Prentice Hall. 92 Li, W. and Davis, E. J. (1995) “Measurement of the thermophoretic force by electrodynamic levitation: microspheres in air,” J. Aerosol Sci. Vol. 26, pp. 1063- 1083. Libbrecht, K.G. (2005) “The physics of snow crystals,” Rep. Prog. Phys. Vol. 68, pp. 855-895 Lide, D. edit. (2005) CRC Handbook of Chemistry and Physics, CRC Press, Boca Raton, FL, USA. Liepmann, H.W. and Roshko, A. (1957) Elements of Gasdynamics, John Wiley & Sons, New York. Liger-Belair, G. and Jeandet, P. (2002) “Turbulence modulation in gas-particle flows: a compariosn of selected models,” Europhysics News, Vol. 33, No. 1. Lightstone & Hodgson (2004) “Effervescence in a glass of champagne: a bubble story,” Canadian Journal of Chemical Engineering, Vol. 82, pp. 209-220. Lim, E.A., Coimbra, C.F.M & Kobayashi, M.H. (2005) "Dynamics of suspended particles in eccentrically rotating flows," J. Fluid Mech. Vol. 535, pp. 101-110. Lin, B. Yu, J. and Rics, S.A. (2000) “Direct measurements of constrained Brownian motion of an isolated sphere between two walls,” Physical Review E, Vol. 62, pp. 3909-3918. Lin, C.J., Perry, J.H. & Scholwater, W.R. (1970) “Simple shear flow round a rigid sphere: inertial effects and suspension rheology,” Journal of Fluid Mechanics, Vol. 44, pp. 1-17. Ling, W. Chung, J.N. Troutt, R.T. and Crowe, C.T. (1998) “Direct numerical simulation of a three-dimensional temporal mixing layer with particle dispersion,” J. Fluid Mech. Vol. 358, pp. 61–85. List, R., Retsch, U.W., Byram, A.C. & Lozowski, E.P. (1973) "On the aerodynamics of spheroidal hailstone models," J. Atmos. Sci. Vol. 30, No. 4, pp. 653-661, List, R., & Schemenauer (1971) "Free-fall behavior of planar snow crystals, conical graupel and small hail," J. Atmos. Sci. Vol. 28, pp. 110-115. Loewenberg, M. (1993) “Stokes resistance, added mass, and Basset force for arbitrarily oriented finite-length cylinder,” Physics of Fluids A, Vol. 5, pp. 765-767. Löhner, R. and Ambrosiano, J. (1990) “A vectorized particle tracer for unstructured grids,” Journal of Computational Physics, Vol. 91, No. 1, November, pp. 22-31. Lomholt, S. & Maxey, M. (2003) “Force-coupling method for particulate two-phase flow: stokes flow,” Journal of Computational Physics, Vol. 184, pp. 381-405. Lopez de Bertodano, M., Lahey, R.T. and Jones, O.C. (1994) "Development of a k- model for two-phase flow," AMSE Journal of Fluids Engineering, Vol. 116, pp. 128- 134. Loth, E. (2000) “Numerical approaches for motion of dispersed particles, bubbles, and droplets,” Progress in Energy and Combustion Sciences, Vol. 26, pp. 161-223. Loth, E. (2001) "A Eulerian turbulent diffusion model for particles and bubbles," International Journal of Multiphase Flow, Vol. 27, pp. 1051-1063. Loth, E. (2008a) "“Review: Quasi-Steady Shape and Drag of Deformable Bubbles and Drops” Intl. J. of Multiphase Flow Vol. 34, pp. 523-546. Loth, E. (2008b) “Review: Lift of a Spherical Particle Subject to Vorticity and/or Spin” AIAA J. Vol. 46, pp. 801-809, April. 93 Loth, E. (2008c) “Drag of Non-Spherical Solid Particles of Regular and Irregular Shape” Powder Technology Vol. 182, pp.342-353, March. Loth, E. (2008d) "Compressibility and rarefaction effects on drag of a spherical particle," AIAA Journal, (to appear). Loth, E., Boris, J. & Emery, M. (1998) "Very large bubble cavitation in a temporally- evolving free shear layer," ASME Summer Fluids Engineering Meeting, Washington, D.C., June. Loth, E. and Dorgan, A.J. (2009) “An equation of motion for particles of finite size and Reynolds number in a liquid," Environmental Fluid Mechanics (to appear). Loth, E., March, J. & K. Krishnan (2008) “Modeling drop-drop collision regimes for variable pressures and viscosities” Aerospace Sciences Meeting, AIAA Paper, 2008- 1022, Reno NV. Loth, E. and Faeth, G.M. (1989) "Structure of under-expanded round air jets submerged in water," Intl. J. of Multiphase Flow, Vol. 15, No. 4, pp. 589-603. Loth, E., O’Brien, T.J., Syamlal, M. and Cantero, M. (2004) "Effective diameter for group motion of polydisperse particle mixtures," Powder Technology, Vol. 142, No. 2-3, pp. 209-218. 30 Apr. Loth, E. and Stedl, J. (1999) "Taylor and Lagrangian correlations in a turbulent free shear layer," Experiments in Fluids, Vol. 26, pp. 1-6. Loth, E, Taebi-Rahni, M. and Tryggvason, G. (1997) "Deformable bubbles in a free shear layer," Intl. J. of Multiphase Flow, Vol. 23, No. 56, pp. 977-1001, September. Lovalenti, P.M. & Brady, J.F. (1993) “The force on a sphere in a uniform flow with small-amplitude oscillations at finite Reynolds number,” J. Fluid Mechanics, Vol. 256, pp. 607-614. Lovalenti, P.M. & Brady, J.F. (1995) “The temporal behavior of the hydrodynamics force on a body in response to an abrupt change in velocity at small but finite Reynolds number,” Journal of Fluid Mechanics, Vol. 293, p. 35-46. Loyalka, S. K. (1992) “Thermophoretic force on a single particle I. Numerical solution of the linearized Boltzmann equation,” J. Aerosol Sci. Vol. 23, pp. 291-300. Luo, H. and Svendsen, H.F. (1996) "Theoretical model for drop and bubble breakup in turbulent dispersions,” Amer. Inst. of Chem Engr. J. Vol. 42, pp. 1225-1233. MacCormack, R.W. (1969) "The effect of viscosity in hypervelocity impact cratering," AIAA Paper, 69-354. MacInnes, J.M. and Bracco, F.V. (1992) "Stochastic particle dispersion modeling and the tracer-particle limit," Phys. Fluids, Vol. 4, No. 12, pp. 2809-2824. Maccoll, J.H. (1928) "Aerodynamics of a spinning sphere," J. Royal Aeronautical Society, Vol. 32, pp. 777-798. Macrossan, N. (2004) “Scaling parameters in rarefied flow and the breakdown of the Navier-Stokes equations,” Technical Report No: 2004/09, Department of Mechanical Engineering, University of Queensland, St Lucia 4072, Australia. Madhav, G.V. & Chhabra, R.P. (1995) "Drag on non-spherical particles in viscous fluids," Int. J. Miner. Process, Vol. 43, pp. 15-29. Maeda, M., Hishida, K. and Furutani, T. (1980) “Optical measurements of local gas and particle velocity in an upward flowing dilute gas–solids suspension,” Polyphase Flow Transport Technology, Century 2-ETC, San Francisco, CA. 94 Magnaudet, J. (2003) “Small inertial effects on a spherical bubble, drop or particle moving near a wall in a time-dependent linear flow,” J. of Fluid Mechanics, Vol. 485, pp. 115-142. Magnaudet, J. and Eames, I. (2000) “The motion of high-Reynolds-number bubbles in inhomogeneous flows,” Annual Reviews of Fluid Mechanics, Vol. 32, pp. 659-708. Magnaudet, J., Takagi, S. & Legendre, D. (2003) “Drag, deformation and lateral migration of a buoyant drop moving near a wall,” J. of Fluid Mechanics, Vol. 476, pp. 115-157. Mainardi, F. & Pironi, P. (1996) “The fractional Langevin equation Brownian motion revisited,” Extracta Mathematicae, Vol. 11, pp. 140-154. Marble, E.F. (1970) “Dynamics of dusty gases,” Ann. Rev. Fluid Mech, Vol. 2, pp. 397- 446. Marie, J.L., Moursali, E., and Tran-Cong, S. (1997) "Similarity law and turbulence intensity profiles in a bubbly boundary layer at low void fractions," Intl. J. of Multiphase Flow, Vol. 23, No. 2, pp. 227-247. Martin, G., Loth, E. & Lankford, D. (2009) “Particle Host-Cell Determination in Unstructured Grids,” Computers & Fluids, Vol. 38, pp. 101-110. Martinez-Bazan, C., Montanes, J.L. and Lasheras, J.C. (1999a) "On the breakup of an air bubble injected into a fully developed turbulent flow. Part 1. Breakup frequency,” J. of Fluid Mechanics, Vol. 401, pp. 157-182. Martinez-Bazan, C., Monatnes, J.L., Lasheras, J.C. (1999b) “On the breakup of an air bubble injected into a fully developed turbulent flow. Part 2. Size PDF of the resulting daughter bubbles,” J. of Fluid Mechanics, Vol. 401, pp. 183-207. Mason, J. (1978) “Physics of a raindrop,” Physics Bull. Vol. 29, pp. 364-369. Maude, A.D. and Whitmore, R.L. (1958) “A generalized theory of sedimentation,” Brit. J. Appl. Phys. Vol. 9, pp. 477-482. Maw, N., Barber, J.R. and Fawcett, J.N. (1976) “The oblique impact of elastic spheres,” Wear. Vol. 38, pp. 101-114. Maxey, M.R. (1987) “The gravitational settling of aerosol particles in homogeneous turbulence and random flow fields,” J. Fluid Mechanics, Vol. 174, pp. 441-465. Maxey, M.R. (1993) "Equation of motion for a small rigid sphere in a non-uniform or unsteady flow," ASME/FED Gas-Solid Flows, Vol. 166, pp. 57-62. Maxey, M.R., Chang, E.J. and Wang, L.-P. (1994) "Simulations of interactions between microbubbles and turbulent flows," Appl. Mech Rev. Vol. 47, No. 6, pp. S70-S74. Maxey, M.R., Chang, E.J. and Wang, L.-P. (1996) “Interactions of particles and microbubbles with turbulence,” Experimental Thermal and Fluid Science, Vol. 12, No. 4, pp. 417-425. Maxey, M.R. & Patel, B.K. (2001) "Localized force representations for particles sedimenting in Stokes flow," International Journal of Multiphase Flow, Vol. 27, pp. 1603-1626. Maxey, M.R. (1987) "Gravitational settling of aerosol particles in homogeneous turbulence and random flow fields," Journal of Fluid Mechanics, Vol. 174, pp. 441- 465. Maxey, M.R., Patel, B.K., Chang, E.J. and Wang, L.-P. (1997) "Simulations of dispersed turbulent multiphase flow," Fluid Dynamics Research, Vol. 20, pp. 143-156. 95 Maxey, M.R. and Riley, J.J. (1983) "Equation of motion for a small rigid sphere in a non- uniform flow," Phys. Fluids, Vol. 26, No. 4, pp. 883-889. McComb, W.D. and Salih, S.M. (1977) “Measurement of normalised radial concentration profiles in a turbulent aerosol jet, using a laser-doppler anemometer,” Journal of Aerosol Science, Vol. 8, No. 3, pp. 171-181. McComb, W.D. and Salih S. M. (1978), “The calculation of three dimensional turbulent flow jets,” Proceedings of the symposium on Turbulent Shear Flows, Penn. State University. McLaughlin, M.H. (1968) “An experimental study of particle-wall collision relating of flow of solid particles in a fluid,” Engineer's degree thesis, California Institute of Technology, Pasadena, California. McLaughlin, J.B. (1991) "Inertial migration of a small sphere in linear shear flows," Journal of Fluid Mechanics, Vol. 224, pp. 261-274. McLaughlin, J.B. (1996) “Numerical simulation of bubble motion in water,” Journal of Colloid and Interface Science, Vol. 184, pp. 614–625. McNamara, S. and Falcon, E. (2005) “Simulations of vibrated granular medium with impact-velocity-dependent restitution coefficient,'' Phys Rev E, Vol. 71. Mei, R. (1992) "An approximate expression for shear lift force on a spherical particle at a finite Reynolds number," Intl. J. of Multiphase Flow, Vol. 18, pp. 145-147. Mei, R. & Adrian, R.J. (1992) "Flow past a sphere with an oscillation in the free-stream and unsteady drag at finite Reynolds number," J. Fluid Mech. Vol. 237, pp. 323-341. Mei, R. and Klausner, J.F. (1992) "Unsteady force on a spherical bubble with finite Reynolds number with small fluctuations in the free-stream velocity," Physics of Fluids A, Vol. 4, No. 1, pp. 63. Mei, R. (1994) "Effect of turbulence on the particle settling velocity in the nonlinear drag regime," Int. J. Multiphase Flow, Vol. 20, pp. 273-284. Mei, R. and Klausner, J.F., (1994) "Shear lift force on spherical bubbles," Int. J. Heat and Fluid Flow, Vol. 15, No. 1, pp. 62-65. Mei, R., Klausner, J.F. & Lawrence, C.J. (1994) "A note on the history force on a spherical bubble at finite Reynolds number," Physics of Fluids, Vol. 6, No. 1, pp. 418-420. Mei, R., Lawrence, C.J., and Adrian, R.J. (1991) "Unsteady drag on a sphere at finite Reynolds number with small fluctuations in the free-stream velocity," J. Fluid Mech. Vol. 233, pp. 613-631. Mendelson, H.D. (1967) "The prediction of bubble terminal velocities from wave theory," AIChE J. Vol. 13, pp. 250-253. Menter, F.R. (1993) “Zonal two-equation k-w turbulence models for aerodynamic flows,” AIAA Paper 93-2906. Mercier, J., Lyrio, A. and Forslund, R. (1973) “Three-dimensional study of the nonrectilinear trajectory of air bubbles rising in water,” Journal of Applied Mechanics, Transactions ASME, Vol. 40, Ser. E, No. 3, pp. 650-654. Merle, A., Legendre, D. & Magnaudet, J. (2005) “Forces on a high–Reynolds–number spherical bubble in a turbulent flow,” J. Fluid Mech. Vol. 532, pp. 53–62. Michaelidis, E. (1997) “Review – the transient equation of motion for particles, bubbles, and droplets,” AMSE Journal of Fluids Engineering, Vol. 119, pp. 233-247. 96 Michaelides, E. (2003) “Hydrodynamic force and heat/mass transfer from particles, bubbles and drops––the Freeman scholar lecture,” J. Fluids Eng. Vol. 125, pp. 209– 238. Michaelidis, E. (2006) Particles, Bubbles, and Droplets: Their Motion, Hear and Mass Transfer, World Scientific Publishing, Hackensack, NJ. Miller, R.S. and Bellan, J. (1999) “Direct numerical simulation of a confined three- dimensional gas mixing layer with one evaporating hydrocarbon-droplet-laden stream,” Journal of Fluid Mechanics, Vol. 384, pp. 293-338. Millies, M. & Mewes, D. (1999) “Interfacial area density in bubbly flow,” Chemical Engineering and Processing, Vol. 38, pp. 307-319. Millikan, R.A. (1911) “The isolation of an ion: a precision measurement of its charge, and the correction of Stokes’s Law,” Physical Review, Vol. 32, p. 349-397. Millikan, R.A. (1923) “Coefficients of slip in gasses and the law of reflection of molecules from the surface of solids and liquids,” Physical Reviews, Vol. 21, pp. 217- 238. Minnaert, M. (1933) “On musical air-bubbles and the sounds of running water,” Phil Mag. Vol. 16, pp. 235-248. Mishin, G.I. (1997) “Experimental investigation of the flight of a sphere in weakly ionized air,” Applied Aerodynamics Conference, AIAA-1997-2298, Atlanta, GA, June 23-25. Mito, Y. & Hanratty, T.J. (2002) “Use of a modified Langevin equation to describe turbulent dispersion of fluid particles in a channel flow,” Flow, Turbulence and Combustion, Vol. 68, pp. 1-26. Miyahara, T. & Takahashi, T. (1985) “Drag coefficient of a single bubble rising through a quiescent liquid,” International Chemical Engineering, Vol. 25, pp. 146-148. Moore, D.W. (1963) “The boundary layer on a spherical gas bubble,” J. Fluid Mechanics, Vol. 16, pp. 161-176. Moore, D.W. (1965) “The velocity rise of distorted gas bubbles in a liquid of small viscosity,” J. Fluid Mechanics, Vol. 23, pp. 749-766. Moorman, R.W. (1955) “Motion of a spherical particle in the accelerated portion of free- fall,” Doctor of Philosophy Dissertation, University of Iowa. Moshaii, A., Sadighi-Bonabi, R. & Taeibi-Rahni, M. (2004), “Effects of bulk viscosity in non-linear bubble dynamics,” J. of Physics: Condensed Matter, Vol. 16, pp. 1687- 1694. Mudde, R.F. and Saito, T. (2001), “Hydrodynamical similarities between bubble column and bubbly pipe flow,” J. Fluid Mech. Vol. 437, pp. 203-228. Mudde, R.F. and Simonin, O. (1999), “Two- and three-dimensional simulations of a bubble plume using a two-fluid model,” Chemical Engineering Science, Vol. 54, pp. 5061-5069. Mugele, R. and Evans, H.D. (1951) “Droplet size distribution in sprays,” Ind. Engrg. Chem. Vol. 43, pp.1317-1324. Murai, Y. & Matsumoto, Y. (1999) “Eulerian analysis of bubbly two-phase flow using CIP scheme,” Computational Fluid Dynamics J. Vol. 8, pp. 26-33. Najjar, F. (2005) University of Illinois at Urbaba-Champaign, personal communication. 97 Narciri, M.A. (1992) “Contribution a l’etude des forces exercees par un liquide sur une bulle de gaz: portance, masse ajoutee et interactions hydrodynamiques,” Ph.D. thesis, L’ecole Centrale De Lyon, Lyons. Nathan, A.M., Hopkins, J., Chong, L. & Kaczmarski, H. (2006) “The effect of spin on the flight of a baseball,” International Sports Engineering Conference, Munich, July. Nelson, E.M., Eppley, K.R., Petillo, J.J. and Levush, B. (2004) “Particle tracking on unstructured grids,” Proceedings of the 2001 Particle Accelerator Conference, Vol. 4, pp. 3057-3059, June. Neve, R.S. & Jaafar, F.B. (1982) “The effect of turbulence and surface roughness on the drag of spheres in thin jets,” The Aeronautical Journal, Vol. 86, pp. 331-336. Neve, R.S. & Shansonga, T. (1989) “The effects of turbulence characteristics on sphere drag,” Intl. J. of Heat & Fluid Flow, Vol. 10, pp. 318-321. Niazmand, H. & Renksizbulut, M. (2003) "Surface effects on transient three-dimensional flows around rotating spheres at moderate Reynolds numbers,” Computers & Fluids, Vol. 32, pp. 1405-1433. Nichols, R.H. and Nelson, C.C. (2003) “Application of hybrid RANS/LES turbulence models,” 41st Aerospace Sciences Meeting and Exhibit, AIAA-2003-083, Reno, NV. Nielsen, P. (2007) "Mean and variance of the velocity of solid particles in turbulence,” Particle Laden Flow, Vol. 11, ERCOFTAC Series, Eds. H. Clercx et al. Springer,. Nihei and Nadaoka (1988) "Improvement of GAL model for convective transport simulation,” International Conference on Multiphase Flow, Lyon, France, June. Niro A/S (2005), personal communication. Niven, R.W., Lott, D.F., Ip, A.Y., Somaratne, K.D. and Kearney, M. (1994) "Development and use of an in vitro system to evaluate inhaler devices," Int. J. Pharmaceutics, Vol. 101, pp. 81-87. Oberbeck, A. (1876) Crelles J.,Vol. 81, pp. 62. O’Brien, J.A. (1990) “Brownian dynamics simulation of very small particles: point source in uniform flow,” Journal of Colloid and Interface Science, Vol. 134, pp. 497- 521. Odar, F. (1962) “Forces acting on a sphere accelerating in a viscous fluid,” Doctor of Philosophy Dissertation, Northwestern University. Odar, F. (1964) “Forces on a sphere accelerating in a viscous fluid,” US Army Cold Regions Research and Engineering Laboratory, Research Rep. 128. Odar, F. & Hamilton, W.S. (1964) “Forces on a sphere accelerating in a viscous fluid,” J. Fluid Mech. Vol. 18, pp. 302-314. Oesterle, B. and Zaichik, L.I. (2004) “On Lagrangian time scales and particle dispersion modeling in equilibrium turbulent shear flows,” Physics of Fluids, Vol. 16, No. 9, pp. 3374-3384. Oguz, H.N. & Prosperetti, A. (1990) “Bubble entrainment by the impact of drops on liquid surfaces,” Journal of Fluid Mechanics, Vol. 219, pp. 143-179. Ohta, M., Haranaka, S., Yoshida, Y. and Sussman, M. (2004), “Three-dimensional numerical simulations of the motion of a gas bubble rising in viscous liquids,” J. Chem. Eng. Japan, Vol. 37, No. 8, pp. 968-975. Olim, Igra, Mond, and Ben-Dor (1995) “A general attenuation law of planar shock waves propagating into dusty gases,” Proceedings of the 16th International Symposium on Shock Waves And Shock Tubes, 1985, pp. 217-225. 98 O’Neill, M.E. (1964) “A slow motion of viscous liquid caused by a slowly moving solid sphere,” Mathematika, Vol. 121, pp. 67-74. Oran, E.S. and Boris, J.P. (1987) Numerical Simulation of Reactive Flow, Elsevier, New York. Orszag, S.A. (1980) “Spectral methods for problems in complex geometries,” J. Comput. Phys. Vol. 37, pp. 70–92. Oseen, C.W. (1910) “Uber die Stokessche Formel und uber die verwandte Aufgabe in der Hydrodynamik,” Arkiv Math. Astron. Fys. Vol. 6, No. 29. Oseen, C.W. (1927) Hydrodynamik, Leipzig: Akademische Verlagsgesellschaft. Osher, S. and Sethian, J.A. (1988) “Fronts propagating with curvature-dependent speed: algorithms based on Hamilton-Jacobi formulations,” Journal of Computational Physics, Vol. 70, pp. 12-49. Osiptsov, A.N. (1988) “Modified Lagrangian method for calculating the particle concentration in dusty-gas flows with intersecting particle trajectories” International Conference on Multiphase Flow, Lyon, France, June, Paper 236. Ounis, H., Ahmadi, G. & McLaughlin, J.B. (1991) “Brownian diffusion of submicrometer particles in viscous sublayer,” Journal of Colloid and Interface Science, Vol. 143, pp. 266-277. Pal, S., Merkle, C.L. and Deutsch, S. (1988) “Bubble characteristics in a microbubble boundary layer,” Physics of Fluids, Vol. 31, pp. 774-751. Pan, F. and Acrivos, A. (1968) "Shape of a drop or bubble at low Reynolds number," I and EC Fund. Vol. 7, pp. 227. Mashayek, F. & Pandya, R.V.R. (2003) “Analytical description of particle/droplet-laden turbulent flows ,” Progress in Energy and Combustion Science, Vol. 29, pp. 329-378 . Parthasarathy, R.N. and Faeth, G.M. (1990) “Turbulence modulation in homogeneous dilute particle-laden flow,” J. Fluid Mech. Vol. 220, pp. 485-537. Park, W.C., Klauner, J.F. & Mei, R. (1995) “Unsteady forces on spherical bubbles,” Experiments in Fluid. Vol. 19, pp. 167-172. Patankar, N.A. and Joseph, D.D. (2001) “Modeling and numerical simulation of particulate flows by the Eulerian-Lagrangian approach,” Int. J. Multiphase Flow, Vol. 27, pp. 1659 – 1684. Patankar and Karki (1999) “Calculation of particle trajectories using complex meshes,” Numerical Heat Transfer, Part B, Vol. 35, pp. 431-437. Peregrine, D.H. and Thais, L. (1996) "The effect of entrained air in violent water wave impacts," Journal of Fluid Mechanics, Vol. 325, pp. 377-397, Oct. Peskin, C.S. (1977) “Numerical analysis of blood flow in the heart,” Journal of Computational Physics, Vol. 25, pp. 220-252. Pettyjohn, E.S. & Christiansen, E.B. (1948) “Effect of particle shape on free settling rates of isometric particles,” Chem. Eng. Progr. Vol. 4, pp. 157-172. Phillips, W.F. (1975) Physics of Fluids, Vol. 18, pp. 1089-1093. Pilch, M. and Erdman, C.A. (1987), “Use of break-up time data to predict the maximum size of stable fragment for acceleration induced breakup of a liquid drop,” Int. J. Multiphase Flow, Vol. 13, pp. 741–757. Piomelli, U. (1997) "Introduction to the modeling of turbulence: large eddy and direct simulation of turbulent flows," von Karman Institute for Fluid Dynamics, Lecture Series 1997-03, March. 99 Plesset and Prosperetti (1977) “Bubble dynamics and cavitation,” Annual Reviews of Fluid Mechanics, Vol. 9, pp. 145-185. Poe, G.G. & Acrivos, A (1975) "Closed-streamline flows past rotating single cylinders and spheres: inertia effects," Journal of Fluid Mechanics, Vol. 72, pp. 605-623. Pope, S.B. (2000) Turbulent Flows, Cambridge University Press, Cambridge, U.K. Poorte, R.E.G. & Biesheuvel, A. (2002) "Experiments on the motion of gas bubbles in turbulence generated by an active grid,” Journal of Fluid Mechanics, Vol. 461, pp. 127-154. Pozorski, J. & Minier, J.-P. (1998) “On the Lagrangian turbulent dispersion models based on the Langevin equation,” International Journal of Multiphase Flow, Vol. 24, pp. 913-945. Prosperetti. A. (1987) "The equation of bubble dynamics in a compressible liquid," Physics of Fluids, Vol. 30, pp. 3626. Prosperetti, A., (1998) "Ensemble averaging techniques for disperse flows," Particulate Flows: Processing and Rheology, Editors: Drew, D.A., Joseph, D.D. and Passman, S.L., Springer, New York, pp. 99-136. Prosperetti, A., (2007) “Averaged equations for multiphase flow,” Computational Methods for Multiphase Flow, Editors: Prosperetti, A. & Tryggvason, G., Cambridge University Press, Cambridge, U.K. Prosperetti, A., & Tryggvason, G., (2007) Computational Methods for Multiphase Flow, Cambridge University Press, Cambridge, U.K. Prosperetti, A., Sundaresan, Pannala, S. & Zhang, D.Z. (2007) “Segregated methods for two-fluid models” Computational Methods for Multiphase Flow, Editors: Prosperetti, A. & Tryggvason, G., Cambridge University Press, Cambridge, U.K. Proudman, I. & Pearson, J.R.A. (1957) "Expansions at small Reynolds number for the flow past a sphere and a cylinder," J. Fluid Mechanics, Vol. 2, pp. 237-262. Proudman, I. (1969) “On the flow past a sphere at low Reynolds number,” J. Fluid Mech. Vol. 37, pp.759-760. Putnam, A. (1961) “Integrable form of droplet drag coefficient,” American Rocket Society Journal, Vol. 31, Oct., pp. 1467-1468. Qian, J. and Law, C.K. (1997) "Regimes of coalescence and separation in droplet collision," Journal of Fluid Mechanics, Vol. 331, pp. 59-80. Qiand, D. (2003) "Bubble motion, deformation, and breakup in stirred tanks," Ph.D. Dissertation, Chemical Engineering, Clarkson University. Ramirez, L.E.S, Lim, E.A., Coimbra, C.F.M & Kobayashi, M.H. (2003), "On the dynamics of a spherical scaffold in rotating bioreactors," Biotechnology and Bioengineering, Vol. 84, pp. 382-389. Rani, S.L. & Balachandar, S. (2003) “Evaluation of the equilibrium Eulerian approach for the evolution of particle concentration in isotropic turbulence,” International Journal of Multiphase Flow, Vol. 29, pp. 1793-1816. Range, K. and Feuillebois, F. (1998) “Influence of surface roughness on liquid drop impact,” Journal of Colloid and Interface Science, Vol. 203, pp. 16-30. Ranz, W.E. & Marshall, W.R. (1952) “Evaporation from drops – part II,” Chemical Engineering Progress, Vol. 48, pp. 141-146. 100 Raymond, F. & Rosant, J.-M. (2000) "A numerical and experimental study of the terminal velocity and shape of bubbles in viscous liquids," Chem. Engr. Sci. Vol. 55, pp. 943-955. Reeks, M.W. (1977) "On the dispersion of small particles suspended in an isotropic turbulent fluid," Journal of Fluid Mechanics, Vol. 83, Part 3, pp. 529-546. Reeks, M.W. (1993) "On the constitutive relations for dispersed particles in non-uniform flows, 1: dispersion in a simple shear flow," Physics of Fluids A, Vol. 5, pp. 750-761. Reeks, M.W. and McKee, S. (1984) "The dispersive effects of Basset history forces on particle motion in a turbulent flow," Physics of Fluids, Vol. 27, pp. 1573-1582. Rein, Martin (1993) “Phenomena of liquid drop impact on solid and liquid surfaces,” Fluid Dynamics Research, Vol. 12, pp.61-93. Reinhart, A. (1964) "Das Verhalten fallender Topfen," Chemie-Ing.-Techn., Vol. 36, pp. 740-746. Reynolds, A.M. and Lo Iacono, G. (2004) “On the simulation of particle trajectories in turbulent flows,” Physics of Fluids, Vol. 16, pp. 4353-4358. Reynolds, W. (1989) “The potential and limitations of direct and large eddy simulations,” Whither Turbulence Conference, Cornell University, Ithaca, NY, March 22-24. Richard, D. & Quere, D. (2000) “Bouncing water drops,” Europhys. Lett. Vol. 50, pp. 769–775. Richardson, J. F. and Zaki, W. N. (1954a) “Sedimentation and fluidization: part I,” Trans. Instn. Chem. Engrs. Vol. 32, pp. 35-53. Richardson, J. F. and Zaki, W. N. (1954b) “The sedimentation of a suspension of uniform spheres under conditions of viscous flow,” Chem. Engr. Sci. Vol. 8, pp. 65-73. Risso, F. (1999) "Experimental investigation of the motion of a bubble in a gradient of turbulence," Physics of Fluids, Vol. 11, pp. 3567-3569. Risso, F. and Fabre, J. (1998) "Oscillations and break-up of a bubble immersed in a turbulent field," Journal of Fluid Mechanics, Vol. 372, pp. 323-355. Rodi, W. (1980) Turbulence Models and Their Applications in Hydraulics, International Association of Hydraulic Research, Delft, Netherlands. Roessler (1982) “Diesel particle mass concentration by optical techniques,” Applied Optics, pp. 4077-4086, Vol. 21, No. 22, November. Rogers, C. (1999) Tufts University, personal communication. Rosin, P. and Rammler, E. (1933) “The laws governing the fineness of powdered coal,” Institute of Fuel, Vol. 7, pp. 29-36. Rowe, P. N. (1987) “A convenient empirical equation for estimation of the Richardson- Zaki exponent,” Chem. Engng Sci. Vol. 42, pp. 2795-2796. Rubinow, S.I and Keller, J.B. (1961) "The transverse force on spinning spheres moving in a viscous liquid," Journal of Fluid Mechanics, Vol. 11, No. 3, pp. 447. Rudinger, G. (1964) "Some properties of shock relaxation in gas flows carrying small particles," Physics of Fluids, Vol. 7, No. 5, pp. 658-663. Russel, W.B. (1981) "Brownian motion of small particles suspended in liquids," Annual Review of Fluid Mechanics, Vol. 13, pp. 425-455. Rybalko, M., Loth, E. & Lankford, D. (2008) "Lagrangian sub-grid particle diffusion for LES/RANS flows,” ASME Fluids Engineering Division Summer Meeting, FEDSM2008-55207, Jacksonville, FL. 101 Ryskin, G. and Leal, L.G. (1984) "Numerical simulation of free-boundary problems in fluid mechanics. Part 3, bubble deformation in an axisymmetric straining flow," Journal of Fluid Mechanics, Vol. 148, No.37. Sabnis, J.S. and de Jong, F.J. (1990) "Calculation of the two-phase flow in an evaporating spray using an Eulerian-Lagrangian analysis," AIAA Aerospace Sciences Meeting, Reno NV, AIAA-90-0447, Jan. Sadhal, S.S. & Johnson, R.E. (1983) “Stokes flow past bubbles and drops partially coated with thin films. Part 1, stagnant cap of surfactant film—exact solution,” Journal of Fluid Mechanics, Vol. 126, pp. 237-250. Saffman, P.G. (1956) “ On rise of small air bubbles in water,” Journal of Fluid Mechanics, Vol. 1, pp. 249-275. Saffman, P.G. (1965) "The lift on a sphere in slow shear flow,” Journal of Fluid Mechanics, Vol. 22, pp. 385-400; Corrigendum, Vol. 31, 1968, pp. 624. Saffman, P.G. and Turner, J.S. (1956) “ On the collision of drops in turbulent clouds,” Journal of Fluid Mechanics, Vol. 1, pp. 16-30. Saito, S. (1913) “On the shape of the nearly spherical drop which falls through a viscous fluid,” Science Rep. Tohuku Imp. Univ. Sendai, Japan, Vol. 2, pp. 179. Salem, M.B. & Osterlé, B. (1998) “A shear flow around a spinning sphere: numerical study at moderate Reynolds numbers,” Intl. J. Multiphase Flow, Vol. 24, pp. 563-585. Sangini, A.S. & Didwania, A.K. (1993) "Dynamic simulations of bubbly liquids at large Reynolds number," Journal of Fluid Mechanics, Vol. 250, pp. 307-337. Sangani, A.S., Zhang, D. Z. and Prosperetti, A. (1991) “The added mass, Basset, and viscous drag coefficients in non-dilute bubbly liquids undergoing small-amplitude oscillatory motion,” Phys. Fluids A. Vol. 3, pp. 2955-2970. Sankaran, V. and Menon, S. (2002) “Vorticity-scalar alignments and small-scale structures in swirling spray combustion,” Proceedings of the Combustion Institute, Vol. 29, No. 1, pp. 577-583 Sankaranarayanan, K., Shan, X., Kevrekidid, I.G & Sundaresan, S. (2003) “Analysis of drag and virtual mass forces in bubble suspension using an implicit formulation of the Lattice Boltzmann Method,” J. Fluid Mechanics, Vol. 452, pp. 61-966. Sano (1981) “Unsteady flow past a sphere at low Reynolds number,” J. Fluid Mechanics, Vol. 112, pp. 443-441. Saurel, R. Forrestier, A., Veyret, D. & Loraud, J.-C. (1994) “A finite volume scheme for two-phase compressible flows,” Intl. J. for Numer. Methods in Fluids, Vol. 18, pp. 803-819. Sawford, B. L. (1991) “Reynolds number effects in Lagrangian stochastic models of turbulent dispersion,” Phys. Fluids A, Vol. 3, pp. 1577–1586. Sawicki, G. S., Hubbard, M. and Stronge, W. (2003) “How to hit home runs: optimum base-ball bat swing parameters for maximum range trajectories,” Am. J. Phys. Vol. 71, pp. 1152-1162. Saxton, R.L. and Ranz, W.E. (1952) “Thermal force on an aerosol particle in a temperature gradient,” J. Appl. Phys. Vol. 23, pp. 917-923. Schaaf, S.A. & Chambre, P.L. (1958) “The flow of rarefied gases,” Fundamentals of Gas Dynamics, Sec. H, H.W. Emmons, Ed. Schadt, C.F. and Cadle, R.D. (1961) “Thermal forces on aerosol particles,” J. Phys. Chem. Vol. 65, pp. 1689-1694. 102 Schiller, L. and Naumann, A.Z. (1933) “Über die grundlegenden Berechungen bei der Schwerkraftaufbereitung,”Ver. Deut. Ing, Vol. 77, ppp. 318-320. Schlichting, H. (1979) Boundary Layer Theory, Springer-Verlag, Berlin Heidleberg, Germany. Schlichting, H. & Gresten, G. (2000) Boundary Layer Theory, ninth ed., Springer-Verlag, New York. Schmitt, K. H. (1959) “Untersuchungen an Schwebstoffteilchen im Temperaturfeld,” Z. Naturforsch Vol. 14a, pp. 870-881. Schneider, F. & Merzkirch, W. (2001) "Distribution of coherent structures in turbulent pipe flow," J. of Flow Visualization and Image Processing, Vol. 8, pp. 253-261. Schuen, J-S., Chen, L-D. & Faeth, G.M. (1983) “Evaluation of a stochastic model of particle dispersion in a turbulent round jet,” AIChE Journal, Vol. 29, pp. 396-404. Sene, K.J., Hunt, J.C.R. and Thomas, N.H. (1994) "The role of coherent structures in bubble transport by turbulent shear flows," Journal of Fluid Mechanics, Vol. 259, pp. 219-240. Sevik, M. & Park, S.H. (1973) “The splitting of drops and bubbles by turbulent fluid flow,” Trans. ASME: J. Fluids Engng, March, pp. 53-60. Shirolkar, J.S., Coimbra, C.F.M and Quirez McQuay, M. (1996) “Fundamental aspects of modeling turbulent particle dispersion in dilute flows,” Progress in Energy and Combustion Science, Vol. 22, pp. 363-399. Shoemaaker, J. (2005), Tulas photography company, http://www.tulsaphotographer.com/pages/indk.html Shotorban, B., Afshari, A., Jaberi, F. & Mashavek, F. (2004) “A particle-tracking algorithm for LES of two-phase flow,” AIAA-2004-332, 42nd AIAA Aerospace Science Meeting and Exhibit, Reno, Nevada. Shrayber, A.A. (1979) “Euler and Lagrange Methods in the theory of two-phase flows with variable particle size of the discrete phase,” Fluid Mechanics - Soviet Research, Vol. 8, No. 5, Sept.-Oct., pp. 79-87. Shyy, W., Thakur, S.S., Ouyang, H., Liu, J. & Blosch, E. (1997) Computational Techniques for Complex Transport Phenomena, Cambridge University Press. Sigurgeirsson, H., Stuart, A.M. & Wan, J. (2001) “Collision detection for particles in a flow,” Journal of Computational Physics, Vol. 172, pp. 766–807. Šikalo, S., Tropea, C. and Ganić, E.N. (2005) “Impact of droplets onto inclined surfaces,” Journal of Colloid and Interface Science, Vol. 286, pp. 661-669. Simonnet, M. Gentric, C., Olmos, E. & Midoux, N. (2007) "Experimental determination of the drag coefficient in a swarm of bubbles," Chemical Engineering Science, Vol. 62, pp. 858-866. Simpkins, P.G. and Bales, E.L. (1972) "Water drop response to sudden accelerations," Journal of Fluid Mechanics, Vol. 55, pp. 629-639. Sirignano, W.A. (1993) "Fluid dynamics of sprays - 1992 Freeman Scholar Lecture," ASME Journal of Fluids Engineering, Vol. 115, pp. 345-378. Sirignano, W.A. (1999) Fluid Dynamics and Transport of Droplets and Sprays, Cambridge University Press. Sivier, S.A., Loth, E., and Baum, J.D. (1996) "Dusty shock flow with unstructured adaptive finite elements and parcels," AIAA Journal, Vol. 34, No. 5, pp. 1078-1080. 103 Sivier, S.A., Loth, E., Baum, J.D. and Lohner, R. (1994) "Unstructured adaptive remeshing finite element method for dusty shock flows," Shock Waves, Intl. J. Vol. 4, No. 1, pp. 31-41. Slater, S.A., Leeming, A.D. & Young, J.B. (2003) “Particle deposition from two- dimensional turbulent gas flows”, International Journal of Multiphase Flow, Vol. 29, pp. 721-750. Slater, S.A. and Young, J.B. (1988) “The calculation of inertial particle transport using an Eulerian formulation”, Proceedings of the 16th ASME Fluids Engineering Summer Mtg, June, Washington, D.C., FEDSM98-5007. Smith, D.A. & Cheung, K.F. (2003) “Settling characteristics of calcareous sand,” Journal of Hydraulic Engineering, June, pp. 479-483. Smolouchowski, M. (1916) “Drei vortrageuber diffusion, brownsche bewegung und koagulation von kolloidteilchen,” Phys. Zeitschr. Vol.17, pp. 557-559. Snider, D.M., P.J. Rourke, and Andrews, M.J. (1998) “Sediment flow in inclined vessels calculated using a multiphase particle-in-cell model for dense particle flows,” International Journal of Multiphase Flow, Vol. 24, pp. 1359-1382. Snyder, W.H. and Lumley, J.L. (1971) "Some measurements of particle velocity autocorrelation in approximately isotropic turbulence," Journal of Fluid Mechanics, Vol. 48, pp. 41. Sommerfeld, M. (1987) “Numerical simulation of supersonic two-phase gas-particle flows,” Proceedings of the 16th International Symposium on Shock Tubes and Waves, pp. 235-241. Sommerfeld, M. (1985) “The unsteadiness of shock waves propagating through gas- particle mixtures,” Experiments in Fluids, Vol. 3, pp.197-206. Sommerfeld, M. (2001) “Validation of a stochastic Lagrangian modeling approach for inter-particle collisions in homogeneous isotropic turbulence,” International. Journal of Multiphase Flow, Vol. 27, pp. 1829-1858. Sommerfeld, M. (2002) “Kinetic simulations for analyzing the wall collision process of non-spherical particles,” Proceedings of ASME Fluids Engineering Division Summer Meeting FEDSM2002-31239, pp. 539-547, ASME, NY. Sommerfeld, M. & Huber, N. (1999) “Experimental analysis and modeling of particle- wall collisions,” Int. J. Multiphase Flow, Vol. 25, pp. 1457-1489. Sone, W. & Aoki, K. (1983) “Forces on a spherical particle in a slightly rarefied gas,” Progress in Aeronautics and Astronautics, Rarefied Gas Dynamics, AIAA, New York, Vol. 51, pp. 417-433. Soo, S.L. (1990) Multiphase Fluid Dynamics, Gower Technical, Aldershot-Brookfield, USA. Spalart, P.R. (1999) “Strategies for turbulence modeling and simulations,” 4th Intl. Symp. On Engr. Turb. Modeling ,and Measurements, Corsica, France, May 24-26. Spalart, P.R. & Allmaras, S.R. (1994) "A one-equation turbulence model for aerodynamic flows," La Recherche Aerospatiale, Vol. 1, p. 521. Spalart, P.R., Jou, W.H., Stretlets, M. & Allmaras, S.R. (1997) "Comments on the feasibility of LES for wings and on a hybrid RANS/LES approach," Advances in DNS/LES, ed. C. Liu and Z. Liu, Greyden Press, Columbus, pp. 137-148. Spelt, P.D.M. and Biesheuvel, A. (1997) "On the motion of gas bubbles in homogeneous isotropic turbulence," Journal of Fluid Mechanics, Vol. 336, pp. 221-244. 104 Speziale, C.G., Sarkar, S. and Gatski, T.B (1991) "Modeling the pressure-strain correlation of turbulence: an invariant dynamical systems approach," Journal of Fluid Mechanics, Vol. 227, pp. 245-272. Squires, K.D. (2007) "Point particle models," Computational Methods for Multiphase Flow, Edited by Prosperetti, A. and Tryggvason and G., Cambridge University Press, Cambridge, U.K., pp. 282-319. Squires, K. D. and Eaton, J. K. (1990) “Particle response and turbulence modification in isotropic turbulence,” Phys. Fluids A. Vol. 2, pp. 1191-1203. Sridhar, G. and Katz, J. (1995) "Drag and lift forces on microscopic bubbles entrained by a vortex," Phys. Fluids, Vol. 7, No. 2, pp. 389-399. Stadler, J.R. & Zurick, V.J. (1951) “Theoretical aerodynamic characteristics of bodies in free-molecule flow field,” NACA TN 2423, July, pp.12. Stickel, J.J. and Powell, R.L. (2005) “Fluid mechanics and rheology of dense suspensions,” Annual Reviews of Fluid Mechanics, Vol. 37, pp. 129-149. Stock, D.E. (1996) "Particle dispersion in flowing gases - 1994 Freeman Scholar Lecture," ASME Journal of Fluids Engineering , Vol. 118, pp. 4-17. Stokes, G.G. (1851) "On the effect of the inertial friction of fluids on the motion of pendulums," Trans. Cambr. Phil. Soc., Vol. 9 (part II), pp. 8-106 Stone, H.A. and Samuel, A.D.T (1996) "Propulsion of microorganisms by surface distortions," Physical Review Letter, Vol. 77, No. 19, pp. 4102- 4104, November 4. Stout, J.E., Arya, S.P. & Genikhovich, E.L. (1995) “Effect of nonlinear drag on the motion and settling of heavy particles,” J. of Atmospheric Sciences, Vol. 52, pp. 3836-3848. Stringham, G.E., Simons, D.B. and Guy H. P. (1969), “The behavior of large particles falling in quiescent liquids,” Prof. Pap. US Geol. Surv. 562-C. Subia, S.R., Ingber, M.S., Mondy, L.A., Altobelli, S.A. and Graham, A.L. (1998) “Modeling of concentrated suspensions using a continuum constitutive equation,” J. Fluid Mech. Vol. 373, pp. 193–219. Sugiyama, K., Takagi, S. & Matsumoto, Y. (2001) “Multi-scale analysis of bubbly flows,” Comput. Methods Appl. Mech. Engng. Vol. 191, pp. 689–704. Sun, T.-Y., and Faeth, G.M. (1986) "Structure of turbulent bubbly jets," Intl. J. Multiphase Flow, Vol. 12, pp. 115. Sundaram, S. & Collins, L.R. (1996) “Numerical considerations in simulating a turbulent suspension of finite-volume particles,” Journal of Computational Physics, Vol. 124, pp. 337-350. Sundaram, S. and Collins, L. (1997) “Collision statistics in an isotropic particle-laden turbulent suspension, part 1: direct numerical simulations,” J. Fluid Mech. Vol. 335, pp. 75-109. Suslick, K.S. and Flannigan (2005), “Plasma formation and temperature measurement during single-bubble cavitation,” Nature, March, Vol. 434, pp. 52-55. Suslick, K. S. & Price, G. J. (1999) “Applications of ultrasound to materials chemistry,” Annual Review of Materials Science, Vol. 29, pp. 295-326. Sussman (2003) "A second-order coupled level set method and volume-of-fluid method for computing growth and collapse of vapor bubbles," Journal of Computational Physics, 187, pp. 110-136. 105 Sussman, M., Almgren, A.S., Bell, J.B., Colella, P., Howell, L.H. and Welcome, M.L. (1999) “An adaptive level set approach for incompressible two-phase flows,” J. Computational Physics, Vol. 148, pp. 81-124. Sussman, M., Smereka, P. and Osher, S.J. (1994) “A Level-Set approach for computing solutions to incompressible two-phase flow,” J. Comput Phys, Vol. 114, pp. 146-159. Suzuki, T., Tobita, Y., Yamano, H., Kondo, S., Morita, K., Matsumoto, T., Akasaka, R. & Fukuda, K. (2003) “Development of multi-component vaporization/condensation model for a reactor safety analysis code SIMMER-III: extended verification using multi-bubble condensation experiment,” Nuclear Engineering and Design, Vol. 220, pp. 240–254. Swamee, P.K. and Ojha, C.S.P. (1991) "Drag coefficient and fall velocity of non- spherical particles," Journal of Hydraulic Engineering, Vol. 117, pp. 660-667. Swamy, N.V.C., Gowda, B.H.L & Lakshminath, V.R. (1979) "Auto-correlation measurements and integral time-scales in three-dimensional turbulent boundary layers," Applied Sc. Res. Vol. 35, pp. 265-316. Syamlal, M., Rogers, W. and O’Brien, T.J. (1993) “MFIX documentation, theory guide,” US Dept. of Energy Technical Note DOE/METC-94/1004. Tadaki, T. & Maeda, S. (1961) ) ”On the shape and velocity of single air bubbles rising in various liquids,” Kagaku Kogaku, Vol. 25, p. 254. Takagi, S. and Matsumoto, Y. (1999) "Force acting on a rising bubbles in a quiescent liquid," ASME Fluids Engineering Conference, FED 236, pp. 575-580. Takagi, S. and Matsumoto, Y. (1999) "Numerical investigations of the lift force acting on bubbles and particles," 3rd JSME/ASME Joint Fluids Engineering Conference, San Francisco, FEDS99-7848, July. Takagi, S., Prosperetti, A. and Matsumoto, Y. (1994) "Drag coefficient of a gas bubble in an axisymmetric shear flow," Physics of Fluids, Vol. 6 , No. 9, September, pp. 3186- 3188. Takata, S., Aoki, K. & Sone, Y. (1994) “Thermophoresis of a sphere with a uniform temperature: numerical analysis of the Boltzmann equation for hard-sphere molecules,” Progress in Astronautics and Aeronautics, Vol. 159, pp. 626_639. Takemura, F. (2004) "Migration velocities of spherical solid particles near a vertical wall for Reynolds number from 0.1 to 5," Physics of Fluids, Vol. 16, pp. 204-207. Takemura, F. & Magnaudet, J. (2003) “The transverse force on clean and contaminated bubbles rising near a vertical wall at moderate Reynolds number,” J. Fluid Mech. Vol. 495, pp. 234–253. Takemura F. & Magnaudet, J. (2004) “The history force on a rapidly shrinking bubble rising at finite Reynolds number,” Physics of Fluids, Vol. 16, pp. 3247–3255. Takemura, F., Takagi, S., Magnaudet, J. & Matsumoto, Y. (2002) "Drag and lift forces on a bubble rising near a vertical wall in a viscous liquid," J. Fluid Mechanics, Vol. 461, pp. 277-300. Talbot, L., Cheng, R.K., Schefer, R.W. & Willis, D.R. (1980) “Thermophoresis of particles in a heated boundary layer,” J. Fluid Mech. Vol. 101, pp. 737-758. Tanaka, T., Yonemura, S. & Tsuji, Y. (1990) “Experiments of fluid forces on a rotating sphere and spheroid,” Proc. 2nd KSME-JSME Fluids Engr. Conf., Vol. 1, p. 366. Tani, I. (1950) “Baseball’s curved balls,” Kagaku., Vol. 20, p. 405-409. 106 Taylor, G. I., (1949) “The shape and acceleration of a drop in a high-speed air stream,” for Advisory Council on Scientific Research and Technical Development, Ministry of Supply, AC 10647/Phys. C69. Taylor, T.D. & Acrivos, A. (1964) “On the deformation and drag of a falling viscous drop at low Reynolds number,” J. of Fluid Mechanics, Vol. 18, pp. 466-476. Ten Cate, A. & Sundaresan, S. (2006) “Analysis of unsteady forces in ordered arrays,” J. of Fluid Mechanics, Vol. 552, pp. 257–287. Tennekes, H. and Lumley, J.L. (1972) A First Course in Turbulence, MIT Press, Cambridge, Massachusetts. Theofanous, T.G. & Sullivan, J. (1982) “Turbulence in two-phase dispersed flow,” Journal Fluid Mech. Vol. 116, pp. 343-362. Theofanous, T.G., Li, G.J. and Dinh, T.N. (2004) "Aerobreakup in rarefied supersonic gas flows," ASME Journal of Fluids Engineering, Vol. 126, pp. 516-527. Thompson, P.A. (1972) Compressible Fluid Dynamics, McGraw Hill, New York. Thompson, T.L. & Clark, N.N. (1991) “A holistic approach to particle drag prediction,” Powder Technology, Vol. 67, pp. 57-66. Thorsen, G., Stordalen, R.M. & Terjesen, S.G. (1968) “On the terminal velocity of circulating and oscillating liquid drops,” Chem. Engr. Sci. Vol. 23, pp. 413. Tio, K.-K., Liñán, A., Lasheras, J.C. and Gañán-Calvo, A.M. (1993) "The dynamics of bubbles in periodic vortex flows," Applied Scientific Research, Vol. 51, pp. 285-290. Tomiyama, A. (1998) "Plenary lecture: struggle with computational bubble dynamics," International Conference on Multiphase Flow, Lyon, France, June. Tomiyama, A., Celata, G.P., Hosokawa, S. and Yoshida, S. (2002a) “Terminal velocity of single bubbles in surface tension force dominant regime,” International Journal of Multiphase Flow, Vol. 28, pp. 1497-1519. Tomiyama, A., Tamai, H., Zun, I. & Hosokawa, S. (2002b) “Transverse migration of single bubbles in simple shear layers,” Chemical Engineering Science, Vol. 57, pp. 1849-1858. Tomiyama, A. & Shimada, N. (2001) “A numerical method for bubbly flow simulation based on a multi-fluid model,” Journal of Pressure Vessel Technology, Vol. 123, pp. 510-516. Torobin, L.B. and Gauvin, W.H. (1960) Canadian Journal of Chemical Engineering, Vol. 38, pp. 142-153. Tran-Cong, S., Gay, M. & Michaelides, E.E. (2004) “Drag coefficients of irregularly shaped particles,” Powder Technology, Vol. 139, pp. 21-32. Tri, B.D., Oesterle, B. and Deneu, F. (1990) “Premiers resultants sur la portance d’une sphere en rotation aux nombres de Reynolds intermediaies,” C.R. Acad. Sci., Ser.II: Mec., Phys., Chim., Sci. Terre Universe 311, pp. 27. Tsao, H.K. and Koch, D.T (1997) "Observations of high Reynolds number bubbles interacting with a rigid wall," Physics of Fluids, Vol. 9, No. 1, January, pp. 44-56. Tsirkunov, Y.M. and Panfilov, S.V. (1998) "Modelling of particle-wall interaction in two-phase flows at moderate and high particle impact velocity," International Conference on Multiphase Flow, Lyon, France, June. Tsouris, C. and Tavlarides, L.L (1994) “Breakage and coalescence models for drops in turbulent dispersions,” AIChE Journal, Vol. 40, No. 395-406. 107 Tsuge, H. & Hibino, S.I. (1977) "The onset of oscillatory motion of single gas bubbles rising in various liquids," J. Chem. Engr. Japan. Vol. 10, pp. 66-68. Tsuji, Y. and Morikawa, Y. (1982) "LDV measurements of an air-solid two-phase flow in a horizontal pipe," J. Fluid Mech., Vol. 120, pp. 385-409. Tsuji, Y., Morikawa, Y. and Mizuno, O. (1985) "Experimental measurements of the Magnus Force on a rotating sphere at low Reynolds numbers," J. Fluids Engr. Vol. 107, No. 9, pp. 484. Tsuji, Y., Morikawa, Y. and Shiomi, H. (1984) “LDV measurements of an air-solid two- phase flow in a vertical pipe,” J. Fluid Mech. Vol. 139, pp. 417-434. Tunstall, E.B. & Houghton, G. (1968) “Retardation of falling spheres by hydrodynamic oscillations,” Chemical Eng Science, Vol. 23, No. 9, pp. 1067-1081. Underwood, B.Y. (1993) “Random-walk modelling of turbulent impaction to a smooth wall,” Intl. Journal of Multiphase Flow, Vol. 19, No. 3, pp. 485-500. Unverdi, S.O. and Tryggvason, G. (1992) “A front-tracking method for viscous incompressible, multi-fluid flows,” Journal of Computational Physics, Vol. 100, No. 1, pp. 25-37. Urbin, G. & Knight, D. (2001) “Large-eddy simulation of a supersonic boundary layer using an unstructured grid,” AIAA J. Vol. 39, July, pp. 1288-1295. Ushijima, T. and Perkins, R.J. (1999) "Evaluation of Co and To in turbulent pipe or channel flows," ASME Summer Fluids Engineering Meeting, FEDSM99-7759, San Francisco, July. Vallance, J.M. (1980) http://www.geo.mtu.edu/volcanoes/hazards/primer/pyro.html. Vandromme, D. (1997) "Introduction to the modeling of turbulence: turbulence modeling for compressible flows and acoustics," von Karman Institute for Fluid Dynamics, Lecture Series 1997-03, March. van der Kooij, F.M., Boek, E.S. & Philipse, A.P. (2001) "Rehology of dilute suspensions of hard platelike colloids," J. of Colloid and Interface Science, Vol. 234, pp. 344-349. van der Geld, C.W.M. (1997) "Measurement and prediction of solid sphere trajectories in accelerated gas flow," Int. J. Multiphase Flow, Vol. 23, No. 2, pp. 357-376. Van der Werff, J.C. & de Kruif, C.G. (1989) “Hard-sphere colloidal dispersions: the scaling of rheological properties with particle size, volume fraction, and shear rate,” Journal of Rheology, Vol. 33, pp. 421-454. van Dyke, M. (1982) An Album of Fluid Motion, Parabolic Press, Stanford, CA. van Wijngaarden, L. (1976.) “Hydrodynamic interaction between gas bubbles in liquid,” J. Fluid Mech. Vol. 77, pp. 27–44. van Wiijngaarden, L. (1998) "On pseudo turbulence," Theoretical and Computational Fluid Dynamics, Vol. 10, pp. 449-458. Vasseur, P. & Cox, R.G. (1977) “The lateral migration of spherical particles sedimenting in a stagnant bounded fluid,” Journal of Fluid Mechanics, Vol. 80, pp. 561. Venkateswaran, S., Lindau, J.W., Kunz, R.F. and Merkle, C.L. (2002) “Computation of multiphase mixture flows with compressibility effects,“ Journal of Computational Physics, Vol. 180, pp. 54-77. Verberg, R., Cohen, E.G.D. and de Schepper, I.M. (1997) “Viscosity of colloidal suspensions,” Phys. Review E, Vol. 55, No. 3, pp. 3143-3158. 108 Viollet, P.L. and Simonin, O. (1994) "Modeling dispersed two-phase flows: closure, validation and software development," part of "Mechanics USA 1994" appeared in Appl. Mech. Rev. Vol. 47, No. 6, June, pp. s80-s84, also ASME Reprint No. AMR146. Volkov, A.N. (2007) “The aerodynamic and heat properties of a spinning spherical particle in transitional flow,” 6th Internal Conference on Multiphase Flow, Leipzig, Germany, S2_Mon_C_6. de Vries, A.W.G., Biesheuvel, A. and van Wijngaarden, L. (2001) “Notes on the path and wake of a gas bubble rising in pure water,” International Journal of Multiphase Flow, Vol. 28, No. 11, pp. 1823-1835 Wadell, H. (1934) “The coefficient of resistance as a function of Reynolds number for solids of various shapes,” Journal of Franklin Institute, Vol. 41, pp. 310-331. Wakaba, L.V. & Balachandar, S. (2005) "History force on a sphere in a weak linear shear flow, " International Journal of Multiphase Flow, Vol. 31, pp. 996-1014, 2005. Waldman, L. (1961) Rarefied Gas Dynamics (ed. L. Talbot), Academic Press, pp. 3232. Wallis, G.B. (1969) One-Dimensional Two-Phase Flow, McGraw-Hill, New York. Wallis. G.B. (1974) “The terminal speed of single drops in an infinite medium,” Intl. Journal of Multiphase Flow, Vol. 1, pp. 491-511. Walter, J.F. and Blanch, H.W. (1986) “Bubble break-up in gas-liquid bioreactors: break- up in turbulent flows,” Chemical Engineering Journal Vol. 321, pp. B7-B17. Walton, O.R. (1993) “Numerical simulation of inelastic, frictional particle-particle interactions,” Particulate Two-Phase Flow, edited M.C. Roco, Butterworth- Heinemann, London, pp. 884-907. Wang, B.Y., Wu, Q.S., Wang, C., Igra, O. & Falcovitz, J. (2001) “Shock wave diffraction by a square cavity filled with dusty gas,” Shock Waves, Vol. 11, pp. 7-14. Wang, C.T. (1972) “Free molecular flow over a rotating sphere,” AIAA J. Vol. 10, pp. 713-714. Wang, D.M., Issa, R.I. & Gosman, A.D. (1994) “Numerical prediction of dispersed bubbly flow in a sudden enlargement,” Numerical Methods in Multiphase Flows, FED-Vol. 185, ASME. Wang, L.-P. and Maxey, M.R. (1993) “Settling velocity and concentration distribution of heavy particles in homogeneous, isotropic turbulence,” J. Fluid Mech. Vol. 256, pp. 27-68. Wang, L.-P. and Stock, D.E. (1993) “Dispersion of heavy particles by turbulent motion,” Journal of the Atmospheric Sciences, Vol. 50, No. 13, pp. 1897-1913. Wang, Q., Squires, K.D., Chen, M. and McLaughlin, J.B. (1997) "On the role of the lift force in turbulence simulations of particle deposition," Int. J. Multiphase Flow, Vol. 23, pp. 749-763. Wang, S.K., Lee, S.J., Jones, O.C. & Lahey, R.T. (1987) “3-D turbulence structure and phase distribution measurements in bubbly two-phase flows,” Interr. J. Multiphase Flow, Vol. 13, No. 3. Warnica, W.D., Renksizbulut, M. & Strong, A.B. (1995) “Drag coefficients of spherical liquid droplet. Part II turbulent gaseous fields,” Experiments in Fluids, Vol. 18, pp. 258. Watts, R.G. & Ferrer, R. (1987) “The lateral force on a spinning sphere: aerodynamics of a curveball,” Am. J. Phys. Vol. 55, pp. 40–44. 109 Wegener, P.P. & Ashkenas, H. (1961) "Wind tunnel measurments of sphere drag at supersonic speeds and low Reynolds numbers,” J. Fluid Mechanics, Vol. 10, pp. 550- 560. Wegener, P.P., Sundell, R.E. and Parlange, J.-Y. (1971) "Spherical-cap bubbles rising in liquids," Z. Flugwissenschaften, Vol. 19, pp.347-352. Welleck, R.M., Agrawal, A.K. & Skelland, A.H.P. (1966) "Shape of liquid drops moving in liquid media," AIChE J. Vol. 12, pp. 854-862. Wells, M.R. and Stock, D.E. (1983) "The effects of crossing trajectories on the dispersion of particles in a turbulent flow," Journal of Fluid Mechanics, Vol. 136, pp. 31-62. Wen, C.Y. and Yu, Y.H. (1966) “Mechanics of fluidization,” Chemical Engineering Progr. Symp. Series, Vol. 62, p. 100. Wen, F., Kamalu, N., Chung, J.N., Crowe, C.T. & Trout, T.R. (1992) "Particle dispersion by vortex structures in plane mixing layers," Journal of Fluids Engineering, Vol. 114, pp. 657-666. Wesseleing, P. & Oosterlee, C.W. (2001) "Geometric multigrid with applications to computational fluid dynamics," J. of Comp. & Appl. Math, Vol. 128, pp. 311-334. Westerweel, J., Draad, A.A., van der Hoeven, J.G.T. & van Oord, J. (1996) "Measurement of fully-developed turbulent pipe flow with digital particle image velocimetry," Experiments in Fluids, Vol. 20, pp. 165-177. White, F.M. (1991) Viscous Fluid Flow, McGraw-Hill, New York. Williams, F.A. (1965) Combustion Theory, Addison-Wessley, Reading, Mass. Williams, J.J.E. & Crane, R.I. (1983) “Particle collision rate in turbulent flow,” Intl. Journal of Multiphase Flow, Vol. 9, pp. 421–435. Willis, K.D. and Orme, M.E. (2000) "Experiments on the dynamics of droplet collisions in a vacuum," Experiments in Fluids, Vol. 29, pp 347-358. Willmarth, W.W., Hawk, N.E. and Harvey, R. L. (1964) “Steady and unsteady motions and wakes of freely falling disks,” Phys. Fluids, Vol. 7, pp. 197–208. Winnikow, S. & Chao, B.T. (1966) “Droplet motion in purified systems,” Physics of Fluids, Vol. 9, pp. 50-61. Wolfrum, B., Mettin, R., Kurz, T. and Lauterborn, W. (2003) “Cavitation induced cell detachment and membrane permeabilization,” 2003 IEEE ULTRASONICS SYMPOSIUM-837. Wostenberg, R. (2005) Personal Communication, Graphic Engineering (Montague, CA). Wu, J.-S. and Faeth, G.M. (1994) "Sphere wakes at moderate Reynolds numbers in a turbulent environment," AIAA Journal, Vol. 32, No. 2, pp. 535-54. Wu, M. & Gharib, M. (2002) “Experimental studies on the shape and path of small air bubbles rising in clean water,” Physics of Fluids, Vol. 14, No. 7, pp. 139-150. Wu, S.Y. and Baeyens (1998) “Segregation by size difference in gas fluidized beds,” Powder Technology, pp. 139-150. Wunsch, D., Fede, P. and Simonin, O. (2008) “Development and validation of a binary collision detection algorithm for a polydispersed particle mixture,” FEDSM2008- 55159, ASME Fluids Engineering Conference, August 10-14, 2008, Jacksonville, Florida, USA. Wygnanski, I. and Fiedler, H. (1969) “Some measurements in the self-preserving jet,” J. Fluid Mech. Vol. 38, pp. 577-612. 110 Wygnanski, I. and Fiedler, H.E. (1970), "The two-dimensional mixing region," Journal of Fluid Mechanics, Vol. 41, No. 2, pp. 327-361. Xie, H.-Y. & Zhang, D.-W. (2001) "Stokes shape factor and its application in the measurement of sphericity of non-spherical particles," Powder Technology, Vol. 114, pp. 102-105. Xu, J., Maxey, M. and Karniadakis, G.E. (2002) "Numerical simulations of turbulent drag reduction using micro-bubbles," J. Fluid Mech. Vol. 468, pp. 271-281. Yamamato, K. and Ishihara, Y. (1988) “Thermophoresis of a spherical particle in a rarefied gas of a transition regime,” Phys. Fluids, Vol. 31, pp. 3618-3624. Yan, F., Lightstone, M.F. & Wood, P.E. (2007) “Numerical study on turbulence modulation in gas-particle flows,” Heat Mass Transfer, Vol. 43, pp. 243-253. Yan, H., Knight, D. & Zheltovodov, A.A. (2002) “Large eddy simulation of supersonic flat plate boundary layer – part I,” AIAA Aerospace Sciences Meeting, Reno, January, AIAA2002-0132. Yang, C.Y. and Lei, U. (1998) “Role of the turbulent scales in the settling velocity of heavy particles in homogeneous isotropic turbulence,” Journal of Fluid Mechanics, Vol. 371, pp. 179-205. Yang, T.S. and Shy, S.S. (2003) “The settling velocity of heavy particles in an aqueous near-isotropic turbulence,” Physics Fluids, Vol. 15, No. 4, pp. 868-880. Yang, T.S. and Shy, S.S. (2005) “Two-way interaction between solid particles and homogenous air turbulence: particle settling rate and turbulence modification measurements,” J. Fluid Mechanics, Vol. 526, pp. 171-216. Yang, Y., Chung, J.N., Troutt, T.R. and Crowe, C.T. (1990) “The influence of particles on the spatial stability of two-phase mixing layers,” Physics Fluids A, Vol. 2, No. 10 pp. 1839-1845. Yarin, L. P. and Hetsroni, G. (1994) “Turbulence intensity in dilute two-phase flows,” Int. J. Multiphase Flow, Vol. 20, pp. 1–44. Yih, C. (1969) Fluid Mechanics, New York, McGraw-Hill. Yoon, J. (2001) “Jet engine types,” http://www.aerospaceweb.org/question/propulsion/q0033.shtml. Young, J. and Leeming, A. (1997) "A theory of particle deposition in a turbulent pipe flow," Journal of Fluid Mechanics, Vol. 340, pp. 129-159. Young, J.B. and Hanratty, T.J. (1991)"Trapping of solid particles at a wall in a turbulent flow," AIChE Journal, October, Vol. 37, No. 10, pp. 1529-1536. Yoshizawa A. & Horiuti K. (1985) “A Statistically-Derived Subgrid-Scale Kinetic Energy Model for the Large-Eddy Simulation of Turbulent Flows,” Journal of the Physical Society of Japan , Vol. 54, pp. 2834-2839. Yuan, H. and Prosperetti, A. (1994) "On the in-line motion of two spherical bubbles in a viscous fluid," Journal of Fluid Mechanics, Vol. 278, pp. 325-349. Yuen, M.C. & Chen, L.W. (1976) „On the drag of evaporating liquid droplets,” Cobust. Sci. Technol. Vol. 14, pp. 147. Zakharov, L.V., Ovchinnikov, A.A. & Nikolayev, N.A. (1993) “Modeling of the effect of turbulent two-phase flow friction decrease under the influence of dispersed phase elements,” International J. Heat and Mass Transfer, Vol. 36, pp. 1981-1991. Zarin, N.A. (1970) “Measurement of non-continuum and turbulence effects on subsonic sphere drag,” NASA CR-1585 June. 111 Zeng, L., Balachandar, S., Fischer, P. and Najjar, F. (2008) “Interactions of a stationary finite-sized particle with wall turbulence,” Journal of Fluid Mechanics, Vol. 594, pp. 271-305. Zenit, R., Joseph, G. and Hunt, M. (1999) "The coefficient of restitution for liquid immersed collisions," ASME Summer Fluids Engineering Meeting, FEDSM99-7793, San Francisco, July. Zenit, R., Koch, D.L. & Sangani, A.S. (2001) “Measurements of the average properties of a suspension of bubbles rising in a vertical channel,” J. Fluid Mech. Vol. 429, pp. 307-342. Zhang, D.Z. and Prosperetti, A. (1994) "Averaged equations for inviscid dispersed two- phase flow," Journal of Fluid Mechanics, Vol. 267, pp. 185-219. Zhang, Y. and Reese, J.M. (2003) “The drag force in two fluid models of gas-solid flows,” Chemical Engineering Science, Vol. 58, pp. 1641-1644. Zhao, R. and Rognebakke, O. F. (2004) “Sloshing and impact loads in membrane tanks,” MARINTEK Review No 1 - April – 2004, pp.2-3. Zheng, F. (2002) “Thermophoresis of spherical and non-spherical particles: a review of theories and experiments,” Advances in Colloid and Interface Science, Vol. 97, pp. 255-278. Zhou, Q. & Leschziner, M.A. (1991) “A time-correlated stochastic model for particle dispersion in anisotropic flows,” 8th Turbulence Shear Flow Symposium, Munich, Sept. Zisselmar, R. and Molerus, O. (1979) “Investigation of solid-liquid pipe flow with regard to turbulence modification,” Chem. Eng. Journal, Vol. 18, pp. 233-239. Zuber, N. (1964) “On the dispersed two-phase flow in the laminar flow regime,” Chemical Engineering Science Vol. 19, pp. 897-917. Zun, I., Kljenek, I. and Serizaw, A. (1992) "Bubble coalescence and transition from wall void peaking to core void peaking in turbulent bubbly flow," In Dynamics of Two- Phase Flows, ed. Jones, O.C. and Michiyoshi, I., pp. 233-239, Boca Raton, FL. 112

DOCUMENT INFO

Shared By:

Categories:

Tags:

Stats:

views: | 13 |

posted: | 5/17/2012 |

language: | English |

pages: | 112 |

How are you planning on using Docstoc?
BUSINESS
PERSONAL

By registering with docstoc.com you agree to our
privacy policy and
terms of service, and to receive content and offer notifications.

Docstoc is the premier online destination to start and grow small businesses. It hosts the best quality and widest selection of professional documents (over 20 million) and resources including expert videos, articles and productivity tools to make every small business better.

Search or Browse for any specific document or resource you need for your business. Or explore our curated resources for Starting a Business, Growing a Business or for Professional Development.

Feel free to Contact Us with any questions you might have.