Erosion under turbulent flow: a CFD-based simulation of near-wall turbulent impacts with experimental validation

ABSTRACT Engineering systems operating under particle-laden turbulent flow regimes, such as slurry pumps, are prone to erosive wear, where the specific conditions of particle impingement are crucial for determining the wear rate of the eroding surface and thus its relevance for the mechanical or hydraulic system. However, the near-wall interactions under turbulent flows impose a challenge for the formulation and validation of experimental and numerical models that describe the particle’s trajectories close to the surface. This work analyzes the statistics of particle impacts obtained from a Large-Eddy Simulation (LES) with the Wall-Adapting Local Eddy-Viscosity (WALE) model of the flow and its consistency with experimental data on the statistical distributions of particle impact angle and impact direction determined in a slurry pot configuration. The resulting distributions of impact angle, impact direction, and particle velocity derived from the implemented CFD formulation produced characteristic values similar to those derived from experimental data. The LES-WALE approach used in this work offers shorter computational times and is shown valid by comparative analysis with direct simulation. The model’s applicability and prospective dynamic simulation of practical engineering cases are discussed.


Introduction
Engineering systems associated with fluid-solid multiphase flows, or slurry flows, are present in many industrial systems featuring hydraulic mechanical machinery, such as mining pipelines, turbines, and pumps (Calderón-Hernández et al., 2020;Liu et al., 2019). In these systems, wear degradation of the affected materials is produced by the synergistic effect of the flow with corrosion, commonly referred to as erosioncorrosion (Aguirre & Walczak, 2017;Burson-Thomas & Wood, 2017;Gnanavelu et al., 2009). Despite the certainty regarding the economic risks associated with slurry flow, the current understanding of the involved degradation mechanisms is insufficient to provide an accurate prediction of wear rate by a model that would be applicable (Kuruvila et al., 2018). The comprehension of erosion-corrosion involves a full description of the interactions between the electrochemical and mechanical phenomena on the surface that are strongly related to turbulent flow. Whereas the electrochemical aspects are relatively well understood, a comprehensive description of the particles' impact conditions is missing. Considering that the global damage originates from the sum of CONTACT Magdalena Walczak mwalczak@ing.puc.cl individual impacts, the understanding of which particle impacts contribute efficiently to the damage is fundamental, especially in a turbulent flow. In this research, we focus on pure erosion, i.e. independent of corrosion, by not considering the electrochemical effects, to understand the near-wall turbulent structure's behavior. For this, we develop a validation methodology applicable to the study of erosion-corrosion or other engineering problems related with fluid structure interaction in the future. The application of computational fluid dynamics (CFD) has enabled the study of the effect of turbulence on corrosion and erosion in hydraulic systems (Bremhorst & Brennan, 2011;Gao et al., 2020;Manzar & Shah, 2009;Parsi et al., 2017;Yan et al., 2021). However, the formulation and validation of a numerical model for near-wall flow that would accurately capture the physics of the process is a remaining challenge, limiting the generalization of wear prediction for arbitrary erosion geometry. The velocity of the impacting particles may not be correctly calculated, leading to erroneous erosion patterns (Zhang et al., 2018), which has been attributed to: (i) limited applicability of the turbulence models for predicting non-physical impacts , (ii) dependence on the choice of the cell size in the nearwall region (Xie et al., 2021), (iii) unaccounted for effects of particle size (Wang et al., 2021b), or (iv) unpredicted alterations of steadiness in the near-wall zone (Wang et al., 2021a). Since turbulence modeling is an essential feature for determining the particle impact characteristic over the eroded wall (Zhang et al., 2018), a method of validating the numerical model applied to fluid surface interaction for wear prediction is necessary.
Wear by erosion canbe considered a consequence of the interaction between an erosive fluid and a solid surface, manifested by appearance of surface scars. The geometry of these scars is related to near-wall fluid dynamics, which in turn results from the global flow behavior. The underlying idea of this work is the use of surface scars produced by particles conveyed by the fluid to analyze the magnitude of impact velocity and impact angles onto the surface.
From the material point of view, formation of the scars is the very erosive degradation of metals, the primary mechanism of which is the plastic deformation caused by particles impacting at angles nearly normal to the exposed surface, accompanied by plowing/cutting produced by particles projected from the turbulent flow at lower angles (Clark & Wong, 1995;Lospa et al., 2020). A recent literature review of the erosion mechanisms affirms that surface erosion is influenced by a wide range of factors that can be categorized as material properties, mode of particle impact, and general operating conditions (Tarodiya & Levy, 2021). Because of the multi-factorial character, the current erosion models are derived from, or related to, erosion wear energy as it relates to the cutting and deformation mechanisms. It is also related to the mechanisms of energy dissipation associated with the normal and tangential components of particle impact velocity (Javaheri et al., 2018). Then, in order to predict the erosion rate, it is crucial to correctly describe the impact conditions of a collective flow, i.e. comprising the normal and tangential components of impact velocities of every particle (Zhao et al., 2020).
In erosion wear studies, a variety of experimental tools is available for assessing the erosion behavior. The slurry pot configuration, for example, is often used due to its low cost and expediency of operation. Also, it has been shown effective for quantitative evaluation of the effect of global parameters of the slurry and its flow characteristics, such as mean velocity, density, concentration, and erodent hardness, on the rate of erosion (Tsai et al., 1981;Valtonen et al., 2019). Further, slurry pot testing, combined with characterization of individual scars produced on a ductile metal under highly turbulent diluted slurry, was shown effective for determining the normal and tangential components of dissipated kinetic energy of the impacting particles (Molina et al., 2020). Consequently, information on the global flow can be conveyed through near-wall turbulent fluctuation into individual indents. For further elaboration, it is therefore necessary to verify the consistency between the experimental observation and predictions available through simulation of turbulent flows.
In order to derive erosive wear from a model of fluidsolid turbulent flow the impingement parameters of particle impact velocity, impact angle, and frequency of impacts need to be quantified. These parameters, in turn, are determined by the flow regime and its conditions (Clark, 2002). In this study, we explore a simple configuration of concentric rotating cylinders (see Figure 1), i.e. a slurry pot, in which the flow conditions can be characterized as a Taylor-Couette (T-C) flow (Taylor, 1923). It has been shown that a T-C flow can be simulated by specific numerical models capable of representing the formation of flow structures and near-wall effects (Grossmann et al., 2016;Poncet et al., 2013) relevant for erosion. Furthermore, under T-C flow featuring a turbulent structure, the particle trajectories can be obtained by a two-phase fluid-dynamic model based on the Eulerian-Lagrangian approach, where the fluid flow is solved in an Eulerian framework and the solid phase is simulated by tracking the trajectories in a Lagrangian framework (Chen et al., 2020). These simulations usually involve complex algorithms, entail intensive computation, and require convoluted validation strategies for the case of turbulent slurry flow (Elghobashi, 1994;Messa & Malavasi, 2018). In order to circumvent these disadvantages, this work exploits a single-phase CFD-based approach to derive the impingement parameters based on the premise that lowinertia particles are conveyed mainly along the streamlines of vortices (Clark, 1992) associated with the particle size.
In paricular, this work proposes a Large-Eddy Simulation (LES) with Wall-Adapting Local Eddy-viscosity (WALE) at subgrid-scale to represent the effects of velocity fluctuations at the relevant turbulent length scales of flow structures near the rotating wall (Christensen et al., 2010;Nicoud & Ducros, 1999). This method also allows for the resolution of critical vortex sizes relevant to the transport of physical particles, as well as the inclusion of smaller vortex sizes relevant to particle rotation applicable to Eulerian-Lagrangian models (Breuer & Almohammed, 2015). The objective of the single-phase CFD-based approach is (i) to determine the statistics of the impact conditions of streamlines capable of producing erosive scars on the target surface located at the inner rotating cylinder set-up; and (ii) to verify its consistency with experimental literature data on angle, direction, and velocity distributions obtained experimentally.

