Dynamics of Flow Structures and Transport Phenomena, 1

Jun 4, 2009 - Experimental and Numerical Techniques for Identification and Energy Content of Flow Structures ... Also, several notable mathematical qu...
18 downloads 8 Views 12MB Size
8244

Ind. Eng. Chem. Res. 2009, 48, 8244–8284

Dynamics of Flow Structures and Transport Phenomena, 1. Experimental and Numerical Techniques for Identification and Energy Content of Flow Structures Jyeshtharaj B. Joshi,* Mandar V. Tabib, Sagar S. Deshpande, and Channamallikarjun S. Mathpati Department of Chemical Engineering, Institute of Chemical Technology Matunga, Mumbai-400019, India

Most chemical engineering equipment is operated in the turbulent regime. The flow patterns in this equipment are complex and are characterized by flow structures of wide range of length and time scales. The accurate quantification of these flow structures is very difficult and, hence, the present design practices are still empirical. Abundant literature is available on understanding of these flow structures, but in very few cases efforts have been made to improve the design procedures with this knowledge. There have been several approaches in the literature to identify and characterize the flow structures qualitatively as well as quantitatively. In the last few decades, several numerical as well as experimental methods have been developed that are complementary to each other with the onset of better computational and experimental facilities. In the present work, the methodologies and applications of various experimental fluid dynamics (EFD) techniques (namely, point measurement techniques such as hot film anemometry, laser Doppler velocimetry, and planar measurement techniques such as particle image velocimetry (PIV), high speed photography, Schlieren shadowgraphy, and the recent volume measurement techniques such as holographic PIV, stereo PIV, etc.), and the computational fluid dynamics (CFD) techniques (such as direct numerical simulation (DNS) and large eddy simulation (LES)) have been discussed. Their chronological developments, relative merits, and demerits have been presented to enable readers to make a judgment as to which experimental/numerical technique to adopt. Also, several notable mathematical quantifiers are reviewed (such as quadrant technique, variable integral time average technique, spectral analysis, proper orthogonal decomposition, discrete and continuous wavelet transform, eddy isolation methodology, hybrid POD-Wavelet technique, etc.). All three of these tools (computational, experimental, and mathematical) have evolved over the past 6-7 decades and have shed light on the physics behind the formation and dynamics of various flow structures. The work ends with addressing the present issues, the existing knowledge gaps, and the path forward in this field. 1. Introduction The present methodology of designing the process equipment is based on empiricism, as well as knowledge accumulated from prior experience. Because of inefficient design, the chemical industries have not been able to utilize the resources (raw materials, hardware, utilities, and energy) effectively, which leads to a higher cost of production, a loss of product quality, and environmental pollution. The fundamental reason for empiricism is the complexity of the fluid dynamics in the process equipment. Hence, much work is ongoing for the understanding of the physics of turbulent flows. Because most of the industrial equipment is operated under the turbulent conditions, wherein the motion of compendium of eddies (flow structures) of different length and time scales contribute toward mixing, momentum transfer, heat transfer, and mass transfer. Hence, a proper understanding of the mechanism related to the formation and the role of these turbulent flow structures in determining the transport phenomena can cause improvement in the scaleup and design procedures. The present work reviews some prominent methodologies and techniques devoted to identifying, characterizing, and studying the dynamics of these flow structures, and their effect on transport phenomena. * Author to whom correspondence should be addressed. Tel.: + 9122-2414 0865. Fax: +91 -22-2414 5614. E-mail: [email protected].

Turbulence is generally defined as the fluctuations around the mean flow. These fluctuations are the result of passage of the deterministic organized flow structures and the random disorganized irrotational motions, which, together, constitute the turbulent flows. The organized deterministic patterns are known as eddies or turbulent flow structures. The flow structures are often hidden among the incoherent turbulent motions. Some of these structures are shown in Figure 1. These organized flow structures have been known from as early as the 15th century, as conceptualized by Leonardo Da Vinci. Brown and Roshko1 observed the roller structure in the mixing layer (Figure 1A). Richardson2 gave the first notion of a cascade process of structure breakup. He suggested that these structures are part of the turbulence (Figure 1B) and the energy is injected at larger length scale structures, and the energy is transmitted to smaller and smaller length scale structures, until it reaches a size, where viscosity is dominant, and the viscous dissipation results into direct conversion to heat. To begin with, we review the historically much debated issue of categorizing flow structures, and make an effort toward rationalization. The coherent flow structures have been defined in the following manners in the literature. (1) Hussain3,4 emphasized that coherent vorticity (i.e., instantaneous phase-correlated and space-correlated vorticity) is a primary identifier for detecting characteristic structures (also known as coherent structures) that possess a distinct boundary

10.1021/ie8012506 CCC: $40.75  2009 American Chemical Society Published on Web 06/04/2009

Ind. Eng. Chem. Res., Vol. 48, No. 17, 2009

Figure 1. Turbulence as depicted by (A) mixing layer flow structures (Roshko and Brown1), (B) Richardson2 cascade (notion of turbulence).

and independent territory. He attempted to describe these structures as a turbulent fluid mass with instantaneously phasecorrelated vorticity over its spatial extent, with spatial extent being comparable to transverse length scale of shear flow. Therefore, although length scales of the order of Kolmogorov scales are known to show intense coherent vorticity, they cannot be called characteristic structures, because of the small spatial extent. (2) While, according to McComb,5 coherent structures are discernible patterns in the flow and occur with sufficient regularity, in space and/or time, to be recognizable as quasiperiodic or near-deterministic. Although characteristic flow structures occur prominently in turbulent flows, their contribution to the total kinetic energy is only a small fraction. Hence, to understand the effect of flow structures on the transport phenomena, all the structures should be taken into account. In the literature, more emphasis is given to study the characteristic flow structures, rather than all the structures present in the flow. We have an opinion that, as far as possible, all the flow structures should be studied. Because all the flow structures, from large scale (of the order of reactor geometry) to small scale (Kolmogorov scale), contribute to mixing and/or heat transfer by advecting the species and/or enthalpy from one place to another, at both the microlevel and the macrolevel. For equipment such as ultrasound reactors, the transient flow structures are very important and have significant impact on the transport phenomena. Thus, the focus should be on characterizing each and every flow structure so that one can study their role in the overall transport phenomena, and perhaps not into their classification/categorization categorizing it. In this review, the flow structures are characterized by their (a) mechanism of formation, and their location in space and time (how, when and where they are formed), (b) spatial and temporal coherence (time and length scale), (c) their topological definition (size and shape), (d) their energy content (helps to determine strength of the flow structure), (e) their properties (coherent vorticity value, second invariant gradient values), (f) their dynamics (startup, merging, breakup, interaction with mean flow and other turbulent structures) as we change the operating and geometric conditions, and, finally (g) the effect of these flow structures on the transport phenomena and their relation with design parameters. The last point forms the basis for the second part of this review. In the present review, fluid dynamics in the process equipment (Figure 2) has been classified as flows with stationary walls

8245

(channel flow), flows with some rotating/moving walls (stirred tank, annular centrifugal extractors), flows with discontinuous moving interfaces (bubble column), flows away from the walls (free shear turbulence, free jets, open channel flows), and ultrasound reactor. These equipment are selected so that extreme flow conditions can be covered, e.g., in channel flow, major turbulence generation is at the walls, where the no-slip condition applies and other extreme is free surface turbulence, where practically complete slip conditions prevail and the generation of flow structures is away from the interface. The other cases fall in these two extremes. The extreme flow conditions can be characterized by the range of the energy dissipation rate (ε). For instance, the average value of ε in pipe flow is in the range of 4-10000 W/kg (for water flowing in a 2-in. pipe between Re ) 5600-100000) whereas it is in the range of 50-1000 W/kg in annular centrifugal contactors. The mechanism of formation of characteristic flow structures is discussed in section 2 and literature pertaining to the study of these structures using the experimental and computational tools is presented in sections 3 and 4, respectively. Mathematical tools for identification and characterization of flow structures are reported in section 5. Structures identified using these techniques in the aforementioned equipment are discussed in section 6. Finally, the effect of flow structures on energy spectra and estimation of turbulence parameters is reported in section 7. 2. Mechanism of Formation Turbulent flow structures are described as the motions of fluid parcels that have a life cycle and undergo generation, development, and breakdown. Generally, the initial formation of a flow structure is always the result of instability, i.e., on the relative contribution of forces acting on the fluid element. For example, in pipe flow, the transition to turbulence occurs when the inertial forces are dominated over the viscous forces and structures form because of shear or velocity gradients in the flow. In some cases, it is the amplification of these instabilities that leads to the formation of flow structures with a definite pattern. Then, the vortex stretching mechanism is responsible for producing higherintensity smaller vortices. The stretching of these vortices, because of velocity gradients, leads to a reduction in the sizes of the vortices, and a subsequent increase in their vorticity, to conserve the angular momentum. In this section, we discuss the mechanism of formation of flow structures at some distinct locations: near the solid/fluid interface, near the free surface, behind bluff and blunt bodies, in free shear region (mixing layer and free jet region), and in the multiphase systems (dispersedphase-induced turbulence). 2.1. Solid/Fluid Interface. Turbulent structures near the wall have been studied comprehensively for the last four decades.5-9 Professor Kline’s group6-8 first discovered these structures using a hydrogen bubble technique in the turbulent boundary layer, and called them “near-wall coherent structures”. They called the life cycle of these structures as a “burst”, and the accompanying phenomena as the “bursting process”. They described “bursting” as a phenomenon comprised of regular ejections (outward motion of fluid from the wall) and sweeps (inward motion of fluid toward the wall). They proposed the following mechanism. The total process of bursting is a continuous chain of events starting from a relatively quiescent wall flow and leading to the formation of relatively large and relatively chaotic fluctuations. The first stage of bursting is the lifting of flow structures away from the wall. The fluid motion in the near-wall region is unsteady, and the flow structures also move in the streamwise and wall-normal directions (and they

8246

Ind. Eng. Chem. Res., Vol. 48, No. 17, 2009

Figure 2. Schematic diagram of all the equipments: (A) channel flow, (B) jet loop reactor, (C) stirred tank, (D) Taylor-Couette flow (annular contactor), (E) ultrasonic reactor, and (F) bubble column.

are also called “low speed streaks”). This movement is caused by the local instability in the flow. Once the low-speed streak reaches a critical distance from the wall (at x+ 1 ≈ 12), it appears to turn sharply outward, away from the wall. They then continue to move downstream (Kim et al.7). This more-rapid outward motion is called “low-speed-streak lifting” and is also characterized by low secondary (streamwise) vorticity. Thus, the relatively rapid outward motion creates an inflectional point in the instantaneous velocity profile as very rapid changes in the instantaneous velocity profile10-12 (see Figure 3A). The instantaneous inflectional profiles lead to the growth of an oscillatory disturbance just downstream of the inflectional zone. The dominant mode of this oscillation is the streamwise vortex motion. This is also the reason for a peak in turbulent kinetic energy profile very close to the wall (x1+ ) 10-20; see Figure 3B). An increase in the total kinetic energy occurs during the bursting times, but most of this added energy is in the frequency range near the frequency of oscillatory growth. The oscillation is terminated by the start of a more-chaotic fluctuation called “breakup”. The beginning of the breakup phase indicates the return of the instantaneous velocity profile to a form qualitatively similar to that of the mean profile, including the vanishing of the inflectional zone. The ejection phase contributes to increased interaction with the bulk and, hence, mixing, whereas the sweep phase contributes to surface renewal and, hence, heat/mass transfer. This entire cycle repeats randomly in space and time. The maximum turbulence production is observed primarily during the ejection phase and somewhat less during the sweep phase of bursting. Systematic simulations of homogeneous turbulence with arbitrary shear at random locations in the bulk have shown similar structures, proving that shear at the interface is a important for the formation of such structures.13 2.2. Free Surface Turbulence. The fluid/fluid interface where the shear at the interface is negligibly small enough is called a free surface. The study of transport phenomena

governed by flow structures at free surfaces is of great interest to design contacting equipment (such as evaporators, condensers, and gas absorbers) and to study the interactions between atmosphere and sea/river for understanding the CO2 uptake.14 Mass transfer at these free interfaces (the relative velocities between two phases can be nonzero) plays a central role in many industrial and environmental processes. Generally, in the design of bubble columns, the contact time for surface renewal is taken approximately as the bubble diameter divided by the bubble rise velocity. However, for the accurate estimation of masstransfer rates, an understanding of the turbulence close to free surfaces (a distance on the order of a few micrometers from the surface) is important. This area is still in a state of development, because it is very difficult to measure and/or compute the flow accurately in a region very close to the interface. Many theories have been proposed in the literature to qualitatively represent the free surface turbulence, namely, film theory, penetration theory, surface renewal models, surface divergence models, etc. Systematic measurements and computations by the research groups of Komori and co-workers15-19 and Banerjee and co-workers20-23 have highlighted many features of free surface turbulence. These studies had been mainly performed in open channel flows where at the top interface slip boundary conditions could be manipulated from complete slip (zero shear) to very close to no slip (very high shear). There are persistent large-scale flow structures at the free surface (Figure 4A), comprised of the “upwellings” (caused either by the impingement of bursts emanating from the bottom wall or by the vortices emanating from the bulk), the “downdrafts” (where the flow from adjacent upwellings meet), and the “whirlpool like” attached vortices that form at the edge of the upwellings. In the case of a complete slip boundary, the upwellings are observed at the interface, apparently because of the bulk turbulence or turbulence at the bottom wall (Figure 4B) (see Rashidi and Banerjee20). The upwellings cause the

Ind. Eng. Chem. Res., Vol. 48, No. 17, 2009

8247

Figure 3. (A) Instantaneous velocity profiles close to wall using the hydrogen bubble technique (see Kim et al.7). (B) Turbulent kinetic energy profile in channel flow in the wall-normal direction at different Reynolds numbers (data from Kawamura and co-workers10-12) (legend: line 1, Re ) 5600; line 2, Re ) 14000; line 3, Re ) 24400; and line 4, Re ) 41400.

surface to diverge and surface renewal occurs. These upwelling vortices are energetic and play a significant role in redistributing heat and mass at the interfaces. The edge of the upwellings and the region between the two upwellings are high-shear regions that cause instabilities, which lead to the formation of these whirlpool-like vortices. These attached vortices tend to pair or merge with each other and dissipate very slowly, provided they are not destroyed by other upwellings. These types of vortices are unique to free surfaces and are not found at the solid interface. With the incorporation of shear at the interface, lowspeed and high-speed streaks start to appear at the gas/liquid interface (see Figure 4C). In addition to upwellings, the streaking phenomena also help in heat and mass transfer at the interface. The mechanism of formation of these streaks is similar to turbulence at solid interfaces and the occurrence of these structures is very well-correlated with shear rate at the interface. There are several papers and reviews (Komori and coworkers15-19 and Nakagawa and Nezu24,25), either explicitly or implicitly touching upon the subject of the present discussion. 2.3. Free Shear Flows. In some applications, a fluid is injected through the nozzles in the pool of same fluid or a completely miscible fluid. In this process, the ambient fluid gets

carried along with the jet by viscous drag at the outer layer of the jet. The rate at which fluid from the jet and from its surroundings become mixed at the interface is given by the entrainment rate of the jet, which defines the rate of propagation of the interface between the chaotic and relatively stagnant fluid. This entrainment rate is controlled by the speed at which the interface contortions with the largest scales move into the surrounding fluid.26 As the jet emerges from an axisymmetric nozzle, a shear layer is formed between the jet stream and the surroundings. In the regions nearer to the jet exit, there is an exponential growth of small perturbations, because of the early linear instability jet regime. In regions away from the jet region and beyond the early development stage, the nonlinear Kelvin-Helmholtz instability regime takes over, leading to transitional shear flow. In this regime, the large-scale vortex rings roll up and merge and demerge. Further downstream from the nozzle, the streamwise vorticity component and the threedimensional azimuthal instabilities dominate, leading to the transition to the turbulent jet regime. In the turbulent regime, the azimuthal instabilities cause the breakdown of the vortex rings. The vorticity dynamics witnessed in this regime causes self-induction, vortex stretching, and vortex reconnection, which

8248

Ind. Eng. Chem. Res., Vol. 48, No. 17, 2009

Figure 4. (A) Free surface turbulence (see Turney and Banerjee23). (B) Schematic of turbulence mechanism in open channel flow (see Rashidi and Banerjee20). (C) Instantaneous velocity profile using the hydrogen bubble technique in the presence of shear at the interface (from Rashidi and Banerjee20).

