Numerical assessment of the vulnerability to impact erosion of a pump as turbine in a water supply system

Water transmission, distribution and drainage systems are generally high energy demanding and scarcely efficient. Some studies classify the water industry as the fourth most energy-intensive industrial sector in the UK (Ainger et al. 2009), with 5 tons of greenhouse gas emissions per year. In California, the energy requested for the pumping of water represents of the total energy, being the most energy demanding industrial activity (Lofman et al. 2002). In the US, the energy demand of the whole water utility represents of the total energy consumed with 45 tons of emissions of greenhouse gas and an expense of 4 billion dollars per year (Tunnicliffe 2009). Furthermore, water distribution and transmission systems are affected by high levels of water leakage. Olsson (2015) assessed the total leakage as 45–88 million m³ of water per day, while Cabrera et al. (2010) stated that 40–60% of the energy used to distribute water is lost due to leakage and head losses. Additionally, the strategy generally pursued to reduce the water leakage consists in the pressure control using pressure-reducing valves (PRVs). Indeed, the relation between the water leakage and the water pressure inside the pipe has been largely demonstrated. Thus, the use of a valve to reduce the water pressure is an effective way to control the amount of leakage. Nevertheless, each valve produces a dissipation of energy that contributes to reduce the efficiency of water systems. In this context, some authors proposed to replace the valves with turbines, to both regulate the excess pressure and generate small-scale hydropower. But, even if the overall recoverable power is large, it is generally distributed in many low-power spots. Thus, the use of conventional turbines is not convenient, due to their high design and production costs (Pugliese et al. 2018).

The use of pumps as turbines (PATs) has been presented as a valid choice to pursue the pressure control and the energy recovery with low investment and maintenance costs for pico- and micro-hydropower plants (Novara et al. 2019). PATs are pumps running in the reverse mode by using their asynchronous motor as an electrical generator. Their application has been proposed in both water supply systems and wastewater treatment plants (Power et al. 2011). The large availability of pumps in the market, their simple mechanics, the low purchasing and maintenance costs, and the ease of installation are great advantages that drastically reduce the payback period of the hydropower plants. Nevertheless, the limitations in their use arise from two main issues: (i) the variability of the hydraulic properties of the system (flow rate, head and power) and (ii) the scarce knowledge about the performance of pumps running in the reverse mode.

Unlike the traditional turbines, a PAT has no regulation systems; thus, many authors proposed different design solutions for effective regulation. In this paper, the PAT plant features a hydraulic circuit and the regulation is pursued by two control valves to either increase the head loss or bypass a part of the discharge (Fecarotta et al. 2018b).

The lack of knowledge about the PATs performance occurs because pump manufacturers do not usually provide the performance curves of their machine in the reverse mode. Some semi-empirical models have been developed to predict the behaviour of a PAT when the performances in the pumping mode are known, but their reliability is low. Thus, even if they can be used for preliminary analysis, a proper design of the plant generally requires experimental validation (Pugliese et al. 2018).

A valid alternative to both semi-empirical models, that are poorly reliable, and experimental tests, which surely require much time and expensive costs, is computational fluid dynamics (CFD). Even if CFD models require the complete CAD 3D geometry of all the hydraulic components of the machine, they allow not only the prediction of the performance of the machine but also the complete characterization of the internal flow field (Frosina et al. 2017). Thus, an extensive analysis of the secondary flows, of the internal head losses and of the generated force and torque can be pursued.

The assessment of the material removal from a PAT produced by the impingements of solid particles against the inner surfaces of the device has not yet been deeply investigated. This phenomenon is termed differently in the literature. When talking about hydraulic turbines and slurry pumps, reference is often made to the more general terms ‘hydro-abrasive erosion’ and ‘erosion wear’, respectively. Strictly speaking, these terms indicate the degradation caused by solid particles in water, which includes both the impact erosion and the abrasive action due to particles sliding against curved walls (Pagalthivarthi et al. 2011; Felix et al. 2016). Although the present numerical investigation is mostly oriented to the assessment of the damage produced by particle-wall impingements, for the sake of simplicity – at the risk of some over-simplification – the terms ‘impact erosion’, ‘hydro-abrasive erosion’, ‘erosion wear’ and ‘solid particle erosion’ will be used indifferently in the paper.

