Modeling erosion in a centrifugal pump in an Eulerian-Lagrangian frame using OpenFOAM®

Abstract Erosion induced by solid particle impingement is a very commonwear mechanism in turbomachinery and Computational Fluid Dynamics is one of the most widely used tools for its prediction. In this article, erosion is modeled in one of the channels of a centrifugal pump using OpenFOAM®,which is an Open Source CFD package. A review of some of the most commonly used erosion models is carried out in an Eulerian-Lagrangian frame along with a comparative study of the erosion rates obtained with each model. Results yielded some disparities between models due to the different factors taken into consideration. The mesh is then deformed to obtain the resulting eroded geometry.


Introduction
Computational Fluid Dynamics is one of the most widely used tools for erosion prediction. The use of this method has been enhanced by the major advances in computer technology in recent decades. These allow the attainment of more accurate solutions as well as decreasing computation times signi cantly. In this work, OpenFOAM's erosion modeling capabilities are assessed. OpenFOAM is an Open Source CFD package that allows modeling ow conditions in turbomachinery and calculating erosion rates induced by particle impingement in any geometrical con guration. In this article, a single channel of a centrifugal pump is subjected to solid particle impingement and the location of erosion is analyzed and results yielded by four di erent methods are compared.
The rst part of the text deals both with the selection of the correct approach for the treatment of the disperse phase and the coupling of the lagrangian library and the incompressible solver (SRFPimpleFoam). A basic description of the implementation of this solver is also outlined. After that, the implemented erosion models are brie y commented and the results obtained for each of them are shown and compared. Finally, a deformation algorithm is applied to one of the models in order to get the ampli ed eroded pro le.

CFD model . Particulate phase treatment
Erosion induced by solid particles impinging on complex surfaces is a very common process in many industrial processes. Two di erent approaches may be taken into account when obtaining numerical solutions to the equations governing uid ows. The Eulerian-Eulerian method is suitable for high particle concentrations. In this case the particles are treated as a continuous phase and the eulerian continuum equations are solved for the two different phases. When the concentration of the particles inside the uid is su ciently low, the continuous phase would be solved in the rst place, while the successive positions and velocities of the particles would be integrated from Newton's equations for motion. This would be the so called Eulerian-Lagrangian approach. The in uence of the particulate phase on the continuous phase can be neglected when the concentration of the disperse phase is adequately low. Calculation of the Particle Mass Loading [1] is a convenient tool for making the correct choice of methodology.
In the case being presented in this article, no signicant two-way particle-uid coupling is to be expected due to the calculated particle mass loading being slightly over 0.015. This indicates that the momentum transfer from the particles to the continuous phase is negligible, as it can be veri ed in the simulations by inspecting the recorded values of this momentum transfer.   .

Domain discretisation
For the present article, a random example of a centrifugal pump impeller has been chosen, and its corresponding volume is presented in Figure 1.
In order to keep the computational cost to a minimum while maintaining a reasonable accuracy, only one channel of the impeller has been taken into consideration, along with the proportional part of the inlet boundary. Figure 2 shows the geometrical con guration of the channel, which has been divided into several parts for an optimized discretization procedure and the location of the different boundary conditions. In the case of the inlet boundary condition, a uniform velocity was chosen, no slip was chosen for the walls of the impeller and pressure outlet for the outlet boundary. As it can be observed in Figure 2, the throatbush was sliced in order to include its proportional part in the simulation. It is for that reason that cyclic boundary conditions were chosen at the sides of the inlet.
In this case, the mesh and procedure for the calculation of erosion have been implemented so that the time needed to obtain the contours of erosion and the resulting deformed geometry is minimized. The mesh chosen is composed almost entirely by hexahedral cells. Given that the turbulence model to be employed is the k-omega, the y + was also calculated. The values obtained indicate the suitability of the mesh for this case, ranging from 1 to 60, which was considered adequate for the model used and for the calculation of the particle variables. Quality parameters for the mesh were ensured by the fact that the mesh is composed almost in its entirety by hexahedral cells.

Steady state solution and Implementation of a Simple Reference Frame solver for an incompressible phase and Lagrangian particles
.