Methods
A single-phase CFD-based simulation is applied for numerical replication of experimental data of particle impact angle, particle impact direction, and statistical distribution of impact velocity determined under turbulent flow in a slurry pot configuration and reported elsewhere (Molina et al., 2020). The slurry pot test geometry consisted of an inner rotating cylinder concentric with an outer non-rotating cylinder, resembling a Taylor-Couette (T-C) flow geometry shown schematically in Figure 1. The model uses Large-Eddy Simulations (LES) with sub-grid-scale of Wall-Adapting Local Eddy viscosity (WALE) (Nicoud & Ducros, 1999) to simulate the turbulent flow on the rig as Taylor-Couette Flow (Grossmann et al., 2016). The statistical analysis of the impact angle, impact direction, and impact velocity of the flow onto the target surface, i.e. the erosive flow; was conducted using the near-wall flow velocities obtained from the fluid dynamics simulations.

Erosion by inertial particles under fully developed turbulent flow
Erosive wear is the degradation of solid surfaces caused by successive impacts of particles. In general, there are five factors identified as most significant (Clark, 2002): (i) particle concentration, (ii) flow speed, (iii) impact angle, (iv) particle size, and (v) particle shape. These factors are necessary for the modeling of the three main energy dissipation mechanism in erosive wear: wall deformation, wall cutting, and particle deformation, in which coefficients of restitution, damping factors, and coefficient of friction are incorporated to determine the energy losses. The kinetics and dynamics of the particles prior to impact are mostly determined by the flow state, thus understanding the flow parameters onto the target surface is critical to modeling erosion (Clark, 1992).
After Clark and Wong (Clark & Wong, 1995), total erosion in a unit of time (E T ) is defined as the sum of deformation (E D ) and shear erosion associated with cutting (E C ) (Equation (1)), both detailed in Equation (2): where M P is the total mass uniformly impacted by particles, V T and V N are the values of the tangential and normal components, respectively, of the velocity of the particles impacting at anangle α. The constants ε and ϕ are empirical and related to energy of erosion deformation and energy of erosive cutting, respectively (Javaheri et al., 2018). The erosive wear is then primarily determined by the vectorial components of velocity, the distribution of which is determined numerically as follows.
The dynamics of the particles are considered to be driven by the flow structure of the liquid phase, where inertial particles are dragged along by the major streamlines and clusters following the patterns of the Very Large-and Large-Scale Motion (VLSM & LSM) patterns (Johnson et al., 2020;Lee et al., 2019;Wang & Richter, 2019). In highly turbulent flows, inertial particles are swept by large-scale motions minimizing their acceleration in the process. As described by Chen et al. (Chen et al., 2006), the tendency of particles to be swept can be evaluated by the relation between the Stokes number and the ratio of scale times response, explained with detail in the following paragraph. In sufficiently diluted particleladen flow, the particles dragged by VLSM and LSM can be assumed to have a negligible effect on the liquid phase dynamics (Chen et al., 2006). Thus, a stand-alone analysis of liquid phase dynamics is sufficient for inferring the dynamics of particle impacts, which justifies the use of a single-phase CFD-based approach to determine the impact angle and velocities of particles impinging onto a solid surface.
The experiment of Molina et al. (2020) was configured as a diluted slurry with spherical silica particles of different size, (density of 2.5 g/cm 2 , Vickers hardness of 6.2 ± 0.5 GPa), tested in the configuration of rotating inner cylinder (Figure 1) under constant tangential velocity of 10 m/s. The flow characteristics of this experiment, relevant for our analysis, derive from particle trajectories and the involved time scales. The particle trajectories depend on the ratio of particle time response (τ p ) and the fluid time scale (τ f ). This ratio is known as Stokes number (St): Here, ρ f , d p and μ f define the fluid density, particle diameter and fluid dynamic viscosity, respectively. Equations (3) and (4) determines the Stokes number regime. For St 1, the particles are expected to follow the fluid trajectory. On the other hand, the flow time scale τ f can be analyzed depending on the multiple time scales for the turbulent flow by means of the integral time scale (τ ) associated with large scale motions, also known as eddy turnover time, or the Kolgomorov time scale (τ λ ) associated with the smaller scales that are defined by viscosity. For τ p τ λ , the particles are expected to follow all turbulence scales, but for τ p τ , the particle motion would be driven mainly by integral scales like VLSM & LSM, unaffected by the Kolmogorov time scale (Loth, 2010). In this work, the integral time scale is defined as follow: where U is the velocity scale associated with the integral length scale, ν is the kinematic velocity of the fluid and τ is the eddy turnover time or integral time scale. Equation (8) defines the eddy turnover time for a T-C flow (Ostilla-Mónico et al., 2014). In this work, the particle time response (τ p ) and integral time scale (τ ) are calculated as 0.0187 s and 0.00275 s, respectively, thus obtaining a St = τ p τ = 6.78, which is the particle Stokes number evaluated using the integral time scale. With this number, it is reasonable to assume that the particles follow the flow trajectories of VLSM and LSM (Loth, 2010).
Although the referenced experiment was conducted for five particle size distributions, in this analysis, we focus on the medium size of 548 μm because the experimental results did not show any statistically significant differences in impact conditions between the evaluated sizes. Also, it allows for the use of a coarser mesh that favors shorter time to solution.