Having said that, hydro-abrasive erosion is a heavy concern whenever a PAT works with untreated water, e.g. when the device is applied in wastewater treatment plants or before the inlet of a water purifying plant. Another possible situation in which a PAT might suffer from hydro-abrasive erosion is when this device is used in an irrigation system (Crespo Chacón et al. 2020), which generally delivers untreated water to the end users. Extending the view to water turbines used in hydropower plants helps to better understand how the presence of particles can severely reduce the lifetime of hydraulic machines, with significant economic implications. Review papers by Padhy & Saini (2008), Neopane et al. (2011), Thapa et al. (2012) and Felix et al. (2016) suggest that different types of turbines – e.g. Pelton, Francis and Kaplan – are affected by wear in several rotating and fixed parts, which were subject of specific investigations. Hydro-abrasive erosion is a concern if either of the two conditions is met, namely a high suspended sediment concentration and the presence of hard mineral particles. Similar issues are commonly found in mountain rivers in the Himalayas, especially due to the heavy rains occurring during the monsoon period (Bajracharya et al. 2008), as well as in the European Alps and the Andes (Felix et al. 2016). Considerable research work has been done on this topic. Several case studies were presented to demonstrate the effect of wear on turbines installed in real hydropower plants (Bajracharya et al. 2008; Abgottspon et al. 2016). Laboratory studies on small-scale prototype turbines were also made and predictive correlations derived (Padhy & Saini 2009, 2011), but the size-up scalability of the predictions is the most limiting factor to their applicability. Being free – at least in principle – from scaling concerns, CFD can be a valid approach to assess the vulnerability to wear of a hydraulic machine. However, modelling this type of phenomena is particularly challenging since the difficulties related to the description of the impact erosion combine with those due to the modelling of the particle-laden flow in a rotating system. The rotation of the impeller, in fact, is expected to produce a considerable impact on the particle dynamics, as suggested by some literature experiments (Deng et al. 2008) on a centrifugal erosion tester. All these difficulties might explain why CFD studies devoted to the erosion assessment in hydraulic machines are so limited in number and generally based on a combination of empiricism and simplifications in the modelling of the particle-induced damage. Krüger et al. (2010) simulated the quartz sand impingement in a centrifugal pump, observing significant erosion areas along the blade and the impeller trailing edge, depending on the particle impingement angle and solid concentration. Lopez et al. (2015) investigated how the erosion prediction in one of the channels of a pump impeller was affected by the modelling approach and some of the parameters of the CFD-based wear estimation model. All the considered model formulations yielded the same locations of the maximum wear areas, placed around the blades, but different predictions in terms of the amount of removed material. Later, other authors (López et al. 2018) focused their attention on the modelling of erosion by accounting for the changes in the particle-laden flow caused by the eroding geometry, and with the aim of demonstrating the capability of the method they implemented in OpenFoam, they briefly presented an application case study regarding the wear damage of the volute of a slurry pump. A considerable amount of research work was done on the prediction of the wear occurring in the volute of slurry pumps by other authors (Pagalthivarthi et al. 2011, 2013, 2017). The authors were mostly interested in assessing the effect of modelling, design and operation parameters on the risk of failure of the volute, and they simulated only the fixed casing by accounting for the impeller rotation only indirectly. In most cases, they inferred the erosion from the fluid dynamic characteristics of the slurry flow field via empirical coefficients. Few other studies focused on the internal erosion of turbines. Among them, the recent works by Leguizamón et al. (2019a, 2019b) provided a significant contribution not only to the specific application but also to the erosion modelling. These authors proposed a multiscale approach to erosion prediction, which combines fluid dynamics and solid mechanics to get closer to the physics of the erosion process. The physically based approach is a direction that needs to be sustained in the future but, in its current formulation, the multiscale approach appears to be too computationally expensive to represent a real effective tool for engineering design and management. At this stage, engineering erosion prediction in complex geometries appears still dependent on empirical correlations, and therefore, much effort has to be devoted to guaranteeing the reliability of the information arising from the CFD simulations, particularly when no field or experimental data are available for model validation.

As a final note, it should be mentioned that in real applications, cavitation and corrosion might act in synergy with impact erosion enhancing the deterioration of a hydraulic machine. The synergistic effect of solid particle erosion and these other sources of damage has been the subject of several studies (Wang et al. 2009; Gohil & Saini 2014); nevertheless, at present, there does not appear to exist a well-established mathematical model for its description, although valuable steps in this direction were made (Rajahram et al. 2009). For this reason, only the material removal produced by impinging solid particles will be considered in the present study.

