Access the full text.
Sign up today, get DeepDyve free for 14 days.
The flow of a mixture of liquid and solid particles at medium and high volume fraction through an expansion in a rectangular duct is considered. In order to improve the modelling of the phenomenon with respect to a previous investigation (Messa and Malavasi, 2013), use is made of a two-fluid model specifically derived for dense flows that we developed and implemented in the PHOENICS code via user-defined subroutines. Due to the lack of experimental data, the two-fluid model was validated in the horizontal pipe case, reporting good agreement with measurements from different authors for fully-suspended flows. A 3D system is simulated in order to account for the effect of side walls. A wider range of the parameters characterizing the mixture (particle size, particle density, and delivered solid volume fraction) is considered. A parametric analysis is performed to investigate the role played by the key physical mechanisms on the development of the two-phase flow for different compositions of the mixture. The main focuses are the distribution of the particles in the system and the pressure recovery. Keywords: Two-fluid model; Two-phase flow; Turbulent liquid-solid flows; Sudden expansion. INTRODUCTION Many engineering applications, of which the slurry pipelines are a significant example, involve the presence of solid particles with medium and high loading ratios in a turbulent liquid internal flow. The flow of solid-liquid mixtures is very complex. Doron and Barnea (1996) identified the flow patterns of slurries flowing in horizontal pipes as the velocity decreases: (a) fullysuspended flow, in which all the particles are suspended; (b) flow with a bed (moving/stationary), in which the particles accumulate at the pipe bottom and form a packed bed either sliding or fixed. We already stressed that the knowledge of these flows is far from being exhaustive (Messa and Malavasi, 2013). The majority of the researches concern straight pipes or rectangular ducts, and rely on the experimental approach (Kaushal and Tomita, 2007; Matousek, 2002; Vlasak and Chara, 2011; Vlasak et al., 2012). Most of the few CFD studies employ the Eulerian-Eulerian approach (Ekambara et al., 2009; Kaushal et al., 2012), the Eulerian-Lagrangian models being not applicable to dense mixtures due to their excessive computational cost; anyway, the Eulerian-Eulerian models, often referred to as twofluid models, appear numerically unstable, computationally expensive and their prediction do not often agree with the experimental evidence. These features, observed even for pipe flows, could complicate and even prevent their application to more complex scenarios like pipeline fittings. The literature material about the case of pipeline fittings is rather poor. The few experimental tests, extremely difficult to perform, are limited to rather low particle loadings and only rarely the experimenters performed measurements of flow parameters, such as mean and RMS velocity of the phases (Founti and Klipfel, 1998). Most of the numerical studies were aimed at predicting the erosion of pipeline fittings by coupling CFD and erosion models. Typically the simulations concerned dilute mixtures, and were performed using Eulerian-Lagrangian models (Habib et al., 2008); two-fluid models have been applied to solid-liquid flows through pipeline fittings in few researches only (Frawley et al., 2010; Mohanarangam and Tu, 2009); even less concern dense flows (Kaushal et al., 2013). In Messa and Malavasi (2013) we preliminarily addressed the flow of mixtures of water and mono-dispersed spherical glass beads through an upward-facing step, focusing on the effect of expansion ratio, particle size and delivered solid volume fraction on several flow parameters. We simulated a 2D slice of the system using the two-fluid model developed by Spalding (1980) and embedded in the PHOENICS commercial code. This model, numerically stable and easy to converge, had been successfully applied to rather dilute gas-solid flows (Marjanovic et al., 1999) and, by cursory comparison with experimental data concerning the straight pipe case, was found suitable to provide preliminarily information also for the dense liquid-solid flows under consideration. In this work, we improved the mathematical description of the flow by developing an original two-fluid model which specifically accounts for the physical mechanisms effective in case of liquid-solid flows at medium and high particle loading. We simulated a 3D system, accounting for the effect of the lateral vertical walls. The domain is sketched in Fig. 1: the heights of the upstream and downstream ducts are h = 26 mm and H = 39 mm respectively, corresponding to a step height D of 13 mm and an expansion ratio H/h of 1.5. The duct width B is equal to H. Compared to our previous investigation, we enlarged the range of variability of the parameters describing the composition of the mixture, focusing on configurations of interest for the applications in the field of mining industry. Fig. 1. Sketch of the expansion. The carrier fluid is water and three different values of particle density p (1450, 2650 and 4100 kg/m3), characteristic of coal, zinc, and iron tailings are considered; the diameter dp of the particles (assumed spherical and mono-dispersed) is 90, 175 and 280 m; the delivered solid volume fraction C ranges from 5 to 30%. The superficial velocity in the upstream duct Vs, equal to 4 m/s, is sufficient to avoid particle accumulation for all the flow conditions considered. The flow Reynolds number R = VsD/ ( is the kinematic viscosity coefficient of the fluid) is 52000. The theoretical framework of the problem is briefly illustrated as inferred from the literature available on this topic (which mainly concern single-phase flows or dilute fluid-solid flows) and on the results of our previous investigation (Messa and Malavasi, 2013). 5D V0 (3) where V0 is the mean velocity in the centerline of the smaller duct. Both the particle response to the mean flow and the turbulent dispersion increase as S decreases. Hardalupas et al. (1992), Fessler and Eaton (1999), and Founti et al. (1999), who dealt with dilute gas-solid and rather dilute liquid-solid flows through expansions along the vertical direction, all agree that the particles enter the recirculation region only if S < 1. Other significant physical mechanisms are gravity and particle-particle interactions. Gravity increases with particle mass, and contributes to keep the particles within the main core region of the flow, preventing them from entering the recirculation region. The importance of particle-particle interactions (either quick collisions or long-lasting contacts) increases with the number of particles, i.e. the delivered solid volume fraction. The remainder of this paper is divided in three sections. The first illustrates the theoretical background of the phenomenon; the second describes the original mathematical model; and the third reports the results of the numerical predictions. MATHEMATICAL MODEL Conservation equations and turbulence modelling The two-phase flow is represented by using the Eulerian approach in which the two phases are treated as interpenetrating continua. The flow is assumed statistically steady and so the mass conservation equation for phase k = f,p takes the following form: Fig. 2. Qualitative sketch of pressure and velocity. When a fluid flows through an asymmetric expansion, the flow separates and recirculation occurs downstream the step lip. The decrease of the fluid velocity is accompanied by an increase in pressure. Figure 2 shows a qualitative sketch of pressure and velocity downstream a two-dimensional expansion. The pressure in the system is typically quantified by the pressure coefficient Cp defined as: ( k k Uk ) = ( k Dk ) (4) Cp = p - pref 2 0.5 f Vref (1) being p the pressure, f the density of the fluid, and pref and Vref some reference values of velocity and pressure, respectively. Several physical mechanisms contribute to the flow of the mixture through the expansion. First, the responses of the particles to the mean fluid flow. Secondly, turbulent dispersion, that is how the particles are affected by the turbulent fluctuations of the fluid, and causes the particles to spread throughout the large duct, entering the recirculation region. An important dimensionless parameter affecting these mechanisms is the integral scale Stokes number S, which is the ratio of the particle response time p to a characteristic time scale of the flow . p is a characteristic timescale of the particle's reaction to changes in the velocity of the fluid, and can be evaluated as: where k is the volume fraction of phase k and D a phase diffusion coefficient, which appears in the phase diffusion term that represents the turbulent flux associated with correlations between fluctuating velocity u 'k and volume fraction 'k . These correlations, which appear in all conservation equations, are modelled in terms of a gradient-diffusion approximation in which the phase-diffusion coefficient D is given by t/, where t is the turbulent kinematic viscosity of the fluid and the turbulent Schmidt number for volume fractions, which may in some sense be interpreted as the ratio of turbulent momentum transport to the turbulent transport of phase mass. The origin of the correlations u 'k 'k in all conservation equations has been clarified elsewhere (Burns et al., 2004) and their modelling by means of a gradient diffusion approximation with diffusivity given by t/ is a well-known approach in literature. The mean global continuity is given by the equation that states that the two volume fractions must sum to unity. The momentum equations for phase k are: p = 4 ( p + Cvm f ) d p 3 f Cd U f - U p (2) ( kk Uk Uk ) = -kP + k (Tk +Tt,k ) + +kk g + Mk + ( k DUkk ) (5) being Up and Uf the mean velocities of particles and fluid, and Cvm and Cd the virtual mass and drag coefficients, which will be discussed later. The fluid time scale is based on an approximate large-eddy passing frequency in the separated shear layer: where: P is the pressure, shared by the phases; Tk and Tt,k are the viscous and turbulent stress tensors respectively; g is the gravitational acceleration; and Mk is the generalized drag per unit volume, which will be discussed later. The viscous stress tensor (present only in the fluid phase) and the turbulent stress tensor (present in both phases) are given by: T f = 2 f D f Tt , k = 2k t Dk (6) where Dk = 0.5[Uk+ (Uk)+] is the deformation tensor (the superscript "+" indicates that the transpose of the dyadic Uk is taken). The generalized drag term Mk accounts for the momentum transfer between the phases. Under the assumption of monodispersed spherical particles on which the model relies, this term is given by: Mf = -Mp = 6p d3 p ( Fd + Fl + Fvm ) (7) in which  is the intrinsic viscosity and pm is the maximum packing volume fraction, which are often treated as calibration parameters. Replacing Rp with Rm in Eq. (11) improves the predictive capacity of the model for dense flows. The asymptotic behavior of m, which tends to infinity as p pm, sets an upper limit to the solid volume fraction, preventing the particles from over-packing. This avoids the need to introduce a collisional pressure term in the momentum equation of the dispersed phase. The definition of the particle Reynolds number with respect to the viscosity of the mixture and the inclusion of lift and virtual mass forces are the main features which distinguish the current model from that of Spalding (1980) used in our previous work (Messa and Malavasi, 2013). The following modified form of the k- RNG model is used for evaluating the turbulent kinematic viscosity of the fluid phase, t: where Fd, Fl, and Fvm are the drag, lift, and virtual mass forces, calculated respectively as follows: 1 d Fd = 2 4 2 p Cd f U p - U f ( U p - U f ) t = C k2 (13) (8) (9) Fl =Cl f d 3 ( U p - U f ) × ( × U f ) p ( f f U f k ) = f f + t k + k + f f ( Pk -) + f Dk f ( f f U f ) = f f + t + f f C1 Pk + k 3 C (1 - / 0 ) -C2 - + f D f 1 +3 (14) 4 d3 p Fvm = Cvm f ( U p U p - U f U f ) 3 8 (10) in which Cl is the lift coefficient, set like Cvm to 0.5 as per the indications reported by Kaushal et al. (2012). The Basset force was not included in the model as it proved ineffective for the flows addressed in this work (Chung and Troutt, 1988). Some authors have argued for the existence of a near-wall lift force to account for the repulsion of particles from the pipe wall observed in some experiments (Kaushal and Tomita, 2007), but this effect is not considered in the present work. A semitheoretical model for this force was derived by Antal et al. (1991) for air-water bubbly flow in the laminar regime, but it proved unsuitable for slurry flows in horizontal pipes, confirming the observations of Ekambara et al. (2009). Wilson et al. (2010) proposed a model for the near-wall force in slurry flows, but the global nature of its formulation precludes its implementation in a CFD code. Anyway, the absence of this force is not expected to affect the overall quality of the solution for the flow conditions considered here except marginally for the two cases in which the particles are larger than 90 m. The drag coefficient is given by the well-known formula of Schiller and Naumann (1935): 24 C d = max (1 + 0.15 R 0.687 ) , 0.44 p R p (11) (15) in which: Pk=2tDf :Uf is the volumetric production rate of k due to the working of the Reynolds stresses against the mean flow; and is equal to Sk/ being = 2 : . The usual values of the model constants are employed, namely k = 0.7194, =0.7194, C=0.09, C1=1.42, C2=1.68, 0 = 4.38, and =0.012. Fig. 3. Computational domain and boundary conditions. Computational domain and boundary conditions The computational domain is shown in Fig. 3; the flow and geometrical symmetry of the phenomenon has been exploited by solving only over one half of the duct section. A fully-developed turbulent flow distribution is specified at the inlet section, with the distribution of axial velocity, turbulent kinetic energy and dissipation rate determined from the Nikuradse's boundary layer theory (Schlichting, 1960) for single-phase flow. No slip is allowed between the phases at the inlet section, thus the same velocity distribution is applied to both phases. The fully-developed volume fraction distribution obtained from the simulation of a straight duct is imposed at the in which Rp is the particle Reynolds number, defined as dp|UpUf|/. In order to account for the presence of multiple particles, Rp is replaced by a mixture-based particle Reynolds number Rm= fdp|Up-Uf|/m, where m is a characteristic viscosity of the mixture obtained from the correlation of Mooney (1951): [ ] p m = f exp 1- / p pm (12) inlet section. At the outlet, the normal gradients of all variables and the value of the pressure are set to zero. The inlet and the outlet boundaries are located 8D and 200D upstream and downstream the expansion, in order to analyze the development of the flow downstream the step. The proper wall boundary condition is still an open research problem. In Messa et al. (2013, 2014) we addressed horizontal pipe flows and found that the application of the equilibrium wall function of Launder and Spalding (1972) to model the velocity of the two phases, the turbulent kinetic energy and its dissipation rate in the near-wall cells procures reliable predictions of the viscous and mechanical contributions to friction (Shook and Bartosik, 1994). In order to generalize this result to a flow in which the equilibrium assumption is not satisfied, in this work we replaced the equilibrium wall function with the non-equilibrium one (Launder and Spalding, 1974). The flow variables are not resolved all the way to the walls, but only up to a distance from the walls comprised between 0.5 and 1 mm, and then modeled in the near-wall region. Governing terms and validation of the model Many terms of the two-fluid model contain adjustable coefficients which are in practice treated as calibration parameters, their values being determined by matching computations with experiments. Due to the lack of experiments regarding expansion flows, these coefficients were determined by reference to the horizontal pipe case, virtually the only scenario in which measured data for dense slurries are available. It was found that the key parameters affecting the numerical solution are the turbulent Schmidt number for volume fractions in the phase diffusion coefficient D and the intrinsic viscosity  and the maximum packing volume fraction pm in the Mooney's formula for the viscosity of the mixture (Eq. 12). Messa (2013) and Messa et al. (2014) performed an extensive comparison with a large dataset of measurements from different authors, and found that the triplet = 0.7,  = 2.5, and pm = 0.70 (within the range of variability usually encountered in literature) produces good agreement with the experiments in case of fully suspended flows except close to the pipe wall for those flows in which the already mentioned near-wall lift plays an important role. The comparison with additional data from Kaushal and Tomita (2007) and Nabil et al. (2013), which are very relevant in the context of the present study, is reported in Appendix A. Starting from these results, we performed a sensitivity analysis and found that lift and virtual mass had to be accounted for in the expansion case despite being ineffective in the simpler horizontal pipe one. Computational methodology numerical solution and consistency of the variables to achieve convergence. Inertial relaxation is applied to the momentum equations with a false-time step of 104 s. A linear relaxation factor of 0.4 is applied to all other flow variables. A Cartesian structured mesh is used to discretize the domain. The grid employed consists of 16 by 39 by 795 cells along the x, y, and z directions respectively; a grid independence study, illustrated in Appendix B, proved the suitability of this mesh. The PHOENICS solver was run until the sum of the absolute residual sources over the whole solution domain are less than 1 per cent of reference quantities based on the total inflow of the variable in question. An additional requirement is that the values of the monitored dependent variables at a selected location do not change by more than 0.1% between subsequent iteration cycles. RESULTS AND DISCUSSION The flow field obtained as output of the 3D simulations is much more complex compared to the qualitative 2D solution depicted in Fig. 2, especially in the recirculation region close to the side walls. Fig. 4 shows the projection of the fluid phase vectors on different planes for a specific flow condition. The main vortex in the streamwise direction, induced by the step, is accompanied by secondary motions which tend to push the particles away from the side walls. The effect of these walls on the flow in the central part of the duct becomes less significant as the duct width B increases; in particular, the 2D model considered in Messa and Malavasi (2013) becomes no longer representative of the flow if the aspect ratio of the large duct, B/H, is below 12 (Messa, 2013). Fig. 4. Recirculation downstream the step for C = 5%, dp= 90 m, and p = 2650 kg/m3. The development of the two-phase flow, and in particular the distribution of the solids, is strongly affected by the configuration of the mixture, in terms of delivered solid volume fraction C, particle size dp, and particle density p. A parametric investigation was performed with reference to the flow conditions summarized in Table 1. The integral scale Stokes number S was calculated as the ratio between p (evaluated by Eq. (2) in which Cd, Uf and Up were taken as some characteristic values of the recirculation region) and the characteristic fluid time scale given by Eq. (3). Figure 5 shows the solid volume fraction profiles along the reference sections depicted in Fig. 6. In particular, Figs. 5(ac) The CFD code PHOENICS was used for the numerical solution of the finite-volume analogue of the mathematical described above. This was done by using the Inter-Phase Slip Algorithm (IPSA) of Spalding (1980), together with userdefined functions and subroutines for implementation of specific closures and boundary conditions. The calculations are performed following the elliptic-staggered formulation in which the scalar variables are evaluated at the cell centers and the velocity components at the cell faces. Central differencing is employed for the diffusion terms, while the convection terms are discretized using the Linear Upwind scheme. The numerical solution procedure requires appropriate relaxation of the field Table 1. Flow conditions considered in the parametric investigation. Case E1 E2 E3 E4 E5 E6 E7 p [kg/m3] 2650 2650 2650 2650 2650 1450 4100 dp [m] 90 90 90 175 280 90 90 C [%] 5 15 30 5 5 5 5 S 0.09 0.09 0.09 0.35 0.89 0.05 0.14 focuses on the effect of increasing C for a given kind of particles (dp = 90 m; p = 2650 kg/m3). The essentially uniform profiles reveal that, whatever the solids loading is, the particles tend to spread throughout the large duct entering the recirculation region. The low value of S, in addition to the negligible effect of gravity, contributes to explain the results. Three values of dp (90, 175, and 280 m) are instead depicted in Figs. 5(df) for a unique set of p (2650 kg/m3) and C (5%). The variability of the profiles reveals the different behavior of the system. Increasing dp increases significantly the effectiveness of gravity, which tends to push the solids downwards. Therefore, the small particles (case E1) are almost uniformly distributed throughout the large duct, whilst the large sized ones (case E5) tend to stay in the core region, resulting in almost zero solid volume fraction within the recirculation region. For the medium particles (cases E4) the counteracting effects of gravity, interphase friction and turbulent dispersion determine an intermediate configuration, in which the volume fraction of particles within the main vortex is greater than zero but lower than the mean inlet value. Since S < 1 for case E5, the absence of particles within the recirculation region may appear inconsistent with the experimental observations of Fessler and Eaton (1999), Founti et al. (1999), and Hardalupas et al. (1992). However, the different orientation of the system removes this apparent contradiction. In fact, these authors investigated sudden expansions in vertical pipes; therefore, no stratification is induced by gravity, resulting in a flat particle distribution in the small duct. This, in turn, produces higher solids concentration along the border of the recirculation region, promoting the entrainment of particles in that region. Figs. 5(gi) refer to the flow conditions E1, E6, and E7, in which p varies assuming values characteristic of characteristic of coal tailings (1450 kg/m3), zinc tailings (2650 kg/m3), and iron tailings (4100 kg/m3). The numerical results indicate that the gravitational stratification becomes more and more effective as p increases, but the influence on these parameters appears minor compared to dp. Fig. 5. Solid volume fraction profiles along the sections Y1 to Y3 depicted in Fig. 6 for the flow conditions summarized in Table 1. Fig. 6. Identification of the reference sections. insensitive to dp, probably as a consequence of the two counteracting effects: on one hand, the bigger the particles are, the higher is the momentum they transfer to the fluid (higher pressure recovery); on the other hand, the bigger the particles are, the less is the number of solids which enter the recirculation region due to the increased importance of gravity (lower pressure recovery). At the limit, if all the particles would accumulate at the bottom of the duct one would expect the pressure recovery to be very close to that of pure liquid. At last, Fig. 7(c) focuses on the influence of p (1450, 2650, and 4100 kg/m3). The pressure recovery increases with p, due to the higher momentum transferred from the particle to the fluid. The greater influence of p compared to dp may be explained considering that, even for the heaviest particles the gravitational stratification appears rather limited (Fig. 5(g-i)); consequently, only the former of the two above-cited effects play a significant role. CONCLUSIONS This work focuses on the flow of liquid-solid slurries with medium and high solids loading through a sudden expansion in a horizontal duct. Starting from a preliminary investigation (Messa and Malavasi, 2013), we improved the mathematical description of the two-phase flow by making use of an original two-fluid model that we developed for dense turbulent liquidsolid flows and implemented in the PHOENICS CFD code. Since no experimental data regarding the flows addressed here are available in literature, the model was validated in the horizontal pipe case, reporting good agreement with the experiments from different authors over a wide range of flow conditions. We simulated a 3D model of the expansion, accounting for the crosswise recirculation induced by the side walls (Fig. 4). We considered a wide range of variability of the parameters regarding the composition of the mixture, in terms of delivered solid volume fraction (from 5 to 30%), particle size (from 90 to 280 m), and particle density (1450 to 4100 kg/m3). A parametric analysis was performed to investigate the role played by the key physical mechanisms (particle response to the mean fluid flow; turbulent dispersion; gravity; particle-particle interactions) on the development of the flow as the composition of the mixture changes. If the particles are sufficiently fine, they are very likely to follow the fluid and, whatever the solids loading is tend to spread throughout the large duct entering the recirculation region (Fig. 5(ac)). Conversely, as the particles size and, to a lesser extent, the particle density increase, gravity becomes more and more effective causing the solids to remain in the lower part of the duct (Figs. 5(di)). The pressure recovery downstream the step was found to be higher than the single-phase case and increase with the solids loading and the particle density, whilst the influence of the particle size seems minor (Fig. 7). An interpretation of these results was given in terms of role played by the above mentioned key physical mechanisms, considering that the additional momentum transferred from the particles to the fluid depends on the size and weight of the particles but also on the number of particles in the recirculation region. Acknowledgements. The authors would like to thank CHAM Ltd for contributing to supporting this research. They would also acknowledge the CINECA and the Regione Lombardia award under the LISA initiative, for the availability of high performance computing resources and support. Fig. 7. Distribution of C p downstream the expansion for the flow conditions E1 to E7 in Table 1. The single phase (SP) case is depicted too. The pressure distribution in the system is now considered. As in single-phase flow (Fig. 2), an increase in pressure occurs as a consequence of the decrease in slurry velocity. A parametric investigation aimed at analyzing the way in which C, dp and p affect the pressure distribution is performed. Reference is made to the pressure coefficient Cp defined in Eq. (1) with pref = 0 and Vref = Vs. Figure 7(a) compares the results for a unique set of particles (p = 2650 kg/m3; dp = 90 m) with C ranging from 5 to 30%. The plot indicates that the pressure recovery in the two-phase case is higher compared to the single-phase case and increases with C, possibly as a consequence of the momentum transferred from the particles to the fluid. Similar observations were reported for gas-solid flows (Tomita et al., 1980). Three values of dp (90, 175, and 280 m) are depicted in Fig. 7(b) for p = 2650 kg/m3 and C = 5%. The pressure recovery seems rather
Journal of Hydrology and Hydromechanics – de Gruyter
Published: Sep 1, 2014
Access the full text.
Sign up today, get DeepDyve free for 14 days.