Document Sample

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 and U are equal, but we will generally use the former to be consistent with most of the equations in the main text. 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 though 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 expression, u is the speed 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 q (qu) q A.4 t This expression is often referred to as the conservation equation if q =0. 1 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 expression is particularly convenient for Eulerian formulations of the continuous- phase flow. The second expression 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 uf 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, which can be expressed in terms of shearing and compressive velocity gradients in terms of both shear and bulk viscosities u u j u K i j f i bulk,f k A.8 x j x i x k 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 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 2 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 k x j x i 3 In this expression, ij is defined as 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 For “inviscid” flow, the viscous stresses are negligible so 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 de 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 expression 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 internal energy as u = K ij pij i + k f T De f A.14 Dt x j In contrast to Eq. A.13, this transport equation is based on the substantial derivative of the internal energy such that the pressure and body force work terms do not appear explicitly (White, 1991). 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 As will be discussed in §B.3.3, conservative forms of the transport equations are preferred for compressible flows. 3 Fluid Properties At STP, several of the fluid properties used in the above PDE (viscosity, density, conductivity, etc.) are given in Table A.1 for of air, 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 expression, the subscript “o” indicates STP and absolute temperatures are used. This is form is at least qualitatively reasonable for other gasses up to about 400 K, but more detailed temperature dependencies for several gasses and liquids are given in App. A of White (1991). 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 such that h cp A.18a T p e c A.18b T In the above expression, the proportionalities are referred to as is the specific heat at constant pressure (cp) and specific heat at constant volume (c). 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 7/5, which is a reasonable approximation for air at STP. At higher temperatures, the gas is generally no longer calorically perfect and is reduced. 4 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 qualitative differences between the pressure and density relationship for a gas and a liquid can be seen in Fig. A.2. In general, liquids are much less compressible so that the gradients (and thus speed of sound) in a liquid are typically much greater than that for a gas. For example, the speed of sound in water is nearly five times that for air (Table A.1). The Mach number is defined as the local speed of the continuous-fluid normalized by the speed of sound u M A.21 a This parameter indicates the relative speeds of convection as compared to acoustic signals. It is useful in characterizing the mathematical character of the continuous-phase PDEs as will be discussed in the next section. 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 force 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 as observed by Dalton with a proportionality constant equal to the universal gas constant (R univ) divided by the molecular weight of the gas (MWg). This proportionality constant is also referred to as the specific gas constant, R g, and the resulting equation of state is known as the perfect gas law: R pg R g univ A.22 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.23 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.24 Further assuming constant entropy, i.e., “isentropic flow”, the changes in pressure can then be related to the changes in density and temperature to the local Mach number such that 1 1 1 T p 1 2 1 M A.25 Ttot tot gas p tot gas 2 5 In these relationships, “tot” refers to the total or stagnation condition where the flow has undergoes an isentropic deceleration to zero velocity. Thus changes in pressure and temperature become more pronounced once the Mach number variations are significant. 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 STP: K p+Bl po l l A.26 (1+ Bl )po lo As an example, the compressibility of liquid water at standard temperature is reasonably predicted with Bl=3000 and Kl=7.15, and Thompson (1972) gives values for these constants for other fluids. The speed of sound for liquid (al) from the Tait equation based on isentropic conditions becomes K l (p Bl p o ) al2 A.27 l Note that this expression reverts to the calorically perfect gas conditions for B=0 and K=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 density. As such, the Mach number is a good indicator of the level of compressibility associated with changes in velocity. In the following section, the magnitude of the Mach number will be correlated with the character of the PDEs for inviscid conditions. 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. When neglecting viscosity for the continuous-flow description, the wall boundaries should be treated with 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). For convenience, gravitational effects can also be neglected (which is reasonable if the flow is not stratified and is considered as a single- phase). Based on these two assumptions, the convection and pressure effects will then dominate. Let us first consider the low Mach number regime where the density will not be significantly altered due to the 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 6 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 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 becomes u 0 A.29 This is often referred to as the continuity equation in Cartesian coordinates. 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. In Cartesian coordinates, the continuity equation is written as u x u y u z 0 A.30 x y z Defining r as the radial coordinal, as the polar angle, and as the azimuthal angle about the axis =0, the continuity equation in spherical coordinates becomes 1 r ur 2 1 u sin u 0 A.31 r 2 r r sin The expression in cylindrical coordinates and general curvilinear orthogonal coordinates is given by Batchelor (1972). 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.32a y -u y A.32b 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 sine 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.33 7 The RHS expression is the vorticity in Cartesian coordinates and its magnitude in can be defined by the separate components or in tensor form: u z u y u x u z u y u x 2 2 2 2 2 2 y z z x x y x y z A.34a 1 u i u j u j u i 2 x j x i x i x j A.34b Two simple flow conditions with vorticity are the linear shear flow and the simple vortex flow. For linear shear flow, the flow lines are parallel and the velocity field has a uniform gradient so that the shear rate is equal to the vorticity magnitude. For example, a 2-D flow with u y u z 0 and a linear gradient in x-velocity is given by: u x shear const. for linear shear flow A.35 y In contrast, the simple vortex flow is based on solid-body rotation of the fluid and a 2-D flow with u r u z 0 with a linearly increasing tangential velocity will have a vorticity given by: u vortex 2 const. for simple vortex flow A.36 r If there is no vorticity flow altogether (=0), such a flow is called “irrotational”. The third descriptor is the velocity potential defined as u A.37 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 expression can be combined with the flow continuity (Eq. A.30) to yield: 2 0 A.38 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.32 and A.34a can be used to show the resulting stream function will also satisfy Laplace’s equation: 2 0 A.39 This property allows linear super-position of solutions, i.e. if 1 and 2 satisfy Eq. A.38, 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.40 This relationship may be combined with a vorticity transport equation to describe two- dimensional flow incompressible flow, independent of whether viscosity is included. As a 8 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 u 2 f g - p A.41 For steady flow, this equation can be integrated to yield an incompressible relationship for the pressure throughout the domain p 12 f u 2 f gz const. A.42 If hydrostatic effects are small over the domain scale D, i.e. gD u 2 , it is often reasonable to D neglect the body force, so that the pressure relationship becomes: p 12 f u 2 pstag A.43 This is called the Bernoulli relationship and indicates that pstag is the highest pressure possible in the system and occurs at the stagnation points. The term 12 f u 2 is know as the “dynamic pressure” for incompressible flow and can be used to scale the stresses of many phenomena which are governed by convective effects. If hydrostatic effects are important, a reference pressure may be defined relative to the hydrostatic pressure as p p f gz A.44 This expression assumes constant density and assumes that gravity acts in the -z direction. This pressure is often convenient in that substitution into the momentum equations removes the gravitational term and similarly need not be included in any pressure boundary conditions. Once this reference pressure distribution is known, then Eq. A.44 can be used to determine the local pressure throughout. 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. These equations can be obtained by rearranging Eqs. A.5a, A.6a and A.12 and neglecting hydrostatic effects and the work to overcome potential energy (reasonable if gD u 2 ) yielding: D 9 f f u j 0 A.45a t x j f u i f u i u j pij 0 A.45b t x j f e tot f u je tot pu j 0 A.45c t x j 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.38. 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.46 x y z where the free-stream Mach number is defined as M∞=ux∞/a∞. This expression can then be evaluated for four different Mach number regimes which characterize the relative importance of the compressibility: M«1, M<1, M~1, and M>1. In addition, we can identify the mathematic character of the governing PDE (Eq. A.46) by considering the slope of the characteristic directions (Anderson et al. 1997): dy dz 1 A.47 dx dx M 1 2 These regimes and their mathematical character in space are discussed below and shown in Fig. A.3. 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.46 reverts to the incompressible form given in Eq. A.38) 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.46 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 10 downstream and Eq. A.46 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. In the supersonic regime, the kinetic energy associated with the fluid velocity is on the order of the internal energy from random molecular motion. 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.48 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.49 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. 11 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.50 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.51a t x i x i k k u i f ,k k for constant density A.51b t x i x i 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. Also, if one considers small particles in a fluid these are also not miscible since theses have an interface. 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.7). 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.52 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 their mass diffusion 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). The ratio of momentum diffusion to thermal diffusion is given by the (laminar) Prandtl number f c p,f Pr A.53 kf where cp is the heat capacity at constant pressure. For gases, the Prandtl number is of order unity (0.69 for a diatomic gas) 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 12 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 adiabatic and non-reacting, 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 for Cartesian coordinates is given by: u u j u i u j K i j f i f f 2u u x j x i x j x i x j A.54 Again assuming incompressibility (Eq. A.29), the momentum equation of Eq. A.6b becomes: Du f p f g f 2u A.55 Dt The Cartesian tensor forms of the continuity and the momentum equations are thus: u i 0 A.56a x i u i u i p 2ui f f u j f g i - f A.56b t x j x i x 2j In this expression, gi is the component of gravitational acceleration in the x i direction. In the limit of quiescent flow (ui=0), this yields f g i = p A.57 This is the hydrostatic equation It is often convenient to consider incompressible flow in spherical coordinates. In this case, the continuity equation is given by Eq. A.31. The spherical components of the viscous stress tensor in spherical coordinates are given by Batchelor (1967) as: u K rr 2f r r A.58a u 1 u r K r f r A.58b r r r u 1 u A.58c K 2f r r r 1 u u r u cot K 2 f r sin r r A.58d 1 u r r u K r 2 f A.58e 2r sin 2 r r sin u 1 u A.58f K 2 f 2r sin 2r sin 13 The spherical momentum equations 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 From this, the incompressible momentum conservation equations are given by u r u r u u r u u r u u 2 2 ur gr t r r r sin r r A.60a 1 p 2 2u r 2 u sin u u r 2 2 f 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.60b 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.60c 1 p 2 2 u r u u f u 2 cot f r sin r sin 2sin The expressions in cylindrical coordinates and general curvilinear orthogonal coordinates are given by Batchelor (1972). The momentum equations can also be non-dimensionalized by defining D and uD as the reference macroscopic continuous-phase length and speed, respectively. For a pipe flow, these values could correspond to the pipe diameter and mean speed. The terms in Eq. A.55 can be normalized by f, f, D, and uD, to yield dimensionless equations u * i 0 A.61a x * i Du* 1 g 1 * p* 2 u* A.61b Dt FrD g Re D In these expression ()* super-script indicates non-dimensional values: u*=u/uD, x*=x/D, t*=tuD/D, p*=p/f(uD)2. This non-dimensional form of the pressure is based on Eq. A.43 and thus is reasonable for regions of the flow where viscous effects. The velocity and length scales used to define these variables are generally based on easily defined aspects of the flow and its geometry. The expression also include the non-dimensional macroscopic flow Froude and Reynolds numbers u2 FrD D A.62a gD Du Re D f D A.62b f 14 The Froude number definition used above is consistent with the definition used by 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.42, 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 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 local skin friction, and can control flow separation characteristics which then affect the entire flow field. In 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,t) 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,t 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 (Figs. A.3a and A.3b). For transitional flows, the structures will exhibit detailed features and become increasingly non-linear, though still remain approximately two-dimensional (Fig. A.4b). 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). 15 For ReD>ReD,t (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.4c, where it can be seen that the large scale structures of Fig. A.4b 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 it is 2,000-2,400 for a round pipe flow and even less for a free shear layer. In either case, fully-developed turbulence generally occurs soon after transition, e.g. ReD,t≈4,000 for a pipe flow as shown in Fig. A.5. The influence of the various Reynolds number regimes is illustrated in Fig. A.6 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). 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.6a). For ReD~1, the convective and viscous diffusion rates will be similar order so that upstream diffusion is significant but limited (Fig. A.6b). 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.6c). 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.6d). 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.6e). In this regime, the lateral mixing of both the mass and momentum is high in a mean sense as compared to that for high Reynolds number laminar flow (Fig. A.6c). These sketches are qualitatively similar to the celebrated experiments conducted by Reynolds in 1883 to study flow stability, i.e., where visualization of a dye in water flowing through a pipe was used to implicitly note the effects on momentum diffusion. The effect between laminar and turbulent diffusion is even more profound in this case because the Schmidt number is very large so that lateral diffusion of the dye mass was much less for laminar flow conditions but similar for turbulent flow conditions. Laminar Flow Momentum Equations The momentum conservation given by Eq. A.56b 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 16 hydrodynamic variations (Eq. A.43). This second assumption is also reasonable if pressure distribution is dominated by viscous forces, i.e. small ReD and small ReD /FrD . Laminar flow will generally encompass a large range of flow Reynolds numbers regimes including ReD0, ReD«1, ReD~1, and 1«ReD<ReD,crit. Each of these regimes may have a unique set of momentum equations. 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.56b yields (neglecting hydrostatic effects): u i p 2ui f u j f for ReD~1 A.63 x j x i x 2 j If hydrostatic forces are important, the reference pressure (Eq. A.44) can be used in this equation so that the gravitational terms are implicitly included. Based on Eq. A.61b, as the Reynolds number tends to zero, 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 PDE becomes: p 2ui f for ReD 0 A.64 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 (a further indication that it should be retained) except 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, the convective terms are small but finite (about 0.01<ReD<0.1) so that the convection velocity may be based approximated with the free-stream velocity while the spatial gradients are exact and based on the local fluid velocity. The resulting PDE becomes: u i p 2ui f u j f for ReD 1 A.65 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 (such as a particle) will directly influence the velocity field far away from the body in all directions since the (downstream) momentum diffusion is slow or negligible. 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 (convective) direction are often much smaller 17 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 x u x p 2u x f u x f u y 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 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 influence of viscosity decreases (i.e., as Reynolds number increases). 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.56) must be used. Transitional and Turbulent Momentum Equations While the governing equations are known for transitional and turbulent flow (e.g. Eq. A.56) 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. They are also three- dimensional, unsteady, inherently unstable, and effectively stochastic (since the detailed 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 18 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: T 1 q(x) q(x, t)dt as T A.67 T0 In this expression, 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 vector (as well as density, temperature, pressure, etc.) can be decomposed as: u (x, t) u (x) + u (x, t) A.69 The differences between the mean and instantaneous velocity fields are shown in Fig. A.7 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.8a. 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 T u (x) uu (x) = u (x, t) dt 2 i,rms i i i as T A.70 T0 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.8b, where the variations are on the order of the rms values. The continuous-phase turbulent kinetic energy (k) is defined as one-half of the combined strength of the turbulence in all three directions k(x) 1 uu (x) 1 u1u1 (x) u u (x) u u (x) 2 i i 2 2 2 3 3 A.71 and thus has units of energy per unit mass. One may average the velocity fluctuations on the RHS by assuming isotropy (the rms components have the same magnitude in all three directions) and then relating the average strength to the LHS: u u rms 2 3 k A.72 19 This can be used to non-dimensionalize the velocity fluctuations seen in Fig. A.8b, 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, we can define an expression for turbulent dissipation 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. Note that Schlichting (1979) and 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 D D though there are also effects of ReD (Hinze, 1975). 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.7 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. The RHS approximation in Eq. A.73 can be used if the turbulence in locally isotropic, a concept discussed in the next section. 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). The term “homogeneous” indicates that the strength of the turbulence is independent of position (no gradients in any direction). The term “isotropic” indicates that the turbulence strength is the same for all three components and also indicates, for homogenous turbulence, that the cross-correlations components are negligible: u u1u1 u u u u 2 k ms 2 2 3 3 3 A.74a uu i j i j 0 A.74b As will be discussed later (in A.5.5), the zero cross-correlation assumption (Eq. A.74b) indicates that there is no turbulent transport of velocity in any direction so that one may also expect k 0 , i.e. the homogenous turbulence assumption. The term “stationary” indicates that the mean statistics do not vary in time. In summary, the rms of the velocity fluctuations in 20 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.9 indicate that far from the wall (y/≈1), the flow is nearly homogenous and isotropic (little variations in different locations and between different components). 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 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 ( u1,rms ≠ u ≠ u ). For example, the mean transverse profiles of the 2,rms 3,rms three fluctuating velocity components near the wall (y/≈0.05) in Fig. A.9 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.75a u (x)uj,rms (x x ) i,rms u (x, t)uj (x x , t ) ij (x, x , τ) im jn i A.75b u (x)u (x x ) m,rms n,rms Note that Eq. A.75a is written with no summation of the indices (NSI), while the more formal version given by Eq. A.75b uses Kronecker deltas defined by Eq. A.10. Bother versions are shown as the text contains other expressions with only NSI for sake of conciseness. This correlation function considers the mean relationship of velocity perturbations with temporal 21 and spatial shifts. With no temporal and spatial shifts, the correlation function becomes unity, i.e. ij (0, 0) 1 . If the spatial shift is taken based on the fluid path (Fig. A.10), the result is the Lagrangian correlation tensor defined as t u (x, t)u j (x uDt, t ) L,,ij τ i t A.76 u uj,rms i,rms For homogeneous turbulence, this decorrelation is independent of location, x, and for isotropic turbulence the diagonal component is the same for all the velocity components: L τ L,,11 τ L,,22 τ L,,33 τ A.77 Note that the off-diagonal terms are zero for HIST flows based on Eq. A.74b so only a scalar is needed to represent the correlation. A typical temporal variation of such a Lagrangian correlation function is shown in Fig. A.11, where the time-separation is normalized by the Lagrangian integral time-scale (to be defined below). The auto-correlation result for no time separation obviously gives L (0) 1 while the correlation for long-times becomes negligible, i.e. L () 0 . For intermediate times, there is generally decay in the correlation function as the turbulent features change due to break-up, merging, distortion, dissipation, etc. The time- scale for this reduction is therefore a useful way to define the average temporal variations in the turbulence. Two other reference frames for the correlation function include Eulerian and “moving- Eulerian”. These reference frames correspond to different paths and associated spatial shifts as shown in Fig. A.10. The Eulerian frame is fixed in space (x=0) and thus the correlation function, E () is consistent with a fixed-position time-dependent measurement. The moving- Eulerian path (also called mean-convection path) is similar to the Lagrangian path except that x is based on the mean velocity u (instead of the instantaneous velocity). This resulting correlation function, mE , is a more convenient measurement since the location of the time- shifted velocity fluctuations is fixed and known in advance. These correlation functions can be used to define the respective “integral time-scales” of the turbulence for the Lagrangian (L), moving Eulerian (mE) and Eulerian (E) reference frames: 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 so that a single- scalar value is reasonable for both: 22 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.12a. If one estimates the correlation function to have an exponential decay (as shown in Fig. A.11), the relationships with the integral time-scales become simply L (τ) exp( / L ) A.80a E (τ) exp( / E ) A.80b mE (τ) exp( / mE ) A.80c 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. An interesting aspect of fluid-tracer correlations is that they can sometimes have a negative-loop for large separation times, which is not predicted by the exponential decay form of Eq. A.80. To incorporate this behavior, other correlation functions have been employed. For example, the exponential-cosine function (Fig. A.11) tends to give better qualitative description in the cases when a negative loop occurs. However, it does not represent a clear improvement in terms of experimental agreement. Furthermore, the magnitude of the integral time-scale is generally more important than the exact shape of the correlation function. As such, the simpler exponential function is more commonly used. 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 ( x ) but with no temporal shift (=0). This length-scale is generally defined by a streamwise spatial shift with zero temporal shift: 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 The integral time-scale based on streamwise component 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.13a) 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.4c. 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). 23 The correlation function and integral length-scale based on lateral fluctuations (normal to the mean velocity) are still defined with a streamwise spatial shift: u (x, t)u (x x , t) (x ) A.84a u ,rms u ,rms (x , 0)dx A.84b 0 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 mean direction as 1 and the two lateral directions as 2 and 3. Next consider, the fluctuations associated with u over the surface given by spatial shifts in the 2 other two directions (x 1 and x 3 ) as shown in Fig. A.13b. Since mass must to be conserved in the mean, any positive u 2 on the plane should be counter-acted somewhere else with a negative u 2 . Therefore, integration over this entire plane will give a zero net mass flux: u (x, t) dx dx 2 1 3 0 A.85 If this integral is weighted by u 2 at a given point on this plane and then averaged, the correlation integral over the plane will also be zero, u (x, t)u (x x , t) dx dx 2 2 1 3 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.13b). 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.12b. 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. 24 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 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.12b). 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.9), integral time- scale (Fig. A.12a) and length-scale (Fig. A.12b) 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.12a 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 25 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 integral scales to Kolmogorov scales is large. Based on Eqs. A.91 and A.92a, the ratio of the length-scales is proportional to the 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 (the latter is illustrated in Fig. A.14). At high Reynolds numbers, there is a wide range of scales associated with turbulence as shown in Fig. A.4c and A.5e. This spectrum is often analyzed in terms of the wave-number (n, which is inversely proportional to length-scale) and turbulent energy per wave number ( k , which is proportional to the total kinetic energy, and sometimes called the “specific turbulent energy”): n 2/l A.95a k k n dn A.95b 0 The conceptual relationship between these two parameters is shown in Fig. A.15. 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. These micro-scales correspond to the highest wave-numbers of the flow, 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, 26 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 A.96a kl ~ l 2/3 A.96b Thus, the turbulent energy content decreases as the wave-number becomes higher and the 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.15). 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 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 is smaller but the two length-scales are of similar size (Poorte & Biesheuvel, 2002): l int er l Taylor 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. 27 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 u y t u , y , t A.101 u fr fr fr 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 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., + 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. For the latter, the special case of incompressible stratified-flow is considered, where the fluid is considered incompressible but has variable density. Examples of such flows include low-speed flows (so that the Mach numbers is very small) where fluids of different densities mix. For example, a gas jet exhausting into another gas of lower density due to being colder or of lower molecular weight. A stratified flow example for liquids is that of oceanic salinity currents, where variations in salt concentration yield different densities. 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). 28 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 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. In some cases, it is important to consider the turbulent transport of a tracer species in the continuous-phase. Application of time-averaging to the species transport of Eq. A.51b yields ui f u A.105 x i x i x i i 29 This equation is similar to the steady-state laminar flow species equation where the LHS is the mean convection. The first term on the RHS represents the laminar diffusion associated with the mean gradients of the species. However, the second term with u represents turbulent i diffusion of the mean species fraction. This term is primarily responsible for the mean spread of a fluid tracer (such as a dye) in a turbulent flow. Turbulent mass diffusion of a species is illustrated in Fig. A.16 where turbulent eddies cause a dye to spatially mix by transverse velocity fluctuations in the vicinity of the dye interface. Note that an eddy structure in this region will tend to cause dye from the high concentrations to move into regions of low concentration. If one refers to the vertical spatial average of the concentration as , velocity fluctuations to the left ( u 0 ) will tend to pull black dye to the left, i.e. u 0 . One can x x also see that the diffusion of black dye will be proportional to the eddy strength and the concentration gradient, and in the direction opposing this gradient. Based on the above discussion, one would expect that a model of turbulent diffusion of a species should be proportional to the species gradients. In particular, a common approach is to employ Fick’s law for molecular diffusion (Eq. A.51b) in a similar manner for turbulence, i.e. u t t A.106 x i Sc t x i i In this expression, we have introduced a turbulent diffusivity (t) which is related to a turbulent viscosity (t) and a turbulent Schmidt number (Sct,m), similar to the laminar relation of Eq. A.52. The turbulent Schmidt number is generally taken as 0.7 based on experiments which compared turbulent mass diffusion to turbulent momentum diffusion (Faeth, 1987). Consistent with the illustration of Fig. A.16, the turbulent viscosity is expected to be related to the strength in turbulence and will be discussed in the following sub-section. The heuristic model of Eq. A.106 is a closure model in the sense that a correlation based on fluctuating values has been modeled with a gradient of the mean concentration field. The resulting turbulent mass transport equation then becomes ui f t A.107 x i x i x i 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.51b). As an example of the difference between diffusion in laminar flows and turbulent flows, consider the diffusion of miscible cream into coffee in a cup. Assume that both liquids have the same density and the initial conditions is a horizontal interface, e.g. with the cream on the top and the coffee below (Fig. A.17). If the flow in the cup was then stirred at a low speed such that the flow remains laminar, then the interface may be largely unbroken and it would take a long time for molecular diffusion (associated with Θf) to cause the two liquids to be homogeneously mixed (LHS of Fig. A.17). However, introduction of higher stirring speeds such that turbulent flow is incurred (like used for stirring cream into coffee), will create a cascade of vortices through which the interface between the two liquids is rapidly stretched, convoluted and distributed throughout the domain (RHS of Fig. A.17). This mixing associated 30 with Θt causes the interfacial area to greatly increase so that the molecular mixing needed for a uniform mixture only has to occur over very short distances. As a result, the homogeneous equilibrium is achieved in time-scales which are orders of magnitude smaller. Therefore, overall mixing in turbulent flow can be considered to be a two stage process: 1) interface eddy-mixing whereby a wide range of vortex structures causes pockets of high and low concentration to be transported stretched, convoluted and distorted (due to u ), so as to yield a large amount of interfacial area with high concentration gradients, i 2) laminar diffusivity whereby the high concentration gradients regions lead to rapid but local molecular mixing (associated with lam) at small scales. In laminar flows, only the second phenomenon takes place so that the overall interface area is much smaller and the concentration gradients are weaker. As a result, overall mixing will be much slower since laminar diffusivity alone will be required to reach over very large length- scales with only one interface. Favre-Averaged Mass and Species Transport As mentioned above, Favre-averaging employs a density-weighting to simply 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.108 f Using this notation, an arbitrary property may be decomposed in terms of mean and fluctuating components as q q q A.109 Favre averaging has the same properties given in Eq. A.103 and also leads to the identities: q f q / 0 A.110a qq A.110b 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.111 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.45a), yields the following if written in terms of the mean variables: f u j uj f x x A.112a j j f u j A.112b 0 x j 31 While both expressions 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. Such simplifications are the primary reason for the Favre-averaging popularity for stratified flows. Often, the numerical solution of a RANS PDE is conducted using a “pseudo-unsteady” formulation whereby the unsteady term is retained, even though it should be zero if there is no macroscopic unsteadiness. For example, the above variable-density mean continuity equation can be written as f f u j 0 A.113 t x j Retention of the unsteady term allows the PDE to be conveniently solved with a time-marching scheme. In particular, the unsteady terms allows for changes of the quantities from a “guessed” initial flow-field. If one employs steady boundary conditions and integrates in time for long periods, the unsteady effects will be swept out of the computational domain resulting in a converged steady-state solution. Furthermore, since time-accuracy is not needed, the time- integration can be accelerated using local time-stepping as will be discussed in §B.3.2. This “pseudo-unsteady” formulation is quite commonly used for numerical solution of both Reynolds-average and Favre-averaged transport equations. Application of Favre-averaging to the concentration transport Eq. of A.51a again yields an averaged PDE that is simpler than its Reynolds-averaged counterpart (which would contain additional correlations based on density fluctuations): f k t x i f k u i f ,k f k x i A.114 As in Eq. A.113, this includes an unsteady term on the LHS which will tend to zero for a steady-state solution. 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 u j p u i u j u u j f i f u i f g i f f i A.115a t x j x i x j x j x i x j f u i f u i u j p u i u j 2 f g i f 3 u k f uu j A.115b t x j x i x j x j x i i In both formulations, the LHS represents mean momentum convection and retains the Eulerian time derivative to allow a “pseudo-unsteady” formulation. 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, e.g. uuj . This term is often i referred to as the “turbulent stress” or the Reynolds stress and is discussed in the following using the example of a boundary layer. 32 For shear flows with gradients in the mean velocity, turbulent transport is typically the dominant mechanism for momentum diffusion. In the case of a uniform density 2-D boundary layer in the x-direction, the largest term in this tensor is the cross-correlation: u u . This x y component is primarily responsible for the spatial growth of the boundary layer (Fig. A.8). In particular, the turbulent fluctuations transport upwards the low-speed fluid from the near-wall region thereby diffusing the mean velocity gradients outward and thickening the boundary layer. This can be qualitatively described by considering a boundary layer with a mean flow (ux) in the x-directions as shown in Fig. A.18. The presence of an upward fluctuation velocity ( u y >0) will tend to entrain low-speed fluid from below so that it will be associated with u x <0 at the new location which has a higher mean speed. On the other hand, a downward fluctuation velocity ( u y <0) will tend to entrain high-speed fluid from above and so will be correlated with increased streamwise velocity ( u x >0) at the new location. In both of these scenarios, the net transfer of momentum will be associated with u u <0 (while flows with no mean velocity x y gradient will not have this correlation, i.e. u u =0). The finite anisotropic correlation will x y cause the boundary layer to thicken, as is also the case for free-shear layers such as jet or wakes. There are two principle ways to model the velocity correlation in terms of mean flowfield properties: turbulent viscosity models and Reynolds stress models. The most common approach is to use turbulent viscosity which was introduced in Eqs. A.106 and A.107 to relate turbulent mass transport to the concentration gradients, based on the analogy with laminar mass diffusion. A similar approach is used for turbulent momentum diffusion, based on the Boussinesq (1877) assumption that it is proportional to the mean velocity gradients based on the analogy with the laminar counterpart (Eq. A.9): u u j 2 K f uuj t i k x j x i 3 f ij A.116a ij i u u j 2 u k 2 K f uu t i k ij i j x j x i 3 ij x k 3 f ij A.116b The turbulence velocity correlation is typically negative as discussed above so that the LHS term is typically positive. For algebraic turbulence models, the kinetic energy is not known so that the last term on the RHS of these two formulations is neglected. However, it has been added in recent turbulence models which allow for a transport PDE of kinetic energy to allow consistency with the identity of kinetic energy (Eqs. A.71 and A.72) where jk is the Kronecker delta defined in Eq. A.10. Note that the kinetic energy in the Favre-formulation must be 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 t i f kij x j x i 3 A.118 t x j x i x j 33 Therefore, the Boussinesq assumption effectively assumes that the laminar and turbulent momentum transport have the same form except with different effective viscosities. This assumption also applies to the variable-density momentum equation is similar to that used for mass diffusion PDE of Eq. A.107. It is important to note that Eqs. A.106 and A.116 are empirical approximations since the Boussinesq assumption does not 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. 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.64 for creeping flow, and indicates gradients in the kinetic energy (non-homogenous turbulence) will yield a mean pressure gradient. This expression 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 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.5) 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. Reynolds Stress Closure 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- 34 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 grid-generated turbulence in a wind tunnel of constant cross section. At a downstream station away from the wall, the mean velocity will be uniform with no mean velocity gradients but there will be turbulence which is slowly decaying. An algebraic model bases its stress model on Eq. A.116 but without the last term (associated with k ij ). As such, this equation doe not predict any turbulent viscosity and has no “information” of upstream turbulence. To avoid this problem, most modern turbulence models use at least one transport equation for a turbulence quantity. 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 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 t f ck 2 / A.120 This form introduces an empirical eddy-viscosity constant (c) generally taken as 0.09. Combining this with Eq. A.106 relates the concentration–velocity correlation to the turbulence properties and the mean concentration gradient: c k 2 u t f A.121 x i Sct x i i For variable-density flows, the turbulent viscosity is approximated similar to Eq. A.120: t f t 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.55 and time-average the result. This provides a Reynolds-averaged transport equation for the turbulent kinetic energy: 35 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 expression 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.15. 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. 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 t ij A.125a t x i f x i x i x i 1.44 u j 2 ui K f t 1.92 x i x i A.125b t x i f k 1.3 x i ij 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 t 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 t 1.92f x i 1.3 x i k 36 These expressions 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 Eqs. A.111 and then modeled by the density gradients using A.121: f f c k 2 f u u i u i u / f t A.127 x i Sct 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/ 4 k3/ 2 / A.128 This 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: c / mix 5cc A.129 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. 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, Sct, 1.44, 1.92) 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.9). 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 expression of Eq. A.116 are avoided altogether. 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 37 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. A.5.6. 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.7a). 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.15, where it can be seen that DNS resolves the full spectrum, LES resolves the spectrum with most of the kinetic energy and RANS does not resolve any pf 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 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 38 (unfiltered) and unresolved (filtered) components. The spatial-averaged value of any flow variable (q) can be defined as follows 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 simply as a sharp cut-off given by 1/ G if |x x'| G / 2 G(x) A.132 0 if |x x'|>G / 2 For this filtering, u 0 so that the spatial-averaging is similar to the time-averaging, i.e. i ui u j ui u j uuj i A.133 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. It should be kept in mind that the LES equations are fundamentally unsteady and three-dimensional and thus can be computationally intensive, whereas the converged RANS solution is steady and can be two-dimensional. 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.115 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.134a x i x i 3 i k k 1 u i u j u i u j c G 2 2 A.134b 2 x j x i x j x i In this expression 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 also be modeled with a one-equation transport like Eq. A.125a (Yoshizawa & Horiuti, 1985), from which the sub-grid stress is similarly given by assuming the filter-scale is on the order of the sub-grid mixing length ck k G 0.05 k G A.135 The turbulent dissipation can be similarly modeled as k 3/ 2 / G A.136 39 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.137 The resolved kinetic energy can be obtained by taking a time-average of the unsteady resolved velocity fluctuations. 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.138 G This expression 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 of Eq. A.132 and dynamic sub-grid models whereby c is no longer a constant, but a function of the mean 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. Ideally the spatial filtering is applied at sufficiently small-scales with high-order spatial discretization so that the filtered turbulence is at or below the inertial range and thus is nearly homogeneous and isotropic so that the sub-grid modeling is not significant to improve accuracy. As such, the key performance aspect of LES is the filter width (G), which is generally related to the discretization length-scale () which is based on the local computational cell volume ( ) as: 1/ 3 A.139 A common prescription for the filter-width is set it equal to the discretization length-scale (Garnier et al. 2002): G A.140 For turbulent dispersion, it is desirable that the grid is sufficiently resolved so that the filter Stokes number is greater than unity (Eq. 4.71). However, this must be balanced by the computational cost associated which rapidly increases with the number of continuous-phase nodes. An estimate of the number of nodes needed as a function of ReD is given in §4.3.1. Hybrid RANS/LES Approaches Several hybrid RANS/LES approaches have emerged which intend to treat the separated flow regions are treated 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 40 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: y wall,DES y min c DES , wall A.141 Thus larger wall distances allow a Smagorinsky-type viscosity whereas smaller wall distances revert to the Spallart-Allmaras (1994) turbulence model. Typical values for cDES are 0.65 or less. 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.7a), 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.7b). An example of a DES computation is shown in Fig. A.19 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 with 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 model produces both resolved and sub-grid (unresolved) turbulence which together mimic the actual turbulence, i.e. k = k res + k hy A.142 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 large-scale turbulent properties denoted as kRANS, RANS, and RANS. For example, the large-scale length scale is: RANS max 6.0 RANS /, k 3/2 /ε RANS RANS A.143 In this expression, is the local mean/resolved vorticity. This length-scale is used along with the sub-grid length-scale (, defined by Eq. A.139) and a clipping function to compute the hybrid turbulent kinetic energy: k hy = k RANS f hy A.144a 1 RANS 4/ 3 2 4/ 3 f hy 1 tanh 4/ 3 RANS 2 4/ 3 A.144b 2 The hybrid eddy viscosity can then be calculated using hy RANS f hy 1 f hy A.145a min c hy k hy , RANS A.145b 41 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.20, which reasonably reproduces Eq. A.142. 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 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. 42 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 hexahedra (“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 and indicial means of numerical communication along gridlines which allows significant numerical efficiency for many numerical methods in terms of CPU required 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 a specific spatial resolution. One can identify a local grid resolution length-scale in various directions and x=y=z for an isotropic 3-D grid, but an average resolution can be based on cell volume (e.g. x3=xyz). The length-scale in each direction should be small-enough such that grid- dependent effects do not contaminate the solution. However, increasing grid resolution (smaller x) will also increase the number of unknowns for a given domain. Furthermore, time-dependent simulations must often use a time-step which is proportional to grid size. Therefore, reducing x by a factor of two in a 3-D domain will approximately increase the number of elements by about eight-fold. Furthermore, the time-step scales with grid length- 43 scale (to be discussed in §B.3.2) so that the number of time-steps operations will increase by a factor of sixteen. As such, there is an important balance between 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 resolution than others (e.g. shock waves) or is geometrically complex, then structured grids can become inefficient in terms of the number of elements required and/or difficulty of generation. Unstructured grids do not use grid lines and instead the domain is discretized one element at a time, typically start from boundaries where grid resolutions are specified. As such, hybrid grids can adapt to complex geometries or employ local refinement regions efficiently, i.e. with a reduced number of computational cells as compared to structured grids. As shown in Fig. B.1, unstructured grids often employ triangular elements in two-dimensions. Three- dimensional unstructured elements can take on all the shapes noted in Fig. B.2. This allows adaptation to a wide variety of domain geometries. Because of these attributes, they are becoming increasingly important in computational fluid dynamics. One downside of unstructured grids is that the lack of grid lines of communication 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. Hybrid grids are unstructured in some regions and structured in others. This characteristic allows them to use structured grids for regions with simple geometries (e.g. flat boundary layer regions) while providing reduced number of cells in regions with complex geometry or with localized flow gradients. However, most unstructured and hybrid grids do not allow implementation of third-order or higher spatial discretization schemes which can be beneficial for flows with complex structures, e.g. direct simulation of turbulence. The order of accuracy of a computational method is related to the manner in which the variables are discretized within a grid element, as discussed in the next sub-section. 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= q i i B.1 i 1 In this expression, 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. 44 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 N B.2a i 1 i 1 for cell conservation i (x j ) ij for nodal identity B.2b In this expression, 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 x1 B.3a x* x 2 x / x 2 x1 B.3b 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 * , 2 1 x * B.4 This result combined with Eq. B.1 simply indicates a linear interpolation between nodal values, e.g. the value of pressure half-way between two nodes is given by N p x*=1/2 = pi i p11 p 2 2 1 2 p1 p 2 B.5 i 1 Similarly, the pressure is recovered at the nodes, i.e. p x*=0 =p1 and p x*=1 =p 2 . To allow second- order polynomial variations within a cell one may use quadratic elements. Such a 1-D element is shown in Fig. B.3b which requires three nodes per element. 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 , 2 4 1 x * x * , 3 2x * 1 x * B.6b 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 45 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 can be computed by replacing the x2 and y2 values in Eq. B.8 with x and y. From this the associated shape function id given by A2/A. For a quadrilateral element, A can be computed by dividing it 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 local coordinate system (x*,y*) by defining the local nodes as x*=±1, 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 left face. Based on such a local coordinate system, the four shape functions are: i = 1 x*1 y* . 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.9 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 y3 z2 B.10 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.11 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.3b, B.7, B.10, or B.11, one may define a length scale associated with a computational grid as B.12a x for 1-D grids A 1/ 2 for 2-D grids B.12b 46 B.12c 1/ 3 for 3-D grids 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. 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 number of 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. This is illustrated for 1-D elements with linear shape functions in Fig. B.5. For example, global shape function D2 is non-zero within element 1 and element 2 but zero within element 3 while global shape function D1 is non-zero only within element 1. As such global shape functions associated with internal nodes have a tent-like shape which peaks at unity at the node itself. In two-dimensions, the global shape functions will peak at unity at the corresponding node and then decrease within the corresponding elements, i.e. linearly decrease to zero at the immediate neighboring nodes for linear shape functions. Beyond these nodes, the global function will remain zero. Based on this definition, Eq. B.1 can be expanded to encompass the entire computational domain as: Nf q= q i Di B.13 i 1 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.55 can be written in terms of the flow variables and space as Du L u, p, x f p f g f 2u B.14 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. 47 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.15 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 first discussed, 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.15 becomes Nf L q qi Di , xi L q, x j 0 B.16 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 locations, 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. Note that 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 the local computational volume () and zero elsewhere. Thus, Eq. B.15 becomes Nf N L q qi Di , xi w jdD L qi i xi dk 0 i1 B.17 i 1 This approach thus requires integration over each computational element (i.e. k=1,N e) for which the flow variables spatial distributions are governed by the local shape function distributions (to be discussed in §B.4). As with all the local weighting function methods, 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 (k=1,Ne) becomes N N jL qi i dk 0 i1 j1 B.18 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 Galerkin and Collocation versions (Table B.1). The shape function and weighting function influences the characteristics of the method and its typical usage (Table B.2). The node-based finite 48 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 must be derivatives are 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.19 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.20 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.21 x i x 2 x i 2 x 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.20 which are neglected in the RHS of Eq. B.21 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.22 x i x In both cases, the truncation error is linearly proportional to the discretization, i.e. first-order accurate. 49 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.22 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.20 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.23 x i 2x The differences between the computational stencils of Eq. B.21 and B.23 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.24a x i 2x u 3u i 4u i-1 u i 2 O (x 2 ) B.24b 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 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.23. Note that fixed boundary conditions which simply specify a constant, e.g. q1=const., can be applied directly without no 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.20 and its counterpart). Doing so, eliminates the first derivative terms and yields an expression involving nodes at i+1, i and i-1: 2q qi+1 2qi qi 1 x i2 4q q 2qi qi 1 2 = ... i+1 O (x 2 ) B.25 x i (x) 12 x i (x) 2 4 2 50 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.26 Therefore, the error of the scheme given in Eq. B.21 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.23 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.21 and B.25). 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.27 With such stretching the second derivative finite-difference expression of Eq. B.25 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.28 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. 51 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.29 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 (e.g., Eqs. A.113, A.114, A.115, etc.) 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 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 the following some simple explicit and implicit schemes are discussed. Explicit Schemes 52 To consider various finite-difference expressions, let us consider Eq. A.56b 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 u B.30a t x u u2 f 2 B.30b t 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 expressions, the simplest temporal integration scheme is to evaluate all terms at the current time-step and a nodal location: u u n n u B.31a t i x i n u 2u n f 2 B.31b t 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.20): n u u n 1 u in t 2 u u n 1 u in n i 2 + O (t 2 ) i + O (t) B.32 t i t 2 t i t In this expression 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.31 can be evaluated with second-order spatial (central) difference schemes about location xi. The resulting discrete equations for Eqs. B.31a and B.31b become: t u in 1 =u in (u in1 2u in u in1 ) O (x 2 , t 2 ) B.33a (x) 2 f t u in 1 =u in u (u in1 u in1 ) O (x 2 , t 2 ) B.33b 2x The resulting one-step finite difference equations indicates that un+1 can be obtained at each location directly as a function of the previous unknowns, i.e. these are “explicit” schemes. This is convenient to implement numerically. Note that these schemes are commonly are first- order accurate in time with respect to acceleration (Eq. B.32) but that this converts to second- order in time with respect to velocity (Eq. B.33) because of the multiplication by t. While these two explicit schemes seem reasonable to construct and could be linearly combined to solve Eq. B.30 they each have problems related to numerical stability. The diffusion-based scheme of Eq. B.33a 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 scheme (Eq. B.33b) in that it is unconditionally unstable and thus 53 never appropriate! In the following, explicit convection schemes are considered which are conditionally stable and thus appropriate. The unconditional instability of Eq. B.33b may be traced to a fatal flaw in the mathematical character of the convection term. In particular, this form leads to upstream transfer of information which is non-physical. To see this consider u∞>0. In this case, the influence of the downstream value u in1 on u in 1 indicates that information in this scheme is moving in the direction of -u∞, whereas only upstream information should influence u in 1 . This instability problem may be conditionally solved by “upwinding” the spatial discretization. In particular, one may use Eq. B.21 or Eq. B.22 depending on the flow direction, such that only physically relevant upstream quantities are employed: u in 1 =u in u (u in1 u in )t / x O (x, t 2 ) for u in 0 B.34a u in 1 =u in u (u in u in1 )t / x O (x, t 2 ) for u in 0 B.34b 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.23 and B.24a 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.35a 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.35b 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 close to boundaries. A scheme which avoids the increased stencil width and also allows higher-order accuracy include the Lax-Wendroff method, which is second-order accurate in time. This scheme differentiates the PDE of Eq. B.30a 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 u u u 2 t t i t 2 i t B.36 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.32 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.37 t i t 2 x i Replacing the new spatial derivative with a second-order accurate central difference and using central differencing for the convective spatial derivative (in Eq. B.33b) yields the one-step Lax-Wendroff method (there is also a two-step version): 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 2 ) n i 1 B.38 x x i i 2 2 54 The increased time-accuracy associated with the Lax-Wendroff method importantly results in a scheme which is conditionally stable, and one that can be extended to include the viscous term discretization of Eq. B.33a. Time-Step Limits The convection-based stability constraints for the upwind schemes (Eqs. B.34a-B.34b or Eqs. B.35a-B.35b) and the Lax-Wendroff scheme (Eq. B.38) 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 Eq. B.31a yields a disturbance decay criterion given by the Courant–Friedrichs-Lewy (CFL) limit as: t conv x / u B.39 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 adding a constant (cCFL) to the RHS of Eq. B.39 which can be adjusted empirically between 0 and 1 to ensure numerical stability while allowing reasonable time-steps. This constant is sometimes called the Courant number. If Eq. B.30a 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 then modified to account for the fastest possible wave-speed t conv x /(u a) B.40 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 of Eq. B.33a can also be obtained from stability analysis which yields: t visc ( x) 2 / 2 f B.41 This indicates that the viscous time-step is even more sensitive to grid spacing. These techniques for the diffusion terms and for the convection terms (the RHS terms of Eqs. B.31a and B.31b) can be integrated together by satisfying both stability limits. For example, one may combine central differencing for the diffusive terms (Eq. B.33a) with upwinding for the convective terms (e.g. Eqs. B.35a and B.35b) 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 55 t min t visc , t conv B.42 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.43a 1 u uy u 1 1 1 t conv x z a x y z x y z 2 2 2 B.43b cCFL t conv i, j,k t i, j,k local time - stepping for steady flows 1 2 / Re i, j,k B.43c The last expression 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.44 This global time-stepping can lead to quite small convection time-steps throughout to accommodate regions with small grid resolution (Eq. B.43a). 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. 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 56 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.45 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.19), e.g. if q represents a spatial derivative then: q n 1/2 1 q n 1 q n O (x 2 ) 2 B.46 To apply the Crank-Nicolson method to Eq. B.31b, the LHS is replaced with Eq. B.45 and the RHS is discretized with Eq. B.46, after which the spatial derivatives are approximated with central differencing in space (Eq. B.28) 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.47 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 f u i 1 , u i , u i 1 n n n 1 B.48 2 x 2 Unlike Eq. B.33a, this scheme is implicit due to the n+1 terms on the LHS and RHS. 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.48 happens to yield a tri-diagonal system which can be solved directly. For example, consider Eq. B,21c in the tri-diagonal form given by: a i q in11 bi q in 1 ci q in1 f i 1 for i=1,N B.49 In this expression fi represents the known values associated with the previous time-step and no ghost nodes are used so that a1=cN=0. 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: bimod bi a i ci 1 / bi 1 for i=2,N B.50a f i mod f i a i di 1 / bi 1 for i=2,N B.50b q N f N / b mod mod N B.50c qi f i mod ci qi 1 / b mod i for i=N-1,1 B.50d 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. 57 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.51 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): u in 1 u in f u i 1, j 2u i, j u i 1, j u i, j1 2u i, j u i, j1 n 1 n 1 n 1 u in 1/ 2 n n n 2 B.52a t t x 2 y 2 u n 2 u in 1 f u i 1, j 2u i, j u i 1, j u i, j1 2u i, j u i, j1 n 1 n 1 n 1 n 2 n 2 n 2 u in 1/ 2 i B.52b t t 2 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. For example, the Jacobi implicit solver for Eq. B.48 evaluates the newest estimate of an unknown based upon the surrounding values: f t n 1,k t n 1,k n 1,k 1 2x 2 u i 1, j u in1,k f 2 u i,1 u i,1 + f n 1, j 2y j n 1,k j u i, j B.53 1 f t / x 2 f t / y 2 In this expression, fn refers to f uin1 , uin , uin1 of Eq. B.48, while the k superscript refers to the iteration number. This scheme must be iterated to some pre-specified tolerance or change in u or iteration count (kmax). The convergence can be accelerated by the Gauss-Seidel method which simply uses the most recently updated values as one sweeps from i=1,Nx and j=1,Ny (where Nx and Ny are the number of nodes in the x and y directions). Further acceleration is possible by applying successive over-relaxation (SOR) with a tunable parameter (cSOR). The combination of these techniques applied to Eq. B.50d yields: f t n 1,k t n 1,k 2x 2 u i 1, j u in1, j f 2 u i,1 u i,1 + f n 1,k 2y j n 1,k j u i,1,k 1 n 1 cSOR u i,1,k B.54 n j 1 f t / x f t / y / cSOR 2 2 j Note that cSOR=1 simply reverts this scheme to the Gauss-Seidel, while the range 1<cSOR<2 will accelerate changes in the next time-step, and optimum values are typically on the order of 1.5- 1.8. The upper limit (cSOR=2) corresponds to an instability limit for linear systems and non- linear systems may require a lower maximum. 58 Convergence can be improved further by sweeping in alternating directions, i.e. the sweep of Eq. B.54 (i=1,Nx and j=1,Ny) is followed by a similar sweep for i=Nx,1 and j=Ny,1 and then the process is repeated. 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.4. However, compressible flows typically use explicit schemes for the convection and pressure terms, as will be discussed next. B.3.3. 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.55 t x For transport equations in conservative form, q can represent , u, etot. An explicit 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 2 ) n+1/2 n+1/2 qin+1 qin B.56 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.57a ui+1/2 1 n+1/2 2 u n+1/2 i+1 uin+1/2 B.57b In this expression, 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 these expressions into Eq. B.56 yields a central-difference scheme, similar to Eq. B.23, 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.45) can be written as q Q x Q y Qz S B.58 t x y z 59 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.59 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 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.56 to a two-dimensional explicit 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.60 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 of Eq. B.33b). One solution is to use the upwinding (similar to Eqs. B.35a and B.35b) 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.61a 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.61b This is referred to as the donor-cell scheme and q denotes a convected quantity determined (Prosperetti et al. 2007). A more accurate solution is to compute these flux terms using inviscid solutions to the discrete local flow as proposed by Gudunov in 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.62a qi+1/2,right qi B.62b 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 60 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 uq i+1/2 and uq i-1/2 , can then be used to obtain the cell values at i, i+1, etc. at the next time n n (as in Eq. B.56). 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 reconstruction 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 for 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.34) can be used for the flux differentiation of Eq. B.56. 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.60) but 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.63 1 q i+1 2 clim,i 1/2 q i+1 q i for u i+1 0 As such, clim=1 reverts the scheme to second-order accurate central difference scheme with no limiting (Eq. B.60), while eliminating the mid-point terms (clim=0) reverts the scheme to first- order accurate upwind difference scheme (Eq. B.21 or B.22). 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), min-mod limiter, 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: 61 csmooth,i 1/ 2 qi qi-1 / qi+1 qi for u i+1 0 B.64 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 scheme limiters 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 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 steep discontinuities but can be related to instabilities in some cases. Application of B.64 and B.61a to Eq. B.59, 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.65 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 schemes) 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 undershoot 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.66 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 discontinuities to be monotonically resolved over 3-4 grid cells, whereas low-order schemes would be distributed 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.67 In this form q and F can be vectors corresponding to conservation equations for mass, momentum and energy. The corresponding MacCormack predictor-corrector scheme (with coordinate node index i) is then 62 q ipred =q in (Qin1 Qin )t / x B.68a q in 1 = q in q ipred (Qipred Qipred )t / x / 2 1 B.68b 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 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.69 A similar equation can be written for the backward differencing of the predictor method. For multi-dimensional problems, the forward-differencing for the first-step is repeated for the other spatial derivatives and backward-differencing is used for all the coordinate spatial derivatives for the corrector step. For even higher-order accuracy, the four-step Runge-Kutta scheme can be used which is: q ipred1 =q in Qin t / 2x B.70a q ipred2 =q in Qipred1 t / 2x B.70b B.70c q ipred3 =q in Qipred2 t / x B.70d q in 1 =q in Qin 2Qipred1 2Qipred2 Qipred3 t / 6x B.70e Qi Qi 1/ 2 Qi 1/ 2 Qi 1 Qi 1 / 2 B.70f This scheme uses three intermediate predictions of the variable q (noted by superscripts p1, p2 and p3) 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 a 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.70f assumes a simple linear variation. A similar assumption is made for the MacCormack and Lax-Wendroff schemes. This scheme requires approximately two-fold the storage of the MacCormack and other two-step schemes. B.3.4. 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.39 and B.41 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.43c in 63 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 and sometimes use fourth-order diffusivity to avoid dispersion errors. Multiple dimensions can be addressed efficiently by using “Alternating Difference Implicit” or “Approximate Factorization”, whereby the unknowns are treated implicitly along one given coordinate direction and explicitly in the other directions for the first time-step, and succeeding time-steps specify the implicit direction as each of the other directions, followed again by the first direction, and so on. This allows the set of equations to be only simultaneous in 1-D for each time-step while ensuring overall implicit behavior for all directions for many time-steps. 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 and pressure variables are each defined on a unique stencil (Fig. B.7b) to avoid 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.56a into Eq. A.56b, 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.56a), and assuming the temporal and spatial derivatives are independent. In tensor notation (where i and j are Cartesian indices), the PISO governing equations become: u i f u i u j p f f 2 u i - f g i B.71a t x j x i 2p u f u i u j 2p f i f 2 u i B.71b x i 2 t x i x i x j If hydrostatic forces are important, the reference pressure (Eq. A.44) 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: 64 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.72a uiu j pred1 u n 2 pred1 2 p pred1 f i f f u i t x i x i x j B.72b 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.72c uiu j pred 2 n 1 f u in 2 pred2 p 2 p 2 pred 2 f f u i t x i x i x j B.72d 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.72e Note that u associated with continuity equation appears as the first term on the RHS of Eqs. B.72b and B.72d. For increased accuracy and stability, the terms in the square brackets can be 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.71a) in written terms of a predictor function fpred and a pressure gradient: u i 1 p f pred u i - B.73 t f x i Thus, fpred represents the RHS of the momentum equation uncorrected for the pressure gradient. To discretize this equation, the projection method represents the LHS using central- differencing in time and employs an initial prediction (upred) for the velocity: n 1/ 2 u i u in 1 u in u n 1 u ipred u ipred u in B.74 t t t t The prediction velocity is based on fpred and thus neglects the pressure gradient. The initial prediction for the velocity can be accomplished with 1st-order accurate schemes if steady-state convergence is the goal. For unsteady problems, higher-order schemes are preferable since temporal accuracy is critical. This can include 4th-order accurate Runge-Kutta scheme given by Eq. B.70. However, a more common approach is the Adams-Bashforth method whose 2nd-order temporal accuracy will be generally consistent the ensuing pressure implicit solver (since implicit schemes are limited to 2nd-order accuracy). The explicit Adams- Bashforth method is based on a Taylor series expansion to describe the function at the n+1 65 time-level in terms of the function at previous time levels using a forward-difference extrapolation: f n 1 f n t f pred n f n t f n pred n f pred1 2f n n f pred1 B.75 t t pred pred pred pred Central-difference temporal schemes can then be applied to Eq. B.73 at n+1/2: u pred u n t n f pred1/ 2 1 n 2 n 3 n 1 n f pred f pred1 f pred f pred1 2 2 B.76 Equating the LHS and RHS terms yields a 2nd-order accurate Adams-Bashforth temporal integration method for the uncorrected velocity vector, i.e. 2 n 1 u ipred u in 1 t 3 f pred,i f pred,i n B.77 This predictor step is conveniently explicit in that only variables at the previous time-step are needed to find the solution. Note that this step requires finite differencing of all the spatial derivatives within the function fpred. This can be accomplished with a higher-order operator, e.g. a 3rd-order upwind scheme, while the viscous stress terms can be evaluated with a central- difference scheme. An alternative predictor step to Eq. B.77 it to treat the viscous terms implicitly with a Crank-Nicolson scheme (Eq. B.48) by separating fpred of Eq. B.73 into convective and viscous portions: 2 n n 1 uipred uin 1 t 3 f conv,i f conv,i 1 t f visc,i f visc,i 2 pred n B.78 The resulting predictor step is implicit and can be solved with an iterative method, such as the ADI scheme used in Eqs. B.52a- B.52b. For the second step, a fully implicit pressure correction can be specified by combining Eq. B.73 and B.74 to yield: u in 1 u ipred 1 pn 1 B.79 t f x i Taking the gradient and assuming incompressibility for the final velocity ( un 1 0 ) yields a semi-discrete equation for pressure: 2 pn 1 f u ipred B.80 x i2 t x i This is an implicit equation for the pressure which requires a Poisson solver. One can again use an SOR technique (Eq. B.73), but this step requires special attention since solution of this pressure equation dominates the overall computational time. This is because the pressure field is often very slow to converge as its influence is generally quite global (owing to the infinite relative acoustic speed for the incompressible flow assumption) in contrast to diffusion effects which are more localized. Therefore, advanced iterative solvers are often employed for the pressure equation including pre-conditioning and multi-grid methods (Wesseling & Oosterlee, 2001). Once the pressure field based is obtained to within a suitable convergence level, the predicted velocity is corrected for the (previously missing) pressure gradient term: 66 t pn 1 u in 1 u ipred B.81 f x i This final step for the MAC method completes the momentum-based velocity field prediction of either Eq. B.76 or B.77 while indirectly satisfying the continuity equation by using a pressure solution from Eq. B.80. A variant of the above technique is to use the old pressure gradient for the prediction and then use the change in pressure (from n to n+1) for the correction. For example, Eq. B.77 can be transformed as 1 pn 2 n 1 u ipred u in 1 t 3 f pred,i f pred,i - n f x i B.82 This expression 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 which can be obtained through a Poisson equation based on Eq. B.80: p p n 1 p n B.83a 2 p u ipred f x i2 t x i B.83b t p u in 1 u ipred B.83c f x i Similar to the conventional MAC approach, the predicted velocity (Eq. B.82) is used to solve a Poisson equation (Eq. B.83b) with an implicit iterative scheme from which the corrected velocity can be computed (Eq. B.83b). This variant is quite similar to the SIMPLE method described earlier. 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 the spatial derivative terms in Eq. B.84 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.4. Finite Volume Method While spatial discretization of the differential form is called the Finite Difference 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 67 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 over each all volumes. To illustrate the FV method, consider the conservative form of the compressible Euler equations (Eq. A.45) as a vector of PDEs: q Q x Q y Q z = B.84 t x y z Where 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.85 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.12). Writing the above PDEs in FV form for a discrete volume (Δ) yields q Q x Q y Qz t x y z d 0 B.86 If we define a volume-averaged of the conservative variables as qΔ and consider an explicit temporal scheme, a first-order accurate forward differencing in time yields t Q Q y Q q1 q q xx y zz d n n B.87 This expression 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 stored at either the cell centers (Fig. B.8a) 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.4b). The RHS integral of Eq. B.87 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.88 x y Assuming a linear variation of the fluxes along the cell edges and using area of Eq. B.8, 68 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.89 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.84 (Chung, 2002). Then the contributions from the K cells surrounding a vertex are averaged together to get the net nodal correction and new value 1 K qin 1 = qin qi,k K k 1 B.90 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.87 to get qn+1 with 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. The 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.18). 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.84. 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.91 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.92 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.93 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.93 is the called the mass matrix, which does not appear in the FD or FV methods (see Eqs. B.33 & B.87). The mass matrix causes the 69 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 with the consistent mass matrix and the lumped mass matrix can be used for Eq. B.93 by changing qn+1 to qpred1 and adding the revised LHS added to the RHS while adding a new LHS is added of (qpred2-qn)A . Initially, qpred1 is set as qn and the process is repeated until convergence (typically three to four iterations are suitable for second-order accurate temporal schemes). The term in brackets on the RHS of Eq. B.93 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), Eqs. B.4-B.5 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.93, 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 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.94 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.95 Comparing this with Eq. B.33b, it can be seen that this one-step forward-difference central- space yields the same scheme 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 or 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 which have second spatial derivatives 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.93 would result in terms such as 2/x2 which becomes zero. This problem can be solved with higher-order elements, but is more generally addressed by employing integration by parts. For example, a RHS terms with the second derivative of velocity would be revised as follows: 70 2u u j u x 2 jd x x d x jdA B.96 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. Another benefit of this form, is that it introduces a boundary surface integral which can be used explicitly impose boundary conditions. Because of this, sometimes integration by parts is 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 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. 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.97 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.98 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, 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.99 i 1 71 This is similar to Eq. B.18 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.100 i 1 x =x j This is similar to Eq. B.16, 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 a direct integral 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 specific particle velocity distribution function for each fluid component. The fluid particles move in discrete time and can collide with each other under prescribed interaction forces. The 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+ei, where ei are the local particle 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 B = i B.101a i=1 B u = iei B.101b i=1 3R g T B (ei u) 2 = i 2 i=1 2 B.101c In this expression, 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, which can be written as: 72 i (x ei x, t t) i (x, t) i ( j ( x, t)) for i=1,N and j=1,N B.102 In this expression, i i ( j (x, t)) is the collision operator which represents the rate of change of each fluid particle density resulting from collisions. It may be expressed as i =- i ieq / B.103 where ieq is local equilibrium distribution function, is the equilibrium relation time-scale to be consistent with pressure and shear stresses, and the collision operator is required to satisfy conservation of total mass and total momentum at each lattice: M M i 0, i=1 e i=1 i i 0. B.104 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 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.17) 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. Based on Bird, 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 (Bird, 1976). 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- 73 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 in 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. (2003) “Numerical simulation of turbulent particle diffusion,” Doctoral Dissertation, University of Illinois at Urbana-Champaign. 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. 78 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 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. 79 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. Chang, Y.C., Hou, T.Y., Merriman, B. & Osher, S. (1996) “A level-set formulation of Euleran interface capturing methods for incompressible fluid flows,” (1992) 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. 80 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. 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 Kleinsteuer, 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. 81 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. 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. 82 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. di Felice (1994) "The voidage function for fluid-particle interaction systems," Int. J. Multiphase Flow, Vol. 20, pp. 153-159. Dorgan, A. and 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. and Loth, E. (2007) "Efficient calculation of the history force at finite Reynolds numbers," Intl. Journal of Multiphase Flow (to appear). Dorgan, A. and Loth, E. (2008) "Dispersion of bubbles in a turbulent boundary layer," Intl. Journal of Multiphase Flow (submitted). Dorgan, A.J., Loth, E., Bocksell, T.L. and Yeung, P.K. (2005) "Boundary layer dispersion of near-wall injected particles of various inertias," AIAA Journal, Vol. 43, pp. 1537- 1548, July. 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. 83 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. 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. 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. 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. 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. 85 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. 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. 86 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. 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. 87 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. 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. 88 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: sessionA10.001. 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. Physicsi, 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. 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. 89 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. 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. 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. 90 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. 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. 91 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. 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. 92 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. 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. (2008) “An equation of motion for particles of finite size and Reynolds number in a liquid," Environmental Fluid Mechanics (to appear). 93 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. 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. 94 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. 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. 95 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. 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. 96 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 suing CIP scheme,” Computational Fluid Dynamics J. Vol. 8, pp. 26-33. Najjar, F. (2005) University of Illinois at Urbaba-Champaign, personal communication. 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. 97 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. 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. 98 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. 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. 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. 99 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. 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. 100 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. 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. 101 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. 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. 102 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. Sikalo, 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. 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 DC, 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. 103 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. (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. 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. (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 104 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. 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. 105 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. 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. 106 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. 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. 107 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. 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. 108 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. 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. 109 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. 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. 110 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. 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. 111 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: | 4 |

posted: | 5/17/2012 |

language: | |

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.