Mathematical description of flow in a slurry pot tester
The geometrical parameters of the slurry pot tester are shown in Figure 1, where, adopting the notation by Grossmann et al. (Grossmann et al., 2016), the inner and the outer cylinder radii are r i and r o , respectively. The annular gap is defined as d = r o − r 1 , and L corresponds to the domain height. The dimensionless parameters for the system can be expressed by the ratio of radii η = r i /r o and the aspect ratio = L/d. The rotation is quantified, in general, by the angular velocities ω i and ω o , corresponding to the inner and outer rotational cylinders of the T-C flow, respectively. The corresponding Reynolds numbers Re i and Re o are defined by Equation (9): where ν is the kinematic viscosity of the fluid domain. By convention, Re i is a positive number, and Re o > 0 and Re o < 0 correspond to co-rotating and counter-rotating outer cylinders, respectively. Alternatively to using the inner and outer Reynolds numbers to characterize the driving of T-C flow, it is possible to use the Taylor number (Ta), represented by the Equation (10): The rotation ratio, a = −ω o /ω i , can be represented as a nondimensional number by the inverse Rossby number: Finally, the Navier-Stokes equation for this system can be formulated in a co-rotating coordinate system as: where f (η) = (1 + η) 3 /(8η 2 ) and ∇p is the pressure gradient. Under this formulation, Ro −1 characterizes the strength of the driving Coriolis force directly. In the spatial case considered here, the outer cylinder does not move, thus Re o = 0, and the Navier-Stokes equation is simplified to Equation (13):

Numerical simulation and calculation
The T-C flow in this configuration is particularly challenging for numerical approaches because of the complexity associated with high rotation rates, the effect of confinement, and singularities in the boundary conditions. Several numerical approaches have been applied to address turbulence in the simulation of T-C flows including direct numerical simulation (DNS) (Ostilla-Mónico et al., 2014), Large-Eddy Simulation (LES) (Poncet et al., 2014), and Reynolds Averaged Navier-Stokes (RANS) (Poncet et al., 2013) simulations. RANS models have shown to produce fair results for the mean and turbulent average fields in a comparison test, but fail to predict the mean tangential velocity profiles. This results in significant discrepancies, even in the core of the flow (Poncet et al., 2013) . On the other hand, the LES approach based on the WALE subgrid-scale modeling has shown to be applicable to study the mean and turbulent flow fields. In particular, thin boundary layers developed along the walls and coherent structures can be represented satisfactorily (Poncet et al., 2014). The WALE strategy has proven capable of producing accurate values. Also, it can save about 12% of computational time compared to other subgrid-scale models, such as the Smagorinsky model (Poncet et al., 2014).

LES-WALE model
The present work employs the WALE subgrid-scale model developed by Nicoud and Ducros (Nicoud & Ducros, 1999). The model is based on the gradient velocity tensor g ij (Equation (14)): to represent the physics of the flow. The velocity fluctuations close to the rotating wall associated with the subgrid-scale characteristic length¯ are thus related to the rotation rate of turbulent structures (Poncet et al., 2014).
In the LES-WALE model, the turbulent eddy-viscosity for eddy scales smaller than the grid is modeled by Equation (15) where $S_ij$ is the deformation tensor of the resolved field and the S d ij operator is defined by Equation (16): The S d ij operator is based on the traceless symmetric part of the square of the gradient velocity tensor g ij . The constant C w considers the same ensemble-average subgrid kinetic energy dissipation as the classical Smagorinsky model, and it can be quantified numerically using a variety of homogeneous isotropic turbulence fields. As suggested by Nicoud & Ducros (1999), C w can be considered a fixed constant specific to a particular condition of flow.
It should be noted that the viscous layer is modeled as a function of both the strain and the rotation rates, which necessarily implies that velocity goes to zero at the wall. For this reason and as explained by (Poncet et al., 2014), no damping function or dynamic method is included to fit the effect of non-slip condition.

Geometry and flow parameters
In the experiment, the fluid is confined between the concentric cylinders of the outer and inner radii r i = 7.5 mm and r o = 35 mm, respectively, and height of L = 105 mm, producing a gap of d = 27.5 mm, corresponding to ratio of radii η = 0.2143 and aspect ratio of = 3.8182. The volumetric domain is modeled as distilled water with a kinematic viscosity of ν = 3.667×10 −6 m 2 s −1 . The flow is modeled as isothermal, as temperature variations on short simulation time do not contribute to convective fluctuation on the domain. In addition, temperature variations have not been reported on experimental data on which this work is based.
The outer rotating cylinder is static, i.e. its angular velocity is ω o = 0 rad/s, and thus Re o = 0. An angular velocity of ω i = 1333.33 rad/s is imposed on the inner rotating cylinder, resulting in a tangential velocity of Vp = 10 m/s. The Reynolds and Taylor numbers evaluated by Equation (9) and Equation (10) are then Re i = 7.5×10 4 and Ta = 9.783×10 16 , respectively.