Indeed, this paper, driven by the above considerations, provides a comprehensive analysis of the vulnerability to wear of a PAT working in variable conditions. Up to the authors' best knowledge, no study regarding this type of hydraulic machine has been previously published, except for a preliminary exploration performed by Fecarotta et al. (2018a), which constitutes the basis of this work. A CFD model was implemented by the University of Naples ‘Federico II’ to predict the behaviour of a single-stage horizontal axis centrifugal PAT. The performance curves have been applied to design a PAT plant for a pressure-reducing station in Southern Italy: the efficiency of the plant has been optimized considering both the flow rate and the head variability. Then, the CFD model has been extended to simulate the particle-laden flow inside the machine and, coupled with the in-house numerical software E-CODE developed by Politecnico di Milano, to assess the dependence between the vulnerability to wear and the flow rate. The overall erosion of the PAT due to its variable working conditions has been discussed. Finally, the influence of the characteristics of the solid phase on the wear behaviour has been explored at the best efficiency point (BEP). The focus of this work was not as much as to quantify material removal. Conversely, it will be shown that the used CFD-based methodology allows assessing how the vulnerability of the PAT to erosion is affected by the working conditions, providing justification of the results in the light of the fluid dynamics characteristics of the particle-laden flow.

Generally, a PRV is used to set the downstream pressure to a selected value, namely the backpressure (BP), independently of the upstream pressure. Thus, the PRV generates an energy dissipation that can be converted into energy recovery/production by the use of a hydropower plant equipped with a PAT. The hydropower plant must not modify the regulation of the water system, ensuring the same service to the end users (Fecarotta & McNabola 2017). Thus, the PAT must be provided with a regulation system in order to set the same BP. In this work, a hydraulic regulation system has been considered, as it ensures high efficiency and effectiveness values and a simple regulation procedure (Fecarotta et al. 2018b), namely the variable operating strategy (VOS). The PAT is inserted in a series–parallel circuit with two regulating valves (Figure 1), running at constant rotational speed, N. Thus, the relation between the turbined flow rate, ⁠, and the produced head drop, ⁠, is described by a strictly increasing function, usually assessed as a second-order polynomial. On the other hand, the available jump, ⁠, is intended as the difference between the upstream pressure head and the BP. Since the total discharge flowing into the system, ⁠, is variable in time, also the available jump is variable and due to the flow resistances, it decreases at increasing discharge. For the smallest values of flow rate, when is smaller than ⁠, the series valve generates an additional dissipation to match the BP value downstream of the plant. On the other hand, when the total flow is high and the turbine would produce a head drop higher than the available one, the parallel valve partially bypasses the discharge (⁠⁠) to limit the turbined head to the required value. A small scheme of the hydraulic circuit and a chart with the effect of the regulation are shown in Figure 1.

CFD is widely applied to numerically reproduce the kinematic and dynamic flows field across turbomachines because it is able to predict secondary flows and physical phenomena such as cavitation, internal losses and fluid–structure interactions. CFD modelling is usually based on the numerical solution of the Reynolds-averaged Navier–Stokes (RANS) equations both in the steady and transient state.

The RANS equations are a time-averaged formulation of the Navier–Stokes equations, in which both the velocity and the pressure p are expressed as a sum of a time-averaged value and a fluctuating one (Reynolds decomposition). The RANS are not a closed system of the equation; therefore, they need to be solved together with a turbulence model.

Among the many turbulence models available in the literature, the standard k–ɛ model (Launder & Spalding 1974) is widely applied to industrial flows because of its ability to provide good accuracy with relatively low computational burden. Due to these advantages, the standard k–ɛ turbulence model has been used in this paper to solve the flow field within the PAT. Although it has been suggested not being the best choice for simulating the flow inside reaction turbomachines, especially at off-design conditions, i.e. for working points far from the BEP of the machine, it has been successfully used in many literature studies (Derakhshan & Nourbakhsh 2008; Frosina et al. 2017). In this study, the validation against the experimental tests was taken as a criterion to justify the choice of the turbulence model. Indeed, as shown below, the experimental results describe only global variables, and therefore, the simulated fluid dynamic field inside the PAT may still be affected by inaccuracies arising from the chosen turbulence model. In order to further validate the results, some calculations have been compared with the results of the SST k–ɛ, revealing certain but not essential differences in terms of pressure and velocity distribution.

With the reference to CFD applications to turbomachines, two main approaches are applied: the frozen rotor and the transient one. The frozen rotor allows steady-state prediction because the relative position of the rotor with respect to the stator is frozen in time. Thus, the flow from the rotor to the stator is treated by moving the reference frame (MRF), whereas the relative position of components is fixed in time. On the contrary, the transient approach allows for the assessment of the transient effects, deriving from the interaction between the rotating part (the rotor) and the fixed one (the stator), to be performed through a sliding mesh (SM) technique. Indeed, the latter is intended as a dynamic mesh, reproducing a rigid motion in a given dynamic mesh zone. For further information about the two approaches, the reader is referred to the ANSYS Fluent User Guide (2013).