leads to an additional production of vorticity. The jet development is determined by the dynamics of large-scale vortex rings and braid vortices, and by strong vortex interactions that result in a more-disorganized flow regime characterized by smallerscale elongated vortex tubes.27 The time-averaged and instantaneous flow structures in jet loop reactor28 are shown in Figures 5A and 5B, respectively. An instantaneous LES snapshot is shown in Figure 5B. The instantaneous flow shows the following observations: (a) The flow showed major convective motion in the jet region. The shear region shows the formation of small circulation cells in the shear region. (b) The structures penetrate until the bottom wall and break. (c) The jet region shows fluctuations of the jet plume that are due to the back pressure exerted by the bottom wall, which is the result of the jet impact on the bottom. This causes the jet instability (JI) in the flow. (d) The bulk region shows the motion of eddies 0.04 m in size along the circulation path, and they disappear after a travel length of ∼0.1-0.15 m. 2.4. Flow Past Solid and Blunt Bodies. In the flow past the solid bodies, the two extremes are (i) parallel flow over a streamlined body29 (Figure 6A) and (ii) a circular or rectangular

disk normal to the flow29 (Figure 6B). The intermediate between these extremes is flow past a cylinder29 (Figure 6C). Whenever any such obstacle is placed in a free stream, at very low Reynolds number (Re < 0.01), flow just creeps on the surface. With an increase in Re, the flow gets separated from the surface at some point. The separation point is dependent on the shape of the surface. Beyond this point, the flow structures get formed in the wake. At low to moderate Re values, the zone of separation contains a rather stable flow structure, in which circulatory motion is maintained through the transmission of shear stress across the dividing streamline. Therefore, the velocities within the flow structure are considerably below those of the surrounding flow. As the Reynolds number increases further, essentially the same relative velocities are maintained, but the flow structure becomes more unstable. Eventually, a point comes where these structures tend to grow and detach themselves from the boundary and pass off into the wake, as new flow structures form to take their place. At very high Re values, this process is extremely rapid and complex. Because the separation occurs where the velocity of the surrounding flow is highest, the pressure at the rear of the blunt body is lower than that at the front. The separation produces a net force in the direction of flow and, hence, increases resistance between

Ind. Eng. Chem. Res., Vol. 48, No. 17, 2009

8249

Figure 5. (A) Mean velocity profile in JLR. (B) Instantaneous flow structures in JLR. (From Deshpande.28)

the body and the fluid. This resistance is form drag. The kinetic energy of the flow structures that pass off into the wake cannot be regained, and their repeated production represents the steady drainage upon the energy of the flow. In the case of a circular disk normal to the flow, the separation begins at the sharp edge of the disk (Figure 6B). Knowledge of the pressure at various boundary points is often essential, because the role of skin drag is negligibly small. The force exerted by a flowing fluid on an immersed body is a combination of tangential or shear stress, normal stress, or pressure. If the body is a thin plate parallel to the flow, the entire force will be tangential and plate at right angles, and the effective force will be normal (when it is parallel the entire force is due to shear stress). However, for morecomplex geometries, integration of the shear and pressure over the entire exposed surface is required. The total drag force is related to the kinetic energy head and the proportionality constant is the drag coefficient. In the case of a cylinder perpendicular to the flow, very interesting structures are formed (Figure 6C). At the beginning of the separation behind a cylindrical body, the separation zone includes two symmetrically arranged eddies. As the Re increases, one eddy or the other has a tendency to grow more rapidly and to pass off into the wake, while the remaining one is still in the process of development. The wake eventually takes the form of alternating series of eddies, the formation of which entails not only a longitudinal force upon the body but also a transverse force that repeatedly changes in direction. Above a critical Re value, a wake behind the cylinder is unstable and oscillations are seen in which velocity is periodic, with regard to time and the downstream distance. The oscillating wake rolls up into two staggered rows of vortices with opposite sense of rotation; these are also known as “Ka`rma`n vortex sheets”. 2.5. Taylor-Couette Flow. Flow between two concentric cylinders, with either or both of them rotating, is termed as the Taylor-Couette flow. When the rotational speed is sufficiently high, the flow becomes turbulent and can be used to disperse one liquid into another. The flow is circumferential due to the

Figure 6. Flow patterns over a bluff body: (A) flow over a streamlined body, (B) flat plate perpendicular to the flow, and (C) flow past a cylinder at Re ) 10000. (From Homsy et al.29)

shearing action of the rotating cylinders. With an increase in rotation speed, centrifugal instability develops. Three-dimensional secondary flows are generated because of this instability. These vortex patterns are of great interest, because they provide very high values of heat- and mass-transfer coefficients. These flows are characterized by the Taylor number (Ta), which is the ratio of the centrifugal forces to the viscous forces. These flows have been extensively investigated since the end of the 19th century. Many investigators30-32 have performed a theoretical analysis of the stability of flow patterns to estimate the critical Taylor number (Tacr). This signifies the transition from Couette flow to Taylor-vortex flow. As the Ta value continues to increase beyond Tacr, for the case of a stationary outer cylinder, the flow structures change to wavy vortex flow (WVF), then to chaotic vortex flow (CVF), and finally to turbulent Taylor vortex flow (TTVF) (see Figure 7).33 The Taylor-Couette flow has been extensively investigated since the end of the 19th century and state-of-the-art reviews have been presented by Koschmieder34 and Vedantam and Joshi.35 2.6. Particle (Dispersed-Phase)-Induced Flow Structure. Here, the term “particle” has been used to represent solid particles or gas bubbles or liquid drops. In equipment such as bubble columns, spray columns, and fluidized beds, the dispersed phase (bubbles, drops or solids) imparst energy to the liquid by means of momentum exchange as they rise/settle. This

8250

Ind. Eng. Chem. Res., Vol. 48, No. 17, 2009

Figure 7. Flow structures in Taylor-Couette flow (from Deshmukh et al.33).

imparted energy ultimately sets the liquid into motion and can be measured as the total kinetic energy of the liquid, which is comprised of the kinetic energy due to mean velocity and the turbulent kinetic energy. Both of these quantities are important, because they dictate the performance of these reactors. Bulk motion is mainly responsible for mixing/dispersion, whereas turbulent motion has a greater influence on interfacial heat and mass transfer. Understanding the sources of the fluctuating velocity components, in terms of frequency (time scale) and size (length scale) of flow structures, is extremely important, because the energy imparted to cause these turbulent fluctuations ultimately dissipates into heat at the smallest scale through viscous dissipation. Joshi et al.36 has discussed the chronological development in the understanding of bubble-induced coherent structures in the bubble column. In the bubble column, the mechanism for the onset of large-scale structures is related to the reversal in the direction of the lift force that is acting on the bubbles.37,38,121 The lift force is responsible for the radial hold-up profile (see Figure 8A) and the density gradients, which cause the circulations within the bubble column (see Figure 8B). The magnitude and direction of lift force is dependent upon the established velocity gradient in the liquid phase and the bubble diameter. The size of the bubble formed at the sparger is determined by the balance of four forces: the gas momentum and the buoyancy forces that push and expand the bubbles in opposition to the liquid drag force, liquid inertial force, and the surface tension force. At lower superficial gas velocities, uniform bubbles are formed at the sparger (3-4 mm in diameter) and there is no further coalescence and breakup. The lift force acts on the bubble toward the wall direction. This results in a homogeneous regime having a flatter holdup profile and results in weaker circulation flow, because of a lower density gradient (the driving force for circulation). Thus, the large-scale structures

are absent and the flow structures are of the size of bubblebubble spacing. As the superficial gas velocity is increased, the size of bubble formed at the sparger increases, as a result of an increase in liquid drag force and inertial force (see Joshi et al.39). The higher drag often results in increased bubble coalescence and larger bubble sizes. These larger bubbles experience a lift force in the direction of center of the column, resulting in the creation of hold-up gradients and the density gradients that drive intense liquid circulation. At the top of the column, because of the reduction in hydrostatic pressure, the bubble sizes are larger than that those at the bottom. This local bubble size distribution causes the setup of local density gradients and local circulation patterns (flow structures). These flow structures may, in turn, trap and disperse many small-sized bubbles. Depending on the local void fraction within the flow structure, they may have their own local density. If this density is relatively higher than that in other regions of the column, then they may move down; otherwise, they may move up in the column. Thus, at higher superficial gas velocity, a heterogeneous regime exists that has a wide bubble size distribution and several vortical structures. These vortical structures mostly lie in the midregion between the central plume region and wall region, and they move up and down the column. To suggest the size, number, and locations of circulation cells formed in the bubble column, Joshi and Sharma40 and Joshi41 and Joshi et al.36 proposed two models: the multiple noninteracting circulation cells model and the interacting cells with a considerable intercirculation model (see Figures 8C and 8D). Similar to the bubble column, the fluidized bed and spray columns may experience distributions in particle and bubble sizes, respectively, leading to local density gradients and the formation of flow structures.

Ind. Eng. Chem. Res., Vol. 48, No. 17, 2009

Figure 8. Flow structures in the bubble column (from Joshi et al.36): (A) hold-up profile in radial direction, (B) gross circulation cells in bubble column, (C) multiple noninteracting cell model, and (D) multiple interacting cell model.

3. Experimental Tools Experimental fluid dynamics (EFD) has played a very important role in identification and characterization of flow structures. These techniques can be classified as intrusive and nonintrusive techniques or as point measurement techniques, planar measurement techniques, and volumetric measurement

8251

techniques. In the present work, a second methodology of classification has been adopted. Few techniques (such as smoke injection, dye injection, hydrogen bubbles and photochromic tracers, and fluorescent particle tracing) are widely used for flow visualization and less for quantitative characterization. These techniques are listed separately. Table 1 gives a brief account of various techniques, as well as their principles, advantages, and limitations. The majority of the studies, to date, have involved the use of point and planar techniques; however, recently, holographic and stereoscopic particle image velocimetry (HPIV and SPIV, respectively) techniques are getting much attention. The performance of different major measurement techniques, in terms of qualitative and quantitative characterization of flow structures, is summarized in Table 2. The application of the aforementioned techniques for different flows under consideration are discussed in this section. 3.1. Solid/Fluid Interface: Channel Flow. Because of the importance of near-wall physics, in terms of understanding the transport phenomena, a large number of efforts were presented in the experimentation and conceptualization of near-wall turbulence. Experimental contributions of Fage and Townend,42 Einstein and Li,43 Popovich and Hummel,44 Kline et al.,6 and Meek and Baer45 are extensively referenced while developing the theories of heat and mass transfer. Fage and Townend42 used ultramicroscopy to study the motion of dust particles in the immediate vicinity of the wall on a turbulent fluid. The axial movements of these particles were jerky, and sometimes they almost came to rest. The jerkiness of the axial motion due to the large fluctuations in axial velocity was determined to be associated with the combined effect of very small changes in wall-normal velocity and the large velocity gradient du2/dx1 at the boundary. They noticed the existence of large lateral displacements of groups of these particles, and they mentioned the fact that their motion can be regarded as practically rectilinear only during the interval between such displacements. These observations were contradictory to the Prandtl’s concept of laminar sublayer with completely rectilinear motion of fluid elements without any lateral movements. Their work provided an experimental support for the renewal mechanism, which was later popularized by Higbie46 and Danckwerts.47 Popovich and Hummel44 used a flash photolysis method to study viscous sublayers in nondisturbing turbulent flow. They took photographs of the near-wall region, introducing a photosensitive fluid and focusing the light from a Xenon flash tube. The tracer was observed in two dimensions. They observed rectilinear motion in laminar sublayer for 38.1% of the measurements. An additional 7.5% of the photographs showed a viscous region without turbulent motion, but with continually changing velocity gradient, and the remaining 54.4% of the photographs indicated the presence of some turbulent motion or threedimensional distortion. The authors also investigated the penetration depth of eddies, because the rates of heat and mass transfer from the wall are greatly influenced by them. The authors found the closest distance to the wall reached by eddies was x1+ ) 2.09 ( 0.2. The observed average thickness of the laminar sublayer corresponded to x1+ ) 6.17, with a most probable thickness of x1+ ) 4.3. Their results indicated that, adjacent to the wall, there is a layer of very small thickness (x1+ ) 1.6 ( 0.4) in which a linear gradient occurs virtually at all times, but the slope of the gradient changes with time. Beyond x1+ ) 34.6, essentially turbulent flow exists. The major breakthrough, in terms of understanding the near-wall flow physics, emerged with the work of Kline et al.6 They used hot wire anemometry, dye injection, and

8252

Ind. Eng. Chem. Res., Vol. 48, No. 17, 2009

Table 1. Experimental Tools for Structure Characterization Sr. No.

Technique

1

marker image tracking methods: (a) smoke injection (b) dye injection (c) hydrogen bubbles and photochromic tracers and fluorescent particle tracing

Remarks Principle: The principle underlying each technique is the measurement of the simultaneous displacements marked fluid particles in consecutive images

Advantages: very good (1) for flow visualization and (2) in correlating probe-type turbulent burst detection techniques with the corresponding visualization data Limitations: (1) extracting accurate quantitative data is very difficult 2

laser Doppler velocimetry, LDV (Figure 9A)

3

hot film anemometry, HFA (Figure 9B)

4

particle image velocimetry, PIV (Figure 9C)

5

laser-induced fluorescence, LIF

6

stereoscopic particle image velocimetry, SPIV (Figure 9D)

7

defocusing particle image velocimetry, DPIV

Point Measurement Techniques Principle: LDV provides the instantaneous velocity components at the point in the flow field where two or more mutually perpendicular laser beams (viz. blue, green and cyan) intersect to form a fringe pattern. As the seeding or entrained particle passes through this fringe pattern, it scatters the incident laser beam and induces a shift in the frequency of scattered beam (known as Doppler shift). The Doppler shift depends on the fringe spacing and the velocity of the particle normal to the fringes. Advantages: (1) nonintrusive technique, (2) factory calibrated, and (3) micrometer-sized seedings are traceable Limitations: (a) unequispaced dataset, (b) Interpolation is needed for equispacing, and (c) high frequency values of energy are anomalous Principle: The HFA is based on convective heat transfer from a heated wire or film element placed in a fluid flow. Any change in the fluid flow condition that affects the heat transfer from the heated element will be detected virtually instantaneously by a constant temperature/ constant current HFA system. Therefore, HFA can be used to provide information related to for example, the velocity and temperature of the flow, concentration changes in gas mixtures, and phase change in multiphase flows. Advantages: (1) high data rate of ∼20 kHz, (2) equispaced data, and (3) evaluation of entire energy spectrum possible Limitations: (1) intrusive technique, (2) requires some prior knowledge to enable calibration, (3) HFA probe is very sensitive and very costly, (4) point datasets, and they cannot give us information on the spatial topology of flow structure, (5) voltage velocity calibration is critical, (6) limited response at high frequencies, (7) voltage stability is crucial, (8) variation in temperature is crucial, and (9) negative velocity cannot be measured and, therefore, the results of mean velocity cannot be reported Planar Measurement Techniques Principle: The PIV technique is based on the following steps: seeding the fluid flow volume under investigation with a few micrometer-sized particles, which are assumed to follow the fluid flow closely, illuminating a slice of the flow field with a pulsing light sheet; recording two images of the fluid flow with a short time interval between them, using a digital CCD camera; processing these two successive images to get the instantaneous velocity field. The entire image is then divided into interrogation areas. An intercorrelation technique is used to evaluate the most probable displacement of the seeding particles within each interrogation area. Advantages: (a) Planar measurements, (b) two and three component measurements possible, (c) evaluation of spatial derivatives is possible, (d) velocity measurement over a wide range is possible, (e) PIV gives much more space information on flow instabilities and about noncoherent and coherent turbulent structures. Limitations: (a) low data rate (∼7-10 Hz), (b) window size limits spatial resolution, (c) for strong velocity gradients, accuracy is reduced Principle: we seed the flow with fluorescent dye. Based on temperature and/or concentration, the intensity of the reflected light changes and local transient variations can be obtained. It is an add-on to PIV. Advantages: Quantification of scalar flux and cross-correlations Limitations: same as those with PIV Volumetric Measurement Techniques Principle: The stereoscopic PIV (SPIV), which is commonly considered a 3D extension of PIV, provides three components of the flow velocities confined in a thin slice of moving fluid medium. It employs the statistical average (correlation) of particles in two separate viewing angles before combining the averaged 2D vectors into 3D vectors, thereby losing information about individual particles. Even if one circumvents this problem by resorting to particle tracking instead of correlation, the information attainable with PIV is only limited in the laser sheet thickness. Advantages: shape estimation of flow structures Limitations: (a) qualitative information, (b) very low data rate, and (c) still applicability is limited to simple flows Principle: three-dimensionality is achieved through defocusing principle Advantages: Same as that for SPIV Limitations: qualitative information

8

holographic particle image velocimetry, HPIV

