International Journal of Heat and Mass Transfer 91 (2015) 885–897
Contents lists available at ScienceDirect
International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt
Large eddy simulation of passive scalar transport in a stirred tank for different diffusivities Hyun Sik Yoon a,⇑, S. Balachandar b, Man Yeong Ha c a
Global Core Research Center for Ships and Offshore Plants, Pusan National University, 2, Busandaehak-ro 63beon-gil, Geumjeong-gu, Busan 609-735, Republic of Korea Department of Mechanical and Aerospace Engineering, University of Florida, MAE-A 231, P.O. Box 116250, Gainesville, FL 32611-6250, USA c School of Mechanical Engineering, Pusan National University, 2, Busandaehak-ro 63beon-gil, Geumjeong-gu, Busan 609-735, Republic of Korea b
a r t i c l e
i n f o
Article history: Received 17 April 2015 Received in revised form 8 August 2015 Accepted 10 August 2015 Available online 27 August 2015 Keywords: Passive scalar transport Molecular diffusivity Concentration Stirred tank Large eddy simulation
a b s t r a c t Large eddy simulation of flow and passive scalar transport in a stirred tank has been carried out for three different Reynolds numbers of 4000, 16,000 and 64,000 and corresponding molecular diffusivities of 1.6 104, 4.3 105 and 1.1 105, respectively, using a spectral multi-domain method. The time sequence of concentration and flow fields are explored to investigate the dependency of the concentration development on the flow structure. The sharp gradients of the instantaneous and mean concentrations are distributed around the center of larger ring vortices at high molecular diffusivity. With decreasing molecular diffusivity and increasing Reynolds number, the distribution of concentration becomes much uniform. The eddy diffusivity and subgrid scale dissipation are widely spread from the impeller jet toward outer wall and also the regions possessing small dissipation becomes diminished with increasing Reynolds number and decreasing the molecular diffusivity. Ó 2015 Elsevier Ltd. All rights reserved.
1. Introduction Transport of a scalar quantity, such as chemical concentration or temperature, is important in many engineering applications and environmental flows. A stirred tank reactor is one of the most commonly used devices in industry for achieving mixing and heat and mass transfer. Even a simple stirred tank is geometrically complex due to the presence of one or more impellers. Furthermore, under most operating conditions the flow is fully turbulent. Usually for the stirred tank, the Reynolds number is defined as Rem ¼ ND2 =m, where N, D and m are the number of the revolutions per second, the diameter of the impeller and the kinematic viscosity, respectively. The fully turbulent flow in the stirred tank is in the range of Rem > 104 . The flow is often turbulent over the entire tank; the spatio-temporal complexity of the flow is not limited to the vicinity of the impellers [1]. As a result, prediction of flow and mixing in a stirred tank is challenging. Experimental investigations before 1990 were reviewed excellently by Ranade and Joshi [2]. In recent years, significant works have been carried out to characterize complex fluid dynamics of stirred tank by using laser Doppler anemometer (LDA) [3,4] and particle image velocimetry (PIV), which is a whole field technique, ⇑ Corresponding author. Tel.: +82 51 510 3685; fax: +82 51 581 3718. E-mail address:
[email protected] (H.S. Yoon). http://dx.doi.org/10.1016/j.ijheatmasstransfer.2015.08.030 0017-9310/Ó 2015 Elsevier Ltd. All rights reserved.
to characterize instantaneous flow structures around rotating impeller blades [5–8]. Computational investigation is becoming increasingly important as a tool to numerically analyze the turbulent flow and mixing phenomena. Previous investigators [2,9–20] carried out the CFD calculation and compared their calculated results with experimental data, showing qualitative agreement between calculated and experimental results. In addition, several attempts have been made to develop computational models of gas–liquid flows in stirred vessels. As a result, the multiphase (gas–liquid) flow in the stirred tans has been handled by the CFD simulations Khopkar et al. [21] and [22], Zhang et al. [23], Jahoda et al. [24] and Zhang et al. [25]. In order to predict mixing times and explore the mixing process in stirred tanks, extensive experimental studies have been made over the last several decades [26–32]. Simultaneously, mixing studies have been carried out widely using the computational fluid dynamics (CFD). Ranade et al. [33] used numerical simulation to give detail on the flow and bulk mixing generated by a down flow-pitched blade turbine in a fully baffled cylindrical vessel. Lunden et al. [34] simulated pulse tracer experiments by solving the material balance in three-dimensional flow field with a Rushton turbine. The CFD software was used to predict the flow pattern and mixing process in a stirred tank equipped with dual-Rushton turbines by Jaworski et al. [35]. Javed et al. [36] showed the CFD predicted mixing times at different locations in the tank as well as the overall mixing time show reasonably good agreement with
886
H.S. Yoon et al. / International Journal of Heat and Mass Transfer 91 (2015) 885–897
the measured data. Hartmann et al. [37] studied a dissolution process of solid particles suspended in a turbulent flow of a Rushton turbine stirred tank by LES including passive scalar transport and particle tracking. Min and Gao [38] compared mixing time of their numerical results and experimental data. Wang et al. [39] showed that the predicted mixing time and mean residence time are in good agreement with the experimental measurements. Zadghaffari et al. [40,41] predicted mixing time and power number using LES and compared with experimental results. Consequently, based on CFD, the satisfactory results of comparisons indicate the potential usefulness of this approach as a computational tool for designing stirred reactors. One of main issues for the computational simulation for a stirred tank is how to set up the boundary conditions as the impeller region. Based on the phase-locked velocity measurements of Hill et al. [8] around the impeller region using stereoscopic PIV, recently Yoon et al. [42] have developed a detailed theoretical model for the phase-averaged impeller induced flow that accounts for the circumferential, jet and vortex components. By utilizing the theoretical model of Yoon et al. [42] for the boundary condition of the impeller swept flow, Yoon et al. [43] successfully carried out LES of flow in an unbaffled stirred tank and gave a good comparison with the previous results. Recently, Yoon et al. [44] considered scaling of flow within an unbaffled stirred tank with increasing Reynolds number using the phase-locked stereoscopic PIV measurements. The relation of parameters in the theoretical model [42] associated with the jet, tip vortex components and Reynolds number has been suggested. Recently, Yoon et al. [45] investigated the effect of Reynolds number on the flow generated by the Rushton turbine in an unbaffled tank by LES where the theoretical model of Yoon et al. [42] is adopted and the parameters appeared in this model for different Reynolds numbers are decided by Yoon et al. [44]’s results. The present LES follows author’s previous study [45] to set up the boundary conditions as the impeller region. As reviewed above, numerous researches have well established mixing time and residence time, resulting in greatly improving the understanding of these processes. However, the detailed analysis of mixing structures for different molecular diffusivities in a stirred tank has not been investigated widely. A good understanding of scaling of flow and mixing within a stirred tank reactor is of great importance. Thus, the effects of the Reynolds number and the molecular diffusivity on the flow fields and the scalar transport, respectively, not only close to impeller region but also over entire tank are observed using an accurate spectral multi-domain method. The present study focuses on the investigation of the concentration development at different molecular diffusivities in a stirred tank which is operated under turbulent conditions. The main objective of the work presented here is to study the large-scale mixing structure for different molecular diffusivities in a stirred tank by using the large eddy simulation. The time sequence of concentration and flow fields are explored to investigate the dependency of the concentration development on the flow structure. The distribution of the eddy diffusivity and the dissipation in a tank are given in the present study.
2. Governing equations and mathematical formulation 2.1. Governing equations of fluid flow Here we describe briefly the mathematical formulation for the dynamic model, which will be employed here for the simulation of stirred tank flow. For the incompressible, constant property flow considered here, the basic governing equation for the large eddy
evolution are the grid filtered Navier–Stokes and continuity equations
@ sij @ui @ui uj @p 1 @ 2 ui ¼ þ þ @xi Re @xj @xj @xj @t @xj
ð1aÞ
@ui ¼0 @xi
ð1bÞ
where overbar denotes the large (or resolved) scale flow obtained from grid-filtering. In Eqs. (1a) and (1b), xi are Cartesian coordinates, ui are the corresponding velocity components, t is the time, p is the pressure, and Re is the Reynolds number. The effect of the subgrid scale motion on the dynamics of the resolved velocity is accounted for by the subgrid scale stress tensor
sij ¼ ui uj ui ui
ð2Þ
which must be modeled in terms of the resolved scale velocity, ui , in order to obtain a closure for Eq. (1). The dynamic eddy viscosity model for the subgrid scale stress has been employed here. The detailed approach of the dynamic eddy viscosity to provide a self-consistent way and internally determines the Smagorinsky coefficient can be referred from author’s previous studies [43,45]. 2.2. Scalar transport governing equations and mathematical formulation After filtering, the equation governing convection and diffusion of the concentration of a non-reacting additive, C, can be written as,
@qj @c @ @2c þ ðuj cÞ ¼ Ds @t @xj @xj @xj @xj
ð3Þ
In Eq. (3), Ds is the scalar diffusivity defined as Ds ¼ 1=ðReScÞ, where Sc is Schmidt number which is fixed as ‘1’. The effect of the unresolved subgrid-scales appears in Eq. (3) through the subgrid flux qj . Following Germano [46], theses subgrid terms can be split into three Galilean invariant contributions.
qj ¼ Cuj Cuj
ð4Þ
By using the classical concepts of subgrid diffusivity [47], the expression of subgrid fluxes then becomes
qj ¼ DT
@C @xj
ð5Þ
The subgrid diffusivity is connected to the local grid scale D and to the local strain rate jSj through the well-known expression
DT ¼ ðC c DÞ2 jSj
ð6Þ
where C c is the dimensionless coefficient to be determined dynamically. D ¼ ðD1 ; D2 ; D3 Þ1=3 is the filter width [48] where D1 , D2 and D3 are the dimensions of the computational cell. The dynamic computation of C c is carried out using Germano ~ ~ the characteristic [46]’s approach. First of all, a test filter G ¼ GG, ~ ¼ 2D, is introduced. By refiltering Eq. (3) with width of which is D ~ the filter G, one obtains the LES governing equations at the test ~ . This equation involve new subgrid term, filter scale D
g ~ Q j ¼ Cuj C u~j
ð7Þ
~ the difference between the subgrid fluxes After filtering qj by G, ~ and D obeys Germano [46]’s obtained at the two cutoff levels D algebraic identity:
H.S. Yoon et al. / International Journal of Heat and Mass Transfer 91 (2015) 885–897
g ~ Lcj ¼ Q j q~j ¼ Cuj C u~j
ð8Þ
~ ~ and C ~, u Assuming that Q j can be expressed in terms of D in the same manner as qj is expressed in terms of D, uj and C, using ~ with the Smagorinsky model in Eqs. (5) and (6) at the test-scale D the same coefficient C c , Lcj can be written as follow
Lcj ¼ ðC c DÞMcj
where Mcj ¼
~ D D
!2
g ~ @C ~ @C jSj jSj @xj @xj
ð9Þ
The unknown coefficient C c are computed through the relation
C 2c ¼
1 hLcj Mcj i D2 hMcj Mcj i
ð10Þ
where hi denotes the averaging over the entire computational volume. To prevent numerical instability caused by negative value of C c , the numerator and denominator are averaged over the entire computational volume. This dynamic procedure has been used, for example, by Akselvoll and Moin [49] to study turbulent mixing of fuel and air by tracking a passive scalar without the effects of swirl, chemical reactions and heat for the geometry like a coaxial jet-combustor, and by Ribault et al. [50] to investigate the passive scalar development in a planar jet and show the evolution of the mean scalar and turbulent scalar intensity well representing those of experiments and direct numerical simulation. 3. Numerical methodology The geometry to be considered is a cylindrical vessel of diameter T ¼ 150 mm filled with water at room temperature to a height equal to the diameter as shown in Fig. 1. A simple six-blade Rushton impeller is placed at the center of the tank, halfway between the top and bottom of the fluid. The diameter of the impeller is D ¼ T=3 ¼ 50 mm and the height of the impeller is 0:2D. For simplicity, no baffles are considered in the tank and top surface is assumed to be flat. The computational geometry is greatly simplified in a frame of reference rotating with the blades. In this frame of reference, the geometry has a six-fold symmetry along the circumferential direction. The present computations assume the geometric periodicity to extend and apply to the flow field as well. Thus the velocity and pressure fields are assumed to be periodic over a 60 sector and therefore the computational geometry is confined to this sector. Here we resort to this assumption to maintain the computational cost to the LES manageable. The governing large eddy equation (Eq. (1)) is appropriately modified to suit non-inertial frame of the reference. Owing to periodicity along the h direction, a Fourier expansion with a corresponding uniform grid is used along the circumferential direction. A spectral multi-domain technique is used along the r z plane and typical domain decomposition is shown in Fig. 1(b). No-slip boundary condition is enforced on the top, bottom and outer boundaries of the tank and along the impeller axis the fluid is assumed to rotate with the axis. These boundary conditions are correspondingly adjusted for the rotation of the frame of the reference. The numerical methodology is versatile that momentum sources can be placed along the internal domain boundaries and adjusted to obtain any desired velocity field across this boundary. Here the outer flow section of the impeller swept volume, indicated by the line marked AB on the r z plane in Fig. 1(b), is chosen where the impeller-induced flow (uimp ) is applied as an internal boundary condition. Pressure changes across this interface representing a local momentum source. Thus this approach is similar in spirit to those employed by Eggles [17] and Revstedt et al. [18].
887
Based on the experimental measurements a theoretical description of the impeller-induced flow uimp has been developed [42]. Preliminary results suggests that the impeller-induced flow can be described as a superposition of a circumferential flow, a circular jet and a pair of tip vortices associated with each impeller blade (see Fig. 1(a)). The impeller-induced flow is strongly threedimensional and as a result a complex theoretical description is necessary to capture all the essential features. The parameters in theoretical model used to represent the impeller-induced flow at the three different Reynolds numbers (Rem ¼ 4000, 16,000 and 64,000) are determined using the methods described in Yoon et al. [42,44]. Eventually, the oscillatory inflow model [43] is adopted to the outer flow section of the impeller swept volume, indicated by the line marked AB on the r z plane in Fig. 1(b) as an internal boundary condition. By utilizing the theoretical model of Yoon et al. [42] for the boundary condition of the impeller swept flow, Yoon et al. [43,45] successfully carried out LES of flow in an unbaffled stirred tank and gave a good comparison with the previous results. The more detailed explanation for the impellerinduced flow uimp can be referred Yoon et al. [45]. In the computations the tip radius of the impeller ðD=2Þ and the blade tip velocity ðpNDÞ are chosen as the length and velocity scales, where N is the number of the revolutions per second. The time and pressure scales are accordingly defined as ð1=2pNÞ and ðq½pND2 Þ. Following the tradition used in the chemical engineering literature the Reynolds number of the flow is defined
as Rem ¼ ND2 =m ¼ 2Re=p. In this study, three different Reynolds number of Rem ¼ 4000, Rem ¼ 16; 000 and Rem ¼ 64; 000 are considered. The scalar concentration is set to ‘1’ along the curved surface with the line marked CD on the r z plane in Fig. 1(b). The Ds in Eq. (3) governing convection and diffusion of the concentration of a non-reacting additive denotes 1=ReSc, where the Schmidt number Sc is ‘1’ as above mentioned. Thus, the diffusivity Ds is only dependent on the Reynolds number. In this study, three different Reynolds numbers are considered as above mentioned. The corresponding values of Ds are Ds ¼ 1:6 104 , Ds: ¼ 4:3 105 and Ds ¼ 1:1 105 . The time step considered in this study is 0.0001 which is nondimensionalized by the blade tip velocity and the tip radius. Therefore, the impeller advances about 0:006 along the circumferential direction in each time step. Typical computation involves 32 grid points along the circumferential direction and 72 subdomains in the r z plane, with each subdomain resolved by up to 13 13 points. The grid system in the r z plane is plotted in Fig. 1(c). The spectral multi-domain technique employs Chebyshev expansion within each subdomain and the grid points are accordingly the Gauss–Lobatto points and are therefore nonuniformly distributed [51]. Even at the present Reynolds numbers, grid resolves only a range of large flow scales. A large eddy simulation with appropriate subgrid model is thus required to account for the unresolved small scale eddies. The spectral methodology provides an accurate representation for the range of resolved scales. Large eddy simulations at other grid resolutions have also been performed yielding similar results. The numerical methodology employed for advancing the flow in time is standard. Operator (or time) split scheme is employed with the advection and diffusion effects taken into account first to obtain an intermediate velocity. The advection terms are treated explicitly using a third order Adams–Bashforth scheme and the viscous terms are treated semi-implicitly using the Crank–Nicolson scheme. Thus solution for each component of the intermediate velocity is decoupled and reduces to solving a three-dimensional Helmholtz equation in cylindrical coordinates. The corresponding pressure field is obtained by solving the pressure Poisson equation,
888
H.S. Yoon et al. / International Journal of Heat and Mass Transfer 91 (2015) 885–897
z (uz)
(uθ)
Impeller Swept Volume
θ
Inflow (into computational domain due to impeller)
Return Flow
T
(ur) r
T/2
Radial Jet
D=T/3
Return Flow
T
+
Computational Domain
Tip vortex pairs
r
(a)
θ
(b)
(c)
Fig. 1. (a) Schematic of the stirred tank with a typical six blade Rushton impeller. The plan view shown at the bottom on the right is viewed up from under the tank. (b) The overall computational geometry divided into sub-domains in the r z plane. The gray shaded area corresponds to the area swept by the blade. The impeller-induced flow is applied as inflow along AB. The scalar concentration is set to ‘1’ along the curved surface with height CD. (c) Grid system of 72 subdomains in the r z plane.
Table 1 Validation of mixing times. Impeller Reynolds number
Shekhar and Jayanti [52] Standard k e model
Low Re k e model
960 1920 4800
36.0 23.0 13.0
88.0 36.0 19.0
Nagata correlation [53]
Present results
73.0 36.5 14.6
86.0 35.0 18.0
with the divergence of the intermediate velocity as the forcing term. The intermediate velocity is projected back onto the divergence-free space with pressure correction step. The overall scheme is second-order accurate in time. The three-dimensional Helmholtz and Poisson equations are Fourier transformed in h to
reduce to a set of two-dimensional Helmholtz equations for each circumferential wavenumber. These two-dimensional Helmholtz equations in the r z plane are solved efficiently using patching technique, which enforces continuity of the function and the first derivative (C 1 continuity) at the sub-domain interfaces. In the dynamic LES, the test filter along the circumferential direction is chosen to be a sharp cutoff filter, which retains only half the total number of circumferential modes (8 modes for the case of 32 circumferential grid points). Along the radial (r) and axial (z) directions a box-average filter, that essentially projects the flow onto half as many points along these directions, is applied. Thus, for the present nonuniform grid, even though D ~ vary over the r z plane, their ratio (D ~ =D) in Eq. (9) can be and D taken to be 2. For the range of modest Reynolds numbers to be
H.S. Yoon et al. / International Journal of Heat and Mass Transfer 91 (2015) 885–897
3 0.01 0.02
2
0.03
0.1
0.04 0.05
1
0.08 0.06 0.09
z
0.07 0.1
0
889
considered here, the above grid along with Gauss–Lobatto clustering of points provides adequate resolution near the boundaries. Thus, complex issues pertaining to wall modeling in the context of LES are avoided in the present work. Eventually, by using the present boundary conditions, the computational domain, the number and distribution of the grid, and the numerical methods, Yoon et al. [43,45] successfully carried out LES of flow in an unbaffled stirred tank and gave a good comparison with the previous results. 4. Validation
-1 0.1
-2 -3 0
(a)
1
2
r
3
3 0.01 0.02 0.03
2
0.1 0.08
0.04 0.05
z
1
0.06 0.1
0
-1
0.05
0.1 0.06
Regarding to the fluid flow, the validation of the present results has been successfully conducted by Yoon et al. [43,45] which gave a good comparison with the previous results. The present study investigates the mixing in the unbaffled stirred tank, needing to validate the scalar mixing characteristics. Thus, we followed Shekhar and Jayanti [52]’s conditions. They considered the unbaffled vessel which is the same with the present study, even the tank equips a mechanically stirred eight-blade paddle impeller which is different type impeller with the present one. However, the present numerical method can cover the Shekhar and Jayanti [52]’s conditions such as the different impeller. In this validation study, we consider higher three Reynolds numbers than various Reynolds numbers considered in Shekhar and Jayanti [52]. In addition, the empirical correlations of Nagata [53] have been utilized to obtain the mixing times for three Reynolds numbers to validate the present results. As shown in Table 1, Shekhar and Jayanti [52] used two different turbulence models and the results obtained with the low Reynolds k e model matched very well with those of Nagata [53] while the results from the standard k–e model were much lower. Consequently, the present results are consistent with the findings of Shekhar and Jayanti [52] and Nagata [53], as shown in Table 1.
-2 5. Results and discussion
-3 0
(b)
1
r
2
3
3 0.01
2
0.02 0.03 0.04 0.1 0.09
1
z
0.1
0 0.02 0.03 0.04
-1
0.05
-2
0.1
0.06 0.07 0.08 0.09
-3
(c)
0
1
r
2
3
Fig. 2. Instantaneous flow and concentration fields in r z plane at h ¼ 30 for different Reynolds numbers and corresponding diffusivities; (a) Rem = 4000 and DS = 1.6 104, (b) Rem = 16,000 and DS = 4.3 105, (c) Rem = 64,000 and DS = 1.1 105. Here, the contour levels of the concentration are in between 0.01 and 0.1 at intervals of 0.01.
Fig. 2(a)–(c) show the instantaneous flow and concentration fields in r z plane at h ¼ 30 for different Reynolds numbers and corresponding diffusivities. At DS = 1.6 104 and Rem = 4000 in Fig. 2(a), the scalar concentration is convected along the path of stronger circulation of upper and lower primary ring vortices below the area of radial flow at the level of impeller blades. This pair of primary ring vortices preserves their shape at every instance over the time. But, the location of centers of these cells moves up and down owing to the flapping of radial jet component of the flow originated from the impeller. Yoon et al. [43] showed the presence of two modes of jet instability (sinuous and varicose instability) originated from the jet component of impeller induced flow. Thus, this pair of primary ring vortices contributes to higher concentration to distribute along the vicinity of the outer wall. At this Reynolds number of 4000, the inertial forces are not enough to propel the flow along the outer wall to reach the top and bottom surfaces of tank. Thus, the regions of secondary flow are built up at the top and bottom right hand corners. The tertiary recirculation zones near the shaft are formed by the incoming flow to the impeller separated near the inner corners. Resultingly, the lowest concentration occurs around shaft. As Reynolds number increases, the regions of secondary flow at the top and bottom right hand corners become smaller, because the inertial forces to accelerate flow along the outer wall are larger. Also, the radial jet from the impeller becomes stronger to push the two ring vortices near the outer wall with increasing Rem : Especially, at Rem ¼ 64; 000, each ring vortex at the upper and lower parts of vessel at lower Reynolds numbers in Fig. 2(a) and (b) is
H.S. Yoon et al. / International Journal of Heat and Mass Transfer 91 (2015) 885–897
3
3
2
2
1
1
0
0
z
z
890
-1
-1
-2
-2 -3
-3 0
1
r
2
0
3
1
r
2
(c) t=90
3
3
2
2
1
1
0
0
z
z
(a) t=30
3
-1
-1
-2
-2 -3
-3 0
1
r
2
3
(b) t=60
0
1
r
2
3
(d) t=120
Fig. 3. Sequence of 4 frames, each separated by 30 nondimensional time unit, showing the h-averaged instantaneous flow at Rem = 4000 and corresponding concentration fields at DS = 1.6 104. Contour levels: 0.01–0.1 at steps of 0.01.
broken up into smaller vortices as depicted in Fig. 2(c). The secondary recirculation zones at the top and bottom right hand corners, and the tertiary recirculation zones near the shaft at lower Reynolds numbers are also suppressed by the small vortices at this high Reynolds number. This can be certified by comparing between Fig. 2(a)–(c) for Rem ¼ 4000, 16,000 and 64,000, respectively. Thus, the two main circulation loops with inner small vortices extend over the entire upper and lower halves of tank. Consequently, the variation of the concentration distribution along the diffusivity follows this change of the flow structures according to Rem . With decreasing DS and increasing Rem , the small scale fluctuations of concentration dominate close to the outer wall, because of increasing small scale fluctuation of flow along Rem within the outer wall boundary layer. The thickness of boundary and mixing layer of concentration becomes thinner with decreasing DS and increasing Reynolds number. The contours of concentration along the impeller jet stream are similar to the profiles of jet flow induced by impeller as shown well in Fig. 2(c). The
impeller induced jet entrains surrounding fluid whose concentration is low and transfer it closer to the outer wall with increasing Reynolds number and decreasing DS, compared to that at DS = 1.6 104 and Rem ¼ 4000. Fig. 3 shows time sequence of 4 frames for the h-averaged instantaneous flow for Rem = 4000 and corresponding concentration fields at DS = 1.6 104. The time interval for each frame is t = 30 corresponding to about 4.8 rotations of the impeller. Here, the contour levels of concentration are the same as used in Fig. 2. In order to give a better view, the number of the grids used to plot the velocity vectors in Fig. 3 is reduced, comparing to Fig. 2. After passing t = 30 from the start of solving the passive scalar equation, the concentration is transferred and convected along the path of upper and lower primary circulations as shown in Fig. 3(a). The dense distribution of concentration contours exists along the shear layers between ring vortices and outer wall and between ring vortices and corner vortices in both upper and lower right corners. The fluids around the center of the ring vortices
H.S. Yoon et al. / International Journal of Heat and Mass Transfer 91 (2015) 885–897
3
2
2
1
1
0
0
z
z
3
891
-1
-1
-2
-2 -3
-3 0
1
r
2
0
3
1
r
2
(a) t=30
(c) t=90
3
3
2
2
1
1
0
0
z
z
3
-1
-1
-2
-2 -3
-3 0
1
r
2
3
(b) t=60
0
1
r
2
3
(d) t=120
Fig. 4. Sequence of 4 frames, each separated by 30 nondimensional time unit, showing the h-averaged instantaneous flow at Rem = 16,000 and corresponding concentration fields at DS = 4.3 105. Contour levels: 0.01–0.1 at steps of 0.01.
remain unmixed. Also, we can observe the unmixed region in the region of impeller induced jet flow along the mid-plane. At t = 60 as shown in Fig. 3(b), the h-averaged instantaneous flow field is very similar with that at t = 30 in Fig. 3(a), except the direction of impeller induced jet oscillation close to the outer wall. Thus, the low concentration contours along the mid-plane lean downward to the same direction as oscillating jet. The region containing the concentration over C = 0.1 is moving in the circulation direction and is extended inward to the tank. However, the unmixed fluids still remain around the center of the ring vortices as shown in Fig. 3(b). Fig. 3(c) shows the h-averaged instantaneous concentration field at t = 90. The low concentration reaches close to the impeller shaft and the unmixed region around the center of ring vortices filled with relatively low concentration fluid as shown in Fig. 3 (c). The mixing layer is still observed along the impeller stream. We can also observe the development of new mixing layer on the upper and lower disc.
At t = 120 in Fig. 3(d), the distribution of h-averaged instantaneous concentration from r 1:4 to impeller shaft shows monotonically decreasing nature except the small mixing layer along the mid-plane. Also, the dense contours along the shear layers between ring vortices and outer wall and between ring vortices and corner vortices in both upper and lower right corners shown at an earlier time are almost disappeared. Fig. 4 shows the h-averaged instantaneous flow at Rem = 16,000 and corresponding concentration fields at DS = 4.3 105 at the same time as Fig. 3. As shown in Fig. 4, the size of secondary vortices at the upper and lower right corners at Rem = 16,000 is smaller than that at Rem = 4000 in Fig. 4, since the inertial force along the outer wall becomes larger with increasing Reynolds number. At t = 30, the contours of concentration move mainly upward and downward along the outer wall and are slightly distorted along the direction of the movement of ring vortices as shown in Fig. 4(a). At t = 60, Fig. 4(b) shows that the low concentration fluids are distributed along the mid-plane due to the effect
H.S. Yoon et al. / International Journal of Heat and Mass Transfer 91 (2015) 885–897
3
3
2
2
1
1
0
0
z
z
892
-1
-1
-2
-2
-3
-3 0
1
r
2
3
0
1
r
2
(c) t=90
3
3
2
2
1
1
0
0
z
z
(a) t=30
-1
-1
-2
-2
-3
3
-3 0
1
r
2
3
(b) t=60
0
1
r
2
3
(d) t=120
Fig. 5. Sequence of 4 frames, each separated by 30 nondimensional time unit, showing the h-averaged instantaneous flow field at Rem = 64,000 and corresponding concentration field at DS = 1.1 105. Contour levels: 0.01–0.1 at steps of 0.01.
of jet induced by impeller and around the center of ring vortices, which represent the inactive regions. The low concentration contours along the mid-plane lean to the same direction of jet flapping at t = 90 as shown in Fig. 4(c). The pattern of concentration contours at DS = 1.6 104 is similar to that at DS = 4.3 105 except the regions of top and bottom right corners to each other. This is due to the similar flow pattern of two cases. Thus, as shown in Fig. 4(d) at t = 120, the locally high concentration gradients still exist around the center of ring vortices. Fig. 5 shows the h -averaged instantaneous flow and corresponding concentration fields at Rem ¼ 64; 000 and DS = 1.1 105 at the same time as Figs. 3 and 4. Because the angle of jet oscillation to the upward and downward at Rem ¼ 64; 000 is much lager than that at lower Reynolds numbers, the large ring vortices observed in lower Reynolds numbers are separated into several smaller vortices as also shown in Fig. 2(c). The axial flow at the upper and lower half part of tank moves along the outer wall and reaches the top and bottom wall. Due to this flow characteristics
at Rem ¼ 64; 000, we can observe the high concentration gradients are along the jet stream. The unmixed region which exists near the center of ring vortices at higher two DS s, disappears at DS = 1.1 105. As shown in Fig. 5(d) at t = 120, the distribution of h-averaged instantaneous concentration at DS = 1.1 105 is much uniform at the outer wall than that at higher two DS s. Fig. 6 shows the time and h-averaged concentration fields at DS = 1.6 104, DS = 4.3 105 and DS = 1.1 105 in r z plane. The results have been symmetrized about the mid-plane (z = 0) and only the top half of the tank is plotted. At DS = 1.6 104, the contours of the concentration distort along the path of stronger circulation of primary ring vortices which preserves their shape at every instance over the time. Resultingly, unmixed region with low concentration forms around the center of the ring vortices as shown in Fig. 6(a). At DS = 4.3 105, this unmixed region becomes smaller and the inward distortion of the concentration contours becomes weaker than the case of DS = 1.6 104, as shown in Fig. 6(b).
893
H.S. Yoon et al. / International Journal of Heat and Mass Transfer 91 (2015) 885–897
0.18
Rem= 4000 Rem= 16000 Rem= 64000
0.17
Cs
0.16 0.15 0.14 0.13 0.12 0.11
0
20
40
60
80
100
120
Dimensionless time, t (a)
(a) 0.28
D S= 1.6E-4 D S= 4.3E-5 D S= 1.1E-5
0.26 0.24
CC
0.22 0.2 0.18 0.16 0.14 0
(b)
20
40
60
80
100
120
Dimensionless time, t (b) Fig. 7. Time histories of the Smagorinsky coefficients of (a) C S for momentum equation for three different Reynolds numbers and (b) C C for the passive scalar equation for three different diffusion coefficients.
(c) Fig. 6. Time and h-averaged concentration fields at (a) DS = 1.6 104, (b) DS = 4.3 105 and (c) DS = 1.1 105 in r z plane.
With decreasing DS , the presence of lower concentration around the center of the ring vortices is suppressed and eventually disappeared at DS = 1.1 105, as shown in Fig. 6(c). This results
from that as the Reynolds number increases, the ring vortex becomes slightly more oblate and increases in height until it reaches the lid of the tank. The distribution of time and h-averaged concentration becomes more uniform. However, near the mid-plane, the gradient of the concentration becomes sharper because the radial jet is strongly concentrated along the center of tank and fluids are sucked into the jet. Fig. 7 shows the time histories of the Smagorinsky coefficients of CS for momentum equation for three different Reynolds numbers and CC for the passive scalar equation for three different diffusion coefficients. These values of CC are calculated by the present dynamic large eddy simulation. The definitions of CC and CS are refereed in Eq. (10) in Yoon et al. [45], respectively. Fig. 7(a) shows the time histories of the Smagorinsky coefficient CS obtained from the dynamic large eddy simulation for Rem ¼ 4000, Rem ¼ 16; 000 and Rem ¼ 64; 000. Significant time variation can be observed in Fig. 7(a) and the time-average values in the range from CS = 0.125 to 0.13 for three Reynolds numbers compares well with the value of CS = 0.1 normally used in the Smagorinsky model.
894
H.S. Yoon et al. / International Journal of Heat and Mass Transfer 91 (2015) 885–897
The instantaneous Smagorinsky coefficients, CC, vary significantly along the time with the time averaged values of 0.167, 0.17 and 0.185 at DS = 1.6 104, DS = 4.3 105 and DS = 1.1 105, respectively, as shown in Fig. 7(b). According to Horiuti [54]’s results, CC is slightly higher than CS, and CC∕CS 1.4. In this study, the ratio of CC to CS is between 1.336 and 1.423, which compares well with the value of 1.4 given by Horiuti [54]. The time histories of volume-averaged eddy diffusivities at different diffusion coefficients are shown in Fig. 8. The eddy diffusivity, DT is defined in Eq. (6). The values of time and volumeaveraged eddy diffusivities are 0.00181, 0.00192 and 0.002 at DS = 1.6 104, DS = 4.3 105 and DS = 1.1 105, respectively. The eddy diffusivities are about 10, 40 and 180 times larger than nondimensional molecular diffusivity (1=ReSc) at DS = 1.6 104, DS = 4.3 105 and DS = 1.1 105, respectively. Fig. 9 shows the distributions of time and h-averaged eddy diffusivities DS = 1.6 104, DS = 4.3 105 and DS = 1.1 105 over the entire r z plane. The time and h-averaged eddy diffusivity is focused along the shear layer originated from the jet near the mid-plane of the tank. The overall distribution of time and
0.006
D S= 1.6E-4 D S= 4.3E-5 D S= 1.1E-5
0.005
0.003 0.002 0.001 0
0
20
40
60
80
100
120
Dimensionless time, t Fig. 8. The volume-averaged subgrid scale eddy diffusivity for three different diffusion coefficients.
3
6 4
3
1
2
2
3
4
2
1
2 1 4
5
9 6
0
-1
9 8 5 7 6 5 4 3 2 1
0.0070 0.0055 0.0045 0.0035 0.0025 0.0015 0.0010 0.0008 0.0006
4 3
2
1
1
0
2
5
4
6
7
2
-1 3
-2
-2
3
1
-3
-3
(a)
3
3
2
2
z
6 4
z
< D T>
0.004
0
1
2
r
3
3 6
(b)
0
1
r
2
3
5 4
2
2
z
5
9 7
0
9 8 5 7 4 6 5 4 3 4 2 1
2
1
1
4
3
1
-1
2 2
-2
4
3
5
-3
(c)
0.0070 0.0055 0.0045 0.0035 0.0025 0.0015 0.0010 0.0008 0.0006
0
1
r
2
3
Fig. 9. Time and h-averaged subgrid scale eddy diffusivity at (a) DS = 1.6 104, (b) DS = 4.3 105 and (c) DS = 1.1 105.
895
H.S. Yoon et al. / International Journal of Heat and Mass Transfer 91 (2015) 885–897
3
3
4 3
2
2 3
3
2
5
1
1
1
5
-1
0 1
-1
4
5
7
2
5
-3
(a)
3
2
4
1
0
1
3
2
r 5
4
5
-3
3
0
3
3
2
-2
3
3
1
1 -2
4
2
3
4
z
7
0
3
4
2 1
1
z
5 4
1
2
2
4 1
r
2
3
(b)
3
2 2
1 1
3
1
7 6 5 4 3 2 1
z
4 3
7 6 5
0
3 -1
1 1 2
2
-2 5
-3 0
4
0.0200 0.0100 0.0050 0.0020 0.0006 0.0002 0.0001
3
1
r
2
3
(c)
Fig. 10. The distribution of time and h-averaged subgrid scale dissipation over the entire r z plane for (a) Rem = 4000, (b) Rem = 16,000 and (c) Rem = 64,000.
h-averaged eddy diffusivity at different DS is similar for all values. Because the ring vortices become smaller and convective motion becomes stronger with increasing Rem , the time and h-averaged eddy diffusivity is widely spread toward outer wall with decreasing DS . Dissipation is a quantity of importance in a stirred tank since it provides a direct measure of energy input into the system at the impeller. The efficiency of mixing can then be gauged by evaluating the degree of mixing in relation to total dissipation. In the context of LES, in addition to the resolved scale dissipation, which accounts for the direct loss of energy from the resolved scales to dissipation, there is also subgrid scaled dissipation, which accounts for the net transfer of energy from the resolved scales to the subgrid scale and eventually to heat. Fig. 10 show the distribution of time and h-averaged subgrid scale dissipation over the entire r z plane. The dissipation is focused along the shear layers associated with
the jet near the mid-plane of the tank. The enhanced strain-field associated with the stagnation point also contributes to dissipation. As the Reynolds number increases, the time and h-averaged subgrid scale dissipation spreads widely from the impeller jet and also the regions possessing small dissipation becomes diminished. In general, the peak dissipation occurs along the mid-plane, thus differently influencing the downstream evolution of jet, regardless of Reynolds number. This region contained high energy dissipation rate which correlate quicker mixing. However, the magnitude of the peak dissipation is larger and the distribution is denser with increasing Reynolds number, which can be certified by comparing Fig. 11(a)–(c) for Rem ¼ 4000, Rem ¼ 16; 000 and Rem ¼ 64; 000, respectively. The distribution of resolved scale dissipation looks qualitatively similar, but it of a smaller magnitude, and therefore not shown in here.
896
H.S. Yoon et al. / International Journal of Heat and Mass Transfer 91 (2015) 885–897
0.1
2 3
2
0
z
6 5 4 3 2 1
4
0.012 0.010 0.008 0.006 0.004 0.002
Conflict of interest None declared.
Acknowledgments
-0.1 3
1.2
1.4
r
1.8
2
3
0.1 0
3
4
5
3
2
-0.1
1.2
z
-0.1
1.4
1.6
r
2
1.8
2
4
0
3
6
5
6 5 4 3 2 1
0.012 0.010 0.008 0.006 0.004 0.002
2
3
2
5
0.1
(a)
This work was supported by a National Research Foundation of Korea (NRF) grant funded by the Korea government (MSIP) through GCRC-SOP (No. 2011-0030013).
References
5
z
1.6
4
3
6 5 4 3 2 1
(b) 0.012 0.010 0.008 0.006 0.004 0.002
4
1.2
1.4
r
1.6
1.8
(c)
Fig. 11. The distribution of time and h-averaged subgrid scale dissipation in the region close to the blade tip for (a) Rem = 4000, (b) Rem = 16,000 and (c) Rem = 64,000.
6. Conclusions Large eddy simulation of flow and concentration fields in an unbaffled stirred tank was carried out for three different Reynolds numbers of 4000, 16,000 and 64,000, and corresponding molecular diffusivities of 1.6 10–4, 4.3 10–5 and 1.1 10–5 using the accurate multi-domain spectral method. The steep gradients of the instantaneous and mean concentrations are distributed around the center of larger ring vortices at two higher molecular diffusivities. The impeller flow plays a key role in transporting the unmixed fluid flow to the mixed fluid flow region. The mixing is diffusion-dominated as in the case where there are low-velocity regions close to the center of large ring vortices and impeller shaft. It is advection-dominated as in the case where there are the series of small vortices along the outer wall. With decreasing molecular diffusivity and increasing Reynolds number, the distribution of concentration from the releasing position becomes much uniform. The present result also shows that the mixing structure posses strong diffusive properties at the lowest turbine agitation speed and strong convective properties (and hence low diffusive characteristics) at the highest agitation speed. The time and h-averaged eddy diffusivity and subgrid scale dissipation is widely spread from the impeller jet toward outer wall and also the regions possessing small dissipation becomes diminished with increasing Reynolds number and decreasing the molecular diffusivity.
[1] M.A. Rao, R.S. Brodkey, Continuous flow stirred tank turbulence parameters in the impeller stream, Chem. Eng. Sci. 27 (1972) 137–156. [2] V.V. Ranade, J.B. Joshi, Flow generated by a disc turbine, Trans. IChemE Part A 68 (1990) 34–50. [3] M. Schafer, M. Hofken, F. Durst, Detailed LDV measurements for visualization of the flow field within a stirred tank reactor equipped with a Rushton turbine, Chem. Eng. Res. Des. (1997) 729–736. [4] K.C. Lee, M. Yianneskis, Turbulent properties of the impeller stream of a Rushton turbine, AIChE J. 44 (1998) 13–24. [5] K.V. Sharp, K.S. Kim, R.J. Adrian, A study of vorticity and dissipation around a Rushton turbine using particle image velocimetry, in: Proc. 9th Int’l Symp. Applications Lasers to Fluid Mechanics, Lisbon, July 13–16, 1998, pp. 14.1.1-10. [6] N.G. Deen, B.H. Hjertager, Multiphase particle image velocimetry measurements in an aerated stirred tank, in: AIChe Meeting Dallas, November 1999, Session 06005. [7] V.V. Ranade, M. Perrard, N. Le Sauze, C. Xuereb, J. Bertrand, Trailing vortices of Rushton turbine: PIV measurements and CFD simulations with snapshot approach, Trans. IChemE Part A 79 (2011) 3–12. [8] D.F. Hill, K.V. Sharp, R.J. Adrian, Stereoscopic particle image velocimetry measurements of the flow around a Rushton turbine, Exp. Fluids 29 (2000) 478–485. [9] B.J. Hutchings, R.J. Weetman, B.R. Patel, Computation of flow fields in mixing tanks with experimental verification, in: paper TN-481, ASME Meeting, San Francisco, 1989. [10] V.V. Ranade, J.B. Joshi, A.G. Marathe, Flow generated by pitched blade turbines II: simulation using k–e model, Chem. Eng. Commun. 81 (1989) 225–248. [11] S.M. Kresta, P.E. Wood, Prediction of the three-dimensional turbulent flow in stirred tanks, AIChE. J. 37 (1991) 448–460. [12] A. Bakker, H.E.A. van den Akker, Single-phase flow in stirred reactors, Trans. IChemE 74 (1994) 583–593. [13] L. Dong, S.T. Johansen, T.A. Engh, Flow induced by an impeller in an unbaffled tank – I. Experimental, Chem. Eng. Sci. 49 (1994) 549–560. [14] L. Dong, S.T. Johansen, T.A. Engh, Flow induced by an impeller in an unbaffled tank – II. Numerical modeling, Chem. Eng. Sci. 49 (1994) 3511–3518. [15] R. Verzicco, G. Iaccarino, M. Fatica, P. Orlandi, Flow in an impeller-stirred tank using an immersed-boundary technique, AIChE J. 50 (2004) 1109–1118. [16] R.M. Jones, A.D. Harvey, S. Acharya, Two-equation turbulence modeling for impeller stirred tanks, J. Fluids Eng. 123 (2001) 640–648. [17] J.M.G. Eggels, Direct and large eddy simulation of turbulent fluid flow using the lattice-Boltzmann scheme, Int. J. Heat Fluid Flow 17 (1996) 307–323. [18] J. Revstedt, L. Fuchs, C. Tragardh, Large eddy simulation of the turbulent flow in a stirred reactor, Chem. Eng. Sci. 53 (1998) 4041–4053. [19] K.A. Pericleous, M.K. Patel, The source-sink approach in the modeling of stirred reactors, PhysioChem. Hydrodyn. 9 (1987) 279–297. [20] J. Darksen, H.E.A. Van den Akker, Large eddy simulation on the flow driven by a Rushton turbine, AIChE J. 45 (1999) 209–221. [21] A.R. Khopkar, J. Aubin, C. Xureb, N. Le Sauze, J. Bertrand, V.V. Ranade, Gas– liquid flow generated by a pitched blade turbine: PIV measurements and CFD simulations, Ind. Eng. Chem. Res. 42 (2003) 5318–5332. [22] A.R. Khopkar, A.R. Rammohan, V.V. Ranade, M.P. Dudukovic, Gas–liquid flow generated by a Rushton turbine in stirred vessel: CARPT/CT measurements and CFD simulations, Chem. Eng. Sci. 60 (2005) 2215–2229. [23] Q. Zhang, Y. Yong, Z.S. Mao, C. Yang, C. Zhao, Experimental determination and numerical simulation of mixing time in a gas–liquid stirred tank, Chem. Eng. Sci. 64 (2009) 2926–2933. [24] M. Jahoda, L. Tomášková, M. Mošteˇk, CFD prediction of liquid homogenisation in a gas–liquid stirred tank, Chem. Eng. Res. Des. 87 (2009) 460–467. [25] Y. Zhang, C. Yang, Z.S. Mao, Large eddy simulation of the gas–liquid flow in a stirred tank, AIChE J. 54 (2008) 1963–1974. [26] R.D. Biggs, Mixing rates in stirred tanks, AIChE J. 9 (1963) 636–640. [27] S.J. Kang, O. Levenspiel, New scale-up and design method for stirrer agitated batch mixing vessels, Chem. Eng. Sci. 31 (1976) 569–577. [28] I. Houcine, H. Vivier, E. Plasari, R. David, J. Villermaux, Planar laser induced fluorescence technique for measurements of concentration fields in continuous stirred tank reactors, Exp. Fluids 22 (1996) 95–102. [29] R. Mezaki, M. Mochizuki, K. Ogawa, Engineering Data on Mixing, Elsevier Science, Amsterdam, 2000. pp. 85–115. [30] M.F.W. Distelho, A.J. Marquis, Scalar mixing in the vicinity of two disk turbines and two pitched blade impellers, Chem. Eng. Sci. 55 (2000) 1905–1920.
H.S. Yoon et al. / International Journal of Heat and Mass Transfer 91 (2015) 885–897 [31] S. Benayad, A. Salem, J. Legrand, Investigation of the interaction of velocity and concentration fields by coupling LDV to conductivity probe in a continuous stirred tank reactor (CSTR), Chem. Eng. J. 84 (2001) 239–245. [32] F. Guillard, C. Tragardh, Mixing in industrial Rushton turbine-agitated reactors under aerated conditions, Chem. Eng. Process. 42 (2003) 373–386. [33] V.V. Ranade, J.R. Bourne, J.B. Joshi, Fluid mechanics and blending in agitated tanks, Chem. Eng. Sci. 46 (1991) 1883–1893. [34] M. Lunden, O. Stenberg, B. Andersson, Evaluation of a method for measuring mixing time using numerical simulation and experimental data, Chem. Eng. Commun. 139 (1995) 115–136. [35] Z. Jaworski, W. Bujalski, N. Otomo, A.W. Nienow, CFD study of homogenization with dual Rushton turbines-comparison with experimental results: Part I: initial studies, Trans. Inst. Chem. Eng. Part A Chem. Eng. Res. Des. 78 (2000) 327–333. [36] K.H. Javed, T. Mahmud, J.M. Zhu, Numerical simulation of turbulent batch mixing in a vessel agitated by a Rushton turbine, Chem. Eng. Process. 45 (2006) 99–112. [37] H. Hartmann, J.J. Derksen, H.E.A. van den Akker, Numerical simulation of a dissolution process in a stirred tank reactor, Chem. Eng. Sci. 61 (2006) 3025– 3032. [38] J. Min, Z. Gao, Large eddy simulations of mixing time in a stirred tank, Chin. J. Chem. Eng. 14 (2006) 1–7. [39] Z. Wang, Z. Mao, X. Shen, Numerical simulation of macroscopic mixing in a rushton impeller stirred tank, Chin. J. Process Eng. 6 (2006) 857–863. [40] R. Zadghaffari, J.S. Moghaddas, J. Revstedt, A study on liquid–liquid mixing in a stirred tank with a 6-blade Rushton turbine, Iran. J. Chem. Eng. 5 (2008) 12–22. [41] R. Zadghaffari, J.S. Moghaddas, J. Revstedt, A mixing study in a double-Rushton stirred tank, Comput. Chem. Eng. 33 (2009) 1240–1246.
897
[42] H.S. Yoon, K.V. Sharp, D.F. Hill, R.J. Adrian, S. Balachandar, M.Y. Ha, K. Kar, Integrated experimental and computational approach to simulation of flow in a stirred tank, Chem. Eng. Sci. 56 (2001) 6635–6649. [43] H.S. Yoon, S. Balachandar, M.Y. Ha, K. Kar, Large eddy simulation of flow in a stirred tank, J. Fluids Eng. 125 (2003) 486–499. [44] H.S. Yoon, D.F. Hill, S. Balachandar, R.J. Adrian, M.Y. Ha, Reynolds number scaling of flow in a Rushton turbine stirred tank. Part I – mean flow, circular jet and tip vortex scaling, Chem. Eng. Sci. 60 (2005) 3169–3183. [45] H.S. Yoon, S. Balachandar, M.Y. Ha, Large eddy simulation of flow in an unbaffled stirred tank for different Reynolds numbers, Phys. Fluids 21 (2009) 085102-1–085102-16. [46] M. Germano, U. Piomelli, P. Moin, W.H. Cabot, A dynamic subgrid-scale eddy viscosity model, Phys. Fluids A3 (1991) 1760–1765. [47] J. Smagorinsky, General circulation experiments with the primitive equations. I. The basic experiment, Mon. Weather Rev. 91 (1963) 99–164. [48] J.W. Deardoff, A numerical study of three-dimensional turbulent channel flow at large Reynolds numbers, J. Fluid Mech. 41 (1970) 453–480. [49] K. Akselvoll, P. Moin, Large-eddy simulation of turbulent confined coannular jets, J. Fluid Mech. 315 (1996) 387–411. [50] C.L. Ribault, S. Sarkar, S.A. Stanley, Large eddy simulation of a plane jet, Phys. Fluids 11 (1999) 3069–3083. [51] C. Canuto, A. Quarteroni, T.A. Zang, Spectral Methods in Fluid Dynamics, Springer verlag, Berlin, 1988. [52] S.M. Shekhar, S. Jayanti, CFD study of power and mixing time for paddle mixing in unbaffled vessels, Trans. I Chem E Part A 80 (2002) 482–498. [53] S. Nagata, Mixing-Principles and Applications, Kodansha Ltd, Tokyo, 1975. [54] K. Horiuti, Assessment of two-equation models of turbulent passive-scalar diffusion in channel flow, J. Fluid Mech. 238 (1992) (1992) 405–433.