An improved LES on dense particle-liquid turbulent flows using integrated Boltzmann equations.

Tang, Xuelin ; Wu, Jing

INTRODUCTION

Extensive research has been conducted to investigate dilute particle-fluid two-phase flows using continuum mechanics (Zhou, 2002). Most studies assumed that particle motions were due to fluid motions and the collisions between particles could be ignored. The two-phase Reynolds-averaged equations were developed using various averaging methods and continuum theories, leading to models such as the particle trajectory model (Zhang et al., 1989), the multi-fluid model k-[epsilon]-[A.sub.p] (Laslandes and Sacre, 1998) and the LES model (Wang and Squires, 1996b; Wu et al., 2001). Most models assumed that the fluid-phase was a continuum and the particle-phase was discrete or assumed that both of the fluid- and the particle-phase were continua.

Although the application of continuum theory to dilute two-phase flows has been very successful, this approach does not consider inter-particle collisions, which is a crucial feature that distinguishes dense particle-liquid two-phase flows from dilute flows. Viscous effects and diffusion effects due to the interparticle collisions must be considered in dense particle-liquid two-phase turbulent flows.

The kinetic theory based on the Boltzmann equation accurately describes the microscopic interaction properties, and has been widely used in the simulation of particle-liquid two-phase flows. Ni et al. (2000) investigated sediment-laden two-phase flows. In their work, the Boltzmann equation of the particle-phase was greatly simplified by ignoring the inter-particle collisions in case of dilute flows, and the velocity distribution of both phases could be easily solved. For dense flows, Ni et al. (2000) employed the BGK model to simplify the inter-particle collision term of the Boltzmann equation, using a single relaxation time. Ishii (1975) and Pai (1977) multiplied the Boltzmann equation of each phase by property parameters and integrated over the velocity space to obtain the continuum and momentum equations of each phase. However, their collision integral is very complex and cannot be integrated analytically. Liu (1993) used a similar analysis to develop a solution for collisions in dilute two-phase flows. By assuming the relationship between temperature and velocity are the same for both the particle and gas phases, Enwald et al. (1996) and Mathiesen et al. (2000a; b) analyzed the flow parameters in a fluidized bed. Because of this assumption, the difference between gas molecules and particles cannot be distinguished. Aidun and Lu (1995) and Filippova and Hanel (1997) directly resolved the Boltzmann equation of two-phase flows at the mesh level.

Parallel to the development of governing equations of particle-liquid two-phase flows based on the Boltzmann kinetic theory are the continued improvements on the techniques for solving Navier-Stokes (N-S) equations. There are three types of computational fluid dynamics (CFD): direct numerical simulation (DNS), Reynolds averaged N-S (RANS) equation model, and large-eddy simulation (LES). Since the LES requires less computer resources than the DNS and its applicability is more universal than the RANS, it has become an important and promising approach for engineering applications (Ferziger, 1996, 2000; Murakami, 1998). LES is a turbulence simulation method, which assumes that the turbulent motion can be separated into grid-scale and subgrid-scale eddies. The separation of the subgrid-scale eddies from the grid-scale eddies does not have a significant effect on the evolution of the grid-scale eddies. The subgrid-scale eddies are independent from flow geometries and boundary conditions. Assuming that the subgrid-scale eddies are isotropic and only contribute to dissipations, the LES obtains the grid-scale field by solving a set of filtered governing equations.

Smagorinsky applied LES to single-phase flows for meteorological applications and established the Smagorinsky model (Smagorinsky, 1963), which has been very successful in modelling turbulent flows. Leith (1990) noticed some drawbacks of the Smagorinsky model. To resolve these drawbacks, Germano et al. (1991) proposed a dynamic subgrid stress (SGS) model to directly compute model parameters on the resolved scale. In other words, the parameters of the dynamic SGS model can be obtained as a function of flow domain and time. The dynamic SGS model is capable of predicting the asymptotic behaviour near a wall correctly, and considers the energy backscattering from small to large scales. Although the dynamic SGS model has been successfully applied to simple flows geometries, improvements are necessary for complex geometries. Lill (1992), Yan et al. (1993), Poimelli and Liu (1995) and Ghosal et al. (1995) analyzed and improved the dynamic SGS model. Sinha and Katz (2000) investigated flows in a pump impeller with diffuser vanes by employing the Smagorinsky model.

In recent years, LES has been adapted to investigate two-phase turbulent flows. Using two-way coupling, Boivin et al. (2000) combined LES for the fluid-phase and a particle trajectory simulation for the particle-phase to simulate dilute gas-particle flows. By comparing their computational results to DNS data, they concluded that the subgrid-scale dissipations estimated by the Smagorinsky model were reasonable. The relatively low values of subgrid-scale dissipations could be interpreted by the argument of scale similarity. Using similar approaches, Wang and Squires (1996b) modelled dilute gas-particle two-phase flows, neglecting inter-particle collisions. The results showed that the simulated turbulent flux of the particle-phase was accurate except for near-wall regions. A dynamic SGS model was employed by Wang and Squires (1996a) to solve the same problem. The results were in good agreement with DNS simulation and experimental data. Jaberi and James (1998) analyzed turbulent combustions using a dynamic similarity model. In the LES discussed above, the SGS was assumed to depend only upon the strain-rate tensor, but not upon the rotation-rate tensor. This assumption entails that all these models would be less appropriate for situations where highly rotational flows are prevalent (Olsson and Fuchs, 1998).

In this investigation, a molecular gas kinetic theory analogy (Chapman and Cowling, 1970; Vincenti and Kruger Jr., 1965) is used to analyze dense particle-liquid two-phase turbulent flows. To derive the continuity equation of a phase, the Boltzmann equation is multiplied by the characteristic mass of the phase and is integrated over the velocity space. Similarly, the momentum equation is obtained to by weighting the Boltzmann equation with particle or molecular momentum and integrating over the velocity space. In essence, our approach belongs to the two-fluid strategy, treating the particle phase as a pseudo-fluid. With the governing equation of the particle phase derived from the microscopic Boltzmann equation, our model provides a balanced description of the particle and the fluid nature of the particle-phase.

To further improve existing LES, a second-order SGS with double dynamic coefficients is used in LES. In this way, our SGS is dependent not only upon the strain-rate tensor, but upon the rotation-rate tensor. Therefore, our simulation is applicable to highly rotational flows such as the dense particle-liquid two-phase turbulent flows, in which the rotation of the particle-phase cannot be ignored. It is important to realize that our approach only considers macroscopic rotations of two-phase flows, not microscopic self-rotations of each particle. Our model is used to simulate the dense particle-liquid two-phase turbulent flows in the geometries of a duct and a centrifugal impeller. The simulated pressure and velocity distributions are in good agreement with experimental results.