Steady state results
The rst step in the present methodology is the attainment of the steady state for the uid phase. This was achieved by selecting the appropriate solver from the available ones in OpenFoam: SRFSimpleFoam. This solver is suitable for rotating domains on a Single Reference Frame (SRF). The discretization schemes used were second order for the velocity and rst order for the k-omega turbulence modeling. Final solution for the steady state was achieved when residuals fell below 1*10 − . The resulting velocity eld is represented in Figure 3, which shows a slice of the channel at it's central part. The rotational velocity of the pump was 2850 revolutions per minute and the ow rate at the inlet was 0.0063 m /h. OpenFOAM's customisability is one of its main incentives, together with the ease of manipulation of the di erent libraries and the linkage possibilities between them. For the case being treated here, the lagrangian-intermediate library is chosen for coupling with the incompressible solver. This library contains a series of templated functions, models and forces that allow calculation of the particle variables once the carrier phase has been previously solved. The intermediate library provides the user with ve di erent classes for representing the particles. As no particle-particle collisions, chemical reactions or temperature e ects are taken into consideration; the Kinematic class is the most suitable library to be coupled to the incompressible solver. The Kinematic library allows two different particle treatments. The rst one is to track the particles individually, while the second one avoids extra computational costs by tracking a set of particles, or is also called computational parcel. In a computational parcel, in order to capture the behaviour of the real particles, some real case properties are de ned. For the case being taken into consideration, the individual particle tracking approach is selected, due to the previously mentioned low concentration of particles in the uid phase.
The selection of the incompressible solver was SRF-PimpleFoam, which is a transient solver for incompressible ow in a single rotating frame.
The procedure starts with the solution of the uid ow. Once the continuous phase solution is calculated, the di erent forces on the particles are computed. Newton's equations for motion are then solved and the position of the particles is obtained by means of integration. The di erence between the Simple Reference Solver and the pimpleFoam incompressible solver lies in the treatment of the velocity. The uid's velocity inside the rotating volume is calculated as relative, and the particles are a ected by an SRF Force (due to the rotating component of the movement) which has to be incorporated into the disperse phase properties. An additional momentum source to take into account the drag force exerted by the uid is included.
Linking the lagrangian library to the incompressible solver means adding the constructors and les needed to represent the particles inside the uid. Then inclusion inside the solver of the required functions that will take the continuous phase solution's results in order to apply them to the particles is performed. These functions are responsible for the calculation of the forces on the particles before the mentioned integration of the acceleration in order to obtain their positions and velocities.

Erosion calculation
The Kinematic library is largely based on templates, which is a powerful feature of C++ programming language. This language requires declaring variables, functions and most kinds of entities using types. Templates are functions or classes that are de ned for one or more types that are not yet speci ed. Due to the fact that a lot of parts of the code look the same for di erent types, templates become a very useful tool and that is the reason why they are very common within OpenFOAM's code. The function that calculates erosion is in fact, a template, which in this case is dened for the Kinematic class. Impingement information, such as impact speed and impact angle, is gathered as the individual particles or computational parcels hit the speci ed walls of the geometry on which erosion needs to be calculated. After this, the recorded information is used to calculate erosion at each particular face of the boundary and the eld formed by all the faces is stored to a le for post processing. The erosion eld is calculated at each face of the boundaries by means of the addition of each particle's calculated induced erosion. At the moment each particle reaches the speci ed boundaries, the face that was hit is recorded and erosion is calculated and added together with the previous contributions and stored in the erosion eld. .

Implemented erosion models
Erosion induced by solid particles carried by a uid is a common deterioration mechanism that takes place in many di erent environments and kinds of equipment. It is because of this that the literature that can be found on the subject is vast. However, the level of agreement between the di erent authors is very low. Meng and Ludema [2] carried out a literature review of more than 5000 articles published in the most important erosion related journals and conference proceedings. Within these, they encountered almost 2000 erosion models, from which they separated 28 di erent ones. Discrepancies between di erent authors may be explained by pointing out the extreme complexity of erosion processes and the number of variables involved in its calculation. This complexity was illustrated by J.A.C.
Brought to you by | The University of Strathclyde Authenticated Download Date | 8/31/15 11:42 AM Humphrey in [3], where he listed all the factors a ecting erosion that had been previously investigated and grouped them in three di erent categories, namely, the ones a ecting the particles, the surfaces involved in the process and the uid ow. In fact, one of the only theories on which there seems to be some kind of agreement is the one differentiating between two di erent mechanisms a ecting the process: Cutting and Deformation Wear. De nition of these mechanisms has been adopted by numerous authors ( [4][5][6][7][8][9]).
This theory states that, on one side, particles may hit the surface and tear material away with them in a cutting manner. This is called cutting wear and it is the predominant mechanism for ductile materials and particles impinging at low angles. On the other side, several particles might impact on the same place transferring some of their kinetic energy to the surface in the form of hardening work. Eventually, a piece of material detaches from the surface being subjected to erosion and is carried away by the uid. This mechanism is called deformation wear and it prevails for high angles of impingement and brittle materials.
In this article, in addition to the built-in erosion model [10], two additional equations have been implemented in OpenFOAM. Given that each equation is the result of a personal approach (in fact, literally thousands of di erent models can be encountered in the literature) only the location of erosion will be compared, paying special attention to the location of the maximum erosion. Calculation of wear rates for a particular case would imply experimental tests to obtain some of the constants used in the models. However, if the location of maximum erosion and the relative di erence between the di erent eroded parts is successfully represented by the model, the wear experienced by the geometry can be obtained by scaling that obtained values to t the available results.
It is worth noting that the four di erent models calculate the wear rates for exactly the same number of impact locations, angles and velocities. This is due to the fact that the information concerning each impact is simultaneously given as input to the four erosion models. Due to the very di erent nature of the equations used for erosion calculations, the most signi cant data we can take from the calculations is the location of the maximum erosion along with the relative di erence in erosion between the di erent areas of the channel. It is for this reason that the labeling only speci es maximum and minimum erosion locations. Once we have the real geometry to compare our results with, if the locations and the relative di erence match up, it is only a matter of scaling the numerical values to t the real ones.