Mesh
For the mesh size selection, the following aspects were considered: In LES models, the domain meshing is directly related to the eddy-scale accuracy required by the problem. Accordingly, the eddy-scale larger than the mesh grid is calculated directly, and subgrid-scales are modeled (Figure 2(a)). In the erosion problem, it is essential to obtain high accuracy in the velocity of eddies that convey the particles in the fluid domain near the eroded wall. In an Euler-Lagrange approach, the particle acquires the instantaneous velocity of the mesh point (Elghobashi, 1994). Swirls smaller than particle size mainly induce rotation of the particle rather than translating it to the next mesh point (Figure 2(b)). Therefore, in this single-phase CFD-based approach, the structured grid was refined near the inner cylinder to capture the details relevant to the turbulence scale associated with Figure 2. Selection of length-scale for the flow modeling: (a) Energy of turbulence relevant to LES in reference to DNS approach, (b) Relation between particle size and the size of grid determining for the minimum energy of turbulence shown in a.
the erosion mechanism (translation of eroding particle), whereas a coarser grid was used in the bulk flow domain, extended to the outer cylinder.
The assumption of independence of the distributions of impact angle and impact direction from the particle size is based on experimental observations in a diluted slurry (Molina et al., 2020). This premise allows adjusting the minimum near-wall grid size to the minimum eroding particle size without the necessity of supposing a change in wear mechanism (Clark, 2002). To avoid high computational effort and considering the particle size for which significant erosion damage is expected, the minimum mesh size of 500 μm was chosen. This value is the length-scale filtering for the integral and Kolmogorov scales, following the premise of considering only the scales of particle trajectories relevant for erosion and Kolmogorov scales relevant for sub-grid diffusion and global energy dissipation. Details of the mesh are shown in Figure 4.

Boundaries and initial condition
The boundary condition for the inner cylinder was set to non-slip rotating wall with ω i angular velocity. The outer cylinder and bottom wall were set as non-slip, and the open-top as a slip wall. The initial domain condition was set as static flow 0 m/s and zero-gradient pressure condition. The angular velocity of the inner cylinder was applied at time zero with no ramping up.

Simulations
The proposed model was implemented in the OpenFoam framework (Weller et al., 1998), version 8. The Pim-pleFoam solver with a Gauss linear divergence scheme was selected as suitable for the turbulent incompressible flow under transient conditions. The solver is based on the PISO (Pressure Implicit with Splitting of Operator) algorithm, and has been shown suitable for inherently unstable problems such as transient simulations of incompressible and turbulent flow of Newtonian fluids (Jasak, 1996).
The numerical model was run using a time step of 10 −5 s to obtain a controled Courant Number that minimizes numerical diffusion and iteration divergence induced by an unsuitable time step. The total simulated time was 5 s, selected to exceed 40 large eddy turnover time lapses (r o − r i )/(r i ω i ) (approximately 0.11 s) and to assure a non transient state, turbulent regime and developed flow conditions (Ostilla-Mónico et al., 2014).

Model validation
The validity of our model was assessed by analyzing the consistency of our results with those of DNS modeling, which in turn was reported as experimentally validated (Grossmann et al., 2016;Ostilla-Mónico et al., 2014). We replicated the simulations of Grossmann et al., 2016 and Ostilla-Mónico et al., 2014 by applying our model to the same geometry and flow conditions as those of (Grossmann et al., 2016;Ostilla-Mónico et al., 2014) . The Pseudo-Nusselt numbers along the relevant spectrum of turbulence and the azimuthal velocity profiles obtained by both models were used for comparison.

Pseudo-Nusselt number
The validation process consists of an evaluation of the numerical model at different Taylor numbers that cover the range of the system's response under laminar, classical, and ultimate turbulent regimes known from experimental data. For each Ta number, the Pseudo-Nusselt number, defined as (Nu ω − 1), was determined. This number can be interpreted as the relation between the energies dissipated and added to the turbulent system, respectively, which permits validation of the numerical model with the physical energy dissipation. Following the procedure used by Ostilla-Mï£¡ï£¡ï£¡ï£¡ï£¡ï£¡ï£¡ï£¡ï£¡ï£¡ï£¡ï£¡ï£¡ï£¡nico et al. (2014) to validate DNS models, each pseudo-Nusselt number was compared with experimental results for the corresponding Ta number.
The system response is the torque required to drive the cylinders to obtain a determined Taylor number in the flow, and it is related to energy dissipation balance into the flow (Eckhardt et al., 2007). Based on the theory proposed by Eckhardt et al. (2007) that compares T-C with Rayleigh-Bérnard (RB) convection, a pseudo-Nusselt number has been proposed as a non-dimensional response of the system to the torque required to drive the cylinders. It is defined as Nu ω = T/T pa , where T is the torque and T pa is the torque required to drive the cylinders when the flow is purely azimuthal (Ostilla-Mónico et al., 2014). The pseudo-Nusselt number has been shown useful by other authors for comparing the T-C response under different Ta conditions, e.g. (Grossmann et al., 2016;Ostilla-Mónico et al., 2014). The change of the pseudo-Nusselt number over a range of Ta is a relevant characteristic that describes the global system response. In particular, it facilitates the comparison of numerically simulated flow behavior with experimental data for the different turbulent regimes, because it is indicative of the change of turbulent energy and its dissipation on the flow. This interpretation was demonstrated by torque measurements of inner and outer cylinders (Grossmann et al., 2016) justifying the rationale behind Figure 3(a).
The LES-WALE model used in this work was evaluated using the following values of Ta: 7.04×10 5 , 4.77×10 7 , 4.5×10 8 , and 4.63×10 9 . For laminar and classical turbulent flow, a coarse mesh and C w = 0.325 were used. In the high Reynolds number regime, the model was implemented using a refined mesh and C w constant adjusted to 0.5 according to Nicoud & Ducros (1999). Figure 3(a) compares the pseudo-Nusselt numbers (Nu ω − 1) obtained from our model and the DNS simulation reported by Ostilla-Mónico et al. (2014) along with its experimental measurements covering the range from laminar through classical to ultimate turbulent regimes. Although the available data does not cover the same Ta range for all the methods, as only a limited range of the Ta number was verified experimentally by Ostilla-Mónico et al. (2014), the similarity of the results is apparent and a general consistency between our results and the DNS model can be stated.
Further, the three flow regimes are correctly captured: the laminar, classical turbulent, and the ultimate turbulent regime. This consistency at lower computational cost is attributed to meshing depending on the scale of the coherent structures. In particular, in the ultimate turbulent regime a refined mesh was used, compared with the coarser mesh used in laminar and classical turbulent regime.