The procedure for estimating the erosion of the PAT was developed starting from the two-stage strategy usually employed in CFD-based wear prediction studies. First, the slurry flow was simulated by an Eulerian–Lagrangian model (Crowe et al. 2011), in which the fluid flow is represented in an Eulerian, cell-based framework, while the solid phase is modelled by following a certain number of particle trajectories. Afterwards, an erosion model (Parsi et al. 2014) was employed to attribute to each particle-wall impingement a certain mass of removed material, ⁠. Erosion models of this type are often referred to as ‘single-particle’ erosion models, and they take the form of algebraic correlations relating the to the particle mass, ⁠, and all other parameters involved in the wear process. As it has been well established for long time (Finnie 1972; Humphrey 1990), these parameters which do not include only the mechanical properties of the abrasive particles and the target material but also variables related with the fluid dynamic behaviour of the particle-laden flow, primarily the particle impact velocity, ⁠, and the particle impact angle, (Figure 2). Although some single-particle erosion models have some theoretical fundamentals (e.g. Finnie 1960), most of them are empirically obtained by fitting results of laboratory tests.

In order to achieve simulation results within acceptable time limits, the particles were treated in a steady fashion, in the sense that the ‘frozen’ steady-state fluid flow field of the single-phase MRF simulation was used to solve the particle equation of motion (Equation (4)) sequentially for each computational particle. As already noted, the MRF simulation comprised an outer zone and an inner zone that rotates with the impeller, both simulated as fixed in time. Since the fluid velocity vectors used to integrate Equation (4) referred to as two different coordinate systems, the computed paths did not correspond to the actual trajectories of the particles; nevertheless, for the sake of simplicity, the terms ‘particle paths’ and ‘particle trajectories’ will be used interchangeably. In the same way, also the Lagrangian integration time was not the physical time. Following the ‘parcel approach’ (Crowe et al. 2011), although each path was calculated from the balance of the forces acting on a single particle, the same path was used to represent the behaviour of a group of particles and it was associated with a certain solid mass flux. Although approximated, the steady-state approach – definitely less time-consuming compared to an unsteady simulation with SM – appeared a reasonable way to estimate the particles' fluid dynamic characteristics at the impact stage, which are required as input in erosion calculation. In fact, the particles recirculate around the impeller and hit it several times before leaving the inner zone of the domain. As a result, even if the impeller was modelled as fixed, the erosion was uniformly distributed over its circumference (Figure 3(a)), as it would be expected in a real behaviour. Note that the situation would have been different if the particles were likely to cross the inner zone quickly and without recirculation. In this case, modelling a fixed impeller would have resulted in a preferential location for particle-wall impingement and, consequently, an unphysical erosion pattern, non-uniform over its circumference (Figure 3(b)). Thus, an important novelty of this study is that the axisymmetry of the impeller wear was used as a criterion for the applicability of the steady-state MRF approach in erosion prediction simulations.

A single-phase CFD model was applied to assess the fluid dynamic field of a centrifugal 6 blade PAT with an impeller of diameter D = 189 mm. The model was implemented by using the ANSYS^® framework, from a 3D scanning procedure of the machine. Two connection pipes, located upstream and downstream of the PAT, were included in the computational domain for consistency with the imposed inlet–outlet boundary conditions. The domain, therefore, consisted of four regions, namely inlet pipe, volute, impeller and outlet pipe (the nomenclature inlet and outlet refer to the PAT operation, so that the inlet pipe and the outlet pipe are, respectively, upstream and downstream of the impeller during turbining operation).

A structured hexahedral mesh was generated with the ANSYS^® ICEM CFD^™ tool for connection pipes, whereas, due to their complex configuration, an unstructured tetrahedral mesh was applied for the impeller and volute components (Figure 4). Specific refinements were also introduced around both the volute tongue region and the impeller blades, aimed at increasing the mesh regularity, and at the interface between the impeller and the volute. Although a grid-independence study was performed at this stage of the work, only a little information about it is given here for the sake of brevity. The convergence analysis of the CFD model was performed by simulating the BEP at different rotational speed N_T = 1,200, 1,500, 2,100 and 2,900 rpm, respectively, for four different meshes (0.6 M, 0.9 M, 1.2 M and 2.5 M elements, respectively). The head drop generated by the PAT was considered as the reference parameter for assessing the model convergence. The mesh chosen to assess the PAT hydraulic behaviour is the one composed of 1.2 M elements. A relative scatter not higher than 4% between the selected mesh resolution (1.2 M elements) and the finer one (2.5 M) was observed during simulations, against computational efforts of about three times shorter. The good agreement with the experimental results further confirmed the reliability of the CFD results. Indeed, at this stage of the work, possible inaccuracies in the calculated flow field could be accepted provided that they do not affect the predicted values of net head and mechanical power, which are the only parameters needed for the application of the VOS and the evaluation of the recovered energy. Conversely, higher accuracy in the simulation of the fluid flow field will be required when evaluating the wear due to the impact erosion. Thus, a specific section will be dedicated to illustrate the grid-independence of the wear estimates.

