versión impresa ISSN 0327-0793
Lat. Am. appl. res. v.34 n.1 Bahía Blanca ene./mar. 2004
Numerical simulation of flow in a hydrocyclone
Juan Romero1, Rubens Sampaio2 and Rogerio Saldanha da Gama3
1 Mechanical Engng. Department, Universidade Federal do Espirito Santo, UFES AV. Fernando Ferrari, S/N - Goiabeiras, 29060 - 970 - Vitória - ES - Brasil
2 Mechanical Engng. Department, Catholic University of Rio de Janeiro, PUC-Rio,RJ, Brasil
3 Laboratorio Nacional de Computação Cientifica, CNPq, Petropolis, Brasil
Abstract A new method is presented to predict the localization of an air core, which is usually found on the internal flow of hydrocyclones. The control and stability of the air core affect the flow split between the products of a hydrocyclone. The flow split is one of the least understood aspects of the hydrocyclone operation. This split is greatly influenced by the air core radius, so an understanding of the air core behavior would enhance the prediction of the flow split. The liquid-air interface is characterized by means of the Young-Laplace jump condition. A steady flow of Newtonian fluid is being considered. It was developed a model considering the air core as having a cylindrical shape. This model that we are presenting considers a radius as optimum, when the final expression for the jump condition of the interface liquid - air is minimum. With the velocity field obtained for an optimal air core radius, we can trace the trajectories of the solid particles, thus, it is possible to simulate the performance of the hydrocyclone. The velocity field, the flow split and the selectivity curve obtained are compared to experimental results. Good concordance is achieved.
Keywords Hydrocyclone. Air Core. Turbulence. Free Surface.
Hydrocyclones are used extensively to separate and classify solid particles in the mineral processing industry. A strong rotational movement takes place inside the equipment by means of the tangential feeding of fluid and solid particles, because of this field, the solid particles suspended in the fluid tend to move towards the walls. Besides, by the high tangential velocity of the fluid in the central part of the device, the pressure decreases to values smaller than the atmospheri pressure. A low pressure region is created causing the formation of an air core about the central line. In spite of the simplicity of the geometry and operation of a hydrocyclone,it is extremely complex to explain with details the mechanisms of the dynamics of fluids. The difficulty of finding the actual flow of the hydrocyclones makes it necessary to specify the shape and localization of the air core surface. In the usual models of flows in a hydrocyclone, the interface that bounds the air core is modeled as a fixed cylindrical surface, that simplifies the problem greatly. This approximation avoids the necessity of calculating an unknown boundary that modifies the domain, where the field equations must be solved. Nevertheless, such simplification can produce bad results.
Barrientos et al. (1993) presented a theoretical model to calculate the air core considering that the liquid-air interface is of the Young-Laplace type. In this model, we must know the radial velocity gradient and the pressure jumps on the interface to calculate the air core. Steffens et al. (1993), analyzed experimentall the relation between the air core diameter and the drop of pressure in a single outlet cylindrical vortex chamber in a wide range of operating conditions. They obtained those relationships based on simple models and empirical correlations for the fluid flow. Dyakowski and Williams (1995) showed a model to calculate the flow splits using, as a boundary condition on the inter face liquid-air, the same proposal used by Barrientos et al. (1993). In a first stage, without the air core, the interface location and shape are marked out by spots in which the pressure is the same as the atmospheric pressure. Davidson (1995) analyzed the air core diameter using a simple model for the fluid flow. He considered a non-viscous fluid in each outlet where a factor modifies the flow considering the viscous ef fects. In each pressure drop it is considered the principle that the air core diameter is adapted to the largest flux. The resulting expression is iteratively applied during the calculation of the flow in the hydrocyclone and the computational grid is adjusted in each iteration. Williams et al. (1995) presented an experimental work on which electrical resistance tomography measurements had been used to visualize in real time the movement of the air core. These measurements have confirmed and quantified the oscillatory nature of the air core. When the pressure in the feed flow increases, the air core diameter in the hydrocyclone also increases. In the present work, the liquid-air interface is considered as being of the Young-Laplace type. Using an iterative process (Romero, 1997), the air core radius is updated in order to minimize the error in the Young-Laplace jump condition. The results obtained for the velocity field were compared with experimental results showing a satisfactory concordance. Further, by balancing the forces acting upon a particle, one can arrive at the slip velocity between the solid particles and the liquid phase, consequently the trajectory of every solid particle within the body of the device, can be determined simply by following the path of the particle from the inlet to either the overflow or the underflow, thus, the particle separation efficiency curves can be computed.
II. A MATHEMATICAL MODEL FOR THE INTERFACE
In order to solve the free boundary problem, it is necessary to characterize it. We will consider that the air core surface that is formed in an hydrocyclone operation is an interface of the Young-Laplace type. It means that there is a jump through the interface in the normal stresses, which is proportional to the interface mean curvature, the proportionality constant being the surface tension. Other conditions are the impermeability (kinematic) and slip (no shear tension) condition.
v · n = 0 (1)
[Tn · n] = -2Hσ (2)
[Tn · t] = 0 3)
where [Φ] = Φ1 - Φ2 is the jump of the property Φ through the interface, n is the unitary normal vector, which goes from (1) to (2), t is a unitary tangent vector and T is the stress tensor. [Tn · n] is the jump of normal tensions through the interface, 2H is the interface mean curvature, and σ the surface tension. As there is no fluid flow through the surface, the normal fluid velocity is null.
If the air occupies the region 0 <= r <= Ra where Ra= Ra(z), equation (2) can be written in this way:
Supposing that the air is an ideal fluid, and neglecting the viscous and dynamical effects due to air rotation (Sampaio and Romero 1994), then the pressure in this region is the atmospherical pressure (p0). Therefore, for the air
the liquid in the hydrocyclone is considered as an incompressible newtonian fluid, therefore the stress tensor is given by:
T = -pI + 2μD (6)
where D is the strain rate tensor, , p is the pressure, μ is the liquid viscocity, I is the unitary tensor of second order and ∇ is the nabla operator. In cilindrical coordinates, the velocity vector is: v = ver + weθ + uez, and
where: B = [er eθ ez]T
For an axissymetric surface parametrized by a curve of the type r = Ra(z), where the axis z is the symmetry axis, the interface mean curvature is:
the unit normal vector is
and the unit tangent vector is
A. Model with a cylindrical shape for the ai core
In the present model it is assumed that the air core surface has a cylindrical shape, characterized by one parameter: the radius Ra.
For a cylindrical surface, R'a = 0 and R''a = 0 so, from (8) and (9) the mean curvature 2H and the unit normal vector in the interface are:
|n = -er||(12)|
from (12), (7) on (6):
therefore, equation (4), for this model, can be simplified to:
III. A MATHEMATICAL MODEL FOR THE FLUID FLOW
The fluid flow in the hydrocyclone is modeled by using the Navier-Stokes equations, (continuity and momentum conservation). Because of the rotational fluid movement, a low pressure region is created in the central part of the device, so an air core of an unknown nature is produced. The air core shape can vary with the time or due to the hydrocyclone parameters and operational conditions. This makes the dynamic problem a free boundary problem. Therefore, beyond the calculation of the velocity and pressure fields, the free boundary shape should also be approximated liquid-air interface).
The general case corresponds to a three-dimensional fluid flow with a free boundary of unknown shape. However, the problem complexity can be reduced by considering some simplifications that do not change the essence of the original physical problem. Since the tangential inlet imposes a strong rotational movement in the hydrocyclone interior, generating a symmetry in the velocity field, thus, the velocity can be considered independent of the polar angle in all domain, the only exception are the regions near the feeding or the vortex finder. As the principal body velocity field determines the particles classification, there is no significant error in supposing an axis-symmetric fluid flow.
Due to operation conditions and geometry of the equipment, the flow is turbulent and anisotropic, thus hindering the modelling of the flow. Applying the Reynolds decomposition in the conservation equations of the momentum (Navier-Stokes) and considering the time average, an additional term is obtained it is called Apparent Stress of Reynolds, being thus necessary to use a turbulence model to close the global model.
For a steady turbulent flow of newtonian fluid, the conservation equations are:
- Linear momentum:
with the following boundary conditions (Fig. 1.b)):
- on ∂Ω1
- , on ∂Ω2. The pressure gradient is assumed to be zero at the wall
- , on ∂Ω3 and ∂Ω4, where n is the unit normal vector to the boundary.
- at free boundary ∂Ω5:
- Cylindrical model
- conditions to be imposed
- Condition to be satisfied:
IV. TURBULENCE MODELLING
One of the most used models in turbulent flow simulations is the well know κ - model (Launder and Spalding, 1972), (Mohamadi and Pironeau, 1994), where κ is the turbulent kinetic energy and is the rate of dis sipation of κ. This model supposes that the Reynolds stresses are proportional to the rate of deformation of the main flow. This concept is known as the hypothesis of Boussinesq,
The proportionality factor is the turbulent viscosity μt, a scalar quantity, that is the same for all the components of and different from the molecular viscosity of the fluid, which is a property of the fluid. The turbulent viscosity depends on several details of the flow in consideration. Thus, the distribution of the turbulent viscosity on the field of the low should be determined by the turbulence model.
The main limitation to the κ - model is its appli cation only for isotropic flows, that means that the length and velocity dimensions scales are the same for all directions. In complex flows, where there is a high rotational velocity, the length and velocity scales can vary according to the direction. It is well known that for this kind of flow, the model κ - is inappropriate Duginns and Frith (1987), Zerbini (1992). In the present work it was used the Differential Reynolds Stress Model (DRSM Launder et al. , 1987), implemented in Fluent (1995), which provides a good alternative to simulate the turbulence. This model calculates the local Reynolds stress, that is obtained from the transport equation
- Disipation rate
- kinetic energy κ
with empirical constant: σκ = 1.0, = 1.3 , C3 = 1.8, C4 = 1.44, = 1.44, = 1.98, Cμ = 0.09 Fluent (1995).
A. Boundary condition
For the inlet ∂Ω1 and for the outlets ∂Ω3 and ∂Ω4 are requested the turbulence intensity (I) and the characteristic lenght (L). These inputs are used to compute κ and values at the inlet and outlets, and are given by
where and Cμ is an empirical constant (0.09) and l is a length scale characteristic of the turbulence at the inlet flow. We calculate the inlet turbulence length l by multiplying the characteristic length L by 0.07 (l = 0.07L). This factor of 0.07 is obtained from the average mixing length in a turbulent pipe flow (Fluent, 1995), where L is the hydraulic radius. We used the computed value of κ to derive the Reynolds stress at the inlet.
where is the Reynolds stress component in the streamwise direction (normal to the flow inlet).
At solids walls we used wall function. This wall function is empiric and it is used in the grid points near the wall to estimate the effect of the wall on the flow. These functions are used to solve the entire turbulent boundary layer. The wall function is based on the as sumption that a fully developed equilibrium exist in the turbulent boundary layer. Therefore, all the relevant flow properties can be obtained from the log law that describes such boundary layers. This fact elimi nates the need to solve the equations of the turbulence model completely in the whole region of a turbulent boundary layer and it allows that the grid points near to the wall can be placed relatively far away from the wall. The approach of the wall function provides big savings in the computer effort required for the simulation of the turbulent low.
In the walls, we compute the values of the Reynolds stress near the wall and e from the wall function. It is used an equilibrium assumption and under this as sumption production and turbulence dissipation near the wall are equal. This assumption together with the law of the wall determines the turbulence production only. In the wall the Reynolds stress takes the null value and the gradient of the kinetic energy is zero. As the dissipation rate is infinite in the walls, it is pre scribed in a point close to the wall, (Fluent, 1995), where u* is the friction velocity defined in the log law, k is the Von Karman constant (042), and y is the distance from the point to the wall.
In the free boundary ∂Ω5, it is used the same condition that at the solid wall, with the exception of the slip condition for the mean velocity (18). This is not an appropriate boundary condition (Romero, 1997) because it is necessary to consider the interface fluctuation.
This closes the system.
V. SOLUTION METHOD TO THE FREE BOUNDARY PROBLEM
The conservation equations are solved using the finite volume method in generalized coordinates. This fact allows a good adaptation of the grid to the geometry of the hydrocyclone.
The algorithm used in the present work, makes an aproximation of the velocity and pressure fields (for a certain geometry); after an evaluation of the error obtained from the Young-Laplace equation along the air core surface is done, the algorithm makes an updating of the air core radius. The process is repeated until the desired approximation is reached.
The algorithm follows the following steps:
- An initial curve ∂Ω'5 is chosen to approximate the free interface ∂Ω'5 will be approximated to Ω'.
- The velocity and the pressure field are approximated in ft so they will satisfy the equations (17) and (18)
- The solution obtained in the step 2) is used to evaluate the error in the Young-Laplace jump condition (equation 19). In this case the error is the integral on ∂Ω5 of residual equation. The position ∂Ω'5 is adjusted to a new curve ∂Ω''5, i order to minimize the error in this third condition.
- If this condition is satisfied with the required accuracy, then the process stops and makes of ∂Ω''5 the final approximation of ∂Ω5. Otherwise, proceed this way: ∂Ω'5 = ∂Ω''5, and go to the step 2)
VI. PREDICTION OF HYDROCYCLONE PERFOMANCE
The classification is in some cases a primordial operation, especially when the wanted product has strict specifications of size. In other cases, it is an auxiliary operation of a process and it is here where it has its most important application in the mineral metallurgical industry. Thus, we have the mill operation in closed circuit, where the classification objectives are to make more efficient the mill and to ensure that the product of the operation is under certain size, going back to the mill the larger particles.
The classification performance of a hydrocyclone is influenced by the design variables, such as the hydrocyclone dimension, and the operating variables, such as the feed pressure and the physical properties of the feed solids and the feed pulp. The separation efficiency curve expresses the relationship between the weight fraction or percentage of each particle size of the feed, and the underflow discharge. The particle size is commonly used to represent the performance of the hydrocyclone.
The effect of the addition of solid particles to the fluid is the change of viscosity. The presence of solid particles increases the viscosity of the suspension significantly in relation to the viscosity of the pure fluid. In many cases the suspension becomes a non-newtonian fluid. However, for diluted suspensions it can still be considered as newtonian, but the effect of the solid particles in the flow through the correction of the viscosity must be considered. We used the following expression for the correction of the viscosity (Hsieh, 1988).
where μm is the viscosity of the suspension, μo the viscosity of the fluid. Equation (27) expresses the representative viscosity of the suspension. Cv is the fraction in volume of the solid particle given for:
being Cw the concentration of the solid particle in weight fraction, ρm the density of the suspension, ρp the density of the solid phase and ρl the density of the liquid.
The method of tracing particle trajectories from the predicted axial, tangential and radial velocity components of the flow are derived by tracing each particle location at successive time intervals using an algebraic slip approach. Before one can actually trace the particle trajectories in the hydrocyclone, one must know the relative velocities between particles and the liquid phase ( particle slip velocities). For diluted suspensions it can be considered that the movement of a particle is not affected by the other particles. Thus, we can analyze this as the movement of a solid body submerged in a fluid. The dynamics of the particle is obtained through the characterization of the forces that act upon it (Svarovsky, 1984). One of these forces is the drag force (FD), of difficult characterization. The conventional form of expressing the drag force in a spherical particle is:
where v is the relative velocity between particle and fluid, ρm is the density of the suspension, A is the projected area of the particle, and CD the drag coefficient. For a mass particle m under the influence of a field of acceleration a, the movement equation is given by:
The particle Reynolds number (Rep) that characterizes the flow around the particle is given by , in the common applications of the hydrocyclone, is usually smaller than 1 and rarely larger than 10. For hydrocyclone applications related to the separation of fine particles, the Reynolds number is very low, smaller than 0.2, as a result, the drag coefficient of a spherical particle (CD) depends on the Reynolds number. For a Reynolds number below 0.3, the flow regime belongs to the creeping flow region, The flow is laminar where the viscous forces prevail and the inertial forces can be neglected. The drag force can be de termined theoretically by the solution of the equations of Navier-Stokes without the inertial terms. This force depends on the number of Reynolds. In this way, we have an expression for the drag coefficient so that for the steady state the equation (31) reduces to:
Since significant forces do not exist acting in the particle in the tangent direction, we can suppose that the particle is moved with the same velocity as that of the fluid in the same direction
wp = wf (34)
where up, vp and wp are the components of the velocity vector of the solid particle and uf, vf and wf the components of the velocity vector of the fluid in the axial, radial and tangential direction respectively, and Δρ = ρp - ρm. When the number of Reynolds of the particle in the solid-fluid motion is bigger than 0.2, it is necessary to use a more general equation for the hydrodynamic force. In this case the concept of drag coefficient is used and equations (32) and (33) can be substituted by:
|wp = wf||(37)|
where and the drag coefficient relationship for a sphere can be approached by
The particle trajectories are obtained by the integration of equations (32) and (33) or (35) and (36).
VII. COMPARISON WIHT EXPERIMENTAL DATA
The model presented in this work was applied to two different hydrocyclones, one was the hydrocyclone AKW-100, located in the velocimeter laser laboratory at the etallurgical Engineering Faculty of the University of Concepción-Chile, operating only with water. In this case, only the flow was simulated approximat ing the velocity and the pressure fields. The second hydrocyclone was the one studied by Hsieh (1988), where the velocity and pressure fields were simulated. With these simulated results, the trajectories of solid particles of different sizes were traced and thus the curve of selectivity was plotted. The numeric simulation was carried out in a work station IBM-RISC system/6000 installed in the department of Mechanical Engineering of the Catholic University of Rio de Janeiro.
Table 1: Values of the geometric parameters (mm)
Figure 1: a) Hydrocyclone geometry b) Domain geometry
Table 2: Operational condition for the hydrocyclones AKW-100 and Hsieh (75 mm where Qd is the under flow and Qf is the feed flow
A study of sensibility was made for the operation parameters, shown in Table 3. We observed that, the choice of the radius can affect the feeding pressure and the split flow notably, being this last parameter important for building the classification curve. The results showed that the final values for the optimal radius is close to the experimental values.
Table 3: Numerical results for the hydrocyclone AKW-100 (Ra: air core radius (mm) * Optimal radius)
Figure 2 shows the error of the equation of Young-Laplace in the interface for different values for the ai core radius. The optimal radius obtained numerically is close to another obtained experimentally. We can also observe that the error varies in the center line this means that we can find an optimal shape for the air core where this error vanishes, indicating that the radius could not be constant, which is known by the experience, but to consider the air core in a cylindrical shape is a good approximation.
Figure 2: Error of the Young-Laplace equation in the interface, case 1
Figure 3 compares the predictions of the model with experimental measures for the tangential and axial velocities, showing good agreement. Moreover, the free vortex is very well predicted for the tangential component that characterizes this kind of flow. It is well-known that the air core fluctuates quickly and the experimental measures are distorted by noise near the free interface. However, some deviations are observed, both for axial and tangential velocities near the free interface. The assumption of a cylindrical air core and that the interface is not fluctuating can be responsible for these discrepancies, because the liquid-air interface is in fact a free surface that fluctuates.
Figure 3: Velocity profile to the case 1 (m/s) a) axial b) tangential
Figure 3 shows the typical profile of the tangential velocity (predicted and experimental). We observed that the experimentally measured maximum value is not reached by the theoretical prediction. An explanation for this fact is that a boundary condition adapted for the turbulence was not considered. The kinetic energy is considered zero. The real phenomenon has an oscillatory movement of the free interface. There fore, the kinetic energy as different from zero should be considered. Figure 4 shows the fluctuation of axial and tangential velocity (rms) and figure 5 shows the cross product . The experimental results show that the fluctuation of velocity near the interface grows considerably and it stays almost constant in the region of free vortex.
The model presented to simulate the flow in a hydrocyclone is used in the hydrocyclone of 75 mm. of Hsieh (1988). The operation parameters are given in table (2). The flow of a dilute limestone slurries were studied, where ρm and μm were calculated using the equations (29) and (27). In a first stage, the flow was simulated, approximating the velocitiy and pressure fields as well as the location of the liquid-air interface. With the approximate velocity field, the trajectories of the solid particles are calculated for several sizes, and with these values the selectivity curve is built, being compared with the experimental results. These results are shown in Fig. 7. Figure 6 shows the streamlines. In the Duginns work (1987) the streamlines do not reproduce the area of recirculation as a close region, which is the expected.
Figure 4: Fluctuation velocity profile (rms) to the AKW-100 (m/s) a) axial, b) tangential
Figure 5: Profile for the cross product of the fluctuation velocity AKW-100 () (m2/s2)
Figure 6: Streamlines to the hydrocyclone of 75 mm
Figure 7: Selectivity curve for the hydrocyclone of 75 mm. (selectivity × particle size)
It is possible to trace the particle trajectories by having the means to compute the particle velocities. Since axial symmetry is assumed, the description of the particle path in the radial and axial directions is sufficient only to determine whether a particle reports to the overflow or underflow. The projected trajectories for three sizes of solid particles injected in the feeding are shown in Fig. (8).
Figure 8: Trajectories of solid particle a) dp = 3.13 μm b) dp = 17.74μm c) dp = 35.50μm
From the results we can conclude that the choice of the air core radius for the hydrocyclone flow simulation affects notably the split flow prediction, an important parameter to build the curve of selectivity and the general flow prediction.
This work helps for an initial understanding of the air core behavior and its influence in the classification process of hydrocyclones. Although there were differences between predicted and experimental results near the free interface, a reasonable global agreement was found.
As it was shown in the work, it is of great impor tance to model the turbulent diffusion for a model that adapts itself to the characteristics of the flow. This was reached using the model RSM, which provides qualitative and quantitative solutions quite close to the experimental ones obtained for the velocity field. The use of models of turbulence of larger complexity, as the Differential Model of the Reynolds stress, still represent difficulties because it involves the simultaneous solution of 12 coupled equations. The convergence of the numeric solution is slow, being requested a high number of iterations, with the use of small sub-relaxation factors.
The presented model is applied to diluted suspensions, considered as a newtonian fluid.
1. Barrientos, A., Sampaio, R. and Concha, F., "Effect of the air core in the performance of a hydrocyclones" , XVII International Mineral Processing Congress, Sydney, Australia, 267-270, (1993). [ Links ]
2. Davidson, M. R., "An adaptative method of predicting the air core diameter for numerical models of hydrocyclone flow", Int. Journal Mineral Processing, 43, 167-177, (1995). [ Links ]
3. Duginns, R. K. and Frith, P. C., "Turbulences effects in hydrocyclones", 3 rd International Conference In Hydrocyclones, Oxford, England, 75-81, (1987). [ Links ]
4. Dyakowsky, T. and Williams, R., "Prediction of air-core size and shape in a hydrocyclone", Int. J. Mineral Process, 43, 1-14, (1995). [ Links ]
5. FLUENT Inc., Manual, version 4.32, (1995). [ Links ]
6. Hsieh, Kuo-Tai, "Phenomenological model of the hydrocyclone", Ph. D., Thesis, The University of Utah, Utah, USA, (1988). [ Links ]
7. Mohammadi, B. and Pironeau, "Analisys of the K-Epsilon turbulence model", John Wiley & Sons, England, (1994). [ Links ]
8. Launder, B. E., Reege, G. J. and Rodi, W., "Progress in the development of a Reynolds stress turbulence closure", J. Fluid Mech., 68, 537-566, (1975). [ Links ]
9. Launder, B. E. and Spalding, D. B., "Mathematical models of turbulence", Academic Press, London, (1972). [ Links ]
10. Romero, J. S., "Modelagem matemática, simulação e ajuste do modelo de um hidrociclone", In Portuguese, Dr. Sci. Thesis, PUC, Rio de Janeiro, Brasil, (1997). [ Links ]
11. Romero, J. and Sampaio, R., "A Numerical Model for Predicition of the Air-core Shape of Hydrocyclone Flow", Mechanics Research Communications Journal, 26-3, 379-384 (1999). [ Links ]
12. Sampaio, R. and Romero, J., "Modelo matematico para la simulación numerica del flujo de fluidos en un hidrociclone", In Spanish, Congreso Cientifico 10 años del CYTED, Cancum, Mexico, 254-255 (1994). [ Links ]
13. Steffens, P. R., Whiten, W. J., Appleby, S. and Hitchins, J., "Predicition of air core diameters for hydrocyclones", Int. J. Mineral Process, 39, 61-74 (1993). [ Links ]
14. Svarovsky, L, "Hydrocyclones", Holt Rinehart and Winston Ltd., London, U. K., (1984). [ Links ]
15. Williams, R. A., Ilyas, O. M. and Dyakowsky, T., "Air core imaging in cyclonic coal separators using electrical resistance tomography", Coal Preparation, 15, 149-163, (1995). [ Links ]
16. Zerbini, E. J., "Simulação numerica do escoamento em câmaras ciclonicas", In Portuguese, IV ENCIT, Rio de Janeiro, Brasil, 609-612 (1992). [ Links ]
Received: April 18, 2001.
Accepted for publication: May 27, 2002.
Recommended by Subject Editor Eduardo Dvorkin.