Velocity profile
The azimuthal velocity profile obtained by the proposed model is compared with DNS results obtained from literature data (Grossmann et al., 2016), for which T-C flow presented coherent structures in the relevant range of Taylor numbers. Figure 3(b) compares the mean azimuthal velocity profile obtained by our model along with the experimental profile obtained by Grossmann et al. (2016) and Ostilla-Mónico et al. (2014). It is possible to visualize that the presented model has the same trend and magnitude in the mean azimuthal velocity profile near to inner wall cylinder, which is the zone of interest for this study. Also, in the radii of interest, the flow velocity obtained from our model is between both DNS curves.

Particle impact parameters
The characterization of erosion impact was conducted considering normal and tangential references with respect to the impacted surface. First, cylindrical unwrapping was conducted to obtain the normal and tangential velocities from the global Cartesian coordinates, followed by decomposition of the velocity field in the cylindrical coordinates.

Cylindrical unwrapping
The cylindrical unwrapping consisted in transforming each point of the specimen surface from its Cartesian coordinates (X, Y, Z) to its tangential coordinates (X', Y', Z'), using the convention shown in Figure 4. Each point of the flow velocity field (u) in the local Cartesian coordinate (x, y, z) was assigned prior totransforming into the global tangential coordinate system, from which the tangential and normal components of each flow impact over the cylindrical surface were determined ( Figure 4). The flow velocity in the unwrapped system (X', Y', Z') was defined as flow tangential velocity (V fluid ), which was then decomposed into the x, y and z components according to Equations (17), (18) and (19), respectively.

Impact velocity decomposition
The impact velocity vector was calculated as a relative velocity between the inner cylinder surface and the simulated velocity field in the near-wall region, where V rel is the relative velocity between fluid velocity (V fluid ) and the tangential velocity of the inner cylinder surface (V wall ). This relative velocity is the approach to the impact velocity of the particles that are carried by the fluid. Figure 5 shows the velocity decomposition to determine the impact angle α, impact direction δ, inclination angle θ , and azimuth angle φ, using Equations (21), (22), (23), and (4), respectively:

Statistical analysis of particle impacts
The statistical analysis was conducted on the data of impact velocity decomposition obtained from the simulation results for the selected surface of interest as

Figure 5.
Coordinate system and definition of characteristic angles. Decomposition of the particle impact velocity vector with respect to a plane tangential to the target surface (grey). The angles α and δ correspond to impact angle and directionality of impact, respectively, whereas θ and φ are the inclination and azimuth angles in the cylindrical coordinates.
follows: Each mesh point around and near the specimen or eroded surface determines the relative velocity between flow streamlines and the surface, necessary to obtain a statistical distribution. The mesh points that were selected corresponded to the near-wall layer located at r i < r ≤ r i + r, where r is the minimum mesh size in radii direction (0.5 mm for this work). The outer direction velocities were filtered to represent the streamline velocity at the time instance just before impact onto the eroded surface. From all the computed data on impact velocity decomposition, data points were selected randomly in a manner such that they be statistically representative for the experimental data available for diluted slurry. In particular, Molina et al. (2020) reported about 300 experimental impact scars on the selected surface area, which were analyzed to obtain the impact velocity from the scars geometry employing erosion models. Then, the calculated and experimental angles and velocity of impact were summarized in a frequency distribution to be compared in terms of the interquartile range, mean, and median parameters. The computed velocity data were filtered to suppress values associated with combinations of impact angles and velocities not detected experimentally, because these impacts velocities and angles cannot produce detectable scars on the target surface. In particular, the experimental data show that the measurable velocity is over 1.93 m/s and that the impact angle is between 1.4°and 19°. This velocity limit and impact angle range was applied as filters to the numerical results in order to enable comparison between distributions. However, filtered points did not make a significant difference in the distributions, as is shown in Figure 8.
The calculated angles and velocity of impacts were visually compared with experimental statistics by plotting the distributions. In addition, statistical test of Kruskal-Wallis rank-sum test (Kruskal & Wallis, 1952), Wilcoxon rank-sum test (Wilcoxon, 1945), Holm method for P-value adjustment (Holm, 1979), and Welch two sample t-test as non-parametric methods (Skovlund & Fenstad, 2001) were used to verify the visual observation.

Results
This section shows the simulation results and the statistics of the particle impacts approach using identical model parameters and grid configuration for the model validation under the same regime to simulate the particular slurry pot geometry.  axial section of the slurry pot simulated at Re i = 7.5 × 10 4 . Both, instantaneous t = 2 s (Figure 6(a)) and time averaged data (Figure 6(b)) are shown. Figure 7(a) summarizes the radial velocity profiles displaying the velocity components individually, analyzed at the height corresponding to the target surface as displayed in Figure 6(b). The vertical profile shows upward and downward flows representing a flow recirculation, which along with the visual appearance of the velocity field is indicative of the presence of coherent flow structures. The mean azimuthal and vertical velocity show a characteristic pseudo-Taylor roll corresponding to breaking coherent structures visible in the instantaneous velocity field. Further, the mean azimuthal velocity profile shown in Figure 7(a) is similar to that observed in Figure 3.

Statistics of impact velocity
The impact velocity expected for a low-inertia noninteracting particle corresponds to the relative velocity between near-wall flow, exemplified in Figure 7(b), and the rotating target surface.
The velocity magnitudes obtained for all points are represented as raincloud graphs with the corresponding distribution envelope in Figures 1 and 8. The distribution (at the left) is compared to the experimental data of Molina et al. (2020), which were obtained from analysis of surface deformation per impact by applying the model of wear deformation by either Cheng et al. (2019) or Huang et al. (2008). The results obtained by flow simulation are generally consistent with those obtained by the Cheng model; whereas the median and third quartile of the numerical results are in the range of first quartile and median of the results obtained by the Huang model. The density of the numerical data points is higher than that of the experimental data points because here we consider the entire surface of the exposed sample, that is 500 mm 2 , whereas data of a representative area of 6.12 mm 2 were reported by Molina et al. (2020).
It is important to note that the numerical results do not produce outliers of values higher than 10 m/s corresponding to the linear velocity of the rotating inner cylinder, whereas such outliers were observed in the experimental data. On the other hand, a cut-off at lower velocity is also a feature of the experimental data because not all particles impacting onto the target surface convey sufficient energy, or impact at an angle favorable to produce a scar, so some of the numerical velocity data were filtered out. Both filtered and unfiltered velocities are shown in Figure 8 revealing that the envelope of the filtered velocity distribution is of a shape more similar to the experimental data: it is more centred about the mean  . Impact velocity distributions obtained by the numerical modeling (green) as compared with data derived from experimental analysis of wear scars using erosion models by either Cheng (red) or Huang (orange). Data adapted from Molina et al. (2020). The numerical results are displayed as either filtered or unfiltered (blue). Blue dots consider velocities associated with impact angle for which erosion is expected. Red dots correspond to velocities for which no erosion is feasible and that were filtered out for further analysis. value but with narrower data scattering. The quantitative comparison of the distributions is shown in Table  1, where the similitude in interquartile range between the numerical model (unfiltered) and values derived from the Cheng model can be observed. The numerical and experimental distribution are not comparable in their outliers, which is discussed with more detail in Section 4.1.