The ANSYS^® Fluent™ code was employed for numerically solving the RANS equations coupled with the standard k–ɛ turbulence model, with scalable wall functions for the near-wall treatment. Second-order upwind spatial discretization was applied to all solved variables, and low relaxation factors were used to limit instabilities during the convergence procedure. Three planar interfaces were introduced to simulate the interaction between inlet pipe and volute, volute and impeller, and impeller and outlet pipe, respectively. As boundary conditions at the inlet and outlet, a uniform velocity distribution and a static pressure value were set, respectively.

A two-step computational procedure was followed to calculate the periodic flow field:

• Steady-state conditions, by applying the MRF approach.

• Unsteady-state conditions, applied once the convergence of the steady-state mode was achieved. The velocity and pressure fields from the steady-state mode were applied to initialize the simulations, by reproducing the relative motion between impeller and volute through the SM technique. The time step duration was evaluated so that a complete impeller revolution was performed in 150 time steps and not less than 5 complete impeller revolutions were simulated, useful to achieve the periodicity of flow parameters.

The head drop generated by the PAT, ⁠, was estimated as the difference between the inlet and the outlet pressure, whereas the overall generated power, ⁠, was calculated as a function of the torque on the impeller, estimated from the integration of the pressure and shear stress distributions on the impeller surfaces. The PAT characteristic curves for different rotational speeds were thus derived, in the range = 600–2,900 rpm and a comparison with experimental results were performed for validation (Pugliese et al. 2016). The experimental facility allowed the measurements of the discharge, ⁠, and the head drop, ⁠, while the mechanical power, ⁠, was not measured.

The head number has been directly measured during the experimental tests, and generally shows good agreement, except for high flow rate numbers, i.e. for the lowest rotational speed. The experimental values of the power, instead, are not available since the test rig was not equipped with a torque meter. Thus, only the produced electric power was measured, which can be obtained multiplying by the efficiencies of the generator and the speed driver, which are unknown.

A pressure-reducing station in Southern Italy has been considered, where a PRV is installed to regulate the pressure to the desired BP value. The dataset considered in this study contains the measurements of flow rate, ⁠, and upstream head, ⁠, with a time interval of 15 min for 4 days, as shown in Figure 6.

The average discharge is equal to ⁠, with a peak coefficient of ⁠. The upstream head is not significantly variable, with an average value of ⁠, varying from to ⁠.

This pressure-reducing station belongs to a drinking water supply system. Even if the experimental data regards clear water, i.e. with no transported sediments, such a case has been studied for the sake of presenting a methodology that can be applied whenever the PAT operates in variable conditions. As already said, this can happen in multiple cases, such as at the inlet of a wastewater treatment plant or in irrigation systems, where the water collected by open channels can be conveyed to the farms through pressurized pipelines.

A 4-poles asynchronous generator has been chosen, with a constant efficiency ⁠. When the PAT is connected to a 4-poles generator, its rotational speed is fairly constant and equal to about⁠. Thus, its BEP, as calculated by the plot of Figure 5, occurs for a discharge ⁠, producing a head drop ⁠, with an efficiency ⁠. The dataset of the presented pressure-reducing station has been used for the sake of illustrating a VOS application of the analysed machine. Thus, the VOS has been applied for different BP values, in order to find the value that maximizes the efficiency of the plant when it is equipped with the presented PAT. Although this can be not the optimal match in the real case, Figure 7 shows that the PAT is best suitable when the BP is equal to ⁠. In such conditions, the efficiency of the hydropower plant, ⁠, calculated by Equation (4), is equal to ⁠, with a daily produced energy and a maximum produced power equal to ⁠.

Figure 8 shows the PAT curve at ⁠, together with the available points obtained for ⁠. According to the VOS, when the available point lies above the PAT curve, the series valve of the hydraulic circuit produces an additional head loss. On the other hand, a part of the discharge is bypassed when the point lies below the PAT curve.

Figure 9 shows the frequency distribution of both the available flow rate (⁠⁠) and the turbined flow rate (⁠⁠). The figure shows that presents three peaks, for values of ⁠, and 28 l/s. The frequency distribution of is similar, but the range of is clearly smaller, due to the presence of the bypass line which is activated for the highest discharge values. The most frequent discharge is higher than the BEP flow rate. This probably occurs because, when the discharge is higher than the BEP value, the decrease of efficiency has a smaller effect on the produced power than the increase in discharge.

