2936
Ind. Eng. Chem. Res. 2007, 46, 2936-2950
Modeling Flow Pattern Induced by Ultrasound: The Influence of Modeling Approach and Turbulence Models T. Kumaresan, Ajay Kumar, Aniruddha B. Pandit,* and Jyeshtharaj B. Joshi Institute of Chemical Technology, UniVersity of Mumbai, Matunga, Mumbai-400 019, India
Mean flow and turbulence parameters have been measured using laser Doppler anemometer (LDA) in an ultrasound reactor at the ultrasonic power of 15 kW m-3. The liquid circulation velocities are dominant in the zone nearer to the source of energy and are substantially low at the walls and at the bottom of the reactor. The levels of turbulence kinetic energy and dissipation rate are high near the horn and decrease rapidly with increasing distance from the horn. Comparisons of LDA measurements and computational fluid dynamics (CFD) predictions have been presented. The effects of the model for the horn-fluid interaction (velocity boundary condition approach and acoustic boundary condition approach) and the turbulence model (standard k- model, modified k- model, standard k-ω model, Reynolds stress model (RSM), and large eddy simulation (LES)) on the predicted flow pattern (mean velocities, turbulent kinetic energy, and turbulent energy dissipation rate) have been investigated. 1. Introduction The application of power ultrasound in the chemical process industries is growing because of a wide range of emerging applications. Power ultrasound is now well-known to have significant effects on the rates of various chemical and physical processes in the chemical industry: deagglomeration, particle size reduction, homogenization, emulsification, filtration, crystallization, extraction, degassing and stripping, free radical initiated reactions, etc. The ultrasound can also provide an increased rate together with the selectivity of many different chemical reactions.1-3 All these effects are attributed to acoustic cavitation. Also, the acoustic streaming is effective in improving the rate of the transport processes occurring on the fluid and solid interface including heat transfer, changes in morphology of the biological cells, and removal of loosely adhering surface layers. When ultrasound passes through the liquid, acoustic jet streaming is generated. Acoustic streaming is defined as a tendency for a steady fluid circulation near the surface of obstacles and vibrating elements, and near bounding walls in a high-intensity sound field. Acoustic jet streaming is induced because of spatial attenuation of a sound (pressure) wave in free space and the momentum transfer to the medium by the vibrating object (horn). The propagation of a high-intensity sound wave in a fluid medium causes the decrease of pressure strong enough to produce a steady bulk fluid flow. In turn, acoustic jet streaming is responsible for convective flow produced from the vibrating horn. The flow generated from the ultrasound horn (22 kHz) when mounted vertically goes into the bulk, gets diverted by the vessel bottom, flows along the walls, and again gets diverted at the free liquid interface and returns back to the horn zone. This circulatory three-dimensional flow is usually turbulent and is sustained by using some of the energy released by the horn. The flow characteristics of the ultrasound reactor strongly depend on the power of ultrasound and the geometry of the reactor. Therefore, it is of interest to understand the fluid flow pattern with respect to the ultrasound power. * To whom correspondence should be addressed. Phone: 91-22414 5616. Fax: 91-22-414 5614. E-mail:
[email protected].
The first theoretical analysis of acoustic streaming phenomena was accomplished by Rayleigh.4 More investigation in the theory was carried out by Schlichting,5 Nyborg,6 and Lighthill,7 where emphasis was laid on the fundamental role of dissipation of the acoustic energy in the form of the changing of the gradients in the momentum flux. Jackson and Nyborg8 have also investigated acoustic streaming induced by sonic longitudinal vibration. However, the articles published in the field of acoustic streaming are restricted to qualitative description of the phenomena, and very few researchers have thrown light on acoustic streaming in a quantitative way. Cadwell and Fogler9 have obtained intense convective patterns in a high-frequency reactor of cylindrical shape (800 kHz). They have estimated fluid velocity by noting the length of a streak made by a moving particle, which was measured by photographic technique. The fluid velocities in the vicinity of ultrasound horn were found to be in the range of 0.3-0.5 m/s. Gondrexon et al.10 have investigated the acoustic streaming in a high-frequency transducer (500 kHz, diameter 0.04 m) with visual qualitative observation in a poly(vinyl chloride) (PVC) cylindrical vessel (diameter 0.10 m, height 0.10 m) equipped with a vibrating stainless steel plate fixed to the bottom (diameter 0.15 m, thickness half a wavelength 5.84 × 10-3 m). The transducer is stuck to the stainless steel plate. They found that the velocity ranged from 0.01 to 0.03 m/s for an electric output of 40 W to 0.05-0.10 m/s at maximum power input (100 W). Chouvellon et al.11 have used laser tomography technique and investigated acoustic streaming fluid velocity in a high-frequency ultrasonic reactor (500 kHz). They have observed the effect of the power, the water height, and the fluid viscosity on acoustic streaming velocity. The velocity was found to increase in a nonlinear way with the electric power. It increases more from 50 to 100 W and shows a lesser increase beyond 100 W of power input. They have found streaming velocity in the range of 0.0066-0.033 m/s with the electrical power range of 10-175 W. They have observed the average water velocity to decrease with an increase in the water depth. They have seen the effect of fluid viscosity and concluded that the average velocity increases with an increase in the liquid viscosity until it reaches a threshold and then decreases. Frenkel et al.12 have investigated the liquid velocities in a rectangular ultrasonic reactor having dimensions of (0.25 × 0.17 × 0.14 m) and an ultrasonic frequency of 3 MHz. In this study, the authors
10.1021/ie060587b CCC: $37.00 © 2007 American Chemical Society Published on Web 11/07/2006
Ind. Eng. Chem. Res., Vol. 46, No. 10, 2007 2937
have used particle image velocimeter (PIV) to measure the fluid velocities in the vessel. They found that the velocity increases with an increase in the ultrasound intensity and have developed the following correlation,
VS ) 0.0992I
(1)
where VS is acoustic streaming velocity in cm/s and I is the irradiation intensity in W/cm2. Kim et al.13 have observed that the effects of ultrasonic vibration on flow behavior are vastly different depending on the heat transfer regime and the amount of the dissolved gas in the liquid. Such a difference in the flow behavior eventually affects the degree of heat transfer enhancement. They have also concluded that the number density of cavitation bubbles plays a critical role in enhancing the heat transfer as well as vigorous bubble motion due to acoustic streaming; in addition to increased bubble number density, the acoustic streaming was also revealed to greatly enhance the heat transfer rate. Very few researchers have investigated quantitatively the streaming velocity in a low-frequency ultrasound reactor. Cadwell and Fogler9 have again estimated acoustic streaming velocity in a low-frequency (20 kHz) reactor to be of the order of 0.7-1 m/s using streak photography technique. Dahlem et al.14 have reported acoustic streaming in low-frequency batch ultrasound reactors (20 kHz), namely, telesonic horn (radial horn) and stepped horn (longitudinal horn). They used PIV technique to measure the velocity field associated with the acoustic streaming in the reactor. The average fluid velocity for the case of the radial horn was found to be ∼0.05 m/s, whereas in the case of longitudinal type of horn, the average velocity was found to be ∼1 order of magnitude higher than the radial horn under otherwise similar conditions of the vessel geometry and the power input. The higher velocity detected (∼0.35 m/s) results from the concentration of all the energy on a much smaller emitting surface. Vichare et al.15 and Kumar et al.16 indirectly estimated the average velocity of acoustic streaming indirectly with the help of mixing time in a low-frequency (20 kHz) ultrasonic horn immersed cylindrical reactor. They have observed streaming velocities in the range of 0.05-0.2 m/s, which are substantially lower than those observed by Cadwell and Fogler9 and may be due to the higher level of power input in the latter’s work, resulting in more energy dissipation in the form of cavitational events. All the previous works reported have used either intrusive measurement techniques or visual observations techniques such as streak photography to estimate the velocity magnitude prevailing in the ultrasonic reactors. Dahlem et al.14 used PIV technique for measuring the mean acoustic streaming velocities. The turbulence properties have not been reported. However, equipment such as PIV can give a good spatial resolution, but turbulence measurements need a temporal resolution in order to quantify the turbulent flows. Because of the large number of measurement techniques, one can clearly observe a very wide difference in the measured velocity (0.05-1 m/s). This was probably because the measurements have not been validated by the establishment of overall energy balance. In our earlier work byKumar et al.,17 flow phenomena of acoustic streaming was measured in the low-frequency (20 kHz; 15, 25, and 35 kW/m3) ultrasonic horn reactor using laser Doppler anemometer (LDA) measurement and good energy balance was established. The computational fluid dynamics (CFD) model in the earlier work required measured boundary conditions at four locations to solve the conservation equation, which is not a very satisfactory situation. In the present work, an effort has been made
Figure 1. Experimental setup for LDA measurements.
to model the acoustic streaming flow phenomena using CFD using different numerical solution schemes. So, 3D unsteadystate CFD simulations were performed using various hornfluid and turbulence models to predict the flow pattern and mixing time. The CFD predictions have been compared with LDA measurements. Relative merits of the various models have been presented so as to obtain better flow prediction of ultrasound reactors. 2. Experimental Technique (Flow Pattern Measurements) Flow field in an ultrasound reactor has been studied by laser Doppler anemometry (LDA). The LDA setup comprised a Dantec 55X modular series along with electronic instrumentation and a personal computer. A two-dimensional LDA was used to measure the velocities in an ultrasound reactor. The light source was an argon-ion laser; the power of the emitted beam could be regulated up to 5 W. This beam is split by a modular optical system into three parallel beams equidistant to the optical axis: a blue (488 nm), a green (515 nm), and a blue-green beam. For the measurement of flow reversals, blue and green beams are provided with a shift of 40 MHz in frequency with respect to the blue-green beam, with an aid of a Bragg cell. The three beams pass through the transmitting and focusing optics to intersect and form a measuring volume within a cylindrical vessel (diameter ) 0.135 m, volume ) 2000 mL). The cylindrical vessel was immersed in a square tank filled with water to minimize optical distortions. The laser probe volume was positioned at the desired point with the help of a traversing system with an accuracy and reproducibility of (0.5 mm. The receiving optics with photo detectors were operated in a forwardscattering manner. A Burst spectrum analyzer (BSA) was used to convert analogue signals from the receiving optics into Doppler frequencies and velocities. The Burst signal analyzer uses frequency domain burst detection to convert the analogue signal to instantaneous velocity measurements. The effects of signal sampling rate and sample size on the reproducibility of velocity measurements were investigated. The average validated sample size of 20 000 was obtained at an average data rate range of 1-1.5 kHz. In all the experiments, tap water was used as a working fluid. No additional tracer particles were added in this experiment. The schematic diagram of the experimental setup used is shown in Figure 1. The measurements of axial and radial components of velocities were carried out in a vertical plane of the vessel, which is perpendicular to the optical axis of beams, whereas the tangential components of the velocities were measured in a vertical center plane passing through the optical axis of beam. The measurements of velocities were made at 13 axial locations from the horn and 8 radial locations up to the wall of vessel (Figure 2). In the present work, the radial profiles at z ) 0.13H, 0.53H, and 0.8H axial locations have been shown for sake of brevity. These locations represent (1) near the horn, (2) the center of the vessel, and (3) near the bottom of vessel, respectively, where
2938
Ind. Eng. Chem. Res., Vol. 46, No. 10, 2007
essentially an inhomogeneous wave equation that is derived by manipulating the continuity and Navier-Stokes equations. The Ffowcs-Williams and Hawkings acoustic equation is given below,
1 ∂2p′ - ∇2p′ ) 2 2 co ∂t ∂2 ∂ (T H(f)) - ((Pxyny + Fux(un - Vn))δ(f)) + ∂x∂y xy ∂x ∂ ((F V + F(un - Vn))δ(f)) (2) ∂t o n
Figure 2. Measurement locations in axial and radial directions in the reactor.
H is the height of the liquid (H ) 0.15 m) in the vessel and z is the axial distance from the center of the horn tip. Submergence of the horn tip in liquid was 0.02 m from the upper surface of the liquid. The effects of ultrasound power (P/V ) 15 kW/m3) on velocity were studied for all the axial and radial locations at constant power density. Further, all the experiments have been performed with a constant diameter of the ultrasonic horn (0.013 m). The ultrasound reactor produces substantially downward flow in the central region of the vessel and upward flow in the annular region near the wall. Since the liquid phase is in the batch mode, a mass balance across any cross section of the vessel should indicate the net flow to be zero. The LDA measurements have satisfied this self-consistency in the experimental data, and hence, in turn, these values are used as boundary conditions for the CFD simulations. 3. CFD Simulations 3.1. Modeling of Horn-Fluid Interaction. 3.1.1. Velocity Boundary Condition Approach. The governing equations18 were solved in an unsteady-state mode through the forced velocity boundary condition approach from LDA measurements at different locations. For this purpose, five different sets of locations were considered, as shown in Figure 3. The action of the vibrating horn was thus simulated by providing the boundary conditions (measured fluid flow characteristics using LDA). The boundary conditions consist of mean velocities (radial, axial, and tangential) and the turbulent kinetic energy (k), the values of which were measured experimentally for the system under consideration. Thus, if the tank size/geometry is modified, it becomes necessary to make these measurements again for the implementation of this type of boundary condition. The applicability of velocity boundary condition approach is severely limited by the availability of the experimental data and, hence, cannot be used to predict the flow generated by different ultrasound sources in different vessel configurations. 3.1.2. Acoustic Boundary Condition (ABC) Approach. The time-dependent interaction of the horn and the fluid is modeled using the acoustic equation proposed by Ffowcs-Williams and Hawkings.19 This approach makes use of the translation movinggrid technique in the axial direction, and the ultrasound propagation is modeled using the acoustic equation, which is
where p′ is the sound pressure in the fluid medium. In the above equation whenever f ) 0, it indicates the mathematical surface to embed the exterior flow in the fluid medium (f > 0), which facilitates the use of generalized function theory and the free space Green function to obtain the solution. The surface (f ) 0) indicates the source or emission surface and simultaneously made coincident with the medium permeable.
Txy ) Fuxuy + Pxy - co2(F - Fo)δxy
(
Pxy ) pδxy - µ
∂ux ∂uy 2 ∂uz + δ ∂x ∂x 3 ∂z xy
(3)
)
(4)
The solution to eq 2 is obtained using the free space Green function. The complete solution consists of the surface integrals and volume integrals. The surface integrals represent the contributions from monopole and dipole acoustic sources and partially from quadrupole (volume) sources. The volume integrals solely represent quodrupole sources in the region outside the horn surface. However, the volume integrals are dropped in the present model because of the importance of sound propagation from the horn surface, and hence, the equation contains the following axial terms:
b,t) + p′L(x b,t) p′(x b,t) ) p′T(x
(5)
When the vibrating fluid medium integration surface coincides with the horn tip, the two terms on the right-hand side of eq 5 are referred to as thickness and loading terms, respectively. When the vibrating fluid medium integration surface does not coincide with the horn tip (cavitation and decoupling), then the pressure of the wave is accounted for. The acoustic moving-grid approach suffers from large computational expenses. The presence of moving grids requires proper care to be taken while handling the grids at the interface between the two regions i.e., mass flow continuity should be satisfied at the interface. Because it relies on solution of full-time varying flow in the whole domain, its computational requirements are greater by order of magnitude than those required by the steadystate simulations. Hence, a sufficient (28% of the total grids) number of grids (1 209 000 (z × r × θ ) 93 × 52 × 250)) were used near the ultrasound source zone in the present work. 3.2. Turbulence Modeling. 3.2.1. Standard k- Model. The standard k- model18 is essentially a high Reynolds number model and assumes the existence of isotropic turbulence and the spectral equilibrium. It proposes the relation for the eddy viscosity in terms of k and with the help of the turbulence parameter Cµ. In addition, the triple-velocity correlations in the transport equation for the energy dissipation rate are modeled with the help of two more constants (C1 and C2), whose values have been derived from the measurement of flow characteristics of simple two-dimensional equilibrium flows. In the k- model,
Ind. Eng. Chem. Res., Vol. 46, No. 10, 2007 2939
Figure 3. Locations of boundary conditions (BC) using k- model: (A) 1BC, 1 at z ) 0.05H; (B) 2BC, 1 at z ) 0.05H, 2 at z ) 0.71H; (C) 3BC, 1 at z ) 0.05H, 2 at z ) 0.32H, 3 at z ) 0.71H; (D) 3BC, 1 at z ) 0.1H, 2 at z ) 0.2H, 3 at z ) 0.3H; (E) 4BC, 1 at z ) 0.05H, 2 at z ) 0.32H, 3 at z ) 0.71H, 4 at z ) 0.94H. The contour shows the mean velocity magnitude in m/s.
the length and the time scales are built up from the turbulent kinetic energy and the dissipation rate using the dimensional arguments:
kxk k and τ0 ∝ l0 ∝
(6)
of dissipation is assumed to depend only on the length and the time scales. This can be written in a functional form as
φ ) fφ(l0,τ0) or φ ) fφ(k,) Dimensional analysis for φ requires that
This yields
φ ) C2
k2 ν t ) Cµ
[
]
2 k
(7)
in which Cµ is an empirical constant. The k and are obtained from the solution of the following transport equations:
∂k ∂k ) + ujk ∂t ∂xk
(11)
( )
∂uji ∂ujk ∂uji ∂ νt ∂k νt + + - + ν∇2k (8) ∂xk ∂xi ∂xk ∂xk σk ∂xk
The production of dissipation is assumed to be proportional to the production of turbulent kinetic energy, the length, and the time scale of turbulence. A simple dimensional analysis for P along with the above assumption leads to
P ) C1 P k
(12)
in which C1 is a dimensionless constant, whereas P is given by
The transport equation for is rewritten as
∂ ∂ + uji ) ν∇2 + P + D - φ ∂t ∂xi
(9)
Here, the P, D, and φ represent the production of dissipation, the turbulence diffusion of dissipation, and the turbulence destruction of dissipation, respectively. The diffusion term can be modeled using the gradient approximation similar to the k-equation:
( )
∂ νt ∂ D ) ∂xk σ ∂xk
P ) -τik
(10)
in which σ is an empirical constant. Analogous to the modeling of in the k-equation for one-equation models, the destruction
∂uji ∂xk
(13)
These yield
∂uji P ) -C1 τik k ∂xk
(14)
It was shown that this eq 14 is strictly valid only when: (1) there is a clear separation of time scales of turbulent fluctuations and the mean velocity field and (2) the level of isotropy is not too high. Thus, the modeled version of the transport equation for turbulent dissipation rate reduces to
2940
Ind. Eng. Chem. Res., Vol. 46, No. 10, 2007
∂ ∂ + uji ) ∂t ∂xi
Generation of turbulence kinetic energy (Gk) is given by
( )
(
)
∂uji ∂ujk ∂uji 2 ∂ νt ∂ + C1 νt + - C2 (15) ν∇ + ∂xk σ ∂xk k ∂xk ∂xi ∂xk k 2
Thus, the standard k- model has five empirical constants, Cµ, σk, σ, C1, and C2, in its formulation. C2 is determined using the experiments on the decay of k behind a grid, and a value of C2 ) 1.92 is suggested. The experiments on the local equilibrium shear layers suggested that Cµ is nearly 0.09. σk, σ, and C1 are fixed using the computer optimization, and their suggested values are 1.0, 1.3, and 1.44, respectively. It should be noted that the above values for five constants are not universal. Even though the standard k- model works reasonably well for most of the engineering problems, it suffers from few modeling difficulties such as (i) the assumption of isotropic turbulence (the mean velocity gradients is proportional to the stress imparted) and (ii) the use of the above-mentioned values for the five model parameters which were, in fact, derived for some selected flows from the literature. As a consequence, this model has a limited applicability in the case of (i) curvature effects, (ii) boundary layer separation, and (iii) irrotational strains. 3.2.2. Modified k- Model. Abujelala and Lilley20 have suggested the use of the standard k- model with modified values for the turbulence parameters (Cµ ) 0.125, C1 ) 1.59, and C2 )1.62) for the swirling and rotating flows. Similar modified turbulence parameters have also been tried in the present work to analyze the flow predictions for the case of a cylindrical ultrasound reactor. 3.2.3. Standard k-ω Model. The standard k-ω model is an empirical model based on model transport equations for the turbulence kinetic energy (k) and the specific dissipation rate (ω), which can also be thought of as the ratio of to k. The transport equations for k and ω are as follows:
( )
∂(Fk) ∂(Fkui) ∂ ∂k ) Γ + Gk - Yk + ∂t ∂xj ∂xj k ∂xj
( )
∂(Fω) ∂(Fωui) ∂ ∂ω ) Γ + Gω - Yω + ∂t ∂xi ∂xj ω ∂xj
(16)
Gk ) -Fu′iu′j
∂uj ∂xi
(22)
The following expressions are used to define and estimate the other turbulence parameters. The generation of specific dissipation rate (Gω) is given by
ω Gω ) R G k k
(23)
where
R)
(
)
R∞ R0 + Ret/Rω R* 1 + Ret/Rω
(24)
where Rω ) 2.95. R* and Ret are given by eq 21. For the dissipation of k, the model constants are as follows:
R/∞ ) 1 Rk ) 6
R∞ ) 0.52 Rω ) 2.95
R0 ) 1/9
β/∞ ) 0.09
βi ) 0.072
Rβ ) 8
σk ) 2.0
σω ) 2.0
{
Yk ) Fβ*f/βkω 1 if χk e 0 2 / fβ ) 1 + 680χk 1 + 400χk2 if χk > 0 χk ) β* ) β/∞
(
1 ∂k ∂ω ω3 ∂xj ∂xj
(25)
)
4/15 + (Ret/Rβ)4 1 + (Ret/Rβ)4
Rβ ) 8β/∞ ) 0.09 (17)
Ret is given by eq 21. Dissipation of ω is given by
Effective diffusivities are given by
Γk ) µ +
µt µt Γ )µ+ σk ω σω
Yω ) Fβfβω2 (18)
fβ )
1 + 70χω 1 + 80χω
where
Fk µt ) R* ω
(19)
The coefficient R* damps the turbulent viscosity, causing a lowReynolds-number correction. It is given by
R* )
R/∞
(
)
R/0 + Ret/Rk 1 + Ret/Rk
(20)
where
Ret )
Fk R ) 6R/0 ) βi/3βi ) 0.072 µω k
(21)
χω ) | Ωij )
ΩijΩjkSki | (β/∞ω)3
(
(26)
)
1 ∂ui ∂uj 2 ∂xj ∂xi
3.2.4. RSM Model. The Reynolds stress model (RSM) makes use of transport equations for all six independent components of the Reynolds stress tensor and an equation for the energy dissipation rate. The equation of turbulent stresses contains the terms, namely, pressure rate of strain (which contains fluctuating pressure and velocity gradients) and the flux of Reynolds stresses. In addition to the modeling of the terms in the turbulent energy dissipation rate equation, these terms need to be modeled
Ind. Eng. Chem. Res., Vol. 46, No. 10, 2007 2941
accurately. Therefore, though RSM takes into consideration the transport of turbulent stresses, its parameters are not universal (i.e., model should be calibrated for different types of flows). Further, the extended system often has difficulties to produce converged solutions, and thus, it is computationally expensive. Equations for turbulent fluctuations are obtained by subtracting the Reynolds averaged equations from the governing Navier-Stokes equations. The relevant equations are shown below. The continuity equation is
∂u′i )0 ∂xi
(27)
∂ui ∂ ∂ [ui - ui] + (uj + u′j) (ui + u′i) - uj ) ∂t ∂xj ∂xj -
∂ 1 ∂ [P - P] + ν∇2(ui - ui) + [τij] (28) F ∂xi ∂xj
which gives
∂ui ∂u′i ∂ui ∂u′i ∂ui ∂u′i + uj + u′j + u′j - uj ) + uj ∂t ∂xj ∂xj ∂xj ∂xj ∂xj -
∂τij 1 ∂P′ + ν∇2u′i + F ∂xi ∂xj
or
(29) ∂τij ∂u′i ∂ui ∂u′i 1 ∂P′ ∂u′i ) -u′j - u′j + ν∇2u′i + + uj ∂t ∂xj ∂xj ∂xj F ∂xi ∂xj
Equations 27 and 29 represent the field equations for the fluctuating components. These equations may be solved to obtain the fluctuating component u′i. However, such a solution is very difficult to obtain, since the solution for u′i depends on the global history of the mean-velocity field with an implicit dependence on its own initial and boundary conditions. The Reynolds stress transport equation is obtained from the second moment of eq 29. Equation 29 for the ith component is multiplied by the fluctuating velocity u′j , and the equation for the jth component is multiplied by the velocity, u′i . These two equations are added, and the resulting equation is time averaged to obtain the transport equation for τij.
(
)
∂τij ∂τij ∂uj ∂ui ) - τik + τjk + uk ∂t ∂xk ∂xk ∂xk
[
] ] ( )
∂u′iu′ju′k 1 ∂ + (P′u′iδjk + P′u′jδik) + ∂xk F ∂xk
[
∂u′i ∂u′j P′ ∂u′j ∂u′i + - 2ν + ν∇2τij (30) F ∂xi ∂xj ∂xk ∂xk ∂τij ∂τij + uk ) ∂t ∂xk
(
- τik
)
∂uj ∂ui ∂ + τjk C + Πij - ij + υ∇2τij (31) ∂xk ∂xk ∂xk ijk
The terms Cijk, Πij, and ij in eq 31 are given as
Πij )
[
]
P′ ∂u′j ∂u′i + ) pressure-strain correlation F ∂xi ∂xj
ij ) 2ν
( )
∂u′i ∂u′j ) dissipation rate correlation ∂xk ∂xk
1 Cijk ) u′iu′ju′k + (P′u′iδjk + F P′u′jδik) ) third-order diffusion correlation 3.2.5. LES Model. Large eddy simulation (LES) provides an alternative approach in which the large eddies are computed in a time-dependent simulation that uses a set of “filtered” equations. These can handle laminar, turbulent, and transitional flows. Filtering is essentially a manipulation of the exact Navier-Stokes equations to remove only eddies that are smaller than the size of the filter, which is usually taken as the mesh size. Like Reynolds averaging, the filtering process creates additional unknown terms that must be modeled in order to achieve closure. Statistics of the mean flow quantities, which are generally of the most engineering interest, are gathered during the time-dependent simulation. One might also argue that it ought to be easier to find a “universal” model for the small scales, which tend to be more isotropic and less affected by the macroscopic flow features than the large eddies. In the LES model, the action of the large eddies in turbulent transport is simulated by spatially filtering equations obtained from the first principles. It requires a model or an approximation of the subgrid scale (SGS) effects (i.e., the interaction between the unresolved and the resolved scales). For the generality of the SGS model, one should use adequate spatial and temporal resolution. The large eddies in the turbulent flows are typically comparable in size to the characteristic length of the mean flow, while the smallest eddies are responsible for the dissipation of the turbulent kinetic energy. In large eddy simulations (LESs), the large eddies are resolved directly and the action of smallscale eddies is modeled. The rationale behind LES can be summarized as follows: 1. Momentum, mass, energy, and other passive scalars are transported mostly by large eddies. 2. Large eddies are more dependent on the physical situation under consideration. They are dictated by the geometries and the boundary conditions of the flow involved. 3. Small eddies are less dependent on the geometries, tend to be more isotropic, and are, consequently, more universal. 4. The chance of finding a universal model is much higher when only small eddies are modeled. Governing equations for LES are obtained by filtering the time-dependent Navier-Stokes equations in either Fourier space or physical space. The filtering process effectively filters out eddies whose scales are smaller than the filter width or grid spacing used in the computations. The resulting equations thus govern the dynamics of large eddies. Filtering the incompressible flow Navier-Stokes equations, we obtain
∂F ∂Fui + )0 ∂t ∂xi
( )
∂ui ∂ ∂ ∂P ∂τij ∂ µ (Fui) + (Fuiuj) ) ∂t ∂xj ∂xj ∂xj ∂xi ∂xj
(32)
(33)
2942
Ind. Eng. Chem. Res., Vol. 46, No. 10, 2007
where τij in eq 33 is the subgrid scale (eddies smaller than grid spacing) stress and is defined by
τij ) uiuj - uiuj
(34)
The major difference between the Reynolds averaged Navier Stokes equations and the above equations (eqs 32 and 33) is that the dependent variables are filtered quantities rather than the mean quantities. The subgrid scale stresses resulting from the filtering operation are unknown and require modeling. The majority of subgrid-scale models in use today are eddy viscosity models of the following form:
1 τij - τkkδij ) -2µtShij 3
(35)
where µt is subgrid-scale turbulent viscosity and Shij is the rate of shear strain tensor for the resolved scale defined by
Shij ≡
(
1 ∂uji ∂ujj + 2 ∂xj ∂xi
)
(36)
The basic subgrid-scale model was proposed by Smagorinsky21 and is commonly referred to as the Smagorinsky-Lilly model. In the Smagorinsky-Lilly model, the eddy viscosity is modeled by
νt ) Ls2|Sh|
(37)
where Ls is the mixing length for the subgrid scales and |Sh| ) x2Shij Shij. The mixing length scale is calculated as
Ls ) min(κdw, CSV1/3)
(38)
where κ ) 0.42, dw is the distance of a grid point closest to the wall, and V is the volume of the computational cell. The value of CS ) 0.1 has been found to yield the best results for a wide range of flows. 4. Results and Discussions 4.1. Effect of Velocity Boundary Conditions (VBCs) on Flow Pattern. Various velocity boundary conditions (Figure 3) from LDA measurements are imposed in the CFD simulation of the flow generated by the ultrasound. The results are compared with the experiments at only three different axial locations (z ) 0.13H, 0.53H, and 0.8H) for the sake of brevity. The three locations practically cover the entire region of the vessel from the horn to the vessel base. The comparisons were made at an equal power consumption level of 15 kW/m3. 4.1.1. Mean Axial Velocity. Radial variations in mean axial velocity at different axial locations are shown in Figure 4. From LDA measurements, it has been observed that the axial flow is in the downward direction up to 0.6R and at all axial levels. After this radial distance (r ) 0.6R), the direction of velocity changes from downward to upward at each axial location. The CFD prediction for 4BC shows a good agreement at all three locations, as shown in Figure 3E. It has been also observed that the radial variations in the mean axial velocity at different boundary conditions need not have the same trends at each axial location. The CFD prediction for 1BC (z ) 0.05H) at all the axial locations underpredicts at the radial locations defined by r/R e 0.18. A similar prediction is still seen even after increasing the boundary condition from 1BC to 2BC (z ) 0.05H and z ) 0.71H). Also, the trend for axial mean velocity seems to be in the positive direction for both 1BC and 2BC at z ) 0.8H and
Figure 4. Influence of velocity boundary conditions on the radial profile of axial mean velocity at various locations of the vessel: ()), LDA; (‚ ‚ ‚), 1BC; (- - -), 2BC; (s), 3BC; ()), 4BC.
r/R e 0.4 (Figure 3). This is evident from the contour plots of averaged velocity magnitude, where the jet streaming hits the bottom of the vessel slightly away from the center (r/R ) 0.3). Because of this, a reversible axial vorticity is found at the center of the vessel and near the vessel base, which is not the characteristic flow pattern, as evident from the LDA flow pattern. However, the strength of such a reversible flow pattern increased with an increase in the magnitudes of the values used in the boundary condition from 1BC (0.03 m/s at r/R ) 0.18) to 2BC (0.12 m/s at r/R ) 0.18). Further, with an increase in the boundary to 3BC (z ) 0.05H, z ) 0.32H, and z ) 0.71H), the mean axial velocity trend seems to be near to the LDA measurements. However, 3BC overpredicted the mean velocity values near the vessel base at z ) 0.8H and r/R e 0.2. 4.1.2. Mean Radial Velocity. Mean radial velocity decreases with an increase in the axial distance from the horn. Radial profiles of the mean radial velocity at different axial locations are shown in Figure 5. It has been observed that the mean radial velocity is very low (7% of mean axial velocity), so it can be said that the ultrasound horn reactor behaves as if it was agitated with an axial impeller or a downward plunging liquid jet where mean axial velocity is more than the mean radial velocity. Near the horn tip (z ) 0.13H), there is very little variation with the increase in the distance from the horn, i.e., boundary condition from 1BC to 4BC. Also, the predicted trend for mean radial velocity seems to be closer to the measurement for the 4BC (z ) 0.53H). Similarly, the predictions are very poor near the vessel base (z ) 0.8H). 4.1.3. Mean Tangential Velocity. The origin of the tangential component of the velocity could be due to a slightly off center of the horn and also due to a minor deviation from vertical. Reversible flow is also the cause for getting the tangential velocity stream in the vessel. In reality, no equipment can give
Ind. Eng. Chem. Res., Vol. 46, No. 10, 2007 2943
Figure 5. Influence of velocity boundary conditions on the radial profile of radial mean velocity at various locations of the vessel: ()), LDA; (‚ ‚ ‚), 1BC; (- - -), 2BC; (s), 3BC; ()), 4BC.
only radial and/or axial component velocity by completely excluding the tangentially directed velocity. Even if the equipment is able to produce perfect axial flow, the generated axial jet has to hit the domain wall to produce radial velocity that in turn forms the tangential flow. In 3D jet flows, a small initiation of the tangential component (which may arise out of nonvertical or marginally off-center horn position) will finally impart a greater deviation on the axial/radial direction velocities and, hence, the broadening of jet occurs. However, the proportion of the strength of the individual components of velocities varies from each process equipment. However, the forced boundary conditions in the present 3D CFD simulations will have severe impact in deriving the tangential velocity. Mean tangential velocity is also very low compared to the mean axial velocity and decreases with an increasing axial distance (Figure 6). It has been observed that the mean tangential velocity is very low (4% of mean axial velocity) near the horn tip (z ) 0.13H) and gradually gains its strength equal to mean radial velocity at z ) 0.53H. With 4BC, the predictions seem to be good near the horn tip (z ) 0.13H). The predictions at z ) 0.53H are 4× lower than the LDA measured mean tangential velocity. However, it can be observed that the mean tangential velocity compared well at z ) 0.8H using 4BC. Also, the trend seems to be in good agreement with the 4BC at all the axial locations of the vessel. 4.1.4. Turbulent Kinetic Energy. In this system, the turbulence is 3D and nonisotropic, and hence, the overall turbulent kinetic energy (TKE ) (V′r2 + V′z2 + V′θ2)/2) was estimated using all three fluctuating velocity components at each measurement location. Turbulent kinetic energy decreases with an increase in the axial distance from the horn tip. Radial profiles of turbulent kinetic energy at different axial locations are shown in Figure 7. Turbulent kinetic energy is higher near the horn
Figure 6. Influence of velocity boundary conditions on the radial profile of tangential mean velocity at various axial locations of the vessel: ()), LDA; (‚ ‚ ‚), 1BC; (- - -), 2BC; (s), 3BC; ()), 4BC.
Figure 7. Influence of velocity boundary conditions on the radial profile of turbulent kinetic energy at various locations of the vessel: ()), LDA; (‚ ‚ ‚), 1BC; (- - -), 2BC; (s), 3BC; ()), 4BC.
2944
Ind. Eng. Chem. Res., Vol. 46, No. 10, 2007
∫02π ∫0H ∫0R r dr dz dθ j ) 2π H R ∫0 ∫0 ∫0 r dr dz dθ
Figure 8. Influence of velocity boundary conditions on the radial profile of kinetic energy dissipation rate at various locations of the vessel: ()), LDA; (‚ ‚ ‚), 1BC; (- - -), 2BC; (s), 3BC; ()), 4BC.
and thereafter shows a sharp decrease till r ) 0.14R. From r ) 0.45R to the vessel wall, the turbulent energy values are almost constant thereafter. The modeled turbulent kinetic energy is found to be in good agreement with measured TKE using 4BC at z ) 0.13H. However, the TKE is overpredicted at the vessel center (z ) 0.53H and r/R e 0.2) and very much underpredicted away from the vessel core region (z ) 0.53H and r/R g 0.2). Similar underpredictions are evident near the vessel base (z ) 0.8H). 4.1.5. Energy Dissipation Rate. Energy dissipation rate decreases with an increase in the axial distance from the horn. Radial profiles of local energy dissipation at different axial locations are shown in Figure 8. The energy dissipation rate is maximum at the center line of the horn, and there is a sharp decrease with an increase in the radial distance till r ) 0.14R; then remains practically constant up to the wall. Most of the turbulent kinetic energy (85%) is dissipated near the horn tip, i.e., in 2% of the total volume of liquid, and the remaining (15%) is dissipated in the bulk of the liquid (98% of liquid volume). The trend in the CFD predictions is very similar for turbulent kinetic energy. The CFD prediction for all the boundary conditions is found to be in good agreement with LDA measurements near the horn tip (z/H ) 0.13), whereas both at the vessel center (z/H ) 0.53) and near the vessel base (z/H ) 0.8), the turbulent kinetic energy dissipation is underpredicted. However, there is no large variation in the predicted turbulent kinetic energy dissipation near the horn tip at z ) 0.13H with the variation in the chosen location of the boundary conditions. Knowing the local values of , the average energy dissipated in the vessel was estimated by a volume average integral of local at all locations, as described by eq 39.
(39)
The volume average integral of the predicted turbulent kinetic energy dissipation is found to be nearly half at 8.26 m2/s3 (actual power input is 15 m2/s3) for 4BC, and the values decreased drastically with the decrease in the chosen boundary conditions (4BC to 1BC), indicating that 4BC appears to be the most reasonable location for choosing the measured boundary condition input. 4.2. Comparison of Turbulence Models (By Using Velocity Boundary Conditions). Various turbulence models, namely, standard k- model, modified k- model, standard k-ω model, RSM, and LES, were used to carry out the CFD simulation for the flow generated by the ultrasound reactors. Since the earlier simulations using the 4BC standard k- model gave a realistic characteristic flow, it was decided to compare the performance of other turbulence models using 4BC. 4.2.1. Static Pressure. Figure 9 shows the radial variation of predicted static pressure with the experimentally measured pressure. For the sake of brevity, only one of the locations (23 mm below the horn) is displayed. The predicted pressure profile away from the vessel core region (r/R g 0.4) gives good agreement with LES. However, RSM, standard k-ω model, and modified k- model underpredicted almost by 25-40% away from the vessel core region. It is also important to note that the qualitative predictions of static pressure by standard k-ω model are poor. From all the turbulence models investigated, the LES prediction of static pressure is found to be nearer to the measured values. 4.2.2. Mean Axial Velocity. Radial variations in the mean axial velocity at different axial distances (z ) 0.13H, 0.53H, and 0.8H) from the horn are shown in Figure 10. Near the horn tip, the CFD predictions beyond r/R g 0.2 are good with all five different turbulence models used (z ) 0.13H). The modified k- model, standard k-ω model, and RSM underpredicted the jet streaming, which occurs in the vessel core region of the vessel (r/R e 0.2). However, the standard k- model predicted well the axial mean velocity in the vessel core region at all three locations (Figure 10). From all the turbulence models investigated, the LES prediction of mean axial velocity is found to be nearer with the measured values. Also, it is important to note that the predicted trend by both the modified k- model and the standard k-ω model was the reverse of that measured at all the axial locations in the vessel. 4.2.3. Mean Radial Velocity. Radial profiles of the mean radial velocity at different axial locations are shown in Figure 11. At all the axial locations, the RSM underpredicted the mean radial velocity. However, near the horn tip (z ) 0.13H), apart from RSM, all other turbulence models predicted in a similar manner. Also, the modified k- model overpredicts the mean radial velocity at all the axial locations of the vessel. Near the vessel base (z ) 0.8H), LES predicted closer to the experimental values of the mean radial velocity. Because of the jet streaming, there occurs a curvature near the vessel base (z ) 0.8H) and LES is able to capture the characteristic flow rather than the other turbulence models. 4.2.4. Mean Tangential Velocity. Radial profiles of the mean tangential velocity at different axial locations are shown in Figure 12. LES is able to predict the characteristic nature of the mean tangential velocity at all the axial locations of the vessel. Even though there is 30% underprediction in the vessel core region (r/R e 0.5 and z ) 0.53H), the LES is able to
Ind. Eng. Chem. Res., Vol. 46, No. 10, 2007 2945
Figure 9. Influence of turbulence models on the radial profile of static pressure at 23 mm below the horn: ()), experiment; (‚ ‚ ‚), k-ω; (- - -), modified k-; (s), RSM; ()), standard k-; (gray line), LES.
Figure 11. Influence of turbulence models on the radial profile of radial mean velocity at various locations of the vessel: ()), LDA; (‚ ‚ ‚), k-ω; (- - -), modified k-; (s), RSM; ()), standard k-; (gray line), LES.
Figure 10. Influence of turbulence models on the radial profile of axial mean velocity at various locations of the vessel: ()), LDA; (‚ ‚ ‚), k-ω; (- - -), modified k-; (s), RSM; ()), standard k-; (gray line), LES.
relatively better predict the measured trend at all the locations in the vessel. Also, the RSM predicted a reverse trend, both at z ) 0.53H and z ) 0.8H. Further, RSM severely underpredicted at z ) 0.13H. Except LES, all other turbulence models underpredicted at least 5 times at the center of the vessel (z ) 0.53H and r/R e 0.5). Even though all the predictions are based on the imposed velocity boundary conditions, LES reveals its credibility on the basis of better trend predictions. 4.2.5. Turbulent Kinetic Energy. Radial profiles of turbulent kinetic energy at different axial locations are shown in Figure 13. Near the horn tip, both the standard k- model and LES predicted well the turbulent kinetic energy of the jet streaming (z ) 0.13H and r/R e 0.2), whereas the modified k- model underpredicted 60% lower and the k-ω model overpredicted by 45%. However, both the standard k- model and LES
Figure 12. Influence of turbulence models on the radial profile of tangential mean velocity at various locations of the vessel: ()), LDA; (‚ ‚ ‚), k-ω; (- - -), modified k-; (s), RSM, ()), standard k-; (gray line), LES.
underpredicted severely beyond r/R g 0.2 and z ) 0.53H. Near the vessel base, none of the turbulence models gave a better prediction.
2946
Ind. Eng. Chem. Res., Vol. 46, No. 10, 2007
Figure 13. Influence of turbulence models on the radial profile of turbulent kinetic energy at various locations of the vessel: ()), LDA; (‚ ‚ ‚), k-ω; (- - -), modified k-; (s), RSM; ()), standard k-; (gray line), LES.
Figure 14. Influence of turbulence models on the radial profile of kinetic energy dissipation rate at various locations of the vessel: ()), LDA; (‚ ‚ ‚), k-ω; (- - -), modified k-; (s), RSM; ()), standard k-; (gray line), LES.
4.2.6. Energy Dissipation Rate. Radial profiles of local energy dissipation at different axial locations are shown in Figure 14. Energy dissipation rate is maximum at the center line of the horn, and there is a sharp decrease with an increase in the radial distance till r ) 0.14R; then remains practically constant up to the wall. However, the local energy dissipation increased near the vessel base (z ) 0.8H) because of the impingement of the streaming jet on the vessel base. Also, the standard k- model captured the energy dissipation rate quantitatively near the vessel base. Even though the LES overpredicted the jet streaming by 30% (z ) 0.13H and r/R e 0.2), the standard k- model worked well near the horn tip. At the vessel center (z ) 0.53H), none of the turbulence models are able to match the experimental results. Knowing the local values of , the average energy dissipated in the vessel was estimated by volume average integral of local at all locations, as described by eq 39. The volume average integral of the predicted turbulent kinetic energy dissipation is found to be 8.26, 5.32, 9.78, 8.05, and 13.57 m2/s3 for standard k- model, modified k- model, standard k-ω model, RSM, and LES, respectively (actual power input is 15 m2/s3), showing the ability of the solution scheme in satisfying the energy balance criteria. 4.3. Comparison of Turbulence Models (By Using Acoustic Horn-Fluid Model Interaction as Boundary Conditions). The predicted flow pattern using velocity boundary conditions and its effect on various turbulence models exclusively depends on the boundary conditions provided to the simulation, which needs to be experimental measured, and hence, such a scheme cannot be considered to be completely predictive. Hence, an attempt has been made to predict the flow pattern by providing the Ffowcs-Williams and Hawkings19 acoustic equation, which is essentially an inhomogeneous wave equation as described earlier.
4.3.1. Mean Axial Velocity. Radial variations in the mean axial velocity at different axial locations are shown in Figure 15. At all the axial locations, LES overpredicted in the region (0.3 e r/R g 0.8) where there is bulk upward axial flow. It is found that the 4BC k- model predicted reasonably well in those zones. However, RSM underpredicted by almost 94%, both in the jet streaming (z ) 0.8H and r/R e 0.5) and axial upward flow region (z ) 0.8H and r/R g 0.5). Also, the mean axial velocity trend predicted by RSM was different, both near the horn tip (z ) 0.13H) and at the center of the vessel (z ) 0.53H). The predicted trends by both the modified k- model and the standard k-ω model are very much similar to the predictions of the standard k- model. However, the predictions by the modified k- model are almost 50% higher than the predictions by the standard k- model in the bulk region of the vessel (0.2 e r/R g 0.8). 4.3.2. Mean Radial Velocity. Radial profiles of the mean radial velocity at different axial locations are shown in Figure 16. At all the axial locations, all the turbulence models overpredicted the mean radial velocity in the jet streaming zone (r/R e 0.3) by almost 75%. Also, the strength of the mean radial velocity in the bulk region (0.4 e r/R g 0.75 at z ) 0.53H) is severely underpredicted (almost by 50%) by all the turbulence models. 4.3.3. Mean Tangential Velocity. Radial profiles of the mean tangential velocity at different axial locations are shown in Figure 17. The actual generation of mean tangential velocity has a magnitude