Statistics of particle impact angles and directions
The statistical distributions of impact and direction angles of the near-wall flow structures over the test surface are shown in Figures 9 and 10, respectively, along with their experimental counterparts. The color of each data point of the raincloud indicates velocity of the respective impact. The corresponding statistics of the distributions are summarized in Table 2.
In the case of impact angle, the numerical and experimental distributions are similar in that the median value differs by less than 2°and both distributions are skewed towards lower angles. Most of the particle impacts in the studied turbulent flow set-up are then producing wear scars at relatively low impact angles ( 10°). Further, the lower impact angles tend to be associated with higher impact velocities, which is apparent in both numerical and experimental data. The experimental distribution has an inherent cut-off of the impact angle because it is obtained from analysis of particle scars, where particles with insufficient kinetic energy cannot generate a measurable scar, resulting in not being represented in the statistics. On the other hand, the statistics derived from numerical simulation comprise all the impact angles and impacts direction which are calculated in the full velocity range of the specimen's near-wall flow structures.
The calculated distributions of impact direction present a higher spread and with a median value slightly distant from zero (5.41°) when compared with the experimental data. The directionality of impacts is indicative of upward and downward translation of particles following the fluctuation of flow structures. This observation is equivalent with no preferred direction observed in the experiment (positive and negative direction angles).  However, the direction is bound in an angle range and this bounded fluctuation has similar behavior in the numerical data, reflecting the formation of at least temporary flow structure. Further, the experimental distribution is more centred, whereas a flatter and slightly bimodal distribution is obtained from the simulations. This observation indicates the presence of one coherent flow structure in front of the experimental target surface, whereas the selection of mesh points in the simulation falls in between two coherent flow structures. Numerical data shows a higher velocity concentration on the extremes of the distribution and lower velocity at the centre, reinforcing the idea of the presence of two coherent flow structures. Similar tendencies are found in the data of azimuth and inclination angles shown in Figure 10 with the corresponding summary of statistics summarized in Table 3. Azimuth and inclination angles are related to impact angle and impact direction, showing a phase shift close to 90°, maintaining the same trends. There are two outliers in the experimental data of the azimuth angle as compared to the numerical azimuth angle and experimental impact directions. These outliers could be attributed to an error in the experimental data handling. Because the azimuth and impact direction angle are related, these outliers are not presented in the experimental impact direction. However, their impact is negligible in the context of the total information.

Discussion
The importance of impact angle and particle velocity in understanding the erosion wear mechanisms is widely accepted, as these parameters feature in most of the erosion models (Lospa et al., 2020). Although these models were validated experimentally, their validity is based on the consistency between expected and observed mass loss produced under given conditions of exposure. Such an approach is necessarily integrative in that information on the cumulative effect of all the impacts is available, without considering possible modes and tendencies of individual impacts. Only recently, the erosion by solid particles determined from single impacts was shown to improve the global prediction of dry erosion, contradicting phenomenological models based on mass loss . In the following, we discuss the possible implications of the statistics of impact velocity and impact angles study for erosion by turbulent slurry flow and the applicability for calibration and validation of numerical models.

Impact velocities projected from turbulent flow
The values of impact velocities projected onto the target surface obtained from the simulation are generally consistent with experimental results. The central values fall between those derived from the experimental data, which were obtained by inverse analysis applying two different erosion models. Unlike the simulation results, the experimental values present outliers of higher deviation. This observation might be caused by choosing the time frame from which the simulated statistics were generated. Whereas the statistic of the simulation features 5 s of exposure, the experimental values result from 60 min exposure during which instabilities, not observed during the simulated time, might have occurred. In particular, flow structures comprising higher speeds relative to the surface might arise, as suggested by the evolution of instantaneous velocity, an example of which is shown in Figure 6(b).
On the other hand, the experimental method is limited to collecting information on particles that impinge at sufficient velocity and impact angle to produce a wear scar, limiting the lower bound of quantifiable velocities. Also, it should be considered that the experimental data rely on the accuracy of the erosion model used to inverse determine the velocity. The models involve elastoplastic deformation caused by particle impacts describing the erosion mechanism as relation between deformation induced in a material of given mechanical properties and the essential parameters of impact: velocity and impact angle. The reliability of the models depends on how representative all the involved parameters are, including the precision Table 4. Similarity analysis of experimental and velocity distributions shown in Figure 8: P-value of pairwise comparisons using Wilcoxon rank sum test (applying Holm method as P value adjustment) for impact velocity distributions. of their validation. For instance, while both Cheng and Huang models are derived from a single impact rationale, the former was validated using single impact but the latter using wear of cumulative impacts. Figure 9 reveals that the higher impact velocities are associated with lower impact angles. On the other hand, the experimental data of impact angle α present a cut-off in the lower bound. Both these effects can be explained by the relative magnitudes of normal and tangential components of the velocity vector and their relationship with wear. For very low impact angles, the normal component of velocity is very low, reducing the probability of erosive impact. For this reason, in a post-hoc analysis, the velocity distribution obtained from numerical simulations has been filtered by the lowest impact angle evidenced in experimental results, which resulted in the distribution envelope of impact velocity to have a shape resembling the experimental results; in particular, those featuring Cheng´s model (Figure 8). Since the experimental method of inverse analysis does not necessarily capture the entire velocity spectrum as discussed in Molina et al. (2020), to improve the accuracy of determining the cut-off impact angle and cut-off impact velocity, the model of erosion would have to be carefully revised to correctly describe the particle-target interaction and the corresponding formation of erosive scars, which is beyond the scope of the present work.
Notwithstanding, sufficient data is available to assess the statistical similarity of the velocity distributions obtained from the experimental and numerical data. The analysis is conducted using a Kruskal-Wallis rank-sum test (Kruskal & Wallis, 1952) complemented by pairwise comparisons using the Wilcoxon rank-sum test, applying the Holm method as P-value adjustment (Holm, 1979;Wilcoxon, 1945). Table 4 summarizes the results for all the velocity distributions shown in Figure 8. According to the P-values of pairwise comparisons, the numerical data can be considered almost identical with the experimental data analyzed using Cheng's model.