The frequency distribution of flow rates can be considered to assess the amount of solid particles carried by the flow. The particle dynamics and the subsequent erosion they cause strictly depend on the flow field inside the machine. Thus, the erosion of the PAT, as shown in the following sections, depends on the discharge flowing through it.

The E-CODE was applied to assess the erosion of the PAT in different working conditions, in terms of flow rate (from 18 to 28 l/s), particle size (from 0.1 to 0.5 mm) and particle shape. The rotational speed of the impeller was kept constant at 1,500 rpm, and so the density of the particles, which is 2,650 kg/m³, representative of quartz. In this analysis, the wear damage was estimated separately on the inlet pipe, the volute, the impeller and the outlet pipe. The erosion was quantified by different global and local parameters, which have been obtained from the geometrical and particle-related data through the application of the Oka erosion model (Equation (5)), as detailed in Messa & Malavasi (2017). These are: (1) the integral erosion rate on each of the four PAT components, ⁠, i.e. the time rate of mass removal; (2) the integral erosion ratio, ⁠, i.e. the mass of material removal per unit mass of abrasive entering the PAT; and (3) the local penetration rate distribution, ⁠, the velocity at which the local erosion depth increases.

Firstly, it was necessary to assess the numerical convergence of the wear estimates with respect to the resolution of the computational mesh and the number of tracked particles. Since, as already noted, each computed path was the representative of the behaviour of a group of particles and corresponded to a certain solid mass flux, it was necessary to ensure that a statistically significant set of particle impingements was obtained. This was even more so because the particle's paths had a random character through the discrete random walk model.

Statistical-based techniques have been proposed to find out the number of Lagrangian particles which need to be tracked to obtained reliable numerical results. One of these methods was recently used by other authors (López et al. 2018) to assess the consistency of wear estimates in abrasive jet impingement tests. In the present study, the complexity of the fluid dynamic problem forced us to make several approximations, tolerating a level of uncertainty that, although higher than the one that can be obtained when simulating simpler benchmarks using more refined CFD models, does not limit the engineering usage of the results. In this perspective, it was decided to follow a simpler approach often utilized in the literature to verify the adequacy of the number of tracked particles, for instance by Chen et al. (2004); such an approach consists of increasing the number of tracked particles until when the variation of an erosion-related parameter falls within an acceptable range. Such an erosion-related parameter was taken as the integral erosion ratio of the impeller, carrying out the analysis for the BEP flow rate at (i.e. ⁠) and ⁠, considering the integral erosion ratio of the impeller as a target parameter. At first, the effect of increasing the number of tracked particles for a given spatial grid was assessed, referring to the mesh referred to as ‘1.2 M’ in Table 1. The data in Figure 10(a) show that the variability of the integral erosion ratio of the impeller is bounded by only about ±2% above 10,000 tracked parcels. Given that and also considering that switching from 10,000 to 25,000 increases by a factor of 3 the calculation time, 10,000 particles were judged enough for providing statistically significant wear predictions. After deciding the amount of released particles, the effect of the number of cells was investigated, always making reference to the integral erosion ratio of the impeller as a target parameter. The results, shown in Figure 10(b), clearly reveal that the converging behaviour of the data points. The method of Celik et al. (2008) was applied to estimate the grid discretization error, yielding about 17.9% and 7.1% for the 2.5 M and the 5 M meshes, respectively. Indeed, a 17.9% grid discretization error is not small, but it is definitely lower than the erosion prediction uncertainties arising from the uncertain sub-models and parameters (Messa & Malavasi 2017). Furthermore, it is not expected to weaken the conclusions of the comparative analysis reported in the following sections. Based on these considerations, and considering that the 5 M resolution requires four times higher computational cost, it was decided to run the remainder of the calculations on the 2.5 M mesh.

The E-CODE was then employed to investigate the vulnerability of the PAT to wear. As it is shown in Figure 9, the machine works with a variable flow rate between from 18 to 28 m/s, and therefore, the erosion calculations were first aimed at understanding how the risk of damage is affected by ⁠.