. . Finnie Erosion model
The present model is one of the rst and still widely used models in the literature [10]. It is divided into two di erent equations, which depend on the angle of impingement. These equations yield the volume of material removed by a single abrasive particle of mass m, travelling at a velocity denoted by V. Two experimental constants are present in the set of equations, namely, the ratio of depth of contact to the depth of cut and the ratio of vertical to horizontal force components. The material's plastic ow stress of the material subjected to erosion and the angle of impingement are the last two variables to complete the formulae.
It is worth underlining that this model yields accurate results for low angles of impingement (which is where maximum erosion usually takes place in ductile metals) but it predicts no erosion at normal incidence. In a second publication [11], Finnie modi ed his form in order to account for some other variables that would allow capturing some erosion for normal incidence. Calculated results for Finnie's model are shown in Figure 4 and the model is shown in Equations (1) and (2).
The equation yields the volume of material, Q removed by a single abrasive grain of mass, m, and velocity V: where p is the plastic ow stress of the material, α is the angle of impingment, ψ is a constant that represents the ratio of depth of contact to the depth of cut and K is the ratio of vertical to horizontal force components on the particle.

. . Grant and Tabako Erosion model
The main assumption made for the development of this equation is that erosion depends on two di erent mechanisms that take place at low and high angles of impingement respectively. A combination of these two mechanisms is used when calculating erosion for intermediate angles of impingement. Obtained erosion with Grant and Tabako 's model [9] can be seen in Figure 5    shown in Equation (3).
where ε = the erosion damage per unit mass of impacting particles K = material contact f (β ) = empirical function of particle impact angle V T = tangential component of incoming particle velocity V T = tangential component of rebounding particle velocity f (V N ) = component of erosion due to the normal component of velocity

. . Nandakumar erosion model
The last model was developed by Nandakumar et al. [12] who took as a starting point other authors work (Finnie [10] and Bitter [5,6]). This model calculated the volume eroded by the impacts of the particles on the material's surface. Again empirical constants are present in this model so that experimental data is needed in order to calculate the eroded volume for a particular case. Figure 6 shows the location and relative di erence of the erosion predicted by Nandakumar's model and it's mathematical expression is detailed in Equation (4).
where C and D are empirical constants, m is the mass of the particle, ρp is the density of the particle, Θ is the angle of impact, V is the impact velocity and dp is the diameter of the particle.

Results discussion
It is clear that the three di erent models predict the location of the maximum erosion at the same location. However, discrepancies arise when quantifying that erosion and the relative erosion around that maximum erosion location. While Finnie's and Grant and Tabako 's's erosion models both yield similar results for the considered pump channel, Nandakumar's model predicts wear in locations where the other two don't. Given that the three models yield similar results for the location of the maximum erosion, it seems reasonable to choose Hashish's model for a rst approach, due to the lack of experimental constants in the formulae being thus, not necessary to perform experimental tests to obtain them.
However for a time dependent evolution of the wear scar, a much more complex approach would have to be considered, applying surface deformation so that the evolution of the velocities and angles of impingement with the development of the eroded surface are properly accounted for.    Figure 7 shows the initial shape of the pump channel while Figure 8 shows the same channel after the deformation algorithm that the authors developed in OpenFoam is applied. Finally, Figure 8 provides con rmation that the deformation of the surface follows the erosion contours, i.e., the higher the erosion rate, the larger the deformation. Deformation of the surface is the rst step towards a more accurate calculation of how erosion develops, which is one of the next steps that the authors intend to follow.

Conclusions and future work
Provided that an enormous number of di erent erosion models can be found in the literature, it seems obvious that one could search for a particular method suitable for the application taken into consideration. However, this search might be time consuming. If only a reasonable initial approach is needed, Hashish's model appears to be the best choice, while for more accurate predictions a literature review should be performed searching for a more suitable method for the each particular application.
The work to be continued in the future includes the study of the ow changes with surface deformation due to erosion, thanks to the development of the deformation algorithm commented here along with the validation of the CFD model for low concentrations by means of Particle Image Velocimetry as well as the implementation of a methodology for erosion prediction that accounts for the surface changes with time.