THE GOVERNING EQUATIONS OF DENSE PARTICLE-LIQUID TWO-PHASE TURBULENT FLOWS

The derivation of the governing equations of dense particle-liquid two-phase flows was briefly discussed in Tang et al. (2002). The starting point is the microscopic Boltzmann equation. By weighting the Boltzmann equation of each phase by property parameters and integrating over the velocity space, the continuity and momentum equations are derived.

The derivation details and glossary of symbols are discussed in the Appendix. The results are summarized here. The governing equations for the liquid-phase are:

[partial derivative]/[partial derivative]t([[rho].sub.f] + [partial derivative]/[partial derivative][x.sub.fl] ([[rho].sub.f][u.sub.fl]) = 0 (A42)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (A43)

where [F.sub.fpci] = [F.sub.pfci].

The governing equations for the particle-phase are given by:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (A44)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (A45)

where [F.sub.ppci] = 2/3(1 + e)[nabla]([C.sub.p][[rho].sub.p]<[c.sup.'.sub.p], and [F.sub.pfci] = [[rho].sub.p]/[[tau].sub.rp]([u.sub.pi]-[u.sub.fi]).

FILTERED AND MASS-WEIGHTED GOVERNING EQUATIONS

In a dynamic LES, two filtering operators, i.e., grid filter G and test filter G, are usually defined, with G denoting G G.

By applying the filter functions G and G to Equations (A42) and (A43), the filtered liquid-phase equations are obtained:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (1)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (2)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (3)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (4)

For the liquid-phase, the relationship among its volumetric concentration [C.sub.f], its density [[rho].sub.f] and the pure liquid density [[[rho].sub.ff] can be expressed as: [C.sub.f] = [[rho].sub.f]/[[rho].sub.ff].

It can be seen from Equations (1) and (3) that continuity equation of the filtered liquid-phase cannot satisfy mass conservation because of density fluctuations. To maintain mass conservation, a mass-weighted approach is applied to the governing equations. Mass-weighted operators can be defined as:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]

With the filtered liquid-phase Equations (1), (2), (3) and (4) weighted by mass, the following equations are obtained, respectively:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (5)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (6)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (7)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (8)

where, [q.sub.fij] and [T.sub.fij] are the SGS and the subtest-scale stress, respectively.

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (9)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (10)

The resolved stress can be defined as:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (11)

Substituting the formulas (9) and (10) into (11), we have:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (12)

According to Cauchy-Helmholtz theorem, the velocity at any infinitesimal point of a particle can be decomposed into the translation velocity of the particle centre, the rotation velocity of the particle as a rigid-body with respect to the point, and the velocity related to the particle deformation. Therefore, a flow is more accurately modelled if both deformation and rotation are considered. In LES, this requires SGS be expressed as a function of both the strain-rate tensor [S.sub.fij] and the rotation-rate tensor [R.sub.fij].

Given [S.sub.fij] and [R.sub.fij], Pope (1975) put forth only five Reynolds-averaged invariants:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]

It is believed that the vortex stretching is the dominant mechanism by which turbulence transfer energy from the grid-scale to the subgrid-scale, and the SGS should be a function of the strain rate tensor and the rotational rate tensor. Wallin and Johansson (2002) proposed an explicit higher-order algebraic subgrid-scale Reynolds stress model including the resolved symmetric strain rate tensor and the resolved anti-symmetric rotation rate tensor, which was able to accurately capture rapidly rotating turbulence.

For a theoretical perspective, the most general expression of the SGS formula [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] should be an infinite tensor polynomial expansion of the resolved [S.sub.fij] and [R.sub.fij], where [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] and [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]. Based on the eddy-viscosity scheme, Canuto and Cheng (1997) determined the Smagorinsky constant according to the five irreducible invariants of [S.sub.fij] and [R.sub.fij]. Based on the symmetry of the SGS, the dimensional consistency and the five irreducible invariants of [S.sub.fij] and [R.sub.fij], Craft et al. (1996) and Ehrhard and Moussiopoulos (2000) proposed a stress-strain non-linear relation:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (13)

The coefficients of Equation (13) are listed in Table 1.

where [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] in Table 1.

To capture the highly rotational feature of the two-phase flows and to reduce the computational complexity entailed by the direct use of Equation (13), we propose a second-order SGS model with double dynamic coefficients in Equation (14), wherein only the second and the fourth terms on the right-hand side of Equation (13) are retained. They are non-cross second-order terms.

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (14)

This SGS model can be used in LES of dense particle-liquid two-phase flows. Coefficients [c.sub.f1] and [c.sub.f2] are dynamic model parameters that can be obtained directly at the resolved scales. [DELTA] is the subgrid filter width.

Because the model is foreign to any filter function, the subtest-scale stress model can be written as:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (15)

where [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] the resolved liquid-phase rotation tensor

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.], is the subtest filter width, which is usually set to [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.].

Substituting Equations (14) and (15) into (11), we have:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (16)

Let [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] and

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] we have:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (17)

Applying a least-square approach:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (17)

We have the expressions for the two model parameters as:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (18)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (19)

where () represents the average over a plane parallel to the wall.

Finally, Equations (5), (6), (14), (18) and (19) form the basis of the LES of the liquid-phase of turbulent dense particle-liquid two-phase flows. The LES uses a dynamic second-order SGS model with double dynamic coefficients.

Applying the same derivation procedure to the equations of the particle-phase, the LES on the particle-phase using a dynamic second-order SGS model with double dynamic coefficients can be developed in a similar fashion.

NUMERICAL METHODS

The finite volume method (FVM) is used to discretize the particle-liquid equations together with a staggered grid system, where the second order difference is employed for the diffusion terms and the upwind difference QUICK is used for the convection terms (Leonard, 1979). The SIMPLEC (Van Doormaal and Raith, 1984) algorithm and the iterative method are used to solve the discretized equations for the liquid-phase and the particle-phase, respectively. Body-fitted coordinates (Kollmann, 1984) are used to simulate flows over complex geometries. The fluid used in the simulations is pure water at the atmospheric pressure.

CASE STUDY I: TWO-PHASE FLOWS IN A DUCT

Duct Geometry and the Coordinate System