Figure 11 depicts the integral erosion rate and the integral erosion ratio as a function of the flow rate. Two considerations can be drawn from the plots. Firstly, the erosion parameters span over two order of magnitudes for the four PAT components, the lowest and highest values occurring for the inlet and outlet pipes, respectively. Nevertheless, higher and values do not necessarily mean higher erosion risk since the two indexes are associated with the speed at which the overall mass removal occurs, whereas the failure due to erosion is due to localized thinning effects, finally resulting in material weakening and hole formation. Thus, two other parameters must be considered when assessing the vulnerability of a PAT component to wear, namely the extension of the surface area affected by the removal of material and the thickness of the element. Although higher in the absolute term, the mass removal from the big, thick volute is likely to create less damage than the mass removal from the small, thin impeller, thereby suggesting that the latter is the most critical element. Secondly, Figure 11 also indicates a different effect of the flow rate on the and different flow components. The data monotonically increase for the inlet pipe and the volute, monotonically decrease for the outlet pipe, and show a minimum with a slight decrease for the impeller. These results can be justified by considering that erosion is the result of three variables related with the particle dynamics, namely, the particle impact velocities, the particle impact angles, and the number and location of the impingements, which should be analysed in a statistical framework. Figure 12 shows, for the different flow rates, the mean values of and (with vertical bars representing their standard deviations), and the number of impingements over the impeller. In order to make the results more general, the impact velocity in Figure 12 has been normalized by the inlet bulk-mean velocity at the BEP, equal to 4.38 m/s. Two curves are depicted in Figure 12(c) and labelled on two different vertical axes: the former represents the fraction of escaped (⁠⁠) over injected (⁠⁠) particles, which slightly increase from 0.85 to 0.95 as the flow rate increases. In principle, in steady-state particle tracking, each particle should leave the domain, but in practice, about 10% of the particles were not able to reach the outlet boundary even after increasing the maximum number of integration time steps. Some of these particles were able to enter the inner impeller zone, but then, they continued hitting the same location of the impeller wall until when the integration of their trajectories was terminated, producing unphysical wear hotspots. Other particles went on circulating inside the volute without being able to enter the impeller zone, sporadically hitting the volute wall at different locations uniformly distributed along its circumference. Incomplete trajectories might indicate flaws either in the numerical scheme (e.g. near-wall turbulence modelling and boundary layer meshing), or in the particle tracking using the fluid flow field derived from a steady-state, RANS-based MRF simulation. In order to avoid unphysical erosion maps and provide a more objective and numerically defined analysis criterion, it was here decided to disregard all incomplete trajectories in erosion calculations. This increased the uncertainty in the wear predictions, but without affecting the overall behaviour of the results and, therefore, the significance of this study. The second curve in Figure 12(c) is the ratio between the total number of impingements and the escaped particles, /⁠, which strongly decreases with increasing the flow rate. The data in Figure 12(c) suggest that the considerable rise in and observed in Figure 11(a) at low can be explained by the higher number of particle-wall impingements: when the machine operates far from its BEP, the velocity vectors at the inlet and at the outlet of the impeller are likely not parallel to the blades; thus, the trajectories followed by the solid particles, having higher inertia than the fluid, impact against the solid wall. Conversely, Figure 12(a) shows that, for ⁠, the standard deviation of the particle impact velocities increases slightly although the mean value remains nearly the same. This suggests that, in spite of a similar mean impact velocity, a higher number of high-velocity impacts occur, probably explaining the mild increase in and observed at high ⁠. The minimum erosion is observed for 24 l/s, where the two effects are balanced.

The influence of the particle characteristics on the vulnerability of the PAT to the impact erosion was also studied. Reference was made to the flow rate at the BEP, i.e. ⁠. Two particles-related characteristics were considered, namely their size and shape.

With regard to the first aspect, the calculations were performed varying the particle size from to ⁠. As a preliminary consideration, it is worth noting that, in this study, the particle tracking calculations were performed under the point-particle approximation (i.e. by treating the particles as point sources of momentum), which requires the particle size to be sufficiently smaller than the cell size. In the impeller and volute regions of the 2.5 M mesh, where the grid density is the highest, the average volume of the computational cells is around ⁠, corresponding to about times the volume of a particle and about times the volume of a particle. Although this strengthens the feasibility of the point-particle approximation, the previous work from one of the authors suggests that evaluating the particles' impact characteristics at the wall could result in an underestimation of the wear damage, especially for the coarsest particles (Messa & Wang 2018). Unfortunately, the computational tools at our disposal do not allow accounting for the finite particle size effect in the simulation of the PAT. This modelling shortcoming is not expected to affect the substantial conclusions reached, but of course, it should be considered when analysing the numerical results. The erosion predictions on the impeller, depicted in Figure 17(a), indicate a drop in the integral erosion ratio followed by a slight increase. In order to provide a justification of these findings, the statistics of the particle-wall impingements were inspected. As shown in Figure 17(b) and 17(c), the particle impact velocities and the particle impact angles tend to decrease as the particle size increases. However, these two variables alone cannot fully explain the trend observed in Figure 17(a), unless the number of impingements is also considered. As it is depicted on the right axis of Figure 17(d), the number of impacts per escaped particle varies non-monotonically with the particle size, similarly to the integral erosion ratio. It is also interesting to note that, as the particle size increases, the fraction of escaped to injected particles decreases, reaching around for the particles. In particular, the increase in the number of incomplete trajectories occurred mostly because a greater number of particles went on circulating in the volute without being able to enter into the inner impeller zone. Incomplete trajectories have already been attributed to flaws in the numerical scheme and in the use of the fluid flow field obtained through a steady-state, MRF, RANS-based simulation to perform particle tracking. Nevertheless, there might also be some physical reasons which, combined with the numerical issues, could explain the ‘big particle effect’. The random fluid velocity component induced by particle turbulent dispersion plays a fundamental role in moving the sediments through the interface between the outer volute zone and the inner impeller zone. Big particles are less affected by turbulent dispersion, and their higher inertia makes it more difficult to curve the trajectories towards the impeller zone; therefore, it is not surprising that the number of escaped to injected particles decreases as the particle size increases. Since, as already said, incomplete trajectories were not considered in erosion calculations, this undesired effect will increase the size of the uncertainty bars in Figure 17, especially for the highest values of particle diameter; nevertheless, it is not likely to affect the overall behaviour of the results.