Principle: The HPIV records the 3D information of a large quantity of particles in a fluid volume on a hologram instantaneously and then reconstructs the particle images in a 3D space. From the reconstructed image field, the 3D positions (as well as size and shape information) of these particles can be retrieved. Furthermore, by finding the 3D displacements of the particles in the image volume between two exposures separated by a short time lapse, the instantaneous 3D velocities of these particles in the volume can also be obtained. Advantages: Same as those fot HPIV Limitations: (a) low laser light intensity, (b) low particle density, (c) low sampling rate, and (d) qualitative information

the hydrogen bubble technique for their study. These combined visual and quantitative techniques, which were employed by the authors, revealed many significant features of turbulent boundary layers. They observed that the laminar sublayer was not two-dimensional and steady, as simple

models had previously perceived. It contained threedimensional unsteady motions. Nakagawa and Nezu25 computed the third-order conditional probability distribution of the Reynolds stress obtained from HFA datasets to gain information on ejections, sweeps, and

Ind. Eng. Chem. Res., Vol. 48, No. 17, 2009

8253

Figure 9. Schematic of experimental techniques: (A) laser Doppler velocimetry (LDV), (B) hot film anemometry (HFA), (C) particle image velocimetry (PIV), and (D) stereoscopic PIV (SPIV). (Courtesy of Dantec Dynamics.) Table 2. Relative Performance of Different Experimental Tools Sr. No.

performance parameter

marker image tracking methods

HFA

PIV

accurately measured

accurately measured

accurately measured

gives good estimate

gives good estimate

gives good estimate

can be estimated using EIM method and dimensional analysis equispacing is required

can be estimated using structure function approach can be estimated using structure function approach can be estimated, but reliability reduces with increase in Re gives good estimate

can be estimated using structure function approach can be estimated using structure function approach can be estimated, but reliability reduces with increase in Re can be accurately estimated

LDV

1

mean velocity

cannot be estimated

2

turbulent kinetic energy energy dissipation rate

cannot be estimated cannot be estimated

obtained using energy spectrum

4

energy spectrum

cannot be estimated

3D energy spectrum gives correct values

5

structure time scale

can be estimated

gives good estimate

gives good estimate

6

shape and size of structures

gives good estimate

cannot be directly estimated

cannot be directly estimated

3

can be estimated, but without directional information gives good estimate

volumetric measurement techniques

interactions. Kreplin and Eckelmann48 computed the root-meansquare (rms) values, skewness and flatness factors, and probability density functions in a fully developed turbulent channel flow at Re ) 7700 using hot-film probes. HFA probe measurements showed that the bursts are associated with a high Reynolds shear stress. The wall layer is composed of elongated regions of high-speed and low-speed streamwise velocity. The large streamwise length of the streaks appears to be due to a sequence of vortices following each other, pumping high-speed fluid toward the wall and low-speed fluid away from the wall. Christensen and Adrian9 performed the PIV measurements at the outer region of turbulent channel flow (x1+ < 100) to determine the average flow field associated with spanwise vortical motions. They observed that the mean structure consists of a series of swirling motions ascending at an angle of 12°-13° from the wall. This is quite similar to the pattern followed by

hairpin vortices. The results proved that the instantaneous structures occur with sufficient frequency, strength, and order to leave an imprint on the mean statistics of the flow. Similarly, Tomkins and Adrian49 investigated the spanwise structure and the growth mechanisms in a turbulent boundary layer, by conducting the PIV measurements in the planar region between the buffer layer and top of the logarithmic region, at Ret ) 1015 and 7705. They observed the hairpin vortices at x1+ < 60 and at regions of local minima of streamwise velocity, and also the organized streamwise vortices as envisaged in the vortex parameter paradigm. The authors proposed that the additional scale growth occurs by the merging of vortex packets on an eddy-by-eddy basis via a vortex reconnection mechanism. Ganapathisubramani et al.50 performed a stereoscopic PIV measurements in streamwisespanwise and inclined cross-stream planes (inclined at angles of 45° and 135° to the principal flow direction) of a turbulent

8254

Ind. Eng. Chem. Res., Vol. 48, No. 17, 2009

boundary layer at moderate turbulent Reynolds number (Ret ) 1100). Their work corroborated the existence of hairpin vortices, similar to their predecessors. The two-point spatial velocity correlations (streamwise-streamwise (Ru2u2) and streamwise-wall-normal (Ru1u2) correlations) computed using PIV revealed higher correlation values for streamwise displacements of more than 1500 wall units, and even the zero-crossing data for the streamwise fluctuating component reveal that streamwise strips between zero crossings of 1500 wall units or longer, occur more frequently for negative u2 than positive u2, suggesting that long streamwise correlations in Ru2u2 are dominated by slower streamwise structures. 3.2. Free Shear Flows: Jet Loop Reactor. Sakakibara et al.51 made PIV measurements at the region of jet exiting the nozzle and at region where it impinged on a perpendicular wall. They used the isovorticity and iso-λ2 surfaces, coupled with critical point theory, to identify the flow structures. They observed that, near the nozzles, the spanwise rollers were evolving out of the shear layer. Furthermore, they observed the “cross ribs”, which extended from the downstream side of each roller to its counterpart across the symmetry plane. The cross ribs were seen to be stretching, because of the diverging flow as the rollers approached the wall and move apart, causing the vorticity within them to intensify. This process continues until the cross ribs reach the wall and merge with “wall ribs”. They observed that the wall ribs sustained themselves by merging and stretching rather than reorientation of the vorticity. Later, Matsuda and Sakakibara52 studied the turbulent vortical structures in a round free jet of water, using the stereo PIV. Using the isosurfaces of the swirling strength λi and the linear stochastic estimation, they observed a group of hairpin-like vortex structures around the rim of the shear layer of the jet. The center of curvature of the head of the hairpin was typically observed around x1/b ) 1.5, and the azimuthal spacing between the legs of the hairpin was roughly 0.9b. The typical spacing between the legs of the estimated hairpin was 0.65b, which is generally constant over the range of Re ) 1500-5000. Haven and Kurosaka53 used LIF and PIV methods to investigate the effect of an exit hole shape in a cross-flow jet on the flow pattern. The holes considered have elliptical, square, and rectangular shapes with the same cross-sectional area. The vorticity around the circumference of the jet was tracked to identify its relative contributions to the nascent streamwise vortices, which evolve eventually into kidney vortices downstream. The distinction between sidewall vorticity and that from the leading and trailing edges, though blurred for a round hole, became clear for a square or a rectangular hole. Deshpande28 performed the PIV and HFA of a jet loop reactor (JLR) (a cylindrical column with an internal diameter of 0.3 m and a height of 0.4 m). The experiments were conducted to understand the variation in mixing pattern and flow structure dynamics with the changes in the nozzle geometry. The jet region showed fluctuations of the jet plume due to the back pressure exerted by the bottom wall, which is the result of the jet impact on the bottom. This causes jet instability (JI) in the flow. The bulk region showed the motion of eddies 0.04 m in size along the circulation path, and they disappear after the travel length of ∼0.1-0.15 m. Similar observations were made near the nozzle where the structures from the bulk interact with the shear region and either disappear or revert back, based on the interaction with high speed flow. The size of these structures was determined to be in the range of 0.004-0.012 m. The comparison of two cases, with variation in nozzle diameter, showed that the small size structures play a more predominant

role as the diameter of the nozzle is reduced. This is clearly seen in the mixing time studies, wherein for a reduction in diameter of the jet, the mixing time performance becomes poorer at the same power per unit volume. 3.3. Flow Past Solid and Blunt Bodies: Stirred Tank. Wei and Smith54 used the hydrogen bubble flow visualization to detect the secondary vortices in the near wake of circular cylinders over a Re range of 1200-11000. Their results suggests that the free shear instability causes the separated cylinder boundary layer to roll up into the secondary vortices, and then these vortices undergo a strong three-dimensional distortion, which could provide the mechanism for the transition from laminar to turbulent Strouhal vortices. Lian and Huang55 used the hydrogen bubble technique to study the vortex behind bluff bodies with a sharp edge (flat plates, circular disks, and hollow hemispheres). Wu et al.56 measured turbulence intensities, autocorrelation functions, turbulence scales, energy spectra, integral scales, and turbulence energy dissipation rates in a baffled Rushton turbine agitated vessel with LDV. They found that 60% of the energy transmitted into the tank via impeller was dissipated in that region, and 40% was dissipated in the bulk of the tank. Glezer and Coles57 applied LDV to measure the turbulence intensity in the region of a thin cored vortex formed by a momentary jet discharge from an orifice in a submerged plate. Derksen and Van den Akker58 performed three-dimensional, angle-resolved LDV measurements of the turbulent flow field (Re ) 2.9 × 104) in the vicinity of a Rushton turbine in a baffled mixing tank. Results on the average flow field, as well as on the complete set of Reynolds stresses, are presented. The anisotropy of the turbulence has been characterized by the invariants of the anisotropy tensor. The trailing vortex structure, which is characteristic for the flow induced by a Rushton turbine, is demonstrated to be associated with strong, anisotropic turbulent activity. Hasal et al.59 obtained the LDV velocity field from a baffled pitch-blade turbine impeller and extracted the macroinstability (MI) from this data using the POD technique and spectral analysis. Water and an 85% aqueous glycerine solution were used as working fluids at three Re values of the impeller: 75000, 1200, and 750. The fluid velocity was recorded in a rectangular region of evenly spaced points close to the stirrer region. They also quantified the relative magnitude of MI (a ratio of the kinetic energy captured by the MI to the total kinetic energy of the flow). Kumaresan and Joshi60 conducted the spectral analysis of LDV dataset for various internals and impeller configurations for analyzing the effect of flow instabilities on mixing time. They suggested that the narrow blade hydrofoils generate MIs and JIs that have amplitudes that are very much smaller than that generated by pitched-blade turbines (PBTDs) and disk turbines (DTs). As a result of this, the hydrofoil gives the minimum mixing time, versus impellers that have more MI strength, because there are not many deviations in the dominant circulation pattern that is causing the mixing. Sheng et al.61 performed PIV measurements in a stirred vessel with PBTD45. They used a planar flow field to estimate the dissipation rate. Escudie´ and Line´62 performed PIV measurements for DT. They used the data for two purposes: to identify and quantify the transfer of kinetic energy between mean flow, periodic flow, and turbulence; and to estimate the dissipation rate of turbulent kinetic energy (TKE) from the balance of TKE. Mavros63 conducted a review of experimental techniques in stirred vessels. He reviewed techniques that were available for the study of the flow patterns induced by the various types of agitators (e.g., classical pressure or velocity measurements with

Ind. Eng. Chem. Res., Vol. 48, No. 17, 2009

Pitot tubes or hot-wire anemometers, and novel ones such as LDV, laser-induced fluorescence (LIF), and PIV). 3.4. Taylor-Couette Flow: Annular Centrifugal Extractor. The first set of experiments in Taylor-Couette flow systems were reported by Couette64 and Mallock.65 During their drag measurement experiments, they reported instability at certain rotational speeds of the cylinders. Coles66 visualized wavy vortices by suspending aluminum particles in silicon oil and reported interesting observations. The wavy nature of the Taylor vortex was found to be dependent on the way in which the rotational speed of the cylinders was varied. The numbers of vortices were also found to be dependent on the methodology of increasing or decreasing the speed. Smith and Townsend67 used HFA and Pitot tubes to measure the velocity and turbulence intensity in turbulent Taylor-Couette flow. They introduced a small amount of axial flow to push the toroidal vortices past the stationary probe. They concluded that the toroidal vortices lose their regularity at very high rotation rates and cannot be clearly distinguished beyond a Ta value of (5 × 105)Tacr. Koschmeider68 reported the wavelength of turbulent Taylor vortices up to 40000Tacr, using visualization experiments. Andereck et al.69 used flow visualization (laser light scattering) to construct flow maps that showed the various regimes from Couette flow to Turbulent Taylor vortex flow. The time dependences of the flows have been studied by measuring the intensity of laser light scattered by polymeric flakes. Later, Lueptow et al.70 presented flow maps in the presence of axial flow. They also measured the wavelengths of vortices for different flow regimes. For highly turbulent regimes, Parker and Merati71 used LDV to measure three components of mean velocity and turbulent intensity at various circumferential planes. For the aspect ratio of 4 and 20, they studied the end effects on the vortices. Wereley and Lueptow72 performed measurements of velocity fields with an imposed pressure-driven axial flow using PIV. For a radius ratio of 0.83 and an aspect ratio of 47, they determined the velocity vector field for nonwavy toroidal vortices and helical vortices. Racina and Kind73 measured the distribution of the local dissipation rate of turbulent kinetic energy in TaylorCouette flow with the help of PIV. They observed that the values of dissipation rate are strongly affected by the spatial resolution of PIV measurements. Fehrenbacher et al.74 used LDV to measure Reynolds stresses, integral and micro time scales, and power spectra over a wide range of turbulence intensities. Measured integral and micro time scales and approximated integral length scales were all observed to decrease with the Reynolds number, possibly associated with a confinement of the largest scales (of the order of the cylinder wall separation distance). Power spectra for the independent directions of velocity fluctuation exhibited slopes of -5/3, which suggests that the flow also has some additional isotropic characteristics and demonstrates the role of the Taylor-Couette apparatus as a novel means for generating turbulence. Campero and Vigil75 studied the hydrodynamic structures and studied the effects of various parameters, including density ratio, viscosity ratio, and feed composition in liquid-liquid TaylorCouette flow, incorporating a weak axial flow. At least three distinct structures were found, which included (1) a translating banded structure with alternating water and organic-rich vortices at low organic-phase volume fractions or sufficiently large rotation rates, (2) a spatially homogeneous emulsion with phase inversion at high organic-phase volume fraction and moderate rotation rates, and (3) an axially translating periodic variation between the banded and homogeneous states at low rotation

8255

rates. From the photographic evidence, they found that the banded flow pattern does not consist of separate aqueous and organic-rich vortices. Instead, the banded appearance is caused by disperse phase droplet migration to vortex cores. A series of experiments with various fluid pairs suggested that density and viscosity differences between fluid phases could not entirely explain the droplet migration to vortex cores. 3.5. Particle (Dispersed-Phase)-Induced Flow Structure. 3.5.1. Flow Past Spheres/Drops. The behavior of the wake behind the sphere at varying Re values has been studied by several researchers.76-86 The earlier flow visualization experiments have been performed by Taneda,76 using a string mounted sphere, which has been constantly moved in a water tank with the help of a motor. He measured the size, the separation angle, and the center of the steady axisymmetric wake behind the sphere and reported that the size of the vortex ring is proportional to the logarithm of Re. He found that the Reynolds number at which the axisymmetric toroidal vortex ring begins to form in the rear end of a sphere is Re ) 24. He also observed a faint periodic motion at the rare end of the vortex ring beginning at Re ) 130. The wakes generated by the liquid drops of carbon tetrachloride and chlorobenzene in water have been studied by Magarvey and Bishop,77 and they classified the wakes, based on the nature of the tail of the vortex and the Re value. Up to Re ) 210, the wake was steady and axisymmetric and was named as a single thread wake. For 210 < Re < 270, the vortex became nonaxisymmetric and has been classified as a double-threaded wake. In the range of 270 < Re < 290, the double-threaded wake becomes unstable and wavy nature of vortex tail has been reported. Above Re ) 290, vortex loops begin to release into the free stream. The formation and structure of vortices due to accelerating liquid drops at different intervals of time, at Re ) 340, has been shown experimentally by Magarvey and MacLatchy.78 They observed that the liquid drops follow a spiral path while settling. As noted by Winnikow and Chao79 and Natarajan and Acrivos,80 these experiments with falling liquid drops in immiscible liquids were compared with the standard solid rigid sphere wakes, because of the presence of the surface active impurities at the liquid-liquid interface, which hold the drops in a spherical shape. Masliyah81 has shown the recirculating wakes behind a sphere and three oblate spheroids, using flow visualization techniques for Re ) 15-100. He analyzed the variation of the wake length and angle of separation of the stable wake, with respect to Re. Achenbach87 studied the fixed sphere wakes for the range of 400 < Re < 5 × 106. With a help of a sketch, he explained the periodic formation and release of vortex loops in the free stream. He also determined the shedding frequencies at Re ) 3000 via the timely release of 50 vortices. The characteristics of the steady wake behind liquid-filled spheres have been studied experimentally by Nakamura,82 using dyed water for flow visualization experiments. From these flow visualization experiments, he observed that a stable and everlasting accumulation of dyed water at the rare end of the sphere begins at Re ) 7.3 and the shape of the wake changes from concave to convex as Re increases. He noted that the tracer, aluminum dust, that was used by Taneda76 is not fine enough to drift along with the slow fluid stream when Re < 30. He reported that the maximum Reynolds number at which the toroidal vortex is steady is ∼190, which is contrary to the observation of Magarvey and Bishop.77 This early instability in the wake at Re ) 190 can be attributed to the liquid-filled spheres, where the mass of fluid in these spherical shells was free to move around and potentially affect the sphere’s motion and the wake development. Sakamoto and