Analysis of impact conditions
The simulated distribution of impact angle has a box plot's central values very similar to the experimental one ( Figure 9). The shape of the distribution's envelope is similar, presenting several high angle outliers and skewness towards lower values. The median value of impact angle is notably low, which is characteristic of the turbulent flow of a diluted slurry in a slurry pot configuration, related to the tangential velocity of the inner cylinder's boundary condition. As discussed earlier, the experimental results do not include shallow impact angles captured by the simulated data because the normal component of impact velocity is not expected to produce a measurable wear scar.
The simulated distribution of impact directionality is centred about a more positive number as compared with the experimental distribution centred closer to zero. This observation, along with the distinctive shape of the distribution, could be explained by the structure of the turbulence. The chaotic nature of turbulent flow allows considering it as a sum of coherent structures in time and space, with the possibility of forming periodic or pseudo-periodic flow patterns such as the average flow shown in Figure 6(b). Since the particles carried by these coherent structures follow the flow pattern, this effect is evidenced in the statistics of the impact directionality. The specific location of the coherent structure persistent in time is very sensitive to initial conditions and boundary conditions so that the precise distribution in an experiment cannot be controled. Also, a structure persistent in time results in average directionality different from zero, whereas averaging over longer time periods cancels the effects to central angle value of zero. Then, the difference between the central angle of the experimental and numerical distributions of directionality could be explained by the difference in the specific location of the coherent structures in the experimental run and a possibility of its shift in the duration of the experiment.
The simulated distributions of azimuth and inclination angles are consistent with observations made for impact angle and impact direction, which can be expected since they represent the same impact condition but observed in a plane rotated by 90°. However, they provide additional insight into the flow structure. The inclination angle corresponds tothe directionality of impacts in the axial direction. In this case, the shift of the central value from zero is indicative of a twisted T-C roll. The inclination angle corresponds to the impact angle observed in the axial direction, so that care must be taken when using these values to associate with erosion damage.
In the simulated data, some particles are found to impact the surface from an opposite direction, which is not the case of experimental data in which the azimuth angle is limited to < 90° (Figure 10). To explain this discrepancy, it should be first recalled that impact velocity is defined in relative terms, i.e. by the velocity difference between the particle and the surface. Further, impacts with relative velocity greater than the tangential velocity of the inner cylinder imply that some bulk flow structures drag particles in the opposite direction of the azimuthal flow, crossing the boundary layer flow in local and brief time events. These effects could not have been detected in the experimental data because the algorithm used for automatic analysis of experimental erosion scars included the assumption of φ < 90°limitation when processing images of the eroded surfaces (Molina et al., 2019(Molina et al., , 2020. Further, the experimental and numerical distributions of the impact, direction, azimuth and inclination angles are compared by means of a Welch two sample t-test, considered superior to other non-parametric methods applicable for this type of data sets (Skovlund & Fenstad, 2001). The null hypothesis for this test is that there are no differences between the distributions being compared. The effect of the statistical sample size was determined by Cohen's d estimate (Cohen, 2013) to obtain the net differences between experimental and numerical distributions. From the results summarized in Table 5, it can be stated significance level of 0.001, that there is no statistical difference in experimental and numeric distributions for impact and azimuth angles, and that their net difference is negligible. On the other hand, there is a significant statistical difference in the experimental vs numerical data of both impact direction and inclination angle. Both observations are consistent with the geometrical relation between the corresponding angles.
It is interesting to remark that experimental data shows no statistically relevant differences for the various particles sizes used in the experiment, neither for impact nor for inclination angles measured. This could be explained by high difference between eddy turnover time and the particles' Stokes number, where independently of the particle size, the inertia effects do not significantly affect in the particle trajectories into the flow. This provides evidence for supporting that the particles follow the flow trajectories of VLSM and LSM mainly. The simulated distributions of particle impacts are richer from the information point of view but they show  sufficient consistency with experimental data to elaborate on linking individual impacts with the generation of wear scars. In order to propose a suitable approach for filtering the data, we examine relevant correlations between impact angle, impact direction and the components of impact velocity. Figure 11 shows the scatter plots of impact angles and impact directions against either the corresponding magnitude of impact velocity or its normal component only. The correlations are quantified by adjusted R-squared, results of which are summarized in Tables 6 and 7. It is possible to observe that for impact direction (δ) there is no correlation with neither total impact velocity nor its normal component. This observation is consistent with the inherent chaotic action of turbulence. The impact angle (α) is weakly negatively correlated with the magnitude of impact velocity but a strong positive correlation is observed against the normal component of impact velocity. In other words, for impact angles near zero the corresponding velocity is mainly tangential to the surface, and increasing impact angle increases the amount of energy conveyed through normal component of impact velocity. This observation can be explained by the viscous layer near the surface because it drifts flow structures with low impact angles into the tangential direction, and the flow structures with higher impact angles are thus reduced. In any case, this correlation between impact angle and normal component of impact velocity implies that the probability of producing a wear scar of given type will depend on the properties of both the particle and the target surface, but also the near-wall structure flow. Since the particular mechanism of erosion may further change under a given combination of impact velocity and impact angle, the design of a suitable filter for identifying potential erosive impacts will have to take into account the mechanics of the involved processes. However, the experimental results present outliers that the numerical model does not reproduce. The origin of these outliers are likely unaccounted for interactions at the interface between fluid (numerical model) and material (experimental data). The outliers in experimental data may arise form experimental error, or they could be part of the physical response. To explore this aspect, it could be interesting to extend the simulation in time and conduct a detailed time-space analysis to check whether the outliers emerge from numerical limitations or whether they would be feature of the experimental data only.

Method of numerical model calibration and validation
This work considers the particles as flow tracers because in high turbulent flow with low particle concentration and determined particle Stokes numbers with the potential to generate measurables scars, it could be possible to dismiss the fluid-particle interaction effect for focusing on the near-wall liquid phase behavior. The proposed technique returned promising results in comparison with experimental data because the statistical study did not establish that the numerical simulation's impact features, such as velocity and angles, are drawn from a separate population of impact data collected experimentally. Then it is reasonable to infer that the experimental and numerical distributions are derived from the same populations, implying that the numerical model properly represents near-wall flow structures. Then the slurry pot test and impact velocity and angles measurement could be suitable for calibration and validation of numerical modeling.

Applicability of the numerical model for erosion study
Although a LES-WALE model does not capture the entire spectrum of turbulent energy, it is still suitable for the study of erosive wear because, as already pointed out, small flow structures do not carry enough energy to convey particles capable of producing wear scars. Since the computational effort focuses on solving the large eddies carrying potentially erosive particles, the computational time becomes attractive for practical applications as compared to the more accurate DNS approach. However, for the prediction of wear rate on larger surfaces and extended periods of time, the stability of flow structures in time must be considered, presenting a computational challenge but also an opportunity to explore configurations not feasible for DNS.
Further, it should be kept in mind that the very mechanics of particle-target interaction that eventually produce erosive wear, pose an intrinsic source of error in that validation of the model relies on the accuracy of characterization of the wear loss. The conventional method of weight loss is not applicable because it does not capture the localized, i.e. 'per particle', wear. On the other hand, analysis of individual wear scars is only as accurate as the employed surface characterization technique. In this work, we rely on literature data obtained by automated processing of optical profilometry data. This data permits inferring the impact features as velocity and impact angles, but their accurateness relies on the mechanical model used for describing the process of materials deformation. The materials science of such models is advanced, but a generalized model is still missing due to multiple size scales and variety of involved deformation mechanisms, restricting their validity for specific materials only.
Although the numerical approach presented in this work provides new insights on the near-wall conditions of flow relevant for erosion, it should be kept in mind that there are some limitations in applicability to erosion studies. The CFD-based approach can only be used for diluted slurry because no particle-particle interaction is taken into account. For higher particle content, such interactions produce particle trajectories that do not follow the streamlines of the liquid phase and would thus have to be computed otherwise or included as a factor in post-processing of the calculation. Also, low-inertia of the particles is a strong assumption. Any fluid-particle interaction would modify the results, the importance of which increases with particle mass and particle shapes different from spherical. However, the increasing complexity of the model to capture the effect of particle inertia and particle-particle interaction would imply increasing computational effort diminishing the practicality of the simulation.
Further, the particle Stokes number is a fundamental number used for describing particles' dynamic response under turbulent flow. In this work, the Stokes number is not 1 required for full consistency of particle and fluid flow, but the statistical comparison between numerical and experimental results and the Stokes number evaluated in the integral time scale permit us to infer that the dynamic behavior of particles is sufficiently close to the fluid dynamics for the purpose of erosion study. The results could thus provide a tool for analyzing the particle near-wall behavior applicable for engineering processes. Our work admittedly does not overcome the problems of limiting Stokes number, but it could be interesting to extend this work in the future in order to provide an effective rather than absolute St limit. The experimental work of Molina et al. (2020) does not present statistical differences in velocity and angles of the impacts produced by particle sizes ranging from 276 μm to 774 μm. The corresponding Stokes number ranges from 0.1722-1.3542, respectively; nevertheless, the exact limitation of applicability should be corroborated in future works.
The particularity of the WALE turbulence model is to provide a physical model of the near-wall viscous layer rather than a numerical approach as is the case of other LES approaches. The results obtained in this work (e.g. no difference in impact angles) support the applicability of WALE and it could be interesting to replicate this methodology in future works to compare different turbulence models and their near-wall behavior.
An LES-WALE model could also be a suitable Eulerian base to be coupled with a Lagrangian model for studying erosion mechanisms in the range of slurry concentration interesting for practical application of hydraulic transport.

Conclusion
The results of the numerical study of particle impact conditions from turbulent flow of diluted slurry onto a rotating surface presented in this work, allow for the following conclusions: • The turbulent flow of a diluted slurry in the set-up of rotating cylinder can be simplified under the Eulerian framework, providing meaningful prediction of impact conditions of individual particles onto a target surface. This simplification enables development of dynamic models of erosive wear based on localized impact conditions distributed in time and space. • The LES-WALE approach describes the near-wall flow structures under a range of turbulent regimes, providing an efficient alternative for simulation of near-wall problems of fluid-material interaction. • The consistency of the calculated statistics of impact angles, impact directions, and velocity with those derived from an experiment presents an alternative approach for the study of erosive flows because it is focused specifically on the impacts that effectively produce wear. In particular, the CFD models of particleladen flows that are conventionally described by the Eulerian-Lagrangian framework could be calibrated. • Although the impact angle of individual particles cannot be controled in a slurry pot set-up, depending on the flow conditions, a fairly narrow distribution can be produced. In the configuration studied here, most of the impact angles were ∼ 10°, which is a regime especially suited for the study of wear in ductile alloys.
The statistical approach of quantifying the flow conditions onto a metallic target presented in this work is a demonstration of a new concept for predicting erosive wear that could be extended to complex geometries. However, there are two major limitations that need to be addressed for such type of simulation to be reliable. First, the range of slurry characteristics in terms of concentration and particles size/shape to which this modeling approach is applicable has to be explored. The inertial effect of the particles might result in a deviation that could be corrected by post-processing at least for a range of characteristic that is of practical interest for hydraulic transportation of slurries. Second, the accuracy of experimental validation is subject to a mechanical understanding of the erosion mechanisms as shown by the discrepancies in experimental data derived from the Cheng and Huang models. Thus, the mechanisms of erosion and the possible contribution of corrosion must be clarified in order to allow for accurate integration of individual impacts to predict the overall wear rate.

Disclosure statement
No potential conflict of interest was reported by the author(s).

Funding
This work was supported by Agencia Nacional de Investigación y Desarrollo: [Grant Number FONDECYT Regular 1201547].