Finally, the shape of the particle is quantified by means of the particle spherical coefficient, ⁠, which is the ratio between the surface area of the volume-equivalent sphere and the actual surface area of the particle. In the CFD model, the effect of particle shape was accounted for via the drag coefficient correlation: up to this point in the paper, the Schiller & Naumann (1935) correlation for spherical particles was applied. The Haider & Levenspiel (1989) correlation, which can be regarded as a generalization of the Schiller and Naumann one to non-spherical particles via the dependence upon the particle spherical coefficient is now employed. In particular, was set to 0.66, 0.76 and 0.86, as characteristic of natural particles with fully sharp, nearly rounded and fully rounded grains, respectively (Messa & Malavasi 2017). Figure 18 is the analogue of Figure 17 for different particle spherical coefficients. Angular particles with low produce higher vulnerability of the impeller to erosion. This is mainly due to the higher number of impingements. Considering that a decrease in either particle size or particle spherical coefficient produces a higher effect of the drag, the results in Figure 18 appear consistent with the findings in Figure 17.

As a final remark, it is worth noting that the effect of varying the solid concentration has not been investigated. This is due to the fact, if the one-way coupling regime assumption remains acceptable, a change in the solid concentration does not produce any variation neither in the integral erosion ratio nor the particle-wall impingement statistics. Based on Messa & Malavasi (2017), this condition is expected to be valid up to solid concentration of about 1% by volume.

In this work, the vulnerability of a pump used as turbine to the mechanical wear produced by the presence of solids within the flow was investigated using CFD. The main goal was to explore how the risk of wear failure of the different PAT components was affected by the working conditions of the device, mainly the flow rate and the geometrical properties of the solid particles. In order to better demonstrate the engineering potential of the CFD approach, the chosen case study concerns the installation of a PAT in an existing control station of a WSS in Southern Italy. The VOS was first applied to estimate the frequency distribution of flow rate across the PAT. Then, a commercial CFD code – calibrated against experimental data in the absence of particles – was employed to predict the particle-laden flow inside the machine and the particle-wall impact conditions. Finally, an in-house code converted the particle-related data into erosion estimates. The possibility to obtain reliable erosion prediction of a rotating device through the steady-state MRF instead of the more computationally expensive unsteady SM technique was assessed based on a simple criterion, based on the expected axisymmetric behaviour of the wear damage to the impeller (Figure 3). After demonstrating the numerical convergence of the calculations and estimating the grid discretization error, the following goals were reached. Firstly, the points at higher risk of wear failure were identified by introducing an equivalent penetration rate to provide a time-averaged estimate for variable flow rates; the impeller resulted in the most susceptible component, even if a higher mass removal was observed for the volute. Secondly, the actual influence of the flow rate on the vulnerability to wear of all PAT components was established and interpreted in the light of the primarily variables related with the particle dynamics, namely particle impact velocities, particle impact angles and number of particle-wall impingements. Finally, an analogous sensitivity analysis was made with respect to the geometrical properties of the particles, namely size and shape, at the BEP flow rate.

This research was partly funded by the ERDF Interreg Atlantic Area Programme 2014–2020, through the REDAWN project EAPA 198/2016. The work of the third author (F.P.) was financially supported by the fund ‘PON Ricerca e Innovazione 2014-2020, Asse I, Investimenti in Capitale Umano, Avviso AIM-Attrazione e Mobilità Internazionale, Linea 1’ (CUPE61G18000530007).