8256

Ind. Eng. Chem. Res., Vol. 48, No. 17, 2009

Haniu83 measured the vortex shedding frequencies of a fixed sphere for Re ) 300-40000, using hot wire anemometry and flow visualization experiments. They observed that the wavy instability at Re ) 130, before the onset of shedding, corresponds to the periodic motion observed by Taneda,76 at the same Re value, and the onset of the hairpin vortex shedding occurs at Re ) 300. The wake behind a fixed sphere from Re ) 30 to Re ) 4000 has been visualized using tracers illuminated by a laser light sheet by Wu and Faeth.84 They also took laser velocimetry measurements for the streamwise velocities. Ormieres and Provansal85 qualitatively showed that the doublethreaded wake of a sphere held by a thin metallic pipe was formed at Re ) 220 and periodic vortex shedding occurred at Re ) 300. They also made quantitative measurements of the free stream velocity in the wind tunnel, using LDV and HFA. Visualizations of the vortex structures and measurements of streamwise velocity of a fixed sphere in the Re range of 270-500 in a uniform flow channel have been made by Schouveiler and Provansal.88 Flow visualization experiments that capture the wake structure behind the rising and falling solid spheres in water have been shown quantitatively by Veldhuis et al.,89 using the Schlieren technique. From these experiments, they concluded that the wake generated by a moving sphere is different from the wake generated by a fixed sphere. They have shown a pair of opposite signed vortices threads, subsequently resulting into the formation of kinks onto these threads. This phenomena eventually results in the hairpin vortices. 3.5.2. Bubble Column. The PIV method has been used to understand bubble wake dynamics and bubble-induced flow structures, either for a single bubble train or in a bubble column setup. Initial work toward developing pulsed-laser image velocimetry-based digital data acquisition and analysis techniques, for measuring two component velocities of two-phase bubbly flow, was initiated by Hassan and Blanchat90 and Hassan et al.91 Joshi et al.36 reviewed the experimental observations on flow structures in the bubble column study. They have reviewed the works of Tzeng et al.,92 Reese and Fan,93 Chen et al.,94 and Lin et al.95 There have been several notable works since then. Lin et al.95 used the PIV system in a two-dimensional (2D) bubble column set up to quantify the two-phase flow conditions (four- and three-region flows) with coherent flow structures. The columns operated in the four-region flow condition comprise descending, vortical, fast bubble, and central plume regions. The fast bubble flow region moves in a wavelike manner and, hence, has been characterized macroscopically in terms of wave properties. They observed that, for columns larger than 0.2 m in width, the transition from the dispersed bubble flow regime to the four-region flows, and then to three-region flows in the coalesced bubble regime occur progressively with gas velocities at 1 and 3 cm/s, respectively. Tokuhiro et al.96 investigated the flow around an oscillating bubble and solid ellipsoid with a flat bottom. They measured both the flow field and the size of the bubble. The velocity of the flow field around the bubble was obtained using a CCD camera digital particle image velocimetry (DPIV) system enhanced by LIF. The shape of the bubble was simultaneously recorded along with the velocity using a second CCD camera and an infrared shadow technique (IST). They used the velocity vector plots of flow around and in the wake of a bubble/solid, supplemented by profiles and contours of the average and rms velocities, vorticity, Reynolds stress, and turbulent kinetic energy, to reveal the differences in the wake flow structure behind a bubble and a solid. They observed that the inherent, oscillatory motion of the bubble leads to vorticity and its subsequent stretching in the near-wake region. This vorticity

stretching was determined to be responsible for distributing the turbulent kinetic energy associated with this flow more uniformly on its wake, in contrast to the wake for a solid. Broo¨der and Sommerfeld97 studied a bubble column with a diameter of 140 mm and a height of 650 mm or 1400 mm (the initial water level). The gas holdup was varied in the range of 0.5%-19%. They used a two-phase pulsed-light velocimetry (PLV) system to evaluate instantaneous flow fields of both rising bubbles and the continuous phase. The measurement of the liquid velocities in the bubble swarm was achieved by adding fluorescing seed particles. Images of bubbles and fluorescing tracer particles were acquired by two CCD cameras. The optical interference filters with a bandwidth corresponding to the emitting wavelength of the fluorescing tracer particles and the wavelength of the applied Nd:YAG pulsed laser was used to separate the images from tracers and bubbles. To improve the phase separation of the system, the CCD cameras were additionally placed in a nonperpendicular arrangement, with respect to the light sheet. The acquired images were evaluated with the minimum-quadratic-difference algorithm. They acquired 1000 image pairs and recorded and evaluated data for each phase. They were able to compute the turbulence properties and characterize the bubble-induced turbulence for various bubble mean diameters and gas holdups, and they were able to determine the average bubble slip velocity within the bubble swarm. Fujiwara et al.98 explored the flow structure in the vicinity of the single bubble (dB ≈ 2-6 mm) in one plane and its deformation in two planes by PIV/LIF and a projection technique for two perpendicular planes, respectively. The second and third CCD cameras were used to detect the bubble’s shape and motion via backlighting from an array of infrared lightemitting diodes (LEDs). They studied the three-dimensional (3D) wake structure from measurements of the 2D vortex structure and approximated 3D shape deformation arranged from two perpendicular bubble images. Liu et al.99 investigated the flow structure induced by a chain of gas bubbles in a rectangular bubble column using PIV. It is observed that the bubble rising trajectory changes from one dimension to three dimensions as the liquid viscosity is reduced. Furthermore, they observed a free vortex, cross-flow, and irregular circular flow in the fluid flow. Sathe et al.100 used the shadowgraphy technique and the PIV/ LIF measurements in combination, along with the fluorescent tracer particles and the blue filter. This was done with the purpose of obtaining the shape, size, velocity, and acceleration of gas bubbles, along with the flow structures and liquid velocity profiles. The measurement was performed at a high local gas holdup (∼10%) with a wide variation of bubble sizes (0.1-15 mm). The liquid velocity field was decomposed into seven scales, using wavelet transform. Scales 5-7 were added to generate the vector field showing small scale structures, while scales 1-4 were added to generate the vector field showing the large-scale structures. This separation was observed to be helpful in determining the slip velocities of bubbles in terms of “local” liquid velocity. 3.6. Recent Advances in Experimental Tools (Volumetric Measurement). Several extensions of the classical PIV system have been proposed to obtain volumetric 3D, threecomponent velocity information. The most common extension is the application of the second CCD camera to acquire the stereoscopic view of the flow and, thus, achieve the out-ofplane component of the velocity on a plane (see Raffel et al.101). Some other well-known extensions of PIV toward three

Ind. Eng. Chem. Res., Vol. 48, No. 17, 2009

dimensionality are multiple plane stereoscopic PIV (a combination of two SPIV systems, i.e., two double-pulsed lasers, four cameras, polarized light, and two planes, as proposed by Kahler102), a defocusing PIV (three dimensionality is achieved through a defocusing principle, according to Willert and Gharib103), and the HPIV (volume illumination and holographic recording procedure). Many different variations of HPIV exist, such as light-in-flight, on-axis, off-axis, hybrid recording and reconstruction, etc. The stereoscopic PIV (SPIV), which is commonly considered a 3D extension of PIV, provides three components of the flow velocities confined in a thin slice of moving fluid medium (from Arroyo and Greated,104 Prasad and Adrian,105 and Lecerf et al.106). It uses the statistical average (correlation) of particles in two separate viewing angles before combining the averaged 2D vectors into 3D vectors, thereby losing information about individual particles. Even if one circumvents this problem by resorting to particle tracking instead of correlation, the information attainable with PIV is only limited in the laser sheet thickness. The HPIV method provides a solution to this limitation. HPIV records the 3D information of a large quantity of particles in a fluid volume on a hologram instantaneously and then reconstructs the particle images in a 3D space. Recent works by Tao et al.107,108 have mostly focused on using the HPIV dataset of a square duct facility to study the structure of the filtered, threedimensional vorticity, strain-rate, sub-grid-scale (SGS) stress tensor distributions and 3D flow structures. They observed that the filtered high-vorticity isosurfaces obtained from HPIV showed structures that are only slightly elongated, as opposed to the long and thin “worms” observed in DNS for unfiltered turbulence at much lower Re values. Furthermore, they calculated a 3D geometric relationship between filtered vorticity, strain rate, and subgrid-scale stress tensors at high Re values. They are working on using such relationships to understand the fundamental complexities of turbulence generally, and for the development of turbulence models for large eddy simulations in particular. Along the same lines, Van der Bos et al.109 used HPIV to make better LES models. He studied the effects of small-scale motions on the inertial range structure of turbulence by considering the dynamics of the velocity gradient tensor filtered at scale ∆. An attempt was made to optimize the mixed model. Svizher and Cohen110 recently used an HPIV system to study the evolution of coherent structures artificially generated in a plane Poiseuille air flow. Initially, they used the hot-wire technique and two-dimensional flow visualization technique (PIV) to determine the generation conditions and dimensions of the coherent structures, their shedding frequency, trajectory, and convection velocity. The HPIV method then was used to obtain the instantaneous topology of the hairpin vortex and its associated 3D distribution of the two (streamwise and spanwise) velocity components, as well as the corresponding wall-normal vorticity. It was observed that the generation of hairpins under various base flow conditions was governed by the shear of the base flow and an initial disturbance with large amplitude. 4. Computational Tools Computational fluid dynamics (CFD) is widely used to understand the flow structures and their effect on the turbulent transport of mass, momentum, and heat in chemical process equipment. The application of CFD in chemical process equipment such as mixers, chemical reactors, packed beds, crystallizing/dissolving equipment, pneumatic conveyors and classifiers, flows in pipes, sprays, membrane separation, dust

8257

Figure 10. Figure explaining the RANS-based model, as well as LES and DNS computer resolution on the energy spectrum.

