626
Ind. Eng. Chem. Res. 1995,34, 626-639
Simulation of Flow in Stirred Vessels with Axial Flow Impellers: Effects of Various Numerical Schemes and Turbulence Model Parameters Aksheya K. Sahu and Jyeshtharaj B. Joshi* Department of Chemical Technology, University of Bombay, Matunga Road, Bombay, 400 019, Zndia
The standard two-equation (&E) turbulence model has been used to numerically simulate the flow in a n agitated baffled cylindrical vessel with axial flow impellers. Three numerical schemes, namely, upwind scheme, hybrid scheme, and power-law scheme, were used to evaluate the competitiveness of the various schemes. The solutions were obtained by using the SIMPLE algorithm. The effects of initial guess values of the flow variables, the underrelaxation parameters, internal iterations, etc., on the rate of convergence were analyzed. The sensitivity of the model parameters on the flow characteristics was investigated thoroughly. Further, the effects of global grid sizes and near-wall grid sizes on the solution were also investigated. The boundary conditions in the impeller region for the velocities and turbulent kinetic energy were provided from the experimental data. However, for the energy dissipation rate, as iterative process was adopted taking the local information into account. Computations were carried out for six different impeller designs. It was found that the predicted values were qualitatively in good agreement with the experimental data. 1. Introduction
In recent years, the numerical simulation of turbulent flow in an agitated baffled cylindrical vessel has attracted the attention of many researchers (Ranade and Joshi, 1990; Kresta and Wood, 1991; Harvey and Greaves, 1982). This is due to the need to have a better understanding of the basic hydrodynamics of the flow pattern, which may lead to a better design of the impeller. This will optimize the power consumption for a given design objective such as mixing time. Smith (1991) reported that the lack of fundamental understanding of the proceses in stirred vessels leads to losses of the order $10 billion per year. Further, a detailed experimental study in an agitated vessel may not be feasible. However, a good mathematical model for a turbulent flow can provide detailed information of the flow characteristics in an entire vessel. The obtained numerical results depend on various factors, such as the value of the model parameters, the numerical scheme used, and the model itself. As has been pointed out by Amano (19841, error propagated by a certain type of numerical scheme may undermine the performance of a turbulence model. Further, as the governing equations for a turbulent flow field are nonlinear, generally, an iterative process is adopted for numerical calculation. The computational cost of an iterative process heavily depends on the algorithm used, the guess values of the dependent variables, the underrelaxation parameters, the solution of the set of algebraic equations which are obtained by discretizing the governing equations, etc. Sometimes a good numerical algorithm may not produce a converged solution due t o a bad combination of the aforesaid conditions. Therefore, before starting the computational process, the above-mentioned points should be taken into consideration. In the simulation of the flow field in agitated vessels most of the authors have used available commercial packages (Kresta and Wood, 1991; Bakker and Van Den Akker, 1991); hence attention has not been paid t o the above-mentioned
* Author t o whom correspondence should be addressed.
points, which seem to be very important information for numerical study of such problems. Generally, two types of impellers are used in an agitated baffled vessel: (i) radial flow impeller and (ii) axial flow impeller. Our study will be confined to the axial flow impellers. It has been found that, for a radial flow impeller, more attention has been paid to the specification of the boundary conditions rather than the numerical approach. This is obvious due t o the complex nature of the flow pattern around a radial flow impeller. However, we will only analyze the efficiency of the numerical algorithms cited by the authors. In the study of radial flow impellers, Harvey and Greaves (1982) used the axisymmetric assumption and they estimated the pressure gradient in the angular momentum equation based on a drag coefficient relation. They used a staggered grid control volume method to solve the governing equations. The hybrid difference technique was used for the discretization of convective and diffusive terms. They did not report the number of grids used, a check for grid independence, the number of iterations required to obtain a converged solution, and the convergence criteria. Further, it is also clear that the number of iterations required t o achieve the convergence does not depend only on the numerical algorithm but also on the initial guess values of the independent variables and the underrelaxation parameters. Therefore, the robustness of their algorithm can not be evaluated. Placek et al. (1986) used a two-scale model to take into account the nonhomogeneity of turbulence near a radial flow impeller. The authors used nonuniform grid size near the solid wall; the computational algorithm is an improved version of the one published by Gossman et al. (1969). The upwind scheme was used in their study and care was taken to minimize the false diffusion which is normally associated with the upwind scheme. At the interface of the control volume, the eddy viscosity coefficients were calculated by using the harmonic mean. They did not mention the convergence criteria, the number of iterations taken t o obtain a converged solution, and the effects of other related parameters on the computational
0888-5885l95/2634-0626$09.~0~00 1995 American Chemical Society
Ind. Eng. Chem. Res., Vol. 34, No. 2, 1995 627 aspects of the problem. Without this information, the efficiency of a numerical algorithm cannot be ascertained. Kresta and Wood (1991) used a commercially available code FLUENT for the numerical treatment of the three-dimensional study of a radial flow impeller. They used a nonuniform grid size of grid density 26 x 30 x 30 (z x r x 0) to take into account the rapid variation of the flow variables in the impeller region. At the wall, the boundary conditions were computed by using wall functions built into the FLUENT package. They reported that the convergence was achieved with 600 iterations for the normalized residuals of all the However, emphasis variables taken less than 5 x was not given t o the study of the effects of the turbulence model parameters on the flow characteristics, initial guess values, and underrelaxation parameters. Three-dimensional turbulent flow in an agitated vessel with a nonisotropic turbulence model for a radial flow impeller was investigated by J u et al. (1990). The authors used the power-law scheme of Patankar (1980) to discretize the convective and diffusive term of the transport equation. The solution was tested for grid independence for two sets of grid sizes, 26 x 30 x 25 and 23 x 30 x 25. It may be pointed out here that a small variation in the grid size in r direction may not display the grid independence character of a solution. The authors reported that a solution required about 650 iterations to converge from a good starting point provided by a two-dimensional solution. The convergence criteria of the solution algorithm were clearly defined. However, the effects of initial guess values and underrelaxation parameters on convergence were not discussed. Bakker and Van Den Akker (1991) simulated the flow field for the pitched blade turbine for a single phase and extended their results t o two-phase flow. The calculations were carried out using the available codes FLUENT and GHOST for single- and two-phase flows, respectively. Ranade and Joshi (1989) simulated the flow in an agitated baffled vessel with an axial flow impeller. They used the power-law scheme for discretizing the convective and diffusive terms. They investigated the performance of various iterative processes such as SIMPLE, SIMPLEC, SIMPLER, and PIS0 and found that SIMPLER is the best. The effects of the model parameters CZ,Cd, Uk, and a, on the flow characteristics were analyzed. They opined that, with the suitable adjustment of the parameters CZ and c d , the difference between the predicted value and the experimental value can be minimized. To check the grid independence of the solution, various combinations of grid sizes were taken. However, the effects of the guess values and the underrelaxation parameters on the convergence of the iterative process were not analyzed. For prescribing the boundary conditions at the impeller region, particularly for E , the authors assumed a linear profile in accordance with the experimental value of the total energy dissipation rate. In the present study, it will be shown that the boundary condition for E is actually not needed once the value of k is prescribed. In a very recent study, Fokema et al. (1994) simulated the flow field generated by a downflow pitched blade turbine using FLOW3D package with the k - E model and observed a good agreement between the experimental values and the predicted values. They modeled the impeller using a thin disk with inlets on both the surfaces. The boundary conditions for velocities and turbulent quantities on both the upper and lower surfaces of the impeller discharge were specified through
-O.' >
4
o"!
80
%
0 V
O*"t
8
i
Y
0.21
J
4
a ul v, Z
0
0o.5 t
~
W
0.6 0
006
0.156
0.234
0.32
DIMENSIONLESS RADIAL COORDINATE,r
Figure 1. Dimensionless axial mean velocity at 30 mm below the impeller center plane for PTDl impeller: 0 , 45" from the baffle; 0, 22l/z0 from the baffle; A, -22l/2" from the baffle.
experimental data. The SIMPLEC (Van Doormal and Raithby, 1984) algorithm was used as the velocity pressure coupling algorithm and the hybrid scheme was used for discretization of the diffusion and convection terms. The authors made a detailed discussion on the effects of boundary conditions at the impeller surfaces, which seems to be very important from the modeling point of view. The influence of the underrelaxation parameters on the convergence of the solution was not reported. The review of the existing literature clearly shows that there is a need to study the effects of various numerical schemes, iterative processes for solution of the algebraic equations, guess values of the variables, underrelaxation parameters, and internal iterations, etc., on the numerical simulation of stirred vessel, which may optimize the computational time and hence reduce the computational cost. Further, a detailed discussion on the effects of the model parameters on the flow variables is essential for optimizing the model parameters. In order to take into account the above-mentioned points for numerical simulation of the flow field in an agitated vessel, it is logical to consider a geometrically simple problem. Unlike the flow in a radial flow impeller, the flow phenomenon in an axial flow impeller can be considered as two dimensional without much loss of generality. This behavior is clearly demonstrated by the experimental data. Figure 1 shows that the axial velocity component variation is almost negligible at various angular positions. In the present study, we have considered a two-dimensional turbulent flow in a cylindrical agitated baffled vessel with an axial flow impeller. Ranade et al. (1992) measured the flow generated by different designs of axial flow impellers. As we mentioned earlier, the present study is concerned with the numerical simulation of the axial flow impellers and the parametric sensitivity of the flow characteristics. The impellers considered here are the following: (i)pitched bladed downflow turbine with two different sizes (D= 167 mm, PTD1; D = 155 mm, PTDB), (ii)curved pitched downflow turbine (D= 167 mm, CURPTD), (iii)modified pitched bladed downflow turbine (D = 167 mm, MODPTD), (iv) propeller (D= 167 mm, PROP), and (v)
628 Ind. Eng. Chem. Res., Vol. 34, No. 2, 1995 Table 1. Source Terms for the Generalized Equation 1
I
k
TOP SURFACE
1
- ----- - - --
0
0
SYMMETRIC ALONG TANK AXIS
I I1 1
1 \
I 4
TANK WALL
‘ -‘
/’--‘
where
G = Cdpt[2[($)” lU
p=-
+ ():” +
(2)”] + (E+ E)”]
and ,ut = k 2 k
Cd@RUtip
I
---I M PELLE R 0
f
Z
HI 3
BOT TOM
SURFACE
L
Figure 2. Schematic presentation of solution domain and a general view of the flow pattern (MODPTD).
multiple pitched bladed downflow turbine with 12 blades (D= 167 mm, MPTD). The predicted values will be compared with those of the experimental data for the above-mentionedsix impeller designs. This comparison will validate the usefulness of the present study. In the next phase of our work, three-dimensional turbulent flow will be considered. 2. Mathematical Formulation
We have considered a steady, two-dimensional, axisymmetric turbulent flow in a cylindrical b a e d vessel. The origin is placed at the center of impeller plane, and a cylindrical coordinate system has been chosen for a numerical study of the flow characteristics (Figure 2). The basic equations are consisted of the time-averaged equation of continuity and the Reynolds equations of motion. The eddy viscosity description was obtained through the k-E turbulence model (for high Reynolds number). With the above-mentioned assumptions, the nondimensional form of the governing equations may be written in a compact form as
where the notation 4 is a general representation for the flow characteristics. The expression for reff and S@are given in Table 1. The nondimensional quantities u , u , p , k , E , r, and z are defined as follows:
u = iilUtip, u = B/Utip, p = j j l ~ U ~ ~k ~=’ &/Uti:, , E
= ZR/Utip3 , r = FIR, z = ZIR
In the present study, the effective Prandtl numbers and a,are taken as constants. The standard values of the parameters c d , CI,CZ,Uk, and a,are taken from the work of Launder and Spalding (1974). However, in the present study, a detailed discussion will be provided
to show the effects of these parameters on the flow characteristics. Boundary Conditions. In order to obtain the numerical solution for the set of eq 1, the necessary boundary conditions should be defined on the boundary of the flow domain.
on the axis of symmetry:
at the vertical wall: r=1, u = u = k = ~ = O . O
(2b)
at the bottom of the vessel: z = 213; u = u = k = c = O . O
(2c)
at the top of the vessel:
_-_ a~ a~
z = -413; u = 0.0, - = - - ” - 0.0
az
(2d)
In the impeller region, the boundary conditions for
u, u, and k were provided from the experimental data. As the velocity components (u,u ) are provided in the impeller vicinity, it may be thought of that the eddy viscosity coefficient is known (implicitly). Hence, given the values of either k or E , the other variable can be determined numerically. Therefore, a guess value of E was provided at the impeller region in the starting of the computation and was adjusted iteratively with the local values of the dissipation rate. In the case of high Reynolds number flows, it is a general practice to modify the no-slip condition by the introduction of the wall function (Launder and Spalding, 1974). The wall function approach has several advantages: (i) it requires less computer time and storage, (ii) it allows the introduction of additional empirical information, and (iii) it produces relatively accurate results with fewer node points in the near-wall region. Therefore, the implementation of the no-slip condition is done by setting the normal velocity to the wall to zero and the tangential velocity t o the wall is specified by an empirical relation. In the present analysis, the logarithmic velocity variation has been used in the nearwall region for y+ > 30. Hence
Uk
(3)
Ind. Eng. Chem. Res., Vol. 34, No. 2, 1995 629 where yp+= ,pyU&i, y is the von Karman constant, and UTis the friction velocity. The subscript p denotes the values of the variables of the first node point near the wall. The value of the von Karman constant y is selected in the range of 0.41 f 0.015 (Pate1 et al., 1986). Ater modifying the no-slip condition for incorporating the wall function, it is also necessary to modify the turbulent kinetic energy k and the dissipation rate G. in the near-wall region. The full implication of the assumption of logarithmic velocity profile can be realized if (i) the shear stress in the near-wall region is uniform and (ii)the generation and dissipation of energy are in balance = Cd1I25 = Uz2). These assumptions demand that at the wall
- _ - 0.0 and ?Y
E
= cd3/4kp312
YY
(4)
3. Numerical Method of Solution To solve the set of elliptic partial differential eq 1 numerically, an efficient numerical algorithm that is numerically stable and accurate should be used. A stable method implies that the discretization of the convective and diffusive terms should be made in such a way that the coefficients of the dependent variables should satisfy the Scarborough criteria (Patankar, 1980). In the present analysis, we have used three such numerically stable schemes: (i)upwind scheme (UPS), (ii) hybrid scheme (HBS), and (iii) power-law scheme (PLS). A detailed discussion on the advantages and disadvantages of the various schemes may be found elsewhere (Patankar, 1980). A staggered grid system of control volume method has been used. The grid arrangements for the flow variables are shown on Figure 3. This grid system has no advantage in solving the velocity field, since the pressure gradients are easy to evaluate and the velocities are conveniently located for the calculation of the convective fluxes. The pressure and mean velocities are solved iteratively using the SIMPLE algorithm. The general form of the discretized eq 1 may be put in the following form: ap4p = a E 4 E
+a
h + as4s + a N h + Sc@
where
A(lPl>= 1,for the upwind scheme "(0, 1 - 0.51PI),for the hybrid scheme
I
"(0,
1- O.llPI5),for the power-law scheme
P is the cell Peclet number. Sc@and SP@ are obtained after linearizing the source terms for the respective variables, which is very crucial for obtaining a converged solution. The turbulence variables have (k and E ) positive values; therefore, the Sc part of the source term
-+r
Figure 3. Distribution of nonuniform staggered grid system in r-z plane: -, u velocity; t, u velocity; 0,( p , K , E ) .
for these variables should always be positive. In K equation this objective is achieved by putting the last term of the source expression (Table 1)in Spk and the generation term G in Sck. However, in the dissipation rate equation, the source term has been linearized with the help of the Newton-Raphson method, which automatically satisfies the above-mentioned criterion. Before discussing the solution technique for the set of algebraic equations for the flow variables, it would be appropriate to discuss the procedure for the inclusion of boundary conditions t o these equations for the different flow variables. The gradient condition at the wall can be incorporated by equating the wall coefficient to zero. The incorporation of the wall function can be made in various ways: (i) augmenting the peff at the wall and taking tangential velocity at the wall to be zero (Dutta and Acharya; 1993), (ii) equating the wall coefficient to be zero and augmenting the source term (Harvey and Greaves; 19821, and (iii) calculating the gradient of tangential velocity a t the wall as a@@ = UJyy (4 = u , u ) (Obi and Peric, 1991; Elkaim, et al., 1993). In the present analysis, it is found that the last approach takes less number of iterations to obtain the converged solution as compared to the other two approaches. Moreover, the results obtained by the last two approaches are almost the same, whereas the results obtained by the first approach deviates slightly from the other two. As we have adopted an iterative process to obtain the solution, the algebraic equations for a given variable can be treated as linear. These linear equations have been solved iteratively using the Gauss-Seidel (GS) pointby-point method and the line-by-line method. The lineby-line method, also is known as alternating direction implicit (ADI),was introduced by Peaceman and Rachford (1955). In the line-by-line method, the equations for the grid points along a chosen line would look one dimensional and are solved by the tridiagonal matrix algorithm (TDMA). The advantages and disadvantages of both the iterative processes for solving the algebraic equations have been discussed in section 4. Due to the iterative nature of the solution algorithm, the algebraic
630 Ind. Eng. Chem. Res., Vol. 34, No. 2, 1995 Table 2. Effects of Internal Iterations on the Rate of Convergencea
0 81
pressure turbulent kinetic dissipation corrctn eq velocity eq energy eq rate of eq iterations 4 3 4 4 510 6 3 4 4 518 4 3 7 4 320 4 3 10 4 320 4 3 4 7 544 4 2 7 4 320
u
L
a
c d
= 0.05, CI = 1.44,cz = 1.92,Uk = 1.0, 0, = 1.3, and y =
1.0.
P
hl
P
I -
z=-o.12
X
a Y,
ul W
-0.2
-0.21
I
1
I
0.16
0
I
I
I
0.32
0.48
0.6b
DIMENSIONLESS
1
I
080
0.96
RADIAL COORDINATE,r
Figure 4. Effects of various numerical schemes on the axial velocity and turbulent kinetic energy distribution a t the different axial positions (z = 0.52,-0.12, -0.56): -, corresponds to HBS; - -, corresponds to PLS; - * -, corresponds to UPS.
0
A
250
-1
QX
U
a
500
ITERATION Figure 5. Effects of various schemes and iterative processes on the rate of convergence: A , convergence of HBS with LTDMA iterative process; 0,convergence of PLS with GS iterative process; x , convergence of PLS with LTDMA iterative process; * , convergence of UPS with LTDMA iterative process.
equations for a given variable need not be taken t o complete convergence. It will be shown later that an excessive amount of effort in solving the algebraic equations that are based on only tentative coefficients is an unnecessary waste of computer time. Further, to achieve a better convergence, the mean velocities, the pressure corrections, and the turbulent quantities need t o be underrelaxed with suitable values of the under-
relaxation parameters. Furthermore, the guess values for the flow variables should also be chosen properly to obtain a converged solution. A detailed discussion on these points is given in section 4, which is also an aspect of the present study. The computation was carried out in T808 transputer with PC 386 computer. The iterative process is terminated when the mass source residue over each control volume is less than and the residuals for momentum, turbulent kinetic energy, and dissipation rate are less than It has been found that when mass residue is less than 5 x the residual criteria for the momentum and turbulence variables are automatically satisfied.
Ind. Eng. Chem. Res., Vol. 34, No. 2, 1995 631 Table 3. Variation of Average Dissipation Rate with the Different Grid Systemsa 35 x 41
a c d
21
= 0.05,
x
41
21 x 44
2.43 381
2.50 512
103 iteration
grid size 21 x 26 28 x 26
2.40 373
2.66 431
2.60 360
ci = 1.44, cz = 1.92, Uk = 1.0, U, = 1.3, and y
28 x 44
28 x 33
35 x 56
2.46 450
2.56 440
2.42 504
= 0.42.
Table 4. Effects of the Underrelaxation Parameters on the Rate of Convergence with a Given Initial Guess Solution" a, a, ak ac a, iterations ____
0.25 0.25 0.25 0.5 0.6 0.6 0.25 a
0.4 0.4 0.4 0.4 0.4 0.5 0.4
0.5 0.15 0.25 0.25 0.25 0.25 0.25
0.1 0.1 0.1 0.1 0.1 0.1 0.15
0.5 0.5 0.75 0.5 0.4 0.4 0.5
483 498 337 300 294 294 570
Parameter values same as in Table 2.
Table 5. Effects of Initial Guess Solution of the Flow Variables (u,u, k, 6) on the Rate of Convergence for a Given Set of Underrelaxation Parametersa U
V
k
E
iterations
0.025 0.0025 0.025 0.025 0.025 0.025 0.025 0.025 0.025
-0.021 -0.021 -0.051 -0.021 -0.021 -0.021 -0.021 -0.021 -0.021
0.015 0.015 0.015 0.01 0.005 0.015 0.0025 0.002 0.005
0.005 0.005 0.005 0.005 0.005 0.05 0.005 0.002 0.05
294 294 294 288 271 408 608 681 diverged
a a p = 0 . 6 , a U = 0 . 4 , a ~ = 0 . 1 , ~ = 0 . 2 5 , a n d a , = 0 .Parameter 4. values same as in Table 2.
7
1
-0.25
Table 6. Comparison between the Predicted and Experimental Values of Average Dissipation Rate x
, . . E
1w
4
impellers
IA
-0.2
PTD1 PTD2 CURPTD MODPTD PROP MPTD expt cd cd
=0.05
=0.09 a
3.16 4.07 5.03
2.80 2.56 2.92
3.74 4.45 5.18
0.87 0.79 0.90
c1 = 1.44, c2 = 1.92, Uk = 1.0, 0, = 1.3, and y
0.67 1.3 1.55
2.28 1.68 1.98
= 0.42.
4. Results and Discussion
Before presenting the numerically predicted values for the flow variables at different axial positions and their comparison with experimental data, it is worth discussing the following few points, which are very important for numerical calculations: (i)comparison of the different numerical schemes, (ii) effects of internal iteration on the rate of convergence and the iterative processes for solving the algebraic equations, (iii)effects of the grid size (global and near-wall region), (iv) effects of the underrelaxation parameters on the rate of convergence, (v) effects of guess values of the variables on the rate of convergence, and (vi) effects of the model parameters on the solution. Comparison of Different Numerical Schemes. The effects of various numerical schemes on the predicted values of the flow characteristics is shown in Figure 4. It is observed that the predictions of the flow characteristics by the hybrid and power-law schemes are almost the same, whereas the solution obtained by the upwind scheme differs substantially from the other two schemes. Further, the rate of convergence of the
I
0
016
032
I
048
I
084
I
080
I
096
DIMENSIONLESS RADIAL COORDINATE,r
Figure 7. Effects of the model parameters c characteristics at the different axial positions: -, c d = 0.09; - * -, c d = 0.13.
d on the flow c d = 0.05; - -
power-law scheme is faster than the hybrid scheme. However, the upwind scheme converges at a faster rate than the power-law scheme (Figure 5). The time taken by the upwind, hybrid, and power-law schemes per iteration under identical conditions is 4.0,4.14, and 4.8 s, respectively. The total time consumed by the upwind, hybrid, and power-law schemes to obtain the converged solution under identical conditions is 24,31, and 33 min, respectively. Although the upwind scheme takes much less time compared to the other two schemes, it may not be accepted for the numerical simulation of an axial flow impeller due to the underprediction of the values of the flow variables. Further, there is no significant difference in the computational time between the powerlaw scheme and the hybrid scheme. Therefore, for the robustness of the algorithm, the power-law scheme has been used. Effects of Internal Iteration on the Rate of Convergence. For an iterative process, it is normal
::I
632 Ind. Eng. Chem. Res., Vol. 34, No. 2, 1995
3.2
C
0%
+/-2
-0.56
-- _ _ _ _--- ---
~
0.0
c
u
-0.1 L -
L
' 1
s I-
-0.31
P
o
2.-0.56
I
I
0
I
1
I
I
0.48 0.64 0.80 DIMENSIONLESS RADIAL COORDINATE, r
0.16
0.32
096
Figure 9. Effects of Cz on the flow characteristics at the different axial positions: -, CZ= 1.72, z = 0.52; - -, CZ= 1.72, z = -0.12;
-.- , CZ= 1.72, 2 = -0.56;
.Z I
0
I
I
I
I
= 0.09, 'i = 1.0; 0, c d = 0.09, y = 0.42; 0.05, )' = 0.42.
X, c d
A , cz = 1.92, Z = 0.52; = -0.56.
X,
cz = 1.92,
1
0.16 0.32 048 0.64 OB0 0.36 DIMENSIONLESS R A D I A L COORDINATE, r
Figure 8. Variation of cd and the von Karman constant y on the axial velocity and turbulent kinetic energy a t the different axial positions: -, c d = 0.13, y = 1.0; - -, cd = 0.13, y =0.42; - * -, c d
= -0.12; 0, Cz = 1.92,
= 0.05, y = 1.0;
A, c d =
practice that two or three iterations are used for solving the algebraic equations for u, u , k,and E and the number of iterations used for pressure correction equations is 4-5 times that of said variables with the line-by-line TDMA method. It may be pointed out here that the above-mentioned practice does not yield a better convergence (Table 2). As the coefficients of the velocity variables are influenced by the eddy viscosity, the solution of the k and E equations should be iterated more times than the pressure correction equation. Further, an excessive number of internal iterations for k and 6 also does not enhance the rate of convergence of the solution. It may be remarked that a more accurate solution for the pressure correction equation slows down the convergence. A similar observation was also made by Patankar (1980). It has been observed that, besides the iterative process, underrelaxation parameters, initial guess values, etc., the eddy viscosity calculations at control volume faces also have a significant effect on the rate of convergence. It is generally accepted that, due to the rapid variation of eddy viscosity in a spatial direction, at control volume faces the values of the variables are calculated by using the harmonic mean instead of the arithmetic mean. However, in the present study, it is found that this practice sometimes yields a divergent solution. This may be due to the small magnitudes of the k and E variables in the upper part of the vessel, which augments the round-off error.
Finally, it is worth mentioning the advantages and the disadvantages of the Gauss-Seidel point-by-point and the line-by-line TDMA (LTDMA)iterative processes. It is found that the Gauss-Seidel method takes a greater number of iterations t o obtain a converged solution as compared to the line-by-line TDMA method when the algebraic equations are obtained by using the powerlaw scheme (Figure 5 ) . Further, under identical conditions, the Gauss-Seidel iterative process produces an oscillating solution. The time taken for completing one iteration by using the Gauss-Seidel method and the line-by-line TDMA method is about 4.28 and 4.8 s, respectively. However, the total time required to obtain a converged solution with the help of the Gauss-Seidel method is higher than that of the line-by-line TDMA method. Hence, the line-by-line TDMA method has been used for all the purposes. The oscillating behavior of the Gauss-Seidel method may be overcome by increasing the internal iterations of the variables; consequently, there is a further increase in the computational time. From these observations, one may conclude that a better iterative process for solving the algebraic equations may reduce the computational cost. Effects of Grid Sizes. The grid sizes taken to show the grid independence of a solution are of two types: the near-wall grid size and the overall grid size. It may be observed from Figure 6 that the introduction of a fine grid system in the near-wall region does have an effect in the near-wall region as well as in the bulk region. The only difficulty is that, with an increase in the grid density near the wall, the computational time increases. An approximate increase of one grid point near the wall region results in an increase in total number of iterations of 50 with the same guess values and underrelaxation parameters. It is very interesting to note that the
Ind. Eng. Chem. Res., Vol. 34,No. 2, 1995 633
0.51 -0.3c L
1
CUAPTO
o.c 0.3
-0.31
0.11 D I M E N S I O N L E S S RADIAL COORDINATE, r
’
0.16 0.32 0.48 0.64 0.80 0.96 DI M E N SlON LE S S R A D I A L COORDINATE, r
0.0
Figure 10. Comparison of the predicted and experimental values of axial velocity a t the axial location z = 0.52 for different impeller designs: 0, experimental data; -, c d = 0.09; - -, c d = 0.05.
Figure 11. Comparison of the predicted and experimental values of axial velocity at the axial location z = 0.36 for different impeller designs: symbols same as in Figure 10.
predicted value of the turbulent kinetic energy near the wall has a sudden jump when the solution is obtained by using a uniform grid size (Figure 6d). This behavior is not observed with the introduction of a fine grid size near the wall region. The behavior of the predicted values of k for a fine grid size is in conformity with the experimental data. It may be pointed out that, below the impeller region, one subdivision of near-wall control volume gives a grid-independent solution (Figure 68,121, whereas in the upper part of the vessel three t o four subdivisions are necessary to obtain a grid-independent solution. This variation may be related to the small magnitude of the turbulent kinetic energy predicted in this region. Does it imply that, due to the low values of the turbulent kinetic energy, the logarithm variation of the tangential velocity profile to the wall is confined to a very thin layer? In the present calculation, three node points have been introduced in the near-wall region. To study the effects of overall grid size on the flow variables, the grid sizes chosen are 28 x 19,14 x 44,28 x 33,and 35 x 41. It is of interest t o note that the lower number of grid points in the axial and radial directions underpredict the radial velocity and the turbulent kinetic energy and overpredict the axial velocity. It may be remarked that the variations in the magnitude of the axial velocity, which was predicted by using the grid size 14 x 44 and 28 x 33, is about 15%.
It is also observed that the difference in the results obtained by the grid size 28 x 33 and 35 x 41 is almost negligible. Hence, to conserve the computational time, in the present study, a grid density of 28 x 33 has been used. Although the flow characteristics are independent of the grid size for the grid density 28 x 33, the numerical prediction of the average dissipation rate (cav) (Table 3) converges to a limiting value for the large grid system, particularly in the axial direction. The variation of the average dissipation rate with various grid densities has been shown in Table 3. An analysis of the numerical data for eav suggests that, for large grid sizes, its value approaches a limiting value. The variation in the value of Eav for the grid system 28 x 33 and 35 x 56 is about 5%) whereas its variation for grid system 28 x 19 and 35 x 56 is approximately 16%. Effects of the Underrelaxation Parameters on the Rate of Convergence. It has been shown by Peric et al. (1987)that a nearly optimum convergenceis found for the SIMPLE algorithm for a, = 1 - a,. However, in the present study, it has been observed that, along with the underrelaxation parameters for the pressure and the velocities, the underrelaxation parameters for k, E , and the eddy viscosity are also vital to obtain a converged solution. With the suitable guess values for the flow variables, it may be seen from Table 4 that,
634 Ind. Eng. Chem. Res., Vol. 34,No. 2, 1995 -0.21
-0.21
> 2 -0.2t
0
,-0.11
-
*
-
-Xa
0.0
0
-
/
_ _ _--___-..---W , o.l, z 0-0.12 VI VI
0
0
_-pc
'I
MPTD
"
W
I 0.ot
-0.2sl-
0
0.0
0.08
0.24
0.40
0-56
0.72
0.88
D I M E N S I O N L E S S RAD1 A L C O O R DI N A T E , r
PROP
0.08
0.24
0.40
DIM ENS IONLES S
0.56
0.72
0.88
R A D l A L COORDINATE , r
Figure 12. Comparison of the predicted and experimental values of axial velocity at the axial location z = -0.12 for different impeller designs: symbols same as in Figure 10.
Figure 13. Comparison of the predicted and experimental values of axial velocity at the axial location z = -0.56 for different impeller designs: symbols same as in Figure 10.
for the values of a, = 0.6, a, = 0.4, ak = 0.1,a, = 0.25, and a,= 0.5, an optimum convergence can be achieved t o obtain solution for the PTD2 impeller, where a,, a,, ak, cc, and a, are the underrelaxation parameters for the pressure, velocities, turbulent kinetic energy, dissipation rate, and eddy viscosity calculations, respectively. Effects of Guess Values of the Variables on the Rate of Convergence. The effects of the initial guess values on the rate of convergence is shown in Table 5 . It is interesting to note that the guess values for u and u , practically, have no influence on the rate of convergence. However, the rate of convergence of the iterative process primarily depends upon the guess values of k and E . It is very difficult to establish any mathematical relation between the guess values of k and c that yield an optimum rate of convergence. However, an analysis of the numerical data of Table 5 suggests that a proper choice for the guess values of E and k may be found based on the information of cay (Table 6). It is found that, in the numerical calculation of the flow characteristics, for all the impellers that have been considered in the present study, the guess value of E is approximately 3-4 times that of average dissipation rate and the k value may be of the same order as E . It is clear from Table 5 that a bad guess value for lz and E
not only slows down the rate of convergence but also yields a divergent solution. It may be remarked here that a good starting value for the k and E can drastically reduce the total number of iterations for obtaining a converged solution, which ultimately reduces the computer time. In the present analysis, to start with, we have provided uniform profiles for all the flow variables. However, nonuniform profiles for the flow variables may result in a rapid convergence, which has not been attempted here. Effect of Model Parameters on Flow Variables. As c d is an important model parameter for the calculation of eddy viscosity, which influences the flow phenomena, its effects on the flow characteristics should be studied in detail. An increase in c d causes a reduction in the predicted circulation below the impeller region (Figure 7a). The magnitudes of the axial velocity (below and above the impeller) and radial velocity (below the impeller) (Figure 7a-c) increase with a decrease in the value of Cd. However, above the impeller, the magnitude of radial velocity increases with an increase in the value of Cd (Figure 7d). It is observed that the turbulent kinetic energy along the radial direction for a fxed axial position (near the bottom of the vessel) decreases with a decrease in the value of c d (Figure 7e), except near the wall. Near the wall region the effect of
Ind. Eng. Chem. Res., Vol. 34, No. 2, 1995 635 0.0L
0.12 CURPTD
O.02t
0
,-0 .oeF
L
0
v) -0.071 0
0.OlC
VI
z 0.d-
0.09
0
I/-\ PROP
0
-0.071 0.0
0.08
0
0.0
0.16
0.32
01M E N SI ONCE S S
0.48
0.64
R A DI A L
0.80
0.9
COOR D I N A 1El r
Figure 14. Comparison of the predicted and experimental values of radial velocity at the axial location2 = 0.52 for different impeller designs: symbols same as in Figure 10. c d is opposite. This implies that all the curves representing the k profile for different values Of c d cross each other at a small distance from the wall. It is observed that, as one moves toward the upper part of the vessel, the crossing point slowly shifts toward the axis. Finally, at a large distance from the bottom of the vessel, a decrease in the value of c d increases the k value in the entire radial direction for a fxed axial position (Figure 7e,f). Further, it is noticed that the distribution of k in the upper part of the vessel in the radial direction is similar to that of fully developed flow. As the wall function calculation depends on the parameter c d and the von Karman constant y , it would be interesting t o discuss the effect of y with the variation of c d (Figure 8). It is observed that below the impeller, for a lower value of c d , change in the von m a n constant affects the velocity profiles near the wall, as well as in the bulk. However, as c d increases, such a variation in the velocity profiles with variation of y is not noticed. This implies the change in the eddy viscosity is almost constant. It is of interest to note that, as one moves toward the upper part of the vessel, the variation in the values of y influences the flow characteristics, particularly the axial velocity profiles, irrespective of the values of c d . It may be remarked that, with an increase in the value of y , the rate of convergence of the solution becomes much faster (Tables 2 and 3). Its effects on
1
1
I
0.2L
O*LO
0.56
I
I
0.72
0.88
O I M E N S I O N L E S S RADIAL COORDINATE
r
Figure 15. Comparison of the predicted and experimental values of radial velocity a t the axial location z = -0.12 for different impeller designs: symbols same as in Figure 10.
turbulent quantities such as k and E are almost negligible except in the near-wall region. It has been observed that a lower value of the parameter CZ causes an increase in the predicted circulation below the impeller, which results into a sharp velocity gradients (Figure 9a). The lower value of CZalso resulted in a decrease in the overall predicted values of turbulent kinetic energy and the dissipation rate (Figure 9b). Further, its effects on the axidradial velocity profiles decrease with an increase in the distance from the bottom of the vessel. Similar observations have also been made by Ranade and Joshi (1989). The effects of the parameter C1 on the predicted flow variables were not reported by Ranade and Joshi (1989). It is interesting to note that an increase in the value of C1 influences the flow variables in a similar way if the value of C2 is decreased (Figure 9). It may be remarked that in the prediction of turbulent flow in an agitated vessel the parameters c d , C1, and CZmay be adjusted to make the predicted results comparable with the experimental data. The effect of the effective Prandtl number 0; is negligible on the velocity as well as on the turbulent characteristics except at the region below the impeller and close t o the axis. It is observed that, with an increase in the value of a,, the predicted circulation below the impeller increases. The effects of ok is also
636 Ind. Eng. Chem. Res., Vol. 34,No. 2, 1995
CURPTD
-0.06
3.6[ 2.0
L
0
0
-0.06
2
*
PTD 2
0
0
x
- -"=t
0.L
n.ntL
-0.05
00 . 0 -0.025
.
- - - - - _0_ _ _ _ _ _ _ _ - - 1- - -
1
0.1
L
0
PROP 0
-0.05
0.08
0.2~
0.40
0.56
0.72
o m
,
D I M E N SIONLESS RADIAL COORDINATE r
Figure 16. Comparison of the predicted and experimental values of radial velocity a t the axial location z = -0.56 for different impeller designs: symbols same as in Figure 10.
insignificant. However, its effect is less pronounced as compared to that of uE. Comparison of the Experimental Data with Numerical Predictions. Ranade et al. (1992)reported the flow generated by the different designs of axial flow impeller using the laser Doppler anemometer. In the present study, we have considered the same impellers (PTD1, PTD2, CURPTD, MODPTD, PROP, MPTD). The validity of the present theoretical analysis can be established by comparing the obtained results with the experimental data for the different flow variables at the different axial and radial positions. For the sake of brevity, we have chosen four axial positions (z = 0.52, 0.36, -0.12, -0.56). Although in our study we have investigated the effects of the various parameters on the flow variables, an optimization of parameter study based on the experimental data will not be attempted here. In the next part of our study, an attempt will be made to optimize the model parameters in accordance with the experimental data. However, experimental data have been compared with the predicted values of the flow characteristics for two different values of Cd (0.05,0.09), keeping all the other parameters constant. The comparison of the predicted and the experimental axial velocity profiles at the different axial locations ( z = 0.52,0.36) (Figures 10 and 11,respectively) along the radial direction reveals that for r > 0.32 the agreement is excellent for both the values O f Cd (0.05,0.09), except
0.0 0.16 0.32 0 . 4 8 0.64 0.80 0.96 D I M E N S I O N L E S S R A D I A L COORDINATE ,r
Figure 17. Comparison of the predicted and experimental values of turbulent kinetic energy at the axial location z = 0.52 for different impeller designs: symbols same as in Figure 10.
near the wall. Near the wall, the results obtained by c d = 0.05 yield a better agreement with the experimental data compared to c d = 0.09. When r .C 0.32, the predicted values of the axial velocity for c d = 0.05 are also in good agreement with the experimental data. At z = -0.12, an excellent agreement between the experimental data and the theoretical predictions for both values of c d (0.05,0.09) has been observed (Figure 12). However, at z = -0.56, the results obtained with c d = 0.09 show a better agreement when r < 0.72, whereas c d = 0.05 yields a better agreement for r > 0.72 (Figure 13). This discrepancy may be due to the underprediction of the turbulent kinetic energy in this region, particularly near the shaft. It is observed that a good agreement between the predicted and the experimental values of turbulent kinetic energy also yields a good agreement in the velocity profiles. This implies that by improving the predicted values of turbulent kinetic energy in the upper portion of the vessel a good agreement between the predicted and experimental values may be achieved. This improvement may be achieved by including the tangential velocity component in the numerical simulation, i.e., considering threedimensional flow or using a low Reynolds number model, since the magnitude of turbulent kinetic energy is very small in the upper part of the vessel. In earlier studies by various authors, the comparison between the
Ind. Eng. Chem. Res., Vol. 34, No. 2, 1995 637
o
o
o
_ - - -- --. .--
0
CURPTD O
0.61
0
O
O
O
0
0
0
0
0
0
A--__-MODPT D 00.34
0
0
0
I-
O
0
O
_ - - --. '.
.*
_--0
MODPTD
0 0 I-
_ _ e -
0.oL
0.08
0.24
0.25 l
1
O*ko
0.56
I
0.72
d
0.88
0.0
D I M E N S I O N L E S S RADIAL COORDINATE, r
Figure 18. Comparison of the predicted and experimental values of turbulent kinetic energy a t the axial location z = -0.12 for different impeller designs: symbols same as in Figure 10.
predicted and experimental values was made in the vicinity of the impeller. Their conclusions regarding the performance of a turbulence model and the choice of model parameters are related t o the information available close to the impeller. However, in the present analysis, it is observed that there is a qualitative agreement between the experimental and predicted values of the flow variables near the impeller region, whereas away from the impeller, particularly in the upper part of the vessel, the differences are large. This observation implies that a good agreement between the predicted values and the experimental data near the impeller region does not imply that a good agreement can be achieved far away from the impeller. It is interesting to note that, at the axial position z = 0.52, a good agreement between the experimental and predicted value of the radial velocity is observed in the region near the axis; remote from the axis the differences increase (Figure 14). The experimental data for all the impellers at z = 0.36 reveal that there is a sharp increase and decrease in the radial velocity profile. The distance at which this rapid variation occurs is different for the different impeller designs. For CURPTD, MODPTD, PROP, and MPTD,it occurs, approximately, at r = 0.32 from the axis. For PTDl and PTD2, it is observed at r = 0.08 and r = 0.4. The theoretical value of the radial velocity fails to predict such a sharp change
0.08
0.26
0.40
0.56
D I M E N S ION L E S S RADIAL
0.72
0.88
,
C O O R O l NATE r
Figure 19. Comparison of the predicted and experimental values of turbulent kinetic energy at the axial location z = -0.56 for different impeller designs: symbols same as in Figure 10.
0
z
0
0.16
0.32
DIMENSIONLESS
0.48
0.66
0.80
0.96
R A D I A L COORDINATE,r
Figure 20. Variation of turbulent kinetic energy and dissipation rate a t the different locations: . , 0,z = 0.52; A, A , z = -0.12;0 , 0 , z = -0.56.
for PTD2, MPTD, and PROP. However, there is a qualitative agreement between the experimental and predicted value near the axis and wall region. Although
638 Ind. Eng. Chem. Res., Vol. 34, No. 2, 1995
Table 7. Effects of Various Parameters on the Average Dissipation Rate cay 0.42 0.42 0.42 0.42 0.42 0.42 0.42 1.00 0
Z sO.5 5 3
c IJl
-0,121 9
DIMENSIONLESS
RADIAL
COORDINATE, r
Figure 21. Comparison of the predicted and experimental values of axial and radial velocities for the radial flow impeller at the axial locations z = 0.553 and z = 0.373: symbols same as in Figure 10.
a qualitative agreement between the experimental and predicted value is observed at the axial position z = -0.12 for c d = 0.09 (Figure 15),a total disagreement is observed at z = -0.56 for both values of the parameter Cd (Figure 16). Such a discrepancy also has been observed by Ranade and Joshi (1989) in the threedimensional analysis of an axial flow impeller. It may be due to the very low values of k and E in this region. It is found that the predicted values of k at the axial positions z = 0.52 in the radial direction are qualitatively in good agreement with the experimental data for the Cd value of 0.05 for all the impellers, except MODPTD and PROP (Figure 17). In the case of MODPTD and PROP impellers, for the c d value of 0.09 the agreement is very good. At the axial position z = -0.12, the pr4edictedvalues of k with c d = 0.09 display a better agreement with the experimental data in comparison with the data obtained at c d = 0.05 (Figure 18). The theoretical values of k at the axial location z = -0.56 for both values of the parameter c d (0.05,0.09) do not yield a satisfactory result near the shaft (Figure 19). As one moves away from the shaft, a qualitative agreement between the experimental and predicted values may be noticed, particularly for Cd = 0.05. It is well-known that, below the impeller, turbulence is anisotropic. In spite of the anisotropic behavior, the standard k--E model has predicted most of the features of the flow pattern in a baffled agitated vessel. Unfortunately, there are no experimental data available for the comparison of dissipation rate at the different axial positions. However, it will be worth showing the variation of the k and E values at the corresponding
1.0 1.0 1.4 1.0 1.0 1.0 1.0 1.0
1.3 1.0 1.3 1.3 1.3 1.3 1.3 1.3
1.44 1.44 1.44 1.44 1.64 1.44 1.44 1.44
1.92 1.92 1.92 1.72 1.92 1.92 1.92 1.92
0.05 0.05 0.05 0.05 0.05 0.13 0.09 0.05
2.56 2.48 2.62 2.36 2.44 3.23 2.92 2.32
locations. Figure 20 demonstrates that, except near the wall, the shape of the profiles for both the variables is similar and the magnitude of the E value is lower than that of the k value. In simulation of the flow field for the radial flow impeller (disk turbine), the two-dimensional assumption may be very restrictive. However, Figure 21 shows that the predicted values of axial and radial velocities agree qualitatively with the experimental data. Finally, it would be worthwhile to discuss the effects of various parameters on the predicted value of average dissipation rate and a comparison between the calculated value and the experimental data for the average dissipation rate. Table 7 shows that with an increase in the values Of c d , CZ,a, and u,, the average dissipation rate increases, while an increase in the value of C1 has an opposite effect. The predicted values of average dissipation rate for MODPTD, PTD1, and PROP are in excellent agreement with the experimental data when c d = 0.09. For PTD2, CURPTD, and MPTD impellers, the agreement is not so good (Table 6).
5. Conclusion The comparison of three numerical schemes shows that the power-law scheme is more robust and accurate as compared t o the hybrid and upwind schemes. With a suitable combination of guess values of the variables and underrelaxation parameters, the rate of convergence of the iterative process can be enhanced. The predicted values of axial velocity at the axial location z = 0.52, z = 0.36, and z = -0.12 agreed qualitatively and quantitatively with the experimental data. Once the values of the turbulent kinetic energy along with the velocities are prescribed near the impeller, the values of the dissipation rate can be adjusted using an iterative process. For an axial flow impeller, the model parameters may be optimized to yield a better prediction which may be comparable with the experimental data. The experimental data for the radial velocity profiles at z = -0.56 shows a random behavior, which is not observed by the predicted values. Experimental data reveal that, very close to the impeller, the turbulence is anisotropic. However, the standard k-E model with two-dimensional assumption has been able to predict most of the flow pattern.
Acknowledgment This research work was supported by a grant from the Department of Biotechnology,Government of India. Nomenclature C1 = constant for the k--E model Cz = constant for the k-6 model Cd = drag coefficient D = impeller diameter, m E = roughness constant in wall function
Ind. Eng. Chem. Res., Vol. 34, No. 2,1995 639
H
= height of the liquid in the tank, m k = turbulent kinetic energy per unit mass, m2 c2 k = dimensionless turbulent kinetic energy P = Peclet number p = pressure, N m-2 p = dimensionless pressure R = radius of the vessel, m i-= radial coordinate, m r = dimensionless radial coordinate $4 = source term for the generalized variable 4 U = mean velocity parallel to the wall at the near-wall node, m s-l v t i p = impeller tip velocity, m s-l U, = friction velocity, m c3-l U, = dimensionless friction velocity ii = mean radial velocity, m s-l u = dimensionless mean radial velocity 8 = mean axial velocity, m s-l v = dimensionless mean axial velocity 9 = local coordinate normal to the wall, m y = dimensionless local coordinate normal to the wall y+ = wall Reynolds number I = axial coordinate, m z = dimensionless axial coordinate
Greek Letters a, = underrelaxation parameter for eddy viscosity ak = underrelaxation parameter for the turbulent kinetic energy a, = underrelaxation parameter for the pressure a, = underrelaxation parameter for the velocities a, = underrelaxation parameter for the dissipation rate reff= effective viscosity for the variable 4 y = von Karman constant ? = rate of turbulent energy dissipation per unit mass, m2 s-3 6 = dimensionless turbulent energy dissipation rate cav = dimensionless average dissipation rate 9 = tangential coordinate p = dimensionless parameter p = molecular viscosity of fluids Pa s ,iit= eddy viscosity Pa s peff= dimensionless effective viscosity (p pt) pt = dimensionless eddy viscosity Q = density of water, kg m-3 ak = Prandtl number for turbulent kinetic energy a, = Prandtl number for turbulent energy dissipation rate 4 = generalized notation for transport variable Tw = wall stress, N m-2 zw = dimensionless wall stress
+
Subscript
p = values of the variables at the node point nearest to the wall
Dutta, S.; Acharya, S. Heat Transfer and Flow Past a Backstep with the Nonlinear k - E Turbulence Model and the Modified k-E Turbulence Model. Numer. Heat Transfer, Part A 1993,23,
281-301. Elkaim, D.; Reggio, M.; Camarero, R. Control Volume FiniteElement Solution of a Confined Turbulent Diffusion Flame. Numer. Heat Transfer, Part A 1993,23,259-279. Fokema, M. D.; Kresta, S. M.; Wood, P. E. Importance of Using the Correct Impeller Boundary Conditions for CFD Simulations of Stirred Tanks. Can. J. Chem. Eng. 1994,72, 177-183. Gossman, A. D.; Pun, W. M.; Runchal, A. K.; Spalding, D. B. Heat and Mass Transfer in Recirculating Flows; Academic Press: London, 1969. Harvey, P. S.;Greaves, M. Turbulent Flow in an Agitated Vessel, Part 11: Numerical Solution and Model Predictions. Trans. Z . Chem. Eng. 1982,60,201-210. Ju, S.Y.;Mulvahill, T. M.; Pike, R. W. Agitated Vessels with a Non-Isotropic Viscosity Turbulence Model. Can. J.Chem. Eng. 1990,68,3-16. Kresta, S. M.; Wood, P. E. Prediction of the Three Dimensional Turbulent Flow in a Stirred Tanks. AZChE J. 1991,37,448-
460. Launder, B. E.; Spalding, D. B. The Numerical Computation of Turbulent Flows. Comput. Methods Appl. Mech. Eng. 1974,3,
269-289. Obi, S.;Peric, M. Second-Moment Calculation Procedure for Turbulent Flows with Collocated Variable Arrangement. AZAA J . 1991,29,585-590. Patankar, S.V. Numerical Heat Transfer and Fluid Flow; Hemisphere: New York, 1980. Patel, V. C.; Rodi, W.; Scheuerer, G. Turbulence Models for NearWall and Low Reynolds Number Flows: A Review. AZAA J.
1986,23,1308-1319. Peaceman, D. W.; Rachford, H. H. The Solution of Parabolic and Elliptic Differential equations. J.SOC.Znd. Appl. Math. 1955,
3,28. Peric, M.; Kessler, R.; Schellerer, G. Comparison of Finite Volume Numerical Methods with Staggered and Non-Staggered Grids. Rept. No. 163/"/87,Lehrst. f. Stromungsmechanik, Univ. Erlange, Nbg, 1987. Placek, J.; Tavlarides, L. L.; Smith, G. W.; Fort, I. Turbulent Flow in Stirred Tank 11: A Two Scale Model of Turbulence. AIChE
J. 1986,32,1771-1786. Ranade, V. V.; Joshi, J. B.; Marathe, A. G. Flow Generated by Pitched Blade Turbine Part 11: Simulation using k - E Model. Chem. Eng. Commun. 1989,81,225-245. Ranade, V. V.; Joshi, J. B. Flow Generated by Disc Turbine, Part I1 Mathematical Modelling and Comparison with Experimental Data. Trans. Z . Chem. Eng. 1990,68,34-50. Ranade, V. V.; Mishra, V. P.; Saraph, V. S.; Despande, G. B.; Joshi, J. B. Comparison of Axial Flow Impeller Using a Laser Doppler hemometer. Znd. Eng. Chem. Res. 1992,31,2370-2379. Smith, J. M. Industrial Needs for Mixing Research. Trans. Znst. Chem. Eng. 1991,68 (Part A), 3-6. Van Doormal, J. P.; Raithby, G. D. Enhancements of the SIMPLE Method for Predicting Incompressible Fluid Flows. Numer. Heat Transfer 1984,7 , 147-163.
Received for review March 23, 1994 Revised manuscript received September 16,1994 Accepted October 4, 1994@
Literature Cited Amano, R. S. Development of a Turbulence Near-Wall Model and its Application to Separated and Reattached Flows. Numer. Heat Transfer 19s4,7 , 59-75. Bakker, A.; Van Den Akker, H. E. A Computational Study on Dispersing Gas in a Stirred Reactor. 7th European Congress on Mixing 1991,Part 1, 199-207.
IE9401840 @
Abstract published in Advance ACS Abstracts, December
15,1994.