In order to validate the accuracy and credibility of these models, turbulent flow through a duct is simulated and is compared to the experimental results of Wang and Chien (1984), (1985). Wang and Chien conducted many sets of experiments in a pressured duct with a rectangular cross-section (0.18 x 0.10 [m.sup.2]) using uniform plastic sand with [d.sub.p] = 0.25 mm and = 1.067 x [10.sup.3] kg/[m.sup.3]. In the experiment, the mixture of the water and solid particles in the fed tank was driven by a slurry pump of 100 HP. The averaged discharge and transport concentration of the slurry were measured by a gauge tank. Velocity distributions of the flow in the duct were measured by a specially designed pitot tube, which prevented the solid particles from clogging the opening of the tube. Samples of the slurry were siphoned out at different elevations for concentration distribution determination. The pressure in the sample tank was carefully adjusted so that the velocity at the inlet of the sampler was about the same as the local flow velocity. The frictional loss was measured by manometer. By injecting the coloured liquid into the flow and by observing the extent of spreading of the coloured liquid in a transparent inspection section, the intensity of turbulence was judged qualitatively (Wang and Chien, 1984; 1985).

The dimensions of the duct and flow parameters of the simulation are identical to those used by Wang and Chien. The cross-section of the duct is 100 mm in height and 180 mm in width. The length of the duct is 4 m and is in the streamwise direction. This is shown in Figure 1. The mean velocities and volumetric concentrations are 2.389 m/s and 0.457, respectively.

The numerical results obtained from the FVM increase in accuracy as the grid spacing is decreased. In general, it is desirable to increase grid resolution in regions where the flow variables exhibit large gradients. Our grid topology is non-uniform where the grid is packed densely near the walls. Several computational trials were run with various grids topologies before the final grid topology was chosen, which is shown in Figure 2. The final grid consists of I x J x K = 400 x 46 x 26 cells, where I is the number of divisions in the X direction; J the number of divisions in the Y direction; and K the number of divisions in the Z direction. The simulation results are grid independent. This is demonstrated in Table 2.

[FIGURE 1 OMITTED]

[FIGURE 2 OMITTED]

Boundary Conditions

Inlet conditions: the velocity is assumed uniform for each phase.

Outlet conditions: the velocity gradients are assumed naught for each phase.

Wall conditions: the slip conditions (Ding and Gidaspow, 1990) are used for the particle-phase. The particle-phase slip velocity