separation, pyrolysis (cracking), and biological systems is increasing at a faster rate. Many success stories can be cited regarding the shortened product-process development cycles, optimization and control of existing processes to improve yield and energy efficiency, efficient design of new processes, novel equipments, and improvements in health, safety, and environment. The efforts are directed toward increasing the reliability of CFD to replace the experiments (see “Technology Roadmap for CFD”111 and the work of Bakker112,113). In addition to velocity and pressure data, CFD provides information on quantifiers of flow structures such as vorticity, second invariant gradients over the entire volume of the reactor, which would otherwise be obtained using multiple advanced measurement techniques. In this work, the role of computational flow modeling in understanding the flow structures and its impact on designing of chemical process equipment has been discussed. The present section highlights the chronological development of different turbulence models, vis-a`-vis an improved understanding of flow structures with increasingly complex hydrodynamics. The application of these computational tools and the flow patterns obtained from various equipments then are reviewed. The governing equation of motion for a Newtonian incompressible fluid can be written as ∂ui ∂F +F )0 ∂t ∂xi F

(1)

( )

∂uiuj ∂ui ∂ui ∂p ∂ +F )+ µ ∂t ∂xj ∂xi ∂xj ∂xj

(2)

Depending on the boundary conditions and the strong nonlinearity of the convective term, different solutions are obtained. The accuracy of predictions is dependent on the numerical resolution of all the turbulent scales present in the flow. Different turbulence models resolve turbulent scales up to different degrees. Figure 10 explains the scales that are resolved and modeled by the DNS, LES, and the Reynoldsaveraged Navier-Stokes (RANS)-based models. These models are explained below, along with their performance in capturing the flow structures in different equipment. 4.1. Direct Numerical Simulation (DNS). It has been already discussed earlier that turbulence contains a wide range of length and time scales, and that the large-scale structures extract the energy from the mean flow and this energy is

8258

Ind. Eng. Chem. Res., Vol. 48, No. 17, 2009

transferred down the cascade toward the smaller scales, until it reaches the smallest Kolmogorov scale, where viscous dissipation directly converts the energy into heat. Now, it is theoretically possible to resolve the entire spectrum of turbulent scales directly and solve the Navier strokes equation (see eqs 1 and 2) at a grid with a size less than that of the Kolmogorov scale. This approach is called as DNS, wherein only the instantaneous Navier-Stokes equations are solved and no turbulence models are used. However, this approach is limited by the speed of computation, because, to resolve all the scales in three dimensions, the number of mesh elements required is proportional to 114 3 Re9/2 ). Hence, the limitation of t Pr (from Mathpati et al. computational resources led to the development of turbulence models such as RANS equations and large-eddy simulation, which enabled the simulations to be performed at a coarser grid size with reasonable accuracy concomitant with the assumptions. 4.2. Reynolds-Averaged Navier-Stokes (RANS)-Based Models. The RANS approach involves solving the timeaveraged Navier-Stokes governing equation to obtain the mean profiles. The instantaneous components of velocities and pressure in the Navier-Stokes equation are split into its mean and fluctuating parts, and then the equation is time-averaged. This approach reduces the computational demand. These equations, with time averaging, are reduced to the following form (eqs 3 and 4): ∂〈ui〉 ∂F +F )0 ∂t ∂xi F

(3)

(

)

∂ 〈ui〉 ∂ 〈ui〉 ∂ 〈pi〉 ∂ 〈ui〉 ∂ + F〈uj〉 )+ µ - F〈ui′uj′〉 (4) ∂t ∂xj ∂xi ∂xj ∂xj

Time averaging of the equation of motion produces an additional term ( τij) that is given by τij ) -F〈ui′uj′〉. These stresses are modeled using gradient diffusion hypothesis to close the system of equations.

(

∂〈uj〉 ∂ 〈ui〉 2 -F 〈ui′uj′〉 ) kFδig - µt + 3 ∂xj ∂xi

)

(5)

To close the system of equations, µt has been formulated in many ways in the literature, in terms of zero-, one-, and twoequation models. The most popular way is the two-equation model, where µt is estimated from turbulent length and time scales. In the case of the standard k-ε model, µt is estimated using turbulent kinetic energy (k) and rate of dissipation of turbulent kinetic energy (ε). µt ) CµF

() k2 ε

(6)

The values of k and ε are obtained by solving their transport equations (see Launder and Spalding115). 4.2.1. Standard k-ε Model. The equations for turbulent kinetic energy and dissipation rate are given below: F

[(

) ]

µt ∂k ∂k ∂k ∂ + F〈uk〉 ) P - Fε + µ+ ∂t ∂xk ∂xk σk ∂xk

(7)

The transport equation for ε is rewritten as F

()

[(

) ]

µt ∂ε ∂ε ε2 ε ∂ε ∂ + F〈uj〉 ) FCε1 P - Cε2F + µ+ ∂t ∂xj k k ∂xj σε ∂xj (8)

where

()

P ) τij

∂〈ui〉 ∂xj

(9)

The aforementioned modeling of k and ε equations result into five turbulence parameters: Cµ, Cε1, Cε2, σk, and σε. This model has been successfully applied in the literature for flow away from walls and regions without sharp velocity gradients. As the wall is approached, the results deviate from the experimental values. To model the effect of contribution of viscous stresses to turbulent stresses, wall functions and low Reynolds number models have been developed in the literature. These are semiempirical formulas and functions that link the solution variables at the near-wall cells and the corresponding quantities on the wall. A critical review of low Reynolds number models is presented by Thakre and Joshi.116 4.2.2. Reynolds Stress Model (RSM). The standard k-ε model inherently fails to predict the return to isotropy after the removal of strain or the isotropization of grid-generated turbulence properly.117 In the k-ε model, the production of turbulence is modeled, whereas accurate prediction of turbulence production can improve the overall flow patterns in anisotropic flows. Furthermore, the modeling of transport equations for k and ε clearly gives a glimpse of its inability to properly account for streamline curvature, rotational strains, and other body-force effects. RSM, in theory, will circumvent all the aforementioned deficiencies and also gives it the ability to predict each individual stress more accurately. In RSM, six transport equations for Reynolds stresses are solved to take into account anisotropic behavior. The equation of turbulent stresses contains the following terms: pressure rate of strain (which contains fluctuating pressure velocity gradients) and the flux of Reynolds stresses.

(

) ( ) 〈 〉

∂τij ∂〈uj〉 ∂〈ui〉 ∂τij ∂ µt ∂τij + 〈uk〉 ) - τik + τjk + ∂t ∂xk ∂xk ∂xk ∂xk σk ∂xk ∂ui′ ∂uj′ + µ∇2τij (10) Πij - 2µ ∂xk ∂xk

〈[

Πij ) p′

∂ui′ ∂uj′ + ∂xi ∂xj

]〉

) Πij,slow + Πij,rapid + Πij,wall (11)

The pressure strain term is the most uncertain term in RSM. This term is responsible for making turbulence isotropic and redistribution of energy between components 〈u′21 〉, 〈u′22 〉, and 〈u′23 〉. The incompressibility condition guarantees that Π11 + Π22 + Π33 ) 0. If individual terms of Πii are nonzero, at least one must be positive and one must be negative. The classical approach is to decompose the pressure strain term into a slow component, a rapid component, and a wall reflection term.118 The RANS-based models (k-ε and RSM), because of time averaging, cannot resolve the flow structure, and they suffer from the following shortcomings: (a) the basic gradient diffusion assumption is questionable; (b) the eddy viscosity has isotropic character; (c) there is inadequate incorporation of viscosity damping effects on the turbulence structure (low-Re models); and (d) there is an inability to mimic the preferentially oriented and geometry dependent effects of pressure reflection and eddyflattening and squeezing mechanisms, because of the proximity of the solid or interphase surface. 4.3. Large Eddy Simulation (LES). In most of the chemical process equipments, large flow structures govern the transport phenomena. If one could manage to resolve the large flow structures and accurately model the small isotropic structures, then the computational demand can be significantly decreases, compared to the DNS. In large eddy simulations (LESs), the

Ind. Eng. Chem. Res., Vol. 48, No. 17, 2009 119

instantaneous velocities are filtered in the spatial domain. The filtering induces a closure problem in which the dynamic effects of the small-scale turbulence on the primary flow features must be represented by introducing a so-called subgrid scale (SGS) model. In LES, large scales are resolved numerically and small subgrid scales are modeled. The final results are dependent on the accuracy of SGS models in transferring energy from large scales at cutoff to lower scales. Much of the LES research is directed toward formulating accurate subgrid closures, evaluating the numerical errors, and identifying and characterizing large-scale structures. The governing equations for LESs are as follows: ∂uji ∂F +F )0 ∂t ∂xi

(

(12)

)

∂uji ∂ ∂ ∂pj ∂ F (uji) + F (uiuj) ) + µ - Fui,sgsuj,sgs ∂t ∂xj ∂xi ∂xj ∂xj

(13)

The major difference between the RANS equations and the aforementioned equation is that the dependent variables are filtered quantities rather than the mean quantities. The SGS stresses that result from the filtering operation are unknown and require modeling. The LES can resolve the large scales, as a result, it can be useful for simulating the flows undergoing transition and wake region, wherein large unsteady turbulent structures are resolved. The passage of flow structures are associated with high Reynolds stress, as a result of strong rotational motions and acceleration/ retardation motions that induces the gradients in flow. 4.4. Application of Turbulence Models for Flow Structure Identification. 4.4.1. Solid/Fluid Interface: Channel Flow. In this section, a discussion on those flows where a majority of turbulence production occurs at a stationary wall is presented. Solid/fluid interfaces offer very high shear rates, because of the flow structures that form near the solid wall, and a maxima of turbulent kinetic energy production and dissipation is observed.120,122 Because of an increased contribution of viscous forces close to the wall, standard RANS models fail to reproduce individual stresses, and even kinetic energy. To understand the transient nature of wall dynamics, simulations must be performed using DNS and LES. The first DNS of wall-bounded flows was performed by Kim et al.123 Subsequent studies have modified the channel configurations to examine the response of wall-bounded turbulence to factors such as rotation,124 mean three dimensionality,125 transverse curvature,126 and heat transfer.127 Some of the insights into the boundary layer structure gained from DNS are summarized below (also see the work of Moin and Spalart128 and Robinson129). The DNS data suggests that the observed bursting event may simply be due to the passage of streamwise vortices past the measuring station. Because of very high computational demand, the DNS could not be used at high Re values. In low-Re computations, nearwall streaks and horseshoe vortices have been observed. However, for better understanding, high-Re simulations are indispensable to remove the empiricism involved in the assumption of similarity between low-Re and high-Re near-wall structures. Many times, the DNS were performed at lower Re values, with the objective being to perform controlled studies that allow better insight, scaling laws, and turbulence models to be developed. To study the flow structures at higher Re values, LES has been used. LES of wall-bounded flows has been conducted for different objectives, such as to understand turbulence,125 the development of SGS models,130-136 the effect

137

8259

of filter shape, and wall treatment when coarse mesh is used,131,138 to study the scalar transport.139 Dong and Lee139 investigated passive heat transfer in turbulent channel flow at Re ) 10000 and Pr ) 0.1-200. They validated their LES results for Pr ) 1 and 10 with the DNS data of Na et al.140 They found that the turbulent Prandtl number to be almost constant, independent of the molecular Pr within the range of study. Similar to the mean velocity profile, they observed a buffer layer, followed by a logarithmic region in the mean temperature profile. For their range of study, they confirmed a Pr1/3 dependency of the Nusselt number (Nu). The contours of instantaneous temperature fluctuations clearly indicated that almost all of the turbulent structures that exist in the viscous sublayer were efficient in the heat-transfer process at Pr ) 1, and only the smallest structures subsisting very close to the wall were involved in the heat-transfer process at Pr ) 100. These simulations accounted for the effect of flow structures on heat transfer indirectly. The instantaneous variation of temperature was averaged to get an accurate estimation of the mean temperature gradient and, hence, heat-transfer coefficient. This transfer coefficient then was related empirically to the rate of surface renewal. 4.4.2. Free Shear Flows: Jet Loop Reactor. Jet loop reactors (JLRs) are frequently used in the process industry, as an alternative to impeller mixers, because of the existence of relatively high convection currents, leading to an enhanced rate of macro-mixing at the same power consumption, compared to stirred tank reactors. The resulting high-speed jet entrains some of the surrounding liquid and creates a circulation pattern within the vessel. CFD has been used on numerous occasions to simulate the flow patterns in the jet flows. A review on application of the k-ε model suggests that various researchers were unable to get accurate results for the free jets; hence, they had to introduce modifications in the model constants of the standard k-ε model to get accuracy.141-146 The LES of jet flow has been performed in the literature with the following objective: to understand the initial development region, vortex interaction in the jet region, the entrainment rate, and, subsequently, the mixing performance. Andersson et al.147 studied the effect of inflow conditions and SGS model on LES of mixing processes in turbulent jets. Grinstein and DeVore148 observed, through LESs, that the initial development of the square jet is characterized by the dynamics of vortex rings and braid vortices. Further downstream, strong vortex interactions lead to the breakdown of the vortices, and to a more disorganized flow regime characterized by smaller scale elongated vortices. Entrainment rates significantly larger than those for round jets are directly related to the enhanced fluid and momentum transport between jet and surroundings determined by the vortex dynamics underlying the axis-rotation of the jet cross-section. The first axis-rotation of the jet cross-section can be directly correlated with self-induced vortex-ring deformation. Mathpati et al.149 performed JLR simulation using a hybrid k-ε model and LES to assess the modeling assumption in RANSbased models. They also simulated scalar transport to study the effect of nozzle diameter on mixing time. They found that, for very small nozzle diameters, JLRs were inferior to stirred tank reactors for mixing applications. 4.4.3. Flow Past Solid and Blunt Bodies: Stirred Tank. The performance of the k-ε model in stirred vessels has been reviewed extensively.150-154 The reviews suggest that the k-ε model has been able to capture the gross circulation pattern quite well, but not the macro-instability and rolloff of vortices behind the impellers and other transient flow structures. The RSM

8260

Ind. Eng. Chem. Res., Vol. 48, No. 17, 2009

model has been used by Bakker and Van den Akker155 and Oshinowo et al.156 They observed that the RSM model has better prediction ability than the k-ε model. However, overall, the view that emerged from reviews was that the LES should be able to capture the instantaneous flow profile. Bartels et al.157 conducted DNS of the flow around a sixbladed Rushton turbine with four baffles for Reynolds numbers ranging from Re ) 0.1 to Re ) 106, and compared the results with those of the k-ε model. They observed that the flow field in a plane at an angle of 15° behind the blades is characterized by two “tip” vortices that are created in a plane ∼5° behind the blades. The location of the vortices in the DNS results was further away from the axis, compared to the RANS results (especially for the lower vortex). The two tip vortices were seen to be dying out at an angle of ∼35° behind the blade. This result is identically produced by the DNS and RANS simulations. The very first study of LES for a stirred vessel was performed by Eggels,158 using a lattice-Boltzmann discretization scheme for filtered Navier-Stokes equations and the adaptive force field technique for the impeller rotation. This study presented the instantaneous and mean flow field, which would help in the initial understanding of the given system. However, there were discrepancies in the quantitative predictions of the mean radial and mean axial velocities, especially in the near-impeller region. Revstedt et al.159 investigated both the mean flow and spectral distribution of flow by power density function at various locations in a baffled stirred tank. In this study, the impeller motion was modeled by specifying the time-dependent momentum source term in a stationary grid. The LES could capture the trailing vortices behind the blades and their detailed movement away from the impeller with increasing distance from the blade. The acceleration of the fluid at the impeller tip, because of the fluid entrainment due to the trailing vortex pair, was well-simulated. Their instantaneous data from LES supported the blade frequency and the existence of Kolmogorov’s -5/3 region in the energy spectrum. These observations were in close agreement with the expensive and time-consuming experiment. They demonstrated that the LES could generate an instantaneous flow field, such as sophisticated flow-measuring techniques and sometimes even better, as far as the flow in the impeller region and near-wall regions are concerned. The mean axial velocity from LES was in good agreement with the experimental data of Costes and Couderc161 at different locations. Derksen and Van den Akker55 advanced the Eggels et al.158 study with a more-refined forcing algorithm. The real potential of LES was explored as the simulations could identify a detailed local flow in the wakes behind each blade, which were determined to be significantly higher than the impeller tip speed (2 × tip speed). Furthermore, the fluctuations become less coherent away from the impeller. The axial profiles of random (turbulent), coherent (pseudo-turbulence), and total kinetic energy are in good agreement with LDV phase-resolved experimental data. Derksen160 applied the previously developed lattice-Boltzmann-based solver for the flow generated by a 4-PBTD45 impeller. He focused on the influence of the spatial resolution and SGS modeling on the flow predictions. It was observed that the spatial resolution has significant impact on the overall average flow field results. No significant effect was observed in the overall and phase-averaged flow field results, because of changes in the SGS model. Only the formation and strength of the tip vortex differed for the two SGS approaches (the vortex formed in the structure function simulations was observed to be weaker). Validation was made only within the impeller region. For the PBTD45 impeller, Roussinova et al.162

performed LES simulations using sliding mesh methodology (SDM). The mean axial velocity in the impeller center plane was in agreement with the experimentally measured values. This extended for the identification of low-frequency macro-instabilities (MIs). Their predicted frequency value of the MI compared well with the LDV measurement. Furthermore, Hartmann and co-workers163,164 used LES to investigate the precessional vortex phenomenon with the help of LES. This study opted for the same numerical techniques as those used by Derksen and Van den Akker.58 Their LES model could capture the vortical structure moving around the tank centerline in the same direction as the impeller. Using LES, they also observed that the strength of the vortex below the impeller was much stronger and more pronounced in size, compared to that prevailing above the impeller, and both vortices move with a mutual phase difference, which was qualitatively in good agreement with the experimental data. It was found that the fluctuations in the impeller region were dominated by the blade passage frequency, whereas the fluctuation levels close to the tank centerline were dominated by the low-frequency motion of a vertical structure. Their LES model could accurately predict the characteristic MI frequency. It could even predict a second frequency peak at f ) 0.092N at Re ) 12500. Using LES data, they observed that the flow was dominated by the low-frequency processional vortex in the bulk flow region, and these were associated with a significant amount of kinetic energy. For the first time, Yeoh et al.165 performed LESs using the sliding deforming technique. The main objective was to assess the capability of the LES modeling technique, coupled with the sliding SDM and relative performance, with the RANS approach. RANS calculations have been performed using the standard k-ε model. The predicted mean radial, tangential, and axial velocities, and the turbulent kinetic energy, were compared with LDV data only in the impeller center plane. For all three velocity components, fairly good qualitative and quantitative agreement was obtained using both the RANS-based and LES models. This work revealed that the equilibrium and isotropy hypotheses hold in most of the tank, except for the impeller discharge stream. Hartmann and co-workers163,164 also made a detailed study to evaluate the predictive capabilities of LES and RANS. The main difference was that this study used the lattice-Boltzmann solver and the RANS-based simulations have been performed with the shear-stress-transport (SST) turbulence closure model. The RANS-based and LES predictions have been compared and assessed with LDV-measured radial, tangential, and turbulent kinetic energy in the impeller zone. The radial velocities of both the models were in better agreement with the experimental data. With regard to the predictions of tangential velocity, the overall performance of LES was quite good, rather than RANS simulations. Their simulations revealed that the Smagorinsky SGS model performed better than the Voke SGS model, especially when it comes to prediction of the trailing vortex pair. It was shown that LES predicted the turbulent kinetic energy better in the impeller discharge flow, whereas RANS invariably underperformed in this area. Using LES data, they found almost-isotropic behavior in the circulation loops. However, in the impeller stream, the boundary layers, and at the separation points, turbulence was determined to be more anisotropic, because of high shear rate. They did not see the relative merits of using the Voke SGS model. Furthermore, they calculated the energy dissipation rate by assuming local equilibrium between production and dissipation at and below the SGS level and compared it with the RANS predictions in the baffle midplane. It was observed that the RANS-based model yielded higher values of the energy dissipation rate than LES

Ind. Eng. Chem. Res., Vol. 48, No. 17, 2009

8261

Figure 11. Mathematical tools for structure identification and characterization: (A) eddy isolation method and (B) decomposition of the velocity signal using the EIM, DWT, and CWT techniques.

and very inhomogeneous distribution throughout the tank. Again, Yeoh et al.166 extended their previous work to study the mixing of an inert tracer using LES. This was the first attempt to study the detailed mixing phenomenon using LES. The study clearly shows that LES provides a very detailed evolution pattern of the concentration field in space and time. It was at the impeller midsection where a substantial amount of scalar has been trapped in a region between the two vertical baffles. In the impeller region, they observed instantaneous concentrations to be as high as three times the equilibrium value. They emphasized that LES can be used to identify stagnant zones, mixing inhomogeneities due to the vessel geometry, and the optimal location of the feed pipe. Furthermore, Alcamo et al.167 presented the predictive capabilities of LES for flow field generated in an unbaffled cylindrical vessel. The LES could capture the existence of a trailing vortice pair in unbaffled tanks. Thus, the LES has been able to reveal some qualitative complex flow features in the flow generated by DT and PBTD45 impellers that have been captured by LES-like sophisticated flow-measuring devices (LDV, PIV, HFA). Figures 16 and 17 show some of the flow structures captured in the stirred tank region. 4.4.4. Taylor-Couette Flow: Annular Centrifugal Extractor. To predict the magnitude of vortex velocities and behavior beyond the critical point, Baier163 used a CFD

approach. Haut et al.168 also performed CFD simulations in the annular region in horizontal co-rotating cylinders, for the wavy vortex flow and the turbulent Taylor vortex flow. They used the standard k-ε model to account for the turbulence parameters. The velocity profiles obtained from the simulations for both regimes were in good agreement with their experimental results. The maximum rotational speed they used was 5 rad/s (Ta ) 5.35 × 105). They have investigated instantaneous and timeaveraged velocity fields experimentally,y using PIV. Hwang and Yang169 performed numerical simulations to study the flow fields and bifurcations related to the Taylor-Couette flow with an imposed axial flow. They found the computed speed of vortex movement due to axial flow to be consistent with the experimental observations of Wereley and Lueptow.72 They found that the torque on the inner cylinder surface is reduced significantly, because of the axial flow, which was consistent with previous established experiments. They furthermore found that the “shift and reflect” symmetry (a symmetry rule that shows identical vortices upon shifting by a period) is applicable for flows without an imposed axial flow. Deshmukh et al.33 performed RSM simulations in an annular centrifugal extractor. They found good agreement with PIV data for the mean components of velocity and the turbulent kinetic energy. They identified wavelengths of the vortices in various regimes over a wide range of speeds and annular gaps.

8262

Ind. Eng. Chem. Res., Vol. 48, No. 17, 2009

Figure 12. Channel flow structures as captured by various methods: (A) DWT scales (Re ) 5600), (B) DWT scales (Re ) 18000), (C) POD modes (Re ) 5600), and (D) hybrid wavelet-POD scales.

Deshmukh et al.170 conducted RTD simulations in an annular centrifugal extractor. They studied the effect of flow ratio, annular gap, and rotational speed on RTD. They observed that the flow in the ACE was similar to back-mixed behavior, because of the presence of counter-rotating vortices. The number of vortices is dependent on the rotational speed and the geometrical parameters of ACE. 4.4.5. Particle (Dispersed-Phase)-Induced Flow Structure. 4.4.5.1. Flow Past Spheres/Drops. A very scant amount of literature is available regarding the direct numerical simulation (DNS) of particles and drops at relatively low Re for a finite number of particles. The major focus was on the identification of the point of bifurcation and the drag experienced by the particles. Kim and Pearlstein171 performed a two-dimensional linear stability analysis by DNS, using a pseudo-spectral method. They predicted that the primary bifurcation point occurs at Re ) 175, which was comparatively lower than the experimental observations of Magarvey and Bishop.77 With the modified linear stability analysis, Natarajan and Acrivos80 performed two-dimensional simula-

tions of spheres and circular disks. They predicted that the primary bifurcation of sphere occurs at Re ) 210, which was in perfect agreement with the experimental observation, and they reported that the secondary bifurcation and shedding of vortices occurs at Re ) 278. Tomboulides172 have performed DNS of a fixed sphere from Re ) 20 to Re ) 1000. They observed initial flow separation at Re ) 20 and steady nonaxisymmetric flow at Re ) 212. The vorticity of the nonaxisymmetric flow field resemble the double-thread wake as observed by Magarvey and Bishop77 at Re ) 210. Both numerical and experimental analyses of the wake structure behind a fixed sphere in the Re range of 20-300 were also studied by Johnson and Patel.173 They reported that the flow separation and formation of a vortex behind the sphere begins at Re ) 20, and they predicted the separation angle, the separation length, and the position of the vortex. They observed that the axisymmetry was broken at Re ) 210 and explained this phenomenon by presenting three-dimensional interactions of stream lines. They reported the variations in the drag and lift forces, which arise after a breakage of symmetry. They reported that the periodic vortex shedding

Ind. Eng. Chem. Res., Vol. 48, No. 17, 2009

8263

Figure 13. JLR structures as captured by various methods: (A) DWT scales (Re ) 17000), (B) DWT scales (Re ) 26000), (C) POD modes (Re ) 17000), and (D) hybrid wavelet-POD scales.

begins at Re ) 270 and explained the same with the help of vorticity diagrams. DNS simulations of the flow past fixed oblate spheroidal bubbles at Re ) 100-1000 have been performed by Magnaudet and Mougin.174 At Re ) 180, they showed a nonaxisymmetric wake, with the help of 3D particle path. They also presented the vorticity isosurfaces at Re ) 180, 300, and 700. The wake instability of a fixed sphere has been studied by Ghidersa and Dusek,175 using the spectral element method, and they reported that the primary bifurcation occur at Re ) 215. Vortex methods for DNS of flow past spheres have been developed by Ploumhans et al.176 They presented the isosurfaces of streamwise vorticity at Re ) 300, 500, and 1000. Hydrodynamic forces

acting on the rigid fixed sphere for the Re range of 10-320 have been simulated by Bouchet et al.,177 using the spectral element method. They presented the vortex shedding frequencies and oscillation amplitude of drag and lift forces from Re ) 280 to Re ) 320. DNS of the wake generated by freely falling particles, by Reddy et al.,178 has shown that, at Re ) 1, the drag due to the pressure was 33% of the total drag and remainder was viscous drag. As flow separation began (at Re ) 25), the pressure contribution in the total drag increased and was equal to the viscous drag (CDP) at Re ≈ 180. They observed that, after Re > 200, the pressure drag dominates over the viscous drag. At Re ) 25, they observed a small and steady wake at the rear end of

8264

Ind. Eng. Chem. Res., Vol. 48, No. 17, 2009

Figure 14. Stirred tank structures as captured by various methods for disk turbine, PBTD45, and HF impellers: (A) DWT scales (DT), (B) DWT scales (PBTD45), (C) DWT scales (HF), (D) POD modes (DT), (E) POD modes (PBTD45), and (F) POD modes (HF).

the sphere. As the Re value increases from 25 to 200, the steady wake grew in size and stretched in the streamwise direction. At Re ) 200, it can be observed that the pressure distribution around the sphere is symmetric. Whereas at Re ) 210, the asymmetry in the pressure distribution can be clearly observed. The reason for this nonaxisymmetry was attributed to the increase in the centrifugal instability in the wake as the sphere falls freely, under the action of gravity. As the instantaneous Reynolds number increases (velocity of the sphere increases with time up to its terminal settling velocity), the rotational or centrifugal instability in the form of pressure gradient increases, because of the fluid rotation within the toroidal vortex. For steady axisymmetric flow, the pressure gradient inside the vortex was balanced by the enforced viscous force of the flow. However, beyond Re ) 210, the viscous force is not sufficient

to maintain the balance, and as a result, the axisymmetric vortex breaks. Because of the lift force, the pair of counter-rotating vortices are formed after breakage of the axisymmetric vortex. 4.4.5.2. Bubble Column. In bubble column reactors and horizontal sparged reactors (see Joshi and Sharma179), bubble size is one of the most important parameters that determines the bubble rise velocity, the gas holdup, and the interfacial area, and, as a result, phenomena such as interphase heat and mass transfer. The bubble size is highly dependent on bubble breakup and coalescence, which are influenced by the interaction with turbulence. Although the fundamental equations governing the behavior of gas and liquids, as well as the conditions at a fluid interface are reasonably well-known, in most practical applications, it is necessary to involve systems where the ratio between

Ind. Eng. Chem. Res., Vol. 48, No. 17, 2009

8265

Figure 15. Taylor-Couette flow structures as captured by various methods: (A) DWT scales (Re ) 143000), (B) DWT scales (Re ) 314140), (C) POD modes (Re ) 143000), and (D) hybrid wavelet-POD scales.

Figure 16. Ultrasound reactor flow structures, as captured by various methods: (A) DWT scales (15 W/m3), (B) DWT scales (35 W/m3), and (C) POD modes (15 W/m3).

the size of the system and the smallest continuum scale is many orders of magnitude. CFD has been used extensively in

multiphase systems to capture the vertical regions in the bubble column. Joshi180 and Sokolichin et al.181 have critically analyzed

8266

Ind. Eng. Chem. Res., Vol. 48, No. 17, 2009

Figure 17. Bubble column flow structures, as captured by various methods for single-hole, sieve-plate, and sintered-plate spargers: (A) DWT scales (singlehole sparger), (B) DWT scales (sieve-plate sparger), (C) DWT scales (sintered-plate sparger), (D) POD modes (single-hole sparger), (E) POD modes (sieveplate sparger), and (F) POD modes (sintered-plate sparger).

the CFD literature in bubble columns. According to them, the correct modeling of interphase forces and turbulence is of prime importance for capturing the physics correctly. A wide range of CFD studies on bubble column reactors have been conducted recently. Mainly, two approachessnamely, Euler-Euler182-185 and Euler-Lagrange186-191shave been adopted to simulate a gas-liquid system operating under different conditions. The k-ε model has been used for both 2D181,183,186,192-199 and 3D183,191,200-205 Euler-Euler CFD simulations. On most of the occasions, the k-ε results gave good agreements with the experimental results, but they could give no indication of instantaneous flow profile like the vertical spiral regime in sieve plate, or the plume dynamics of the bubble column. In this work, we cover some LES and DNS of the bubble column. Lakehal et al.206 were the first to employ the LES model for a bubbly flow. Their system involved a vertical bubbly shear flow at a very low void fraction (1.9%). From their study, they suggested

that, to obtain better results, the optimum cut-off filter should be 1.5 times the bubble diameter. They observed that the dynamic approach of Germano et al.207 does not perform better than the Smagorinsky model, which they felt could be due to the inadequate dimensions of their computational domain. This study suggested that the LES approach can be promising for predicting the phase velocities and the void fraction, and it emphasized the need to perform more LES simulations on complex systems such as buoyancy-driven flow operating at higher void fractions. In this context, Deen et al.208 were the first to apply the LES model to simulate a bubble column, and they reported observing a better resolved flow using LES rather than the k-ε model. Bombardelli et al.209 used the same geometry as Deen et al.208 to study the influence of different numerical schemes, of different drag models, and of initial flow conditions on LES performance. They observed that the use of a second-order FCT scheme for LES simulation enhanced its

Ind. Eng. Chem. Res., Vol. 48, No. 17, 2009

performance. However, they suggested the need to conduct more tests on LES with higher-order schemes and finer grid resolutions before identifying the best numerical scheme for LES simulations. It was observed that the LES results were very sensitive to the initial boundary conditions. This work also suggested that there is a need to pay attention toward the nearwall region, because the SGS models do not account for nearwall region processes, which can lead to the erroneous prediction of frictional stresses at the wall. Bombardelli et al.209 used a finite-element-based LES code to simulate and analyze the phenomena of wandering bubble plumes, as in the experiments that were conducted by Becker et al.210 They observed that LES was been able to capture the instability associated with wandering more vividly, compared to k-ε simulation. They analyzed the plume for coherent structures that arose from the interplay of eddies and bubbles by studying the vorticity contours. Their work showed the evolution of several eddies and their inter-relation with the bubble plume. Zhang et al.211 advanced the work of Deen et al.208 by investigating the effect of the Smagorinsky constant and the interfacial closures for drag, lift, and virtual mass force in two columns that had different aspect ratios (height/diameter ratios of H/DT ) 3 and 6). They observed that, by increasing the value of the Smagorinsky constant, the bubble plume dynamics dampens, and, consequently, a steep mean velocity profile is observed. They obtained good agreement with experimental results when the value of the Smagorinsky constant was used in the range of 0.08-0.10. They also suggested that, for taller columns (H/DT ) 6), the interfacial force closures proposed by Tomiyama gave better agreement with the experiments. As a future work, they have suggested the need for simulating the dynamic free liquid surface at the top and observing its effect on the velocity profiles in the top region. Tabib et al.212 used the LES to study the instantaneous phenomena and flow structures in a cylindrical bubble column arising out of a change in the sparger design. In this work, the LES captured the flow structures and the flow regions of the vortical-spiral flow regime in the sieve plate column (as experimentally reported by Chen et al.94), and the bubble plume dynamics along with vortical structures for the single-hole sparger. However, for sintered-plate columns, the LES showed almost-homogeneous flow conditions with an absence of dynamic eddies. Thus, LES was very effective in capturing the dynamics in the continuous phase, as a result of bubble-imparted turbulence. Compared to LES simulations, the DNS of multiphase flows has been limited to some finite number of bubbles, but they have given some excellent information on bubble dynamics. The DNS of multiphase flows (see Tryggvason et al.213) has been used both to (i) generate insight and understanding of the basic behavior of multiphase flow (such as the forces on a single bubble or a drop), how bubbles and drops affect the flow, and how many bubbles and drops interact in dense disperse flows, as well as (ii) provide data for the generation of closure models for engineering simulations of the averaged flow field. The DNS of bubbly flows, where the flow around each bubble is fully resolved and viscosity, inertia, and surface tension are fully taken into account, have already been used to study the behavior of systems that contain up to O(100) bubbles. As the number of bubbles in each period is increased, the bubbles generally rise unsteadily, repeatedly undergoing close interactions with other bubbles. Figure 18 shows some of the flow structures captured in the stirred tank region.

8267

Figure 18. Energy spectra evaluation using various transform methods in channel flow at x+ 1 ) 20 using LES monitor points: (s) FFT, (0) EIM, (+) DWT, (O) CWT.

5. Mathematical Tools From DNS, LES, and advanced experimental techniques, much information has been obtained regarding flow structures. The datasets generated from numerous experimental and computational techniques carry a plethora of information within them. In this heap of information, many times, the most relevant information becomes indiscernible. To extract the information pertinent to flow structures, many statistical and mathematical techniques have been proposed in the literature. These techniques help to identify, characterize, and quantify the flow structures in terms of their size, shape, and energy content and, finally, their relationship with the design parameters. These

8268

Ind. Eng. Chem. Res., Vol. 48, No. 17, 2009

techniques have evolved in terms of accuracy and complexity for last four decades, and they vary in terms of the amount of information that they require to identify a flow structure. These techniques make use of point data (only streamwise velocity and its time series), or a planar dataset or a complex 3D flow field of velocity/strain rate/vorticity to extract the flow structures. However, the problem lies in the fact that different techniques may detect different flow structures. In this section, these techniques are described in brief, and their relative strength and limitations are brought out. An attempt has been made to draw a correspondence between different techniques by obtaining the threshold levels at which they strike some similarity in their identification of flow structures. In the present work, these techniques are used to extract flow structure information in channel flow, where the top surface is subjected to different levels of shear boundary conditions with complete and no-slip as extreme cases. In the next section, we have covered the flow structures that have been identified in different equipment, using these techniques. 5.1. Techniques Based on A Priori Known Features of Flow Structures. These techniques are based on a central theme that the turbulent signal contains some key features, which form their footprints and they may be distinguished from the raw signal by means of a criterion that uses some physics of development of footprints. The raw data used to detect the flow structures are either instantaneous or fluctuating velocity (or pressure, or temperature, or concentration of tracers, or vorticity, or velocity gradients) in space or in time, etc. The planar dataset (PIV) can provide us with the spatial patterns and length scales of the structures; however, PIV datasets have a low data rate (up to 10 Hz), so we may not know the dynamics of these flow structures. Although the methods based on the point dataset (HFA and LDV) are able to capture the age (or time scale) of flow structures, because of the high temporal resolution of the dataset, but they are not able to predict the shape of the structures. These techniques were mainly developed to identify flow structures in turbulent boundary layers. All these techniques214-225 have been summarized in Table 3. A critical review of most of these techniques has been conducted by Yuan and Mokhtarzadeh-Dehghan.216 In the present review, their methodology (with their own abbreviations) has been used to compare various theories of heat and mass transfer. 5.2. Eddy Isolation Method. To identify the eddies, Luk and Lee226 developed an approach based on zero crossing in the fluctuating velocity. This approach is briefly discussed here. As the eddy passes through the reference point, the velocity fluctuates about a mean value. Subtracting the mean value, one gets the fluctuating part of the velocity. The fluctuating velocity has positive and negative values, and it crosses the time axis at many locations (see Figure 11A). Thus, as the fluctuating part changes from positive to negative or vice versa, it indicates that all the fluid elements (moving in a radial or axial direction) that are taken into account during the measurement have the same averaged kinetic properties. This can be considered to define the existence of an eddy. Therefore, in the plot of fluctuating velocity (u2′) versus time, as the velocity vector crosses the time axis, one eddy ends and the other enters the point of measurement. It has also been proposed that the time gap (tE) between the two successive crossings represents one eddy. Thus, different eddies show different time gaps (tE,1, tE,2, tE,3, tE,4, and tE,5), which are taken as the eddy life. For any one eddy, the velocity increases or decreases from zero, attains maxima (or minima), and then approaches zero. To estimate the characteristic eddy energy, all the fluctuating velocities

within an eddy were squared and then added together. Luk and Lee226 assumed that each eddy followed a sinusoidal shape. This approach was modified by Kulkarni and Joshi227 by considering the shapes of the actual velocity-time diagram. The modified methodology has been applied in the present work.

kE,i )

∑ 21 u′ . ∑ 2 2

i

(14)

i

The characteristic eddy velocity was then estimated from the eddy energy as uE,i ) √2kE,i

(15)

Eddy length scale (lE) was calculated from the eddy lifetime (tE) and the characteristic eddy velocity (uE). For the aforementioned calculations, an in-house code was developed. All eddies were sorted according to their age. An age distribution function (φ) was then estimated, where φ(t)∆t is the fraction of eddies that have a lifetime between t and t + ∆t. 5.3. Proper Orthogonal Decomposition (POD) Methodology. Lumley228 defined the coherent structure as a representation of the flow that has the largest projection onto the flow. Mathematically, this leads to the coherent structure being an eigenmode of a two-point correlation matrix built from a database of snapshots. These correlations can be built-up in two ways: either perform a spatial correlation or perform a temporal correlation. Lumley228 used the classical proper orthogonal decomposition technique (POD) to obtain the most optimally correlated feature from the spatial correlation of a given ensemble of flow realization. In the present work, the POD has been implemented using the temporal correlation because it reduces the computational efforts. This procedure is known as the Sirovich229 method of snapshots. The variables that have been used in this study are the out-of-plane vorticity (ω3) planar datasets and the 2D velocity (u1, u2) planar datasets. The ω3 modes are more efficient than the corresponding velocity modes for identification of the flow structures. A higher value of enstrophy, which is the square of vorticity, indicates the presence of organized vortices; however, the same may or may not be indicated by the presence of higher values of kinetic energy (square of velocity). Hence, enstrophy has been an important quantifier in regard to identifying the flow structures. In this work, the enstrophy used is the pseudo-enstrophy, because only one component of the vorticity is considered. The velocity modes have also been computed in this work to obtain the kinetic energy associated with the flow structures. These quantities are defined as k ) (ui′, ui′) ) ui′2.

(16)

ζ ) (ω3, ω3) ) ω32

(17)

The implementation of snapshot POD methodology is discussed below. Given a database having M number of instantaneous spatial snapshots (snapshot refers to the velocity or vorticity planar dataset), with each snapshot having N data points (N ) Nx2 × Nx1, where Nx1 is the number of data points in the x1 direction and Nx2 is the number of data points in the x2 direction). The application of the 2D snapshot POD method on these planar datasets enables the computation of spatial eigenmodes Ψn(x1) and accompanying temporal eigenmodes Ψn(t) that best fit the most energetic flow patterns contributing to the flow. The spatial eigenmodes represent a typical dynamical structure. The tem-

Ind. Eng. Chem. Res., Vol. 48, No. 17, 2009 Table 3. Mathematical Tools for Structure Characterization

8269

8270

Ind. Eng. Chem. Res., Vol. 48, No. 17, 2009

poral eigenmodes represents the variation in time of each spatial eigenmodes. The eigenvalue corresponding to a spatial mode represents the fraction of turbulent kinetic energy (in the case of the velocity dataset) or the fraction of enstrophy (in the case of the vorticity dataset) associated with the mode. 5.4. Discrete Wavelet Transform (DWT). The DWT separates the information content in the data from a fine scale to a coarser scale systematically, by isolating the fine-scale variability, in terms of wavelet coefficients that represent the details from the corresponding coarser-scale coefficients that depict the smoothness. This procedure can be repeated iteratively over smooth scales to obtain scale-wise detailed decomposition in a multiresolution framework. The wavelet coefficient matrix Wa,b(x) may be obtained by the convolution of the 2D snapshot data, say u(x1,x2), with a suitably chosen mother wavelet function ψa,b(x) as W (x) ) a,b

∫ u(x)ψ

a,b

(x) dx

for de-noising purposes. Therefore, de-noised 2D data may be obtained by performing an IWT after setting WT coefficients in the scales higher than the threshold scale to zero. Integrating back from the derivative domain then gives noise-reduced data uid,(a) that advantageously retains the true dynamical features of the system that otherwise would have been obscured. (4) De-noised uid,(a), obtained in the aforementioned manner, are collated in matrices for i ) 1, 2 and are subjected to 1D WT, similar to the methodology suggested by Camussi,231 i.e., u1 has been subjected to convolution in the x2 direction, while u2 has been subjected to convolution in the x1 direction to obtain the coefficient matrix Wid,(a,b)(x). (5) The structures can be identified by setting a threshold cutoff to theWid,(a,b)(x), obtained in step 4, by analyzing the local energy distribution in a scalewise fashion; this is defined as E(x) ) a

(18)

where

|[

Wd,(a,b) (x)Wd,(a,b) (x) 1 2 maxa〈Wd,(a,b) (x)Wd,(a,b) (x)〉b 1 2

E(x)a > 〈[E(x)a - 〈E(x)a〉]2〉1 / 2

ψ (x) ) 2 ψ(2 (x) - b) (19) In this case, x represents the spatial domain of the 2D planar data, with a and b defining the scaling and translation parameters. When the wavelet functions ψa,b(x) are chosen to be orthogonal to their dyadic dilations by 2-a and their translations by discrete steps 2-ab (for b ) 1, ..., 2-a, it allows a multiresolution analysis in wavelet scales a ) Ji-1, ..., 2, 1, 0 that can be performed in the 2D domain. Here, Ji refers to the number of scales in the direction of x and is dependent on the data resolution. For clarity in notation, we refer to the scale S1 as representing the smooth scale for a ) 0 and D1, D2, ... representing the detailed scales for a ) 1, 2, ... etc., with the higher scales studying the high wavenumber information in the data. An analysis of the length scale distributions in the data is then performed as follows: (1) Obtain wavelet coefficients Wa,b(x), as shown in eq 18, from the two components of velocity ui(x), where i ) 1 and 2, representing the radial and axial components. (2) Wa,b(x) may be used to obtain scalewise reconstructions uia, by performing the inverse wavelet transform (IWT). This is performed scalewise by choosing WT coefficients that correspond to a particular scale for an IWT and setting the coefficients at all other WT scales equal to zero. Note that the reconstructions uia lie in the spatial domain x and, therefore, this is assumed in the notation for u. (3) A de-noising algorithm is then applied to uia by studying the energy spectra at the WT scales. The procedure used is similar to that used by Roy et al.,230 except that it has been appropriately formulated here for application to 2D planar data, rather than for the 1D situation conducted in the work of Roy et al.230 The de-noising method uses the derivative of the scalewise data uia for the WT with a chosen mother wavelet function in eq 19. Thus, we obtain Wi′(a,b)(x) in the derivative domain. By doing this, the deterministic component (true data) gets shifted to lower scales and the noise component to higher scales.230 An automatic threshold for identifying scales with a noise component can then be obtained by studying the power spectra Pi′a vs a at different scales, a,b

Pi′P(a) )

a/2

∑ (W

a

(a,b) (x))2 i

(20)

b

and locating an inflection in the plot of Pi′a vs a. Scales higher than the inflection scale have the noisy component. This thresholding scale may be identified in an automated fashion

]|

(21) (22)

a

Here, E(x) represents the local energy value at WT scale a and 〈 · 〉 is the spatial average. As shown in eq 21, energy E(x)a has been normalized using the maximum of the average energy values obtained at all WT scales. The normalization suppresses the effects of local high-frequency perturbations. The right-hand side of eq 22 brings out the spatial energy dominance at every scale and, therefore, is useful as a suitable criteria for identification of structures in a local region. This local energy measure (LEM) methodology is similar to evaluating the local intermittency measure (LIM).227,231-233 However, the present analysis takes care to contain the effects of LIM shooting to high values in the higher scales via improved normalization over all scales rather than within scales. (6) By studying the vorticity patterns in the data, it is possible to obtain the length scale distribution of structures. For this purpose, we show the advantage of using the scalewise vorticity ωa(x), which is obtained by taking the curl of the velocity components u1d,(a) and u2d,(a) as

|

× ud,(a) ωa(x) ) ud,(a) 1 2

|

(23)

and evaluating LEM of this data (i.e., steps 4 and 5). The scalewise LEM values from vorticity data (ωa(x)) may be used to identify the location, magnitude, and shape of the structures from contour maps. Noted that, although LEM of the velocity data identifies structures, the LEM of vorticity data help to quantify the properties of these structures. The area of the contours may be used as an indirect measure to evaluate the equivalent diameter and length scale distribution characteristics of the structures in the data. (7) The scalewise percentage energy content may be calculated as E(x) ) a

∑ u ∑ ∑ (u x

a

d,(a)2 i

x

d,(a)2 ) i

× 100

(24)

The aforementioned methodology has been used for identification and length scale distribution of velocity and vorticity information for each of the equipment. 5.5. Continuous Wavelet Transform (CWT). In the case of CWT, an appropriate wavelet basis function (ψ) is chosen to decompose the 1D data u2(t) into elementary time-scale components by means of translations b and dilations a.

Ind. Eng. Chem. Res., Vol. 48, No. 17, 2009

T a,b i )

∫ u (t)ψ

a,b

2

(t) dt

The wavelet function ψ can be selected as the nth derivative of the Gaussian (DoG) function: 1 dn exp - t2 (25) n 2 dt The organization of scales is made such that a ) Jt - 1, ..., 2, 1, 0 and b ) It - 1, .., 2, 1, 0 for each scale index a. Thus, the resolution of the basis function may be made finer by increasing the value of the index It at each scale 2-a. Using the Gaussian function, the deterministic trends up to the (n - 1)th order may be eliminated by a local polynomial fit, thereby extracting the masked singularity in the data. Furthermore, only | are chosen and grouped to obtain the local maxima of |T(a,b) i the loci of maxima as a function of t across the scales a. These maxima lines represent the locations of singularities. The lines have been traced by starting at lower scales (low wavenumber region) and then following the continuous trend of locus toward the higher scales. The starting point of a singularity locus is defined as the extreme point in the loci toward the lower scale. In other words, this point indicates an evolution of the eddy at that scale which further breaks down the scale. The start-up distribution is obtained by finding all such points over a long dataset. The group is separated from the neighboring set of loci by fixing the maximum threshold limit, which is varied scalewise, i.e., at the lower scale, the limit for neighboring point identification is kept broader and it is reduced exponentially toward the higher scales. The end of the locus line is decided either by the highest scale location or by the nonavailability of neighborhood tracer point within the threshold. The selected groups have been further studied by the symbolic analysis. The symbolic analysis involves the characterization of selected groups based on the various criteria, such as the start of the singularity at the lowermost scale point, other singularities lying within the bigger singularity, breakup occurring at intermediate scales, the energy and length scale distribution because of breakup, the time span of each singularity, the time span between the two similar scale singularities (age distribution), etc. The time and scalewise singularity start-up points ts and as for each group is extrapolated in time axis on both the sides. The consecutive singularity curves have been traced and the time corresponding to them is stored as t< and t>, in the reverse and forward directions, respectively. The length scale fraction of the singularity breakup is obtained as ψ(n)(t) ≡

[ ()]

8271

5.6. Hybrid POD-Wavelet Technique. The snapshot POD and 2D wavelet transform techniques have been used extensively for educting the flow structures; however, both of these techniques suffer from certain limitations. To reduce the intensity of limitations and gain the advantage of both techniques (namely, POD and wavelet transform), a hybrid technique has been proposed by Tabib et al.234 In the hybrid POD and wavelet technique, the 2D wavelet transform technique has been applied on the spatial POD modes. This gave details of space scale information captured by a POD mode, and its scale selectivity. Wavelet transform gives the advantage of local energy decomposition. The wavelet coefficients span all the scales in the flow for a specific time. The 2D snapshot POD method allows the computation of spatial eigenmodes ψn(r) and the accompanying temporal eigenmodes φn(t) that best-fit the most-energetic flow patterns that contribute to the flow. Consider a database having M instantaneous spatial snapshots, and with each snapshot having N data points (Nx × Ny). This data is then reshaped into a data matrix A, whose column represents the pixels of a given snapshot. Thus, if we are considering a POD of the vorticity field, then the column is stacked by the N data points of vorticity that are in a given snapshot, and, thus, the size of matrix A is then N × M. The time correlation matrix C (of size M × M) is then obtained:

∫ u(r, t)u(r, t′) dr

c(t, t′) )

Σ

(27)

where Σ represents the space domain. Now, applying singular value decomposition on the time correlation matrix C will give us a set of temporal eigenmodes of size N × M and their corresponding eigenvalues λn (having a size of M × 1). The corresponding eigenvalue represents the fraction of enstrophy associated with it. The temporal eigenmodes, when projected on the original dataset, gives us the transformed spatial eigenmodes as coefficients: Ψn(r) )

1 µn

∫ φ (t)u(r, t) dt T

n

(28)

where T represents integration over the length of the acquisition sequence. The spatial eigenmodes represents a typical dynamical structure or a flow event, depending on the type of flow. The temporal eigenmodes represent the variation in time of each of the spatial eigenmodes. A linear combination of the appropriate selected K number of spatial eigenmodes, with its associated temporal eigenmodes as coefficients, can help us study the dynamics of coherent structure. K

u(r, t) )

∑ µ φ (t)Ψ (r) n n

n

(29)

n)1

where λn ) (µn)2 235

The fraction of time scale is obtained at each group of singularity, and the structure break-up distribution based on the time/length scales is obtained. Thus, a criteria for eddy passages is decided and a corresponding surface renewal has been identified at various scales. In view of studying the age distribution, span of these singularities have been evaluated. In addition, the phenomena of an eddy within an eddy is also considered and additional possible effects have been evaluated. Figure 11B shows an analysis of the raw velocity data using EIM, DWT, and CWT.

(30)

Kostas et al. suggested that the vorticity modes have a tendency to highlight vorticity structures more than the equivalent velocity modes. Hence, we use vorticity as a variable for identification of coherent structures. To reveal the space scale structure of POD spatial modes ψn(r), we subject them to scalewise decomposition through the application of DWT. The DWT separates the information from a fine scale to a coarser scale by extracting information that describes the fine-scale variability (detail coefficients or wavelet coefficients) and the coarser-scale smoothness (the smooth coefficients or mother function coefficients). The resultant wavelet coefficient matrix can be written as

8272

Ind. Eng. Chem. Res., Vol. 48, No. 17, 2009

Wa,b )

∫ Ψ (r)ψ

dr

(31)

ψa,b ) 2a / 2ψ(2a / 2(r) - b)

(32)

Σ

n

a,b(r)

Here, ψa,b is the wavelet function. In this case, a and b are scaling and translation coefficients; a moves over the scale in a diadic manner, whereas b shifts the basis function in an orthogonal manner. The POD vorticity modes capture the flow events based on their dominance over a period of time. The fraction of total enstrophy captured by the POD vorticity mode (expressed as a percentage) represents the significance or contribution of the flow event of the mode to the overall flow pattern over a period of time. The flow event may comprise of flow structures of different length scales. A scalewise decomposition of these vorticity modes then will give us the contribution of these different length scales to the flow event captured by the mode. This indicates which localized length scale is dominating the flow events captured by a mode, and it provides information about the scale selectivity of modes. Thus, we are able to relate the most-energetic flow events over a period of time (as obtained in spatial modes of snapshot POD) with the localized dominant scales that are contributing to it. 6. Dynamics and Characterization of Flow Structures Through the use of computational, experimental, and data analysis techniques, the characteristics of flow structures (such as their topology (shape and size), their energy, and their dynamics) have been uncovered. In this section, we summarize the flow structures in various equipment, which again are categorized into the following areas: near the boundary layer and walls (channel flow), flow structures observed in rotating equipment (stirred tank, annular centrifugal contactor), dispersedphase-induced flow structures (bubble column), and flows away from the walls (free shear turbulence, free jets, open channel flows). The flow structures are classified according to their energy content as obtained by DWT, POD, and hybrid wavelet-POD techniques.234,236,237 6.1. Solid/Fluid Interface: Channel Flow. DWT analysis of the channel flow at Re ) 5000 revealed scale D1 as the most energy-containing scale (50.7%) (see Figure 12A). The energy content was determined to be reduced toward scale D5. On the other hand, for Re ) 18000, scale D1 showed 21.5% of the energy, while scale D2 shows a maximum energy content that corresponds to 43.6% (see Figure 12B). From scale D2 onward, high vorticity regions were observed, even in the near vicinity of the wall. From DWT analysis, it can be concluded that as the Re value increases, the velocity gradient increases and, hence, the possibility of structures penetrating into the film region also increases. Streaks of various length scales were seen ascending away from wall. Most of the streaks followed each other, and the interactions between structures are observed to be minimal. The modewise decomposition using POD revealed that many vortices are ascending away from the wall. At higher Reynolds numbers (Re ≈ 18000), there was an increase in the frequency of occurrence of streaking phenomena, and the sizes of structures were observed to be reduced, as shown by DWT. Streaking phenomena was observed in mode 2 (18% energy), which is a significant contribution to turbulence production (see Figure 12C). Energy reconstruction with the first 10 modes contributed to 91% of the energy. An application of wavelet transform over the POD modes (i.e., hybrid wavelet-POD technique) revealed that the scalewise decomposition of Mode 1 showed selectivity toward higher length scales (scales 3 and 4). Wavelet transform of the intermediate modes (from mode 2

to mode 15) reveals their selectivity toward streaking phenomena, and these modes were selected to reconstruct the flow field, which is comprised of streaking and bursting events (see Figure 12D). 6.2. Free Shear Flows: Jet Loop Reactor. In the jet loop reactor (JLR), measurements were made at two different Reynolds numbers (Re ) 17000 and 26000). At Re ) 17000, DWT analysis showed that scale D5 was the most energetic (42.6%) and the energy content was smaller at lower scales (see Figure 13A). On the other hand, for Re ) 26000, scale D5 captured 21.5% of the energy, whereas scale D4 showed the maximum energy content (i.e., 41.6%). The energy content was reduced to 18.4% and 4.5% for scales D3 and D2, respectively (see Figure 13B). Thus, for a variation in flow rate and nozzle diameter, smaller-sized structures were observed to become more energetic. When the scalewise shapes of the structures were compared, at scale D1, large vortex tubes were observed along the periphery of the shear region with large circulation loops visible. At scale D2, leading edge vortices were observed in the shear region. The coexistence of smaller structures was observed in scales D2-D4. In POD analysis, Mode 1 (7.4% energy) captured large vortex tubes along the periphery (see Figure 13C). Mode 2 showed leading-edge vortices that are present at the periphery of the vortex tube region. As the Re value was increased from 17000 to 26000, the energetic contributions of the jet vortex tube spatial pattern to the overall flow pattern increased. Reconstruction with the first 50 modes contributed to 90% of the energy. The scalewise decomposition of Mode 1 (hybrid technique) reconfirmed that Mode 1 was indeed capturing the flow structure of the vortex tube. Scale D4 contributed ∼82% of the energy (see Figure 13D). Similarly, all the modes were decomposed scalewise. Near Mode 20, the energy of the modes was reduced and smaller vortices or some leading-edge vortices were captured that have a lesser contribution to the overall flow. 6.3. Stirred Tank (Disk Turbine, PBTD45, and HF). For the analysis of flow structures in stirred tanks, two axial (downflow pitched blade turbine (PBTD 45°) and hydrofoil (HF)) and one radial flow impeller (disk turbine, DT) were selected. In the case of DT, DWT analysis revealed that scales D1 and D2 contained 24% and 14% energy, while scale D3 has the maximum energy with 33% (see Figure 14A). For PBTD45, scale D2 shows a maximum energy of 35% (see Figure 14B). Similarly, for HF, energy was distributed prominently in scales D2 (30%) and D3 (25%) (see Figure 14C). In the case of DT, trailing vortices were observed at scales D1 and D2. Also, the structures were seen to collapse toward the side wall in the impeller plane. These structures lost energy and moved upward and downward. At the shaft of the impeller, a small structure was seen at scales D2, D3, and D4. In the case of PBTD45, a broad jet is formed inclined at an angle of 35° from the vertical axis and in the lower part of the tank, because of the thrust from the impeller. The structures were carried upward in the near-wall region and formed a large circulation loop. In the case of HF, at scale D2, similar structures were observed but with an angle of ∼25°. The energetic structures were mainly observed in the portion below the impeller, and these are carried upward in the near-wall region. Unlike DT and PBTD, no dominant structure was observed at the top surface near the shaft in HF. For all three impellers, POD Mode 1 captured the respective dominant flow pattern (see Figures 14D-F). The second mode (3% energy) represented the vortices in the shear region along the periphery of the respective impeller discharge stream of three impellers. Higher

Ind. Eng. Chem. Res., Vol. 48, No. 17, 2009

modes (>8% energy) showed vortices near the blade of the impeller and fragments of these vortices, as a result of rolloff. 6.4. Annular Centrifugal Extractor (ACE). The annular centrifugal extractor (ACE) was operated at two different Re values (143000 and 314140). At Re ) 143000, DWT analysis showed that scales D2 and D3 had most of the energy (42.4% and 41.5%, respectively; see Figure 15A). However, at Re ) 314140, scales D3 and D4 were the most dominant with energies of 49.4% and 31%, respectively (see Figure 15B). The total energy content for the higher-Re case was ∼4 times higher. The lowest scale (i.e., smooth scale (S1)) showed circulation cells that had a width and height that were almost equal to the annular width. At scale D4, irregularly shaped, slender vortices near the wall and in the central region in the axial direction were observed. Scales D2 and D3 were determined to have the dominant energy content with smaller vortices lying in the periphery of larger structures identified in scale D1. POD analysis revealed that Mode 1 (8.3% energy) had the most irregular shape, and near the wall region, one can observe few vertical, slender, high-intensity vorticity patches (see Figure 15C). Using the hybrid technique, for Mode 1, the decomposition showed that the length scales represented by scales 3 and 4 are contributing more energy, and as the modes increase, the energetic contribution of scale 4 reduces and scale 3 remains dominant in most cases, while the energetic contribution of scale 2 increases (Figure 15D). 6.5. Ultrasonic Reactor. An ultrasound reactor was operated at a power input of 15 W/m3. DWT analysis revealed that scales D4 and D5 contained 40% and 38% of the energy, respectively, and these correspond to a high energy content (see Figure 16A). On the other hand, scales D2 and D3 contribute only 12% and 8%, respectively. At scale D2, two high-intensity counterrotating vortices near the ultrasound source are observed. For a higher power input (35 W/kg), the length of these vortices became elongated (see Figure 16B). Similarly, at scale D3, the use of higher power was seen to distribute the structures more in the bulk region. Sharp energy structures of small sizes were observed in scales D4 and D5 below the horn. The sharp gradients were observed in this region are due to high-frequency contraction and rarefaction of energy packets. Mode 1 in POD analysis had 20% of the energy and showed two high-intensity counter-rotating vortices near the ultrasound source (see Figure 16C). Higher modes showed very small vortices scattered near the source. In all the modes, the region away from the source did not show any high-vorticity region. 6.6. Bubble Column. In bubble columns, the hydrodynamics is significantly affected by the design of sparges. Hence, three different spragers were selected for the study of flow structures, namely, single hole, sieve plate, and sintered plate, respectively. DWT analysis was performed for all three sparger designs. In the case of the single-hole sparger, scale D3 contained 44.2% of the energy, while scales D2 and D4 contained 17.8% and 28.9% of the energy, respectively (see Figure 17A). For the sieve-plate and sintered-plate spargers, scale D4 showed a maximum (i.e., 43% of the energy) while scales D2 and D3 contributed 12% and 31%, respectively (see Figures 17 and 17C). Although the relative energy content was similar for sieveplate and sintered-plate spargers, the spatial distribution of dominant structures was determined to vary in each of the cases. For the single-hole sparger, elongated counter-rotating trailing vortices were observed in scales D1 and D2 just above the sparger, because of the plume formation at the gas entry. In scales D3 and D4, plume-developing regions that correspond to the gas passage were observed in the snapshots. These

8273

structures were determined to oscillate in time. In the case of the sieve-plate sparger, structures were observed over the entire tank in scales D2 and D3. The nature of these vortices was determined to be elongated with irregularities in shape. An interaction of structures was observed at various locations and was observed in a snapshot animation over time. Small structures found in scale D5 were mainly dominant near the wall. In the case of the sintered-plate sparger, symmetric structures of larger size were observed at regular intervals in the axial direction. The structures were elongated and the dominant structures are observed at the near-wall region for scale D3. This was mainly because of the unidirectional motion of bubbles of equal size where significant gradients occur in the near-wall region. In POD analysis, Mode 1 (23.36%) showed that two counter-rotating vortex pairs were seen adjacent to the plume origin. Mode 2 (9.03%) showed a high-energy vortex tubelike structure within the plume region (see Figure 17D). Reconstruction with the first 50 modes contributed 88% of the energy. Compared to the single-hole sparger, there were no prominent dominant flow structures in the sieve-plate columns. The first mode captured few large structures (contributing ∼18% of the energy), while the higher modes show a large number of small scale structures (see Figure 17E). In the case of a sintered-plate bubble column, no dominant vortices in any of the vorticity modes were observed. Also, the energy associated with the vorticity modes was significantly smaller than those in the single-hole and sieveplate case (see Figure 17F). The scalewise decomposition of Modes 1 and 2 (hybrid technique) showed that these modes were selective toward scale 2, where plume-related oscillations and vortices were observed. It did not capture any other flow patterns such as vortices between the plume and the wall. Similarly, all the modes had been decomposed scalewise. At higher modes, the vortices between the plume and the wall got resolved. The wavelet transform at intermediate modes such as Mode 20 reveals their selectivity toward scale 3, where vortices within the plume and the wall are captured, and also shows selectivity for scale 2, where smaller vortices within the plume region are captured. 7. Energy Spectra and Flow Structures The energy spectrum possesses a storehouse of relevant information related to three distinct regions: the low-frequency flow structures region, an inertial region representing the cascade phenomena, and the high-frequency dissipation region. It becomes important that one should select the suitable technique (FFT, EIM, DWT, CWT, and POD) that can help in accurately evaluating the various regions of a spectra in a greater detail. Accurate information will lead to a better understanding of the turbulence cascade phenomena, and in regard to evaluating the local turbulence parameters (such as turbulent kinetic energy, rate of energy dissipation and tubulent viscosity, etc.). The wavelet-based approaches show poorer frequency resolution, because it gives averaged energy density information within the frequency scales considered, but it does offer the advantage of keeping the temporal information intact; thus, it can provide information in developing our understanding of the turbulence cascade by revealing the eddy breakup in time. The snapshot POD gives an uneven wavenumber resolution, because this energy spectrum is based on the computation of length scale and energy scale of the flow structures observed in the modes. These techniques are compared based on (a) the frequency resolution offered by the technique, (b) the phenomena and assumptions considered by the technique while evaluating the

8274

Ind. Eng. Chem. Res., Vol. 48, No. 17, 2009

Table 4. Methods Used for Estimation of Energy Spectra Sr. No.

technique

1 2

FFT EIM

3

DWT

4

CWT

5 6

POD SRLIM

expression for energy spectrum

general remarksa

E(κ) ) ∑κ21+κ22)κ2 ui(κ1,κ2)*ui(κ1,κ2) kE,i ) ∑i 1/2u′22 /∑i where uE,i ) (2kE,i)1/2 and κ ) 2π(∆tE)/uE Eii(κ) ) ∑b Wia,b(x)/nb where κ2 ) κ21 + κ22 and κi ) 2π2a/∆xi Eii(κ) ) ∑b Tia,b(x)/N(x) where κ2 ) κ21 + κ22 and κi ) 2π2a/∆xi Ψn(x) ) (1/µn) ∫Tφn(t)u(x,t) dt SRLIM(a,b) ) LIM(a,b)∑b W(a,b)2/(∑a∑b W(a,b)2)

1, 3, 4, 7 2 1, 3, 5, 7 1, 3, 5 6 8

a General remarks are identified as follows: (1) FFT offers the finest frequency resolution (0.0015 Hz) followed by CWT and DWT. (2) EIM, as it cannot identify an eddy within an eddy, cannot fully capture the lower subeddies comprising the higher-frequency-dissipation region, and the higher-length-scale low-frequency region. But it captures inertial region of the spectra very well. (3) FFT, CWT and DWT implicitly consider the effect of various eddies and subeddies as being spread across various scales or frequencies. (4) The FFT-based energy spectrum is the most accurate, as a result of its finer frequency resolution; however, the lack of temporal information is its major limitation to study eddy breakup. (5) Wavelet-based approaches (DWT and CWT) thus can provide information in developing our understanding of the turbulence cascade by revealing the eddy breakup in time. (6) The snapshot POD gives an uneven wavenumber resolution as this energy spectrum is based on the computation of length scale and energy scale of the flow structures observed in the modes. (7) The inertial range is captured nicely by DWT and FFT. (8) SRLIM involves multiplying the LIM of a wavelet coefficient of a given scale with the relative power of that scale (power of that scale/total power). Hence, it accounts for the relative strength of the scales.

spectra, and (c) the accuracy of the technique over a given range. These techniques are listed in Table 4, along with their relative merits. In the literature, there have been several attempts to accurately model the irregular nature of energy dissipation in an energy cascade (known as intermittency) by means of energy spectrum analysis. An energy spectrum can give vital information regarding the small scale dissipative range responsible for most of the energy dissipation. The fundamental basis of these cascade models is that the energy transfer occurring in the inertial range may be conceptualized as a cascade process, wherein each eddy of size l subdivides into nl pieces of size l/nl in such a way that the energy flux per unit mass is redistributed equally/unequally among the nl subeddies. The unequal split of the energy flux is modeled as the source of intermittency. Most of these intermittency theories have developed on the seminal work of Kolmogorov’s theory of isotropic turbulence,238,239 whose hypothesis suggested a power-law relationship between the spectral density of energy in inertial range to its associated wavenumber. (33) E(κ) ) Cε2 / 3κ-5 / 3 where κ is the wavenumber in the inertial range and C is a universal constant with a value of 0.5. These theories assume that the conservation of energy flux holds in one dimension, as it does in three dimensions. Later, a correction to this Kolmogorov theory was proposed by Landau and Lifshitz, on the pretext that instabilities do not allow the transfer process to be inertial and some amount of it is continuously dissipated. In view of the argument, Kolmogorov239,240 and Obukhov241 proposed a refinement on the basis of nonlinear interaction between eddy velocity and vorticity by incorporating an addition term:

()

E(κ) ) Cε2 / 3κ-5 / 3 ln

κ κ0

β

(34)

where κ0 is the wavenumber at the beginning of the inertial range. The revised relationship indicates that, at every scale, only a fraction of the energy of the earlier stage is transferred to the next scale. Various intermittency thoeries available in the literature are listed in Table 5. Kulkarni et al.246 proposed a strategy based on the synergistic combination of energy spectrum and the intermittency theories for revealing different stages in the turbulent cascades. They used these scales to estimate turbulent viscosity, and the estimates of viscosity obtained using different intermittency models were compared

with each other and with the estimates obtained using the k-ε turbulence model. They found that the random β model and the P model gave good agreement of the turbulent viscosity with the k-ε turbulence model. 7.1. One-Dimensional Energy Spectrum. The FFT-based energy spectrum has been used more often in literature to depict the -5/3 slope obtained in the inertial range of the Kolmogorov theory. Kolmogorov238,239 suggested the equal breakup of energy based on the second-order structure function with a -5/3 law over the energy spectrum. Although Kolmogorov’s K41 theory was able to explain the experimentally observed slope in the inertial range, it is not a realistic model, because it does not consider the intermittent behavior of energy dissipation. Just as K41 has limitations, the FFT method also has limitations, in that it does not give information on temporal occurrence of a frequency; hence, it will not be able to provide information on break-up and eddy-size information on the basis of temporal dynamics. Whereas the CWT-based technique saves the temporal information and provides good frequency resolution. These intermittency theories involve some adjustable parameters to fit the cascade. The relative performance of various techniques in channel flow is shown in Figure 18. As can be seen in the spectra (see Figure 18), the eddy isolation method (EIM) cannot fully capture the lower subeddies that comprise the higherfrequency dissipation region and the higher-length-scale lowfrequency region. The EIM methodology has been able to capture only the inertial region of the spectra. However, in the case of FFT, the velocity signal is approximated by sinusoidal functions of varying frequency that have infinite length. The passage of a dominant eddy is represented by the frequency with a high energetic contribution. Furthermore, the waveletbased approaches approximate the eddies in terms of wavelets that have localized time information. These methods help to implicitly consider the effect of various eddies and subeddies as being spread across various scales or frequencies. In the time series, the singularities are seen to be occurring at different times, and the wavelet approach preserves the temporal information of their occurrences, along with their effect on various scales. The resultant effect can be seen in the WT-based energy spectra (see Figure 18), where it captures all three ranges, unlike the EIM procedure. The FFT-based energy spectra is probably the most accurate, as far as the finer frequency resolution is concerned; however, the method lacks the temporal information related to eddy dynamics and, hence, does not provide much insight on when the eddy breaks up and how this break-up

Ind. Eng. Chem. Res., Vol. 48, No. 17, 2009

8275

Table 5. Various Intermittency Theories in the Literature Sr. No.

theory

reference 242

1

β model

Frisch

2

random β model

Benzi et al.243

3

P model

Meneveau and Sreenivasan244

4

L model

Sreenivasan and Stolovitzky245

5

PL model

Sreenivasan and Stolovitzky245

remarks The energy transferred at each cascade stage is only a fraction (1 - β) of the energy at that stage. This clearly implies that at each stage, fraction β is dissipated in a self-similar manner. This model considers varying values of β at different stages in the cascade. Here, the value of β changes continuously with a definite probability density characterized by the system under consideration. The suggested distribution appears as p(β) ) xδ(β - 0.5) + (1 - x)δ(β - 1) and the distribution of β can be estimated for different x values and a range for β. In the P model, the energy distribution over the cascade occurred through a binary break-up mechanism, such that in the inertial range eddies break into two equal sizes with a fixed proportion of eddy energy (P1:P2, such that P1 + P2 ) 1.0). Different sets of P yield a distinct intermittent nature. The authors suggested P1:P2 values of 0.7:0.3 for a more realistic turbulent flow. In the L model approach, the breakup was assumed to follow an equal energy distribution, such that the length scales of daughter eddy maintain a fixed proportion (L1:L2) over the cascade. PL model envisages an unequal energetic and length wise distribution of eddy break up. This is the more realistic scenario. The PL model (with parameters P1 ) 0.7, P2 ) 0.3, L1 ) 0.7, L2 ) 0.3) are seen to be showing better performance in capturing the cascade as compared to the β and the log-normal model (Tabib et al.212).

information finds its way into the energy spectra. However, the wavelet-based approaches show poorer frequency resolution, because it gives averaged energy density information within the frequency scales considered, but it does offer the advantage of keeping the temporal information intact and, thus, can provide information in regard to developing our understanding of the turbulence cascade by revealing the eddy breakup in time. Among the wavelet-based approaches (DWT and CWT), the CWT offers finer frequency resolution, because of the frequency scale overlap; however, this feature also results in the slight overestimation of energy in the spectrum (see Figure 18). In this section, the flow structures in various equipments have been identified based on two concepts: (a) the frequencies with high energies, whicha re observed as peaks in the spectral analysis, are considered as the flow structures, and (b) given the intermittency-based concept, the flow structure is viewed as an intermittent event with high energy. 7.1.1. Channel Flow. Figures 18A-C show energy spectra obtained in the near-wall regions of channel flow. Near the wall, a peak at a frequency of