on wall can be written as [u.sub.s,w] = [k.sub.0][[partial derivative]u.sub.p/ [[partial derivative]n.sub.w], where [n.sub.w] is the coodinate normal to wall, [u.sub.p] is the velocity in the streamwise direction, [K.sub.0] = [C.sub.p]/[square root of (2[pi][n.sub.p][d.sup.2.sub.p] and [n.sub.p] is the local particle number per unit volume. The non-slip and non-penetration conditions are employed for the liquid phase.

Boundary conditions of pressures: the Neumann condition is used for the liquid-phase. The pressure of the particle-phase is assumed equal to that of the liquid-phase. In order to guarantee computational stability, a reference pressure point is specified.

In order to guarantee the continuity of the liquid-phase in the computational domain, the velocity at the outlet is adjusted according to the difference between the flow rate at the inlet and the outlet.

Velocity Distributions at Various Y-Z Cross-Sections

From Figure 3, it can be observed that the velocity distributions across the Y-Z sections are rectifying themselves and reach steady states, i.e., the fully developed two-phase turbulent flow, at a distance about X = 2.0 m from the inlet. In the regions where the turbulent flow is fully developed, the simulation results are in good agreement with the experimental data. Discrepancies are observed between the measurements and the simulations in under developed zones (near the entrance or near the wall) because our model is developed for simulating turbulent flows. It cannot account for laminar flows or the flows in transition from laminar to turbulent, which takes place near the entrance or near the wall.

Velocity Distributions at Various X-Z Cross Sections

Figures 4 and 5 suggest that the velocity distributions at the various X-Z sections are rectifying themselves and reach steady states at a distance about X = 2.0 m from the inlet. Slip velocities of the particle-phase on walls are also simulated.

[FIGURE 3 OMITTED]

[FIGURE 4 OMITTED]

Volumetric Concentrations

Figure 6 suggests that the distributions of the volumetric concentrations across the Y-Z section are rectifying themselves and form a flat profile at a distance about X = 3.0 m from the inlet. The flat concentration profile was also experimentally observed in literature under the same experimental conditions (Wang and Chien, 1984; 1985).

Pressure Gradients

Figure 7 suggest that the pressure gradients in the streamwise direction are in good agreement with experimental data.

In Figure 7, [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.], where [lambda] is the particle's linear concentration, [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.], D the distance between centres of two particles, s the distance between the spherical surface of two particles, [C.sub.pm] the maximum of volumetric concentration possible, and U the mean velocity on the cross-section of duct in the streamwise direction.

CASE STUDY II: TWO-PHASE FLOWS IN CENTRIFUGAL PUMP

Governing Equations of Dense Two-Phase Flows in a Rotating Coordinate System

Centrifugal slurry pumps are widely used in the chemical industry. Particle abrasion is a severe problem for pumps (Walker and Bodkin, 2000; Wilson et al., 1996), which causes damages to pump impellers. These abrasions are influenced by many factors such as particle shapes, sizes, distributions and process parameters. It is of great importance to carry out investigations on the behaviour of two-phase turbulent flows in impellers to understand particle abrasion damages to impeller so that pump performance and efficiency can be improved and pump life can be extended (Attride, 2002; Roco, 1990). Baumgarten and Muller (2000) improved pump performance by employing CFD commercial software to simulate the two-phase flows in impeller. The kinetic and dynamic LES models developed in the preceding sections are applied to solve for the dense two-phase flows in the geometry of a centrifugal impeller.

[FIGURE 5 OMITTED]

Using a Cartesian coordinate system (x, y, z) fixed on the impeller and rotating around z-axis with a constant angular velocity, the continuity and momentum equations for the twophase flow can be written in the following Cartesian tensor form.

The governing equations for liquid phase are:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (20)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (21)

where [u.sub.fi] is the relative velocity of the liquid-phase, r the vertical distance to the rotating z-axis, [omega] the angular velocity vector, [[epsilon].sub.jki] the permutation tensor, [[omega].sub.j] the angular velocity, and

p = P -1/2[[omega].sup.2][r.sup.2].

[FIGURE 6 OMITTED]

[FIGURE 7 OMITTED]

The governing equations for the particle-phase are:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (22)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (23)

Equations (20) to (23) are solved using the same techniques as those presented in the previous section.

Simulation Parameters

In order to validate the simulation, the computational results are compared to the experimental data of Wu et al. (1998). In their experiment, due to the difficulties in measuring the two-phase velocity distributions in a centrifugal pump, only the pressure distributions were measured using semiconductor pressure transducers with high frequency responses installed in the shroud. The particles, which are an admixture of polypropylene and precipitated barium sulphate, have an averaged diameter of 2.0 mm, a density of 1930 of kg/[m.sup.3], and a volumetric concentration of 5%. Other parameters of the centrifugal pump are listed in Table 3.

In our simulation, the same parameters as those of the experiment are used.

The Computational Domain and Mesh

It is assumed that the flows in impeller are periodic so that the computational domain is restricted to a blade-to-blade region. In order to simulate the real two-phase flow between blades, the computational domain between blades is extended to the pump inlet and the outlet. Therefore, the complete computational domain consists of a passage between blades and extensions to the inlet and the outlet. The numerical results obtained from the FVM increase in accuracy as the grid spacing is decreased. The grid topology is non-uniform. Higher grid resolutions are used for near-wall regions. Several computational trials were run with various grid resolutions in order to insure the choice of final grid resolutions. The non-uniform grid topology is shown in Figure 8. The grid consists of I x J x K = 128 x 87 x 51 cells, where I is the number of divisions in the streamwise direction, from the inlet section to the outlet section; J the number of divisions in the pitchwise direction, from the pressure side to the suction side; and K the number of divisions in the spanwise direction, from hub to front cover of the impeller. The simulation results are grid independent. This is demonstrated in Table 4.

Boundary Conditions

The inlet, outlet, wall and pressure conditions are identical to those used for the duct. In addition, periodic boundary conditions for the extensions to the inlet and the outlet are applied:

[[PHI].sub.left] = [[PHI].sub.right] ([PHI] = [u.sub.kr],[u.sub.k[theta]],[w.sub.k],[p.sub.k]

where [u.sub.kr] is the radial velocity, [u.sub.k[theta]] the peripheral velocity, and [W.sub.k] the axis velocity.

[FIGURE 8 OMITTED]

[FIGURE 9 OMITTED]

[FIGURE 10 OMITTED]

Pressure Distributions

The pressure coefficient is defined as:

[psi] = p/2[[rho].sub.m][u.sup.2.sub.2]

where, [[rho].sub.m] is the mixture density, and [u.sub.2] the rotating velocity at the outlet of the impeller.

Figure 9 shows that the pressure gradually increases from the inlet to the outlet and the acting force due to the pressure gradient between the pressure side and the suction side is along the rotating direction of the impeller. The simulated pressure distribution is in good agreement with the experimental data.

The Distribution of Particle Number Density

The particle number density is shown in Figure 10. It decreases from the inlet to the outlet. This is consistent with the fact that the impeller passage is divergent. The simulation results suggest an almost uniform distribution in the passage, although the number density is slightly higher on the pressure side than that on the suction side.

Velocity Distributions from Hub to Shroud

Figure 11 and Figure 12 show the relative velocity vectors of the liquid- and the particle-phase, respectively. It can be observed that the relative velocities of the liquid- and particle-phase on the suction side are larger than those on the pressure side. This is a consequence of blade rotation. The relative velocities of the liquid- and particle-phase at the inlet are larger than those at the outlet. This is due to the divergence of the impeller passage. The relative velocities of the liquid- and particle-phase around the middle of passage are larger than those near the hub and shroud. This is due to the effect of wall damping, which is especially prominent for the particle-phase. In addition, at a same position, relative velocities of the particle-phase are a slightly larger than those of the liquid-phase, which indicates that the motions of particles are motivated by the ambient liquid. The flow of each phase around the leading edge of the blade is also simulated, as well as the slip velocities of the particle-phase at walls.

[FIGURE 11 OMITTED]

[FIGURE 12 OMITTED]

CONCLUSIONS

In this work, the k-phase flow properties of dense particle-liquid two-phase turbulent flows are modelled by LES. The governing equations are derived from the kinetic theory of molecular gas. By multiplying the k-phase Boltzmann equation by appropriate property parameters and integrating over the velocity space, the continuity and momentum equations are obtained. The interparticle collision terms for dense two-phase flows are derived from kinetic theory. Each term in the derived momentum and continuity equations of the particle- and liquid-phase has clear physical significance.

This simulation is characterized by the following features.

1. The filtered governing equations are simplified by a massweighted approach in order to maintain mass conservation and continuity.

2. A second-order SGS model with double dynamic coefficients is used in LES, accounting for the macroscopic rotational characteristics of the dense particle-liquid flows, but not at a great expense of computational resources.

The simulation results (i.e., the two-phase pressure and velocity distributions) are in good agreement with the experimental data, suggesting that our kinetic model and the LES with the second-order SGS are applicable to turbulent dense particle-liquid flows. The simulation strategies (e.g. the FVM, the SIMPLEC algorithm and the body-fitted coordinates) are capable of handling two-phase flows over the complex flow geometries. Our model and simulation strategies provide a basis for the design of pump impellers and other equipment handling dense particle-liquid two-phase flows. Future work will focus on other complex flow geometries and optimization of equipment design.

APPENDIX: THE DERIVATION OF GOVERNING EQUATIONS OF DENSE PARTICLE-LIQUID TWO-PHASE TURBULENT FLOWS

Microscopic and Macroscopic Variables

From a theoretical perspective, a particle-liquid two-phase flow can be regarded dense when the ratio between the mean particle response time and the particle collision time is much greater than unity. This criterion implies that inter-particle collisions cannot be neglected. In engineering practice, a particle-liquid two-phase flow can be considered dense when the volumetric concentration of particles is greater than 5%.

A statistical approach can be applied to the particle-phase when the volume unit d[R.sub.k] ([x.sub.1],[x.sub.2],[x.sub.3]) is much smaller than flow volume, but much larger than a particle volume, where [R.sub.k] is a vector in geometric space. Employing molecular gas kinetic theory to a particle-liquid two-phase flow, the k-phase molecular or particle velocity distribution, [f.sub.k], can be determined from the Boltzmann equation.

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (A1)

where [c.sub.k] is the k-phase molecular or particle velocity, which is given by d[R.sub.k] * [dc.sub.k] is the k-phase acceleration, which is given by [F.sub.k], where [F.sub.k] is the k-phase molecular or particle mass forces, and [m.sub.k] the k-phase molecular or particle mass. [[theta].sub.k] is the k-phase molecular or particle temperature, and [r.sub.k] the k-phase molecular or particle radius, [f.sub.k] the number of molecules or particles at spatial position [R.sub.k], at time t with velocities between [c.sub.k] and [c.sub.k] + [dc.sub.k], temperatures between [[theta].sub.k] and [[theta].sub.k] + d[[theta].sub.k], and particle radii between [r.sub.k] and [r.sub.k] + [dr.sub.k]. The number density is [f.sub.k][dc.sub.k]d[[theta].sub.k][dr.sub.k]. Variations in this quantity will lead to changes in macroscopic properties.

In the derivation that follows, the subscript k will be replaced by f or p to represent the liquid- or the particle-phase, respectively. For example, [r.sub.p] and [r.sub.f] stand for particle radius and molecular radius, respectively. The symbol k' stands for -k. Therefore, ([partial derivative][f.sub.k]/[partial derivative]t) is the collision term between particles in the same phase, while [([partial derivative][f.sub.k]/[partial derivative]t).sub.kc] is the collision term between particles in different phases, with the subscript c denotes collision.

The Boltzmann equation conserves the equilibrium of the molecules or particles per unit volume. Therefore, the total variation of the k-phase distribution function [f.sub.k] with respect to time t equals the sum of the molecule or particle number variations with respect to all variables.

A statistical average of each property parameter [PHI] of the k-phase is given by:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (A2)

Transport Equations for the Particle- and the Liquid-Phase

In a phase-space, [PHI] stands for a k-phase parameter, which is a function of velocity [c.sub.k] and radius [r.sub.k], but not a function of [R.sub.k] and t. Thus, [c.sub.k], [R.sub.k] and t are independent variables. If the particles have the same radius, [r.sub.k] and temperature, [[theta].sub.k], the fourth and fifth terms on the left side of Equation (A1) are zero. Therefore, Equation (A1) can be simplified to:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (A3)

Multiplying this equation by the characteristic parameters [PHI] and integrating over the velocity space gives:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (A4)

Each term on the LHS of Equation (A4) can be further expressed as follows:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]

As [c.sub.k] [right arrow] [+ or -][infinity], [f.sub.k] [right arrow] 0, therefore:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]

Since [F.sub.ki] and [c.sub.ki] are independent from each other, we have:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]

As a result of the statistical averaging, the transport equation can be written as:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (A5)

Continuity Equations for the Particle- and Liquid-Phase

Substituting [PHI] by [m.sub.k] in Equation (A5) leads to the k-phase continuity equation:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (A6)

Assuming that there are no changes in the number of particles and molecules, -[([S.sub.k]).sub.kc]-[([S.sub.k]).sub.k'c] = 0. Since the k-phase mass density is given by ?k = nkmk and the macroscopic velocity [u.sub.ki] = <[c.sub.ki]> the continuity equation can be written as:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (A7)

Momentum Equations for the Particle- and Liquid-Phase

Substituting [PHI] by [m.sub.k][c.sub.ki] in Equation (A5) leads to the k-phase momentum equation:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (A8)

where <[c.sub.ki][c.sub.kj]> = <([u.sub.ki]+[c.sup.'.sub.ki])([u.sub.kj]+[c.sup.'.sub.kj])> = [u.sub.ki][u.sub.kj] + [T.sub.kij],[F.sub.ki] = [m.sub.k][g.sub.ki], where [g.sub.ki] is the k-phase mass force. The k-phase stress tensor is given by [T.sub.kij] = <[c.sup'.sub.ki][c.sup.'.sub.kj]> = [p.sub.k][[delta].sub.ij] -[[tau].sub.kij], and the k-phase partial pressure is given by [p.sub.k] = 1/3 ([T.sub.k11] + [T.sub.k22] + [T.sub.k33]), where [[tau].sub.kij] is the viscosity stress tensor. With these notations, Equation (A8) can be written as:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (A9)

The k-Phase Pressure

For dense two-phase flows, the k-phase pressure is assumed proportional to its volumetric concentration [p.sub.k] = [C.sub.k]P, where P is the total pressure of the two-phase flow, and [C.sub.k] is the k-phase volumetric concentration, with k=f or k=p for the liquid- or the particle-phase volumetric concentrations, respectively. Moreover, [C.sub.p] and [C.sub.f] satisfy the following equality:

[C.sub.p] + [C.sub.f] = 1 (A10)

The Collision Integrals

The terms on the right-hand side of Equations (A8) and (A9) are collision integrals. While there are many similarities in behaviour between liquid and gas molecules, there are few between gas molecules and solid particles. For example, the kinetic energy of gas molecules is directly related to their temperature, whereas the kinetic energy of particles is not. Additionally, gas molecules are modelled as elastic spheres with no energy losses during collisions, whereas particles suffer energy losses during collisions. Finally, each particle volume is often several orders of magnitude larger than the volume of a gas molecule, and particle volumes cannot be ignored during a simulation. These features complicate the evaluation of collision integrals when the particle-phase is involved. In the following derivations, the subscript k in Equations (A8) and (A9) will be replaced by f or p to represent the liquid- or the particle-phase, respectively.

Liquid-liquid collision integrals

Liquid molecules can be treated as elastic spheres, experiencing no energy losses during collisions. Since the energy is conserved, the collision integral for the liquid phase is zero:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (A11)

Liquid-particle collision integrals

The second term on the right-hand side of Equations (A8) and (A9) is the integral for the collisions between liquid molecules and particles. They are adequately described by the respective theories developed for dilute two-phase flows. The force exerted by the liquid-phase upon the particle-phase is

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.], which is further expressed in

Equation (A12):

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (A12)

where [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] is the averaged relaxation time of particles, [d.sub.p] the particle radius, [[rho].sub.p] the particle density, and is the liquid viscosity. The particle relative Reynolds number is given by [Re.sub.p] = |[u.sub.p] - [u.sub.f]|dp/v, where v is the liquid kinematic viscosity, [u.sub.p] the particle velocity, and [u.sub.f] the liquid velocity. The force of the particle-phase acting on the liquid-phase is given by [F.sub.fpcj] = -[F.sub.pfcj].

Particle-particle collision integrals

When Particles 1 and 2 collide at the contact point [R.sub.p], the instantaneous positions of Particles 1 and 2 are [R.sub.p]--[r.sub.p.sup.h] and [R.sub.p] + [r.sub.p]h, respectively, with h being the unit vector drawn between the centres of Particles 1 and 2. The velocities of Particles 1 and 2 are represented by [c.sub.p1] and [c.sub.p2] before the collision and [c.sup.'.sub.p1] and [c.sup.'.sub.p2] after the collision. The velocity of Particle 1 relative to Particle 2 is represented by [c.sub.p12] = [c.sub.p1]--[c.sub.p2] before the collision and by [c.sup.'.sub.p12] = [c.sup.'.sub.p1]--[c.sup.'.sub.p2] after the collision. The collision is inelastic with a particle elastic recovery coefficient of e. [cp.sub.12] and [c.sup.'.sub.p2] should satisfy the following relationship:

[c.sup.'.sub.p12] * h = -[ec.sub.p12] * h (A13)

The momentum of Particles 1 and 2 should also satisfy the momentum conservation equation:

[m.sub.p] ([c.sub.p1] + [c.sub.p2]) = [m.sub.p]([c.sup'.sub.p1] + [c.sup.'.sub.p2]) = 2[m.sub.p]G (A14)

where G is the velocity of the centre of mass of Particles 1 and 2. The following relationship can be easily obtained:

[c.sub.p1] = G + 1/2 [c.sub.p12] (A15)

[c.sub.p2] = G + 1/2 [c.sub.p12] (A16)

[c.sup.'sub.p1] = G + 1/2 [c.sub.p12] (A17)

[c.sup.'sub.p1] = G - 1/2 [c.sub.p12] (A18)

From Equations (A13) to (A18), we have:

[c.sup.'.sub.p1] - [c.sub.p1] = 1/2(c.sup.'.sub.p12] - [c.sub.p12]) = -1/2(1 + e) [(h * [c.sub.p12]).sup.h] (A19)

The Jacobian is given by:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (A20)

and the relationship between [dc.sub.p1][dc.sub.p2] and dG[dc.sub.p12] is given by:

dG[dc.sub.p12] = |J|[dc.sub.p1][dc.sub.p2] = [dc.sub.p1][dc.sub.p2] (A21)

Similarly, we have:

dG[dc.sup.'.sub.p12] = [dc.sup'.sub.p1][dc.sup'.sub.p2] (A22)

An area element dS can be envisioned at point [R.sub.p] with h being its unit norm, with Particle 2 on the positive side and Particle 1 on the negative side of dS. When Particles 1 and 2 collide, the centres of Particles 2 must be within a cylinder whose height is along the h direction, with dS being its base and [2r.sub.p] being its height. The volume of the cylinder is [2r.sub.p](h * n)dS, and the number of Particles 2 whose velocities are within [[c.sub.p2], [c.sub.p2] + [dc.sub.p2]] is given by:

[2r.sub.p][f.sub.p2]([R.sub.p] + [r.sub.p]h)(h * h)dS (A23)

Particle 1 is located on the surface of a sphere with a radius of [2r.sub.p], whose centre is the centre of the corresponding Particle 2. With [c.sub.p12] as the polar axis, [beta] as the polar angle between h and [c.sub.p12], and [phi] as the azimuthal angle, the number of Particles 1 whose velocities are within [[c.sub.p1], [c.sub.p1] + [dc.sub.p1]] is given by:

[4r.sup.2.sub.p] ([c.sub.p12] * h)[f.sub.p1]([R.sub.p] - [r.sub.p]h sin[beta]d[beta]d[phi][dc.sub.p1] (A24)

Therefore, the collision probability within a unit time is given by:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (A25)

where [f.sub.p1] and [f.sub.p2] are the velocity distribution functions of Particles 1 and 2.

Each collision transfers the quantity ([PHI]-[PHI]) of the property [PHI] from the negative side to the positive side of dS, wherein [PHI] is the property before the collision and [PHI]' is the property after the collision. The vector flux of [PHI] can be written as:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (A26)

The momentum variation per unit time and per unit volume, i.e., the collision force between Particles 1 and 2 can be expressed as:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (A27)

To further express J, the zero-order terms of the Taylor expansions of [f.sub.p1] and [f.sub.p2] at [R.sub.p] is used:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (A28)

With [PHI] = [m.sub.p][c.sub.p1], the momentum transfer of Particles 1 is given by:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (A29)

Using the equation [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] we have:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (A30)

where I is the second order unit tensor.

When the solid particles are at local equilibrium, [f.sub.p] can be assumed to obey the Maxwell distribution:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (A31)

where [T.sub.p] is the fluctuation kinetic energy of the particles. Therefore, J can be further expressed as:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (A32)

Since [m.sub.p] ([c.sup.2.sub.p1] + [c.sup.2.sub.p2]) =[2m.sub.p] ([G.sup.2] + 1/4[c.sup.2.sub.p12] can be replaced by dG[dc.sub.p12] we have:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (A33)

Integrating over all directions of [c.sub.p12] and G, we have:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (A34)

With the help of the following equation:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.], when b is even when b is odd

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (A35)

The relationship between the fluctuation kinetic energy [T.sub.p] and the fluctuation velocity [c.sup.'.sub.p] (c.sup.'.sub.p] = [c.sub.p] - <[c.sub.p]>) of the particles can be obtained by using the approaches of the kinetic theory of molecular gas. Assuming the distribution function ([f.sub.p]) of [c.sup.'.sub.p] is Maxwellian, we have:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (A36)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (A37)

Integrating over all directions of [c.sup.'.sub.p], we have:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (A38)

Finally, J can be expressed as:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (A39)

and the collision force between particles, i.e., the momentum variation per unit time and per unit volume, can be expressed as:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (A40)

For a dilute two-phase flow, the analysis usually assumes that there are no collisions between particles. Therefore, the collision term does affect neither the pressure nor the viscosity of the particle-phase and the collision integral can be ignored.

The particle-phase viscous coefficient can be obtained from molecular kinetic theory:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (A41)

The Governing Equations for the Particle- and Liquid-Phase

In summary, the governing equations for the liquid-phase are given by:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (A42)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (A43)

where [F.sub.fpci] = -[F.sub.pfci].

The governing equations for the particle-phase are given by:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (A44)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (A45)

where [F.sub.ppci] = 2/3(1 + e)[nabla]([C.sub.p[[rho].sub.p]<[c.sup.'.sub.p]2>), and [F.sub.pfci] = [[rho].sub.p]/[[tau].sub.rp] ([u.sub.pi] - [u.sub.fi]).

NOMENCLATURE

[c.sub.k] k-phase molecular or particle velocity (m/s)
[C.sub.k] k-phase volumetric concentration (g/[m.sup.3])
[d.sub.p] particle radius (m)
d[R.sub.k]
 ([x.sub.1],
 [x.sub.2],
 [x.sub.3]) a volume unit in geometric space ([m.sup.3])
[f.sub.k] k-phase molecular or particle velocity distribution
[F.sub.k] k-phase molecular or particle mass forces (N)
g gravitational acceleration (m/[s.sup.2])
G grid filter operator used in LES
G test filter operator used in LES
G GG
[L.sub.fij] fluid phase resolved stress
[m.sub.k] k-phase molecular or particle mass (kg)
[p.sub.k] k-phase pressure (Pa)
[q.sub.fij] fluid phase SGS stress (Pa)
[R.sub.fij] fluid phase rotation-rate tensor (1/s)
[r.sub.k] k-phase molecular or particle radius (m)
[R.sub.k] vector in geometric space
[Re.sub.p] particle relative Reynolds number
[S.sub.fij] fluid phase strain-rate tensor (1/s)
[T.sub.fij] fluid phase subtest-scale stress (Pa)
[u.sub.p] particle velocity (m/s)
[u.sub.f] liquid velocity (m/s)

Greek Symbols

[DELTA] subgrid filter width
[[micro].sub.f] or [micro] fluid viscosity (Pa s)
[v.sub.f] or v fluid kinematic viscosity ([m.sup.2]/s)
[[micro].sub.p] particle phase viscous coefficient (Pa s)
[[theta].sub.k] k-phase molecular or particle temperature (K)
[[rho].sub.k] particle density (1/[m.sup.3])
[[tau].sub.fij] fluid phase stress tensor (Pa)
[[tau].sub.rp] averaged relaxation time of particles (s)

Subscripts

c: collision
f: physical quantity of the fluid phase
i, j: tensor
k: physical quantity of a k-phase
p: physical quantity of the particle phase

Manuscript received November 19, 2005; revised manuscript received December 10, 2006; accepted for publication January 3, 2007.

REFERENCES

Aidun, C. K. and Y. Lu, "Lattice Boltzmann Simulation of Particles Suspended in Fluid," J. Statistical Phys. 81, 49-61 (1995).

Attride, M., "Enhancing Performance and Reducing Operating Costs for Slurry Pump," World Pumps 433, 33-35 (2002).

Baumgarten, S. and T. Muller, "Computer Aided Impeller Optimization," World Pumps 2000 (402), 28-31 (2000).

Boivin, M., O. Simonin and K. D. Squires, "On the Prediction of Gas-Particle Flows with Two-Way Coupling Using Large Eddy Simulation," Phys. Fluids 12, 2080-2090 (2000).

Canuto, V. M. and Y. Cheng, "Determination of the Smagorinsky-Lilly Constant," Phys. Fluids 9, 1368-1378 (1997).

Chapman, S. and T. G. Cowling, "The Mathematical Theory of Non-Uniform Gases," Cambridge University Press, Cambridge (1970).

Craft, T. J., B. E. Launder and K. Suga, "Development and Application of a Cubic Eddy-Viscosity Model of Turbulence," Int. J. Heat and Fluid Flow 17, 108-115 (1996).

Ding, J. X. and D. Gidaspow, "A Bubbling Fluidization Model Using Kinetic Theory of Granular Flow," AICHE J. 36, 523-538 (1990).

Ehrhard, J. and N. Moussiopoulos, "On a New Nonlinear Turbulence Model for Simulating Flows around Building-Shaped Structure," J. Wind Eng. Ind.l Aerodynamics 88, 91-99 (2000).

Enwald, H., E. Peirano and A. E. Almstedt, "Eulerian Two-Phase Flow Theory Applied to Fluidization," Int. J. Multiphase Flow Suppl. 22, 21-66 (1996).

Ferziger, J. H., "Large Eddy Simulation: An Introduction and Perspective," Springer-Verlag, New York (1996).

Ferziger, J. H., "Direct and Large Eddy Simulation of Turbulence," Trans. Japan Society Mech. Eng. (B) 66, 2-11 (2000).

Filippova, O. and D. Hanel, "Lattice-Boltzmann Simulation of Gas-Particle Flow in Filters," Comput. Fluids 26, 697-712 (1997).

Germano, M., U. Piomelli, P. Moin and W. H. Cabot, "A Dynamic Subgrid-Scale Eddy Viscosity Model," Phys. Fluids A 3, 1760-1765 (1991).

Ghosal, S., T. S. Lund, P. Moin and K. Akselvoll, "A Dynamic Localization Model for Large-Eddy Simulation of Turbulent Flows," J. Fluid Mech. 286, 229-255 (1995).

Ishii, M., "Thermo-Fluid Dynamic Theory of Two-Phase Flow," Eyrolles, Paris (1975).

Jaberi, F. A. and S. James, "A Dynamic Similarity Model for Large Eddy Simulation of Turbulent Combustion," Phys. Fluids 10, 1775-1777 (1998).

Kollmann, W., "Computational Fluid Dynamics," Hemisphere Publishing Corporation, New York (1984).

Laslandes, S. and C. Sacre, "Transport of Particles by a Turbulent Flow Around an Obstacle - A Numerical and a Wind Tunnel Approach," J. Wind Eng. Ind. Aerodynamics 74-76, 577-587 (1998).

Leith, C. E., "Stochastic Backscatter in a Subgrid-Scale Model: Plane Shear Mixing Layer," Phys. Fluids 2, 297-299 (1990).

Leonard, B. P., "A Stable and Accurate Convective Modeling Procedure Based on Quadratic Upstream Interpolation," Comput. Meth. Applied Mech. Eng. 19, 59-98 (1979).

Lilly, D. K., "A Proposed Modification of the Germano Subgrid-Scale Closure Method," Phys. Fluids A 4, 633-635 (1992).

Liu, D. Y., "Fluid Dynamics of Two-Phase System," China Higher Education Press, Beijing, China (1993).

Mathiesen, V., T. Solberg and B. H. Hjertager, "An Experimental and Computational Study of Multiphase Flow Behavior in a Circulating Fluidized Bed," Int. J. Multiphase Flow 26, 387-419 (2000a).

Mathiesen, V., T. Solberg and B. H. Hjertager, "Prediction of Gas/Particle Flow with an Eulerian Model Including a Realistic Particle Size Distribution," Powder Technol. 112, 34-45 (2000b).

Murakami, S., "Overview of Turbulent Models Applied in CWE," in "Proc. of the Second European and African Conference on Wind Engineering," G. Solari Ed., Genova, Italy, SGE, Padova, Vol. 1 (1997), pp. 3-24.

Ni, J. R., G. Q. Wang and A. G. L. Borthwick, "Kinetic Theory for Particles in Dilute and Dense Particle-Liquid Flows," J. Hydraulic Eng., ASCE 126, 893-903 (2000).

Olsson, M. and L. Fuchs, "Large Eddy Simulation of a Forced Semiconfined Circular Impinging Jet," Phys. Fluids 10, 476-486 (1998).

Pai, S.-I., "Two-Phase Flow," Vieweg, Braunschioeig (1977).

Poimelli, U. and J. H. Liu, "Large-Eddy Simulation of Rotating Channel Flows Using a Localized Dynamic Model," Phys. Fluids 7, 839-848 (1995).

Pope, S. B., "A More General Effective Viscosity Hypothesis," J. Fluid Mech. 72, 331-340 (1975).

Roco, M. C., "Wear Mechanisms in Centrifugal Slurry Pumps," Corrosion 46, 424-431 (1990).

Sinha, M. and J. Katz, "Quantitative Visualization of the Flow in a Centrifugal Pump with Diffuser Vanes-II: Addressing Passage Averaged and Large-Eddy Simulation Modeling Issues in Turbomachinery Flows," ASME J. Fluids Eng. 122, 108-116 (2000).

Smagorinsky, J., "General Circulation Experiments with the Primitive Equations," Mon. Weather Review 91, 99-165 (1963).

Tang, X., Y. Xu and Y. Wu, "Kinetic Model for Slit-Laden Solid-Liquid Two-Phase Flow," Acta Mechanica Sinica 34, 956-962 (2002).

Van Doormaal, J. P. and G. D. Raith, "Enhancement of the SIMPLE Method for Predicting Incompressible Fluid Flows," Numerical Heat Transfer 7, 147-163 (1984).

Vincenti, W. G. and C. H. Kruger Jr, "Introduction to Physical Gas Dynamics," John Wiley, New York (1965).

Walker, C. I. and G. C. Bodkin, "Empirical Wear Relationships for Centrifugal Slurry Pumps Part 1: Side-Liners," Wear 242, 140-146 (2000).

Wallin, S. and A. V. Johansson, "Modeling Streamline Curvature Effects in Explicit Algebraic Reynolds Stress Turbulence Models," Int. J. Heat Fluid Flow 23, 721-730 (2002).

Wang, Q. and K. D. Squires, "Large Eddy Simulation of Particle-Laden Turbulent Channel Flow," Phys. Fluids 8, 1207-1223 (1996a).

Wang, Q. and K. D. Squires, "Large Eddy Simulation of Particle Deposition in a Vertical Turbulent Channel Flow," Int. J. Multiphase Flow 22, 667-683 (1996b).

Wang, Z. Y. and N. Chien, "Experimental Study of Two-Phase Turbulent Flow with Hyperconcentration of Coarse Particles," Science In China, Ser. A. 27, 1317-1327 (1984).

Wang, Z. Y. and N. Chien, "Experimental Study of Laminated Load," Science In China, Ser. A. 28, 102-112 (1985).

Wilson, K. C., G. R. Addie and A. Sellgren, "Slurry Transport Using Centrifugal Pumps," Kluwer Academic Publishers, Dordrecht, Netherlands (1996).

Wu, Y. L., L. Ge and N. X. Chen, "Large Eddy Simulation of Silt-Liquid Two-Phase Flow Through a Centrifugal Pump Impeller," J. Tsinghua Univ. (Sci. Tech.) 41, 93-96 (2001).

Wu, Y. L., S. L. Cao, L. Ge, L. Wang, H. Tsukamoto, X. M.

Wang and Z. X. Sun, "Solid-Liquid Two-Phase Turbulent Flow Calculation and Measurement Through Centrifugal Slurry Pump Impellers," J. Tsinghua Univ. (Sci. Tech.) 38, 71-74 (1998).

Yan, Z., R. L. Street and J. R. Koseff, "A Dynamic Mixed Subgrid-Scale Model and Its Application to Turbulent Recirculating Flows," Phys. Fluids A 5, 3186-3196 (1993).

Zhang, Y., J. D. Wilson and E. P. Lozowski, "A Trajectory Simulation Model for Heavy Particle Motion in Turbulent Flow," ASME J. Fluids Eng. 111, 492-494 (1989).

Zhou, L. X., "Dynamics of Multiphase Turbulent Reacting Fluid Flows," National Defence Industry Press, Beijing, China (2002).

* Author to whom correspondence may be addressed. E-mail address: jingwu@adm.njit.edu

Xuelin Tang [1] and Jing Wu [2] *

[1.] College of Water Conservancy and Civil Engineering, China Agricultural University, Beijing, China 100083

[2.] Otto H. York Department of Chemical Engineering, New Jersey Institute of Technology, Newark, NJ U.S.A. 07102

Table 1. Various coefficients in Equation (13)

[C.sub.f1] [C.sub.f2] [C.sub.f3] [C.sub.f4] [C.sub.f5]
-0.2 0.4 [MATHEMATICAL 32 [C.sup.2. 0
 EXPRESSION NOT sub.f[micro]]
 REPRODUCIBLE
 IN ASCII.]

[C.sub.f1] [C.sub.f2] [C.sub.f6] [C.sub.f7]
-0.2 0.4 16 [C.sup.2. 16 [C.sup.2.
 sub.f[micro] sub.f[micro]

Table 2. Grid independence of the computational result for the duct
case; the mean experimental velocity is 2.389 m/s.

Grid Computational Computational Relative errors
 time (h) mean velocity of the
 (m/s) computational
 mean velocity (%)
350 X 35 X 20 44 2.225 6.86

400 X 40 X 26 58 2.375 0.59

450 X 45 X 30 87 2.395 0.25

Table 3. Parameters of the centrifugal pump

Flow rate Blade number Impeller outlet Specifi c speed
(m.sup.3]/min) diameter (m) (rpm, m, [m.sup.3]/
 min)

0.52 3 0.26 146

Table 4. Grid independence of the computational result for pump (the
experimental fl ow rate at the impeller outlet is 0.52 [m.sup.3]/s)

Computational Computational Computational The relative errors of
grid time (h) flow rate Computational flow
 ([m.sup.3]/s) rate (%)

120 x 71 x 41 60 0.485 6.73%

128 x 87 x 51 81 0.514 1.15%

134 x 95 x 61 111 0.517 0.58%