Computational analysis of particle-laden-airflow erosion and experimental verification

Computational analysis of particle-laden-airflow erosion can help engineers have a better understanding of the erosion process, maintenance and protection of turbomachinery components. We present an integrated method for this class of computational analysis. The main components of the method are the residual-based Variational Multiscale (VMS) method, a finite element particle-cloud tracking (PCT) method with ellipsoidal clouds, an erosion model based on two time scales, and the Solid-Extension Mesh Moving Technique (SEMMT). The turbulent-flow nature of the analysis is addressed with the VMS, the particle-cloud trajectories are calculated based on the time-averaged computed flow field and closure models defined for the turbulent dispersion of particles, and one-way dependence is assumed between the flow and particle dynamics. Because the target-geometry update due to the erosion has a very long time scale compared to the fluid–particle dynamics, the update takes place in a sequence of “evolution steps” representing the impact of the erosion. A scale-up factor, calculated based on the update threshold criterion, relates the erosions and particle counts in the evolution steps to those in the PCT computation. As the target geometry evolves, the mesh is updated with the SEMMT. We present a computation designed to match the sand-erosion experiment we conducted with an aluminum-alloy target. We show that, despite the problem complexities and model assumptions involved, we have a reasonably good agreement between the computed and experimental data.


Introduction
Turbomachinery blades are quite often subjected to particleladen-flow erosion. If the blade does not have sufficient surface protection, the erosion can damage it to the point of altering its aerodynamics, degrading the performance of turbomachinery system. Computational analysis of blade erosion is becoming inevitable in reducing the prototyping B cost. Reliable computational analysis can help the designers find new blade protection strategies. It can also help them to plan effective maintenance schedules. Earlier, closely-related studies [1,2] focused on the prediction of rain erosion patterns for a full-span wind-turbine blade. In [3,4], the focus was on the effect of the erosion on the aerodynamic performance of the blade, highlighting the correlation between the erosion patterns and geometry evolution. 1 The typical approach we see in the literature does not account for the changes in the flow field due to the geometry change. In most studies the flow computation is performed only on the original geometry, and then the computed flow field is used for predicting the particle transport and the particle erosion/deposit on the blade surface. It was shown in [5] that different blade geometries, at the same operational point and for the same power output, lead to different local aerodynamic fields with very different erosion patterns. Therefore we expect that even a minor modification on the critical parts of a blade section can make a noticeable difference in the flow field. Furthermore, the aerodynamics affects the particle motion and thus, the erosion itself. To account for this effect and to increase the fidelity of the blade erosion computational analysis, an integrated method was developed in [3] for predicting the time evolution of the interaction between the fluid-particle dynamics and blade erosion and geometry change.
This class of applications are characterized by small particle concentration. Therefore one-way dependence is assumed between the flow and particle dynamics, that is, particle (and cloud) motion is driven by the flow but the flow is not influenced by the particles. The concept of one-way dependence has been used in other computational engineering analyses. For example, in [6], the concept is used for computing the aerodynamic forces acting on the suspension lines of spacecraft parachutes, where the suspension lines are assumed to have no influence on the flow. In [7], the same assumption is used to study the particle-shock interaction. In [8][9][10], the assumption is used in flow-driven string dynamics in turbomachinery, where the strings are assumed to have no influence on the flow.
The combination of the SUPG and PSPG stabilizations is called "SUPS" [17]. In [3], the turbulent-flow nature of the analysis was addressed with a Reynolds-Averaged Navier-Stokes (RANS) model, with the Navier-Stokes equations discretized based on the SUPS stabilization and the two RANS closure equations discretized based on the SUPG stabilization. Here, we address the turbulent-flow nature of the analysis with the residual-based variational multiscale (VMS) method [17][18][19][20][21]. The VMS has two additional sta- 6 Department of Engineering, Lancaster University, Gillow Avenue, Lancaster LA1 4YW, UK bilization terms beyond those in the SUPS and therefore subsumes the SUPS. It has good turbulence modeling features. We compute the unsteady aerodynamic field until we reach a fully-developed flow, and then calculate from that the time-averaged flow variables. Methods like the SUPS and VMS involve stabilization parameters. There are various ways of calculating these parameters (see, for example, [2,3,7,). We will describe in a later section what stabilization parameters we use here.
The PCT method was first formulated by Baxter [44], and then further developed [1,2,[45][46][47] and improved to obtain statistically-independent results [48]. The trajectory of the particle-cloud center is calculated with finite element discretization. Of the elements of that "particle mesh," we use the ones inside the cloud, which has a trajectory-dependent size. The tracking method accounts for the drift-velocity gradient in the near-wall regions [46,49]. In [3], the clouds were assumed to be of spherical shape. Here, we generalize the cloud shape to an ellipsoid with independent semi-axis lengths. The particle-cloud trajectories are calculated based on the time-averaged flow field and a closure model for the turbulent dispersion of particles. The closure model in [3] was based on the turbulence closure variables. Here, we propose a closure model based on the diagonal components of the Reynolds stress tensor.
Tabakoff and coworkers [50][51][52] were among the earliest researchers working on numerical prediction of erosion in turbomachinery. Their methods and experiments served as a foundation for other studies in the field. Nowadays several empirical and semi-empirical erosion models can be found in the literature. The model used in [2] for rain erosion was a modified version of the Keegan model [53]. The models used in [3] were defined in terms of the mass eroded per unit mass of the impacting particles. For the rain erosion, the model from Springer et al. [54] was used, and for the sand erosion, the model from Oka et al. [55]. Here, we use the same sand erosion model. We determine the empirical parameters associated with the target material by curve fitting to the data from the experiment.
The geometry update due to the erosion has a very long time scale compared to the fluid-particle dynamics [4,56,56]. Therefore a single time-marching computation with the typical time-step size of the flow computation is not practical. Instead, we use a sequence of "evolution steps" to represent the impact of the erosion. We time-discretize the geometry evolution not in terms of a standard time step, but in terms of a threshold erosion value that we expect to alter the blade aerodynamics from its current operation pattern or a threshold operation period that we expect to be long enough to alter the blade aerodynamics.
The computation associated with an evolution step gives us the erosion distribution for a specific particle size, for a specific number of particles injected to the computational domain. We have two approaches for scaling-up the erosion distribution at each evolution step. A. We scale-up the erosion distribution by a factor that raises its maximum value to the threshold erosion value. We use the same factor to scale-up the number of particles to the actual number of particles for that evolution step. B. We scale-up the number of particles by a factor that raises it to the actual number of particles for the threshold operation period. We use the same factor to scale-up the erosion distribution.
As the target geometry evolves, we update the mesh with the SEMMT. The core mesh moving method in the SEMMT is the Jacobian-based stiffening method [57][58][59][60]. In the core method, the motion of the internal nodes is determined by solving the equations of linear elasticity. The mesh deformation is dealt with selectively based on the sizes of the elements. The selective treatment is based on the way we account for the Jacobian of the transformation from the element domain to the physical domain. The objective is to stiffen the smaller elements, which are typically placed near solid surfaces, more than the larger ones. When the method was introduced in [57][58][59], it consisted of simply dropping the Jacobian from the finite element formulation of the mesh moving (elasticity) equations. This results in the smaller elements being stiffened more than the larger ones. In the SEMMT, the thin layers of elements placed near solid surfaces are treated almost like an extension of the solid structure. In solving the equations of elasticity governing the motion of the fluid nodes, higher stiffness is assigned to the thin layers of elements compared to the other fluid elements. Two ways of accomplishing this were proposed in [13]: solving the elasticity equations for the nodes connected to the thin layers of elements separately from the elasticity equations for the other nodes, or together.
To show how the integrated method works, we perform a computation designed to match the sand-erosion experiment we conducted with an aluminum-alloy target.
In Sect. 2, we provide an overview of the integrated method. Section 3 is an overview of the mathematical model, including the PCT model. The residual-based VMS is given in Sect. 4. In Sect. 5, we describe the discretized particle equations, including the turbulent dispersion of particles. The erosion models and erosion thickness computation, including how the scale-up factors are calculated, are described in Sect. 6. The experiments are described in Sect. 7. The computation is presented in Sect. 8, and the concluding remarks are given in Sect. 9.

An integrated method
We presented in [3] an integrated method to simulate the long-term erosion of a wind-turbine blade. The method is made of a multiphase flow solver coupled with a geometry update method. We use here much of the same method, altering only some of the sub-steps. Therefore, some parts of this section, included for completeness, are from [3].
The time scales associated with the unsteady aerodynamics and turbulence, particle transport and dynamics, erosion of the target material and the change in the geometry that produces a significant variation in the average flow field and particle trajectories are different. The geometry update due to the erosion has a very long time scale compared to the fluid-particle dynamics, making a single time-marching simulation with the typical time-step size of the flow computation not practical. The basic idea is to have a sequence of evolution steps to represent the impact of the erosion. We time-discretize the evolution of the geometry not in terms of a standard time step, but in terms of a threshold erosion value that we expect to alter the target-geometry aerodynamics from its current operation pattern or a threshold operation period that we expect to be long enough to alter the targetgeometry aerodynamics. The alteration of the flow patterns leads to the alteration of the particle dynamics, which in turn alters the erosion patterns.
An evolution step is composed of the sub-steps listed below in the order they are taken.
1. Compute the flow field with the residual-based VMS. 2. Compute the particle-cloud trajectories with the PCT method [3], updated as presented in Sect. 5. 3. Compute the erosion distribution over the target surface. 4. Scale-up the computed erosion distribution using the threshold erosion value or operation period that triggers a geometry alteration. 5. Update the geometry and fluid mechanics mesh to the next configuration. The mesh update is done with the SEMMT.
The computation associated with an evolution step gives us the distribution of the erosion thickness e for a specific particle size, for a specific number of particles injected to the computational domain. Depending on the application, we have two approaches for scaling-up e at each evolution step.
A. We scale-up e by a factor that raises its maximum value e max to the threshold erosion value e thr . We use the same factor to scale-up the number of particles to the actual number of particles for that evolution step. At the end of all the evolution steps, we obtain a correlation map from the damaged configurations to the amount of particles needed to produce those configurations. This map can later be used to estimate by interpolation the geometrical configuration resulting from the amount of particles in a specific application.
B. We scale-up the number of particles by a factor that raises it to the actual number of particles for the threshold operation period. We use the same factor to scale-up e. At the end of all the evolution steps, we obtain a correlation map from the actual number of particles to the damaged configurations. This map can later be used to directly estimate the damaged configurations from a specific rate of exposure to particles during a long, specific period (e.g., 25 years) that we want to know the damage for.

Navier-Stokes equations of incompressible flows
Let Ω ⊂ R nsd be the spatial domain with boundary Γ , and (0, T ) be the time domain. The Navier-Stokes equations of incompressible flows can be written on Ω and ∀t ∈ (0, T ) as where ρ and u are the density and velocity. The stress tensor σ = −p I + 2με(u), where p is the pressure, I the unit tensor, μ = ρν is the viscosity, ν the kinematic viscosity, and the strain rate ε(u) = (∇ ∇ ∇u + ∇ ∇ ∇u T )/2. The essential and natural boundary conditions associated with Eq. (1) are represented as u = g on Γ g and n · σ = h on Γ h , where Γ g and Γ h are the subsets of the boundary Γ , n is the unit outward normal vector, and g and h are given functions.

Dispersed-phase model
Some parts of this section, included for completeness, are from [3]. Particle trajectories are simulated in a Lagrangian reference frame. Since particle concentration in this class of applications is very small (less than 10 −6 in the particle volume fraction), a one-way dependence approach can be used [61]. We use the PCT model [44] to simulate a large number of particles without tracking them individually. The PCT approach was used in turbulent particle dispersion [45,48,[62][63][64] and validated in turbomachinery and biomass furnaces [65,66]. In the PCT model, each trajectory is not for a particle, but for a group ("cloud") of particles, thus represents the evolution of the cloud position as a function of time: Here, subscript c refers to the cloud, v c is the velocity of the cloud, and (x c ) 0 is the initial position of the cloud, which is the inflow boundary in our computations.
The equation of motion for the cloud is given by the Basset-Boussinesq-Oseen formulation, which, with oneway dependence hypothesis according to Armenio and Fiorotto [67], reads as where denotes ensemble average (defined later), ρ p is the particle material density, a GRAV is the gravitational acceleration, and τ R is the particle relaxation time, which, for spherical particles, is Here, d p is the particle diameter and C D is the drag coefficient based on the particle Reynolds number , first introduced in [68]. The Stokes number is defined as where U is the free-stream velocity and L is the flow length scale. The initial value of v c is given as The ensemble average for the dispersed phase within the cloud is defined according to the hypothesis of independent statistical events, and for any ensemble-averaged quantity θ it reads as Here, Ω c is the cloud domain and P DF(x, t) is the probability density function.
In the PCT approach, the particle position distribution within a cloud is assumed to be Gaussian, and the cloud size varies in time, depending on the flow behaviour. In [3], the cloud was assumed to be of spherical shape. The PDF of the particle distribution within the cloud was given 1 as where σ is the square root of the variance of particle position, and the cloud radius is 3σ .
Here, we generalize the cloud shape to an ellipsoid with nsd independent semi-axis lengths. Then the particledistribution PDF becomes where S is the covariance matrix of particle position. We will define it in Sect. 5. The semi-axis lengths of the ellipsoid are 3σ 1 , 3σ 2 and 3σ 3 , where σ 2 1 , σ 2 2 and σ 2 3 are the eigenvalues of S, from the largest to the smallest. The eigenvectors are the corresponding axis directions. Because the semi-axis lengths are independent, the ellipsoid shape can change in time.

Formulation
The residual-based VMS formulation [17] of Eqs. (1) and (2) is given as where are the residuals of the momentum equation and incompressibility constraint, and τ SU P S and ν L SI C are the stabilization parameters (see Sect. 4.2). The test functions associated with u and p are w and q. A superscript "h" indicates that the function is coming from a finite-dimensional space. The superscript "e" is the element counter, and n el is the number of elements.

Stabilization parameters
This section, included for completeness, is mostly from [3]. We first define the element length [69] in the advectiondominated limit: where s = u u , n en is the number of element nodes, and N a is the basis function associated with node a. In the diffusiondominated limit, the element length [23] is defined as where r is the unit vector in the direction of the solution gradient: The components of τ SU P S corresponding to the advection-, transient-, and diffusion-dominated limits were defined in [23,70] as From these, τ SU P S is defined as and ν L SI C is defined [23] as Some new options for the stabilization parameters used with the SUPS and VMS were proposed in [24,30,31,34,36,39,42,71]. These include a fourth τ component, "τ SU G N4 " [36], which was introduced for the VMS, considering one of the two extra stabilization terms the VMS has compared to the SUPS. They also include stabilization parameters [36] for the thermal-transport part of the VMS for the coupled incompressible-flow and thermal-transport aligns. The element lengths and stabilization parameters introduced in [39,71] target isogeometric discretization but are also applicable to finite element discretization. The stabilization parameters given in this subsection have node-numbering invariance for all element types, but they do not extend to isogeometric discretization in an ideal way. The element length expression introduced in [42] has node-numbering invariance for all element types, including simplex elements.

Discretized particle equations
The early parts of this section, included for completeness, are from [3]. In the discretized particle equations, ensemble averaging is carried out over the discretized cloud domain where Ω e c is the cloud element and n elc is the number of elements. The cloud elements come from a fixed mesh, which we call "particle mesh," and consist of the elements of the fixed mesh within a radius of 3σ . With that, the discretized version of ensemble averaging is written as where the element-level integration is performed by Gaussian quadrature.

Trajectory calculation
Spatially-discretized version of Eq. (10) is written as where Time discretization of Eq. (24) is performed with a predictormulticorrector algorithm. Predictor stage: Multicorrector stage: Here the subscript n is the time level, and the superscript i is the counter for the multiple corrections. We stop the corrections when At each time step, the PCT model requires the computation of the cloud mean position and size, and identification of the elements contained within the cloud volume. This is done with the search algorithm described in [49].

Parameters of the turbulent particle dispersion
For spherical clouds, in [3] one of the choices was where τ L = max(τ SU P S , τ R ). For ellipsoidal clouds, to define S, we first define u F = u h − u h , where the overbar indicates time-averaging after a developed solution is reached. Then we write the components of S as where with C μ = 0.09 and ε D i j = 2ν

Erosion models and erosion thickness calculation
Some parts of this section, included for completeness, are from [3]. In Sect. 2 we described two scale-up approaches that drive the sequence of evolution steps. In the "threshold erosion value" approach, we specify e thr , and we assume that when the scaled-up e max reaches e thr , the erosion is at a level to alter the aerodynamics from its current operation pattern. The nominal target geometry plays a role in determining e thr . The e computed in a simulation associated with an evolution step depends on the current target geometry and the size and spatial distribution of the particles. We assume that all these remain unchanged during an evolution step, justifying the scale-up of e to the erosion distribution for that evolution step. In the "threshold operation period" approach, we specify an operation period that we expect to be long enough to alter the aerodynamics, and that becomes the duration associated with each evolution step. We again assume that the current target geometry and the size and spatial distribution of the particles remain unchanged during an evolution step. The scale-up factor becomes the ratio of the number of particles in that duration to the number of particles used in the simulation. More details on the scale-up factors will be given in Sect. 6.3.

Erosion thickness computation
The erosion thickness, calculated at the element level, is expressed as where m e is the eroded mass per unit area, computed in the simulation associated with the evolution step, and ρ t is the density of the target material. Following the notation in [2,3], we obtain the eroded mass by summing up the mass eroded in each time step Δt: where, after a threshold particle impact count is reached, Here E is the erosion rate, Δn p is the particle impact count per unit area in Δt, and m p is the particle mass. In the finite element PCT computation, Δn p is calculated as where C elem is the particle concentration in the element and v i,n,elem is the normal component of the particle impact velocity. Both are evaluated at the element center.

Erosion rate calculation
For E used in Eq. (35), we adopt the expression given in [73]: (2) air supply pipe; (3) pressure control valve; (4) mixing chamber; (5) electric motor; (6) particle tank; (7) cochlea; (8) particle inlet chamber; (9) injection pipe; (10) erosion chamber; (11) target; (12) particle collection chamber where K and n are empirical coefficients that depend on the impact angle (see Fig. 1) and coating material. They are determined from the experiment (see Sect. 7). We note that this model has a double-dependency on the impact angle, in determining the value of K and n, and in calculating the impact count, through the parameter v i,n,elem .

Erosion scale-up
In [3], it was assumed that the erosion is zero until a threshold impact count is reached. Here we use a simpler model by assuming that erosion occurs at any impact count. Then the scaled-up erosion thickness at an evolution step i ≥ 1 becomes with the scale-up factor F (i) AC T determined from one of the two scale-up approaches described in Sect. 2. In the approach where we specify a threshold erosion value: In the approach where we specify a threshold operation period: Here n (i) R is the number of particles per unit area entering the PCT domain in the threshold operation period for the evolution step i, and n T OT is the total number of particles per unit area entering the PCT simulation domain.

Setup
The experimental setup is composed of four main sections: air supply system, particle supply system, mixing chamber, and erosion chamber (see Fig. 2). The air supply system is composed of a reciprocating compressor, an air supply pipe, and a pressure control valve that controls the air velocity. The particle supply system is composed of an electric motor, a particle tank, a cochlea, and a particle inlet chamber. The electric motor drives the cochlea that brings particles to the inlet chamber. The inlet chamber and pressure control valve bring particles and air, respectively, to the mixing chamber. Then the particles are injected though the injection pipe (see Fig. 3) in the erosion chamber (see Fig. 4). In this chamber, particles impact the target and are then collected to the particle collection chamber.

Particles
Red-brown corundum particles are used in the experiment. They are essentially composed of Al 2 O 3 and are quite suitable for erosion applications because of their ability to quickly abrade most materials. Moreover, because of their high toughness and impact resistance, they can be reused several times. Physical properties of the particles are given in "Appendix A" (Table 4).

Target material
The target is made of Peraluman EN AW 5754, an Al-Mg alloy. The physical properties are given in "Appendix A" (Table 5).

Test description
The injection velocity is 15 m/s and the impingement angle is 80 • . The impingement location is 3 cm from the injectionpipe outlet and 0.8 cm off-center, closer to the right edge of the target. The test was carried on for 40 min. Figure 5 shows the time evolution of the erosion zone in the first 12 min of the test, and Table 1 shows the time evolution of the erosion rate in the same period. The eroded material was measured by weighing the target. Figure 6 shows the laser scan of the target at the end of the 40-min test.
We compare our computed data to the data from this test. In addition, we repeat the 40-min test for a set of impingement angles and injection velocities. We use the E values at the end of those tests to calculate, by curve fitting, K and n values (see Eq. (37)) for each impingement angle.

Problem setup
The computational setup reproduces the experimental setup. The injection velocity is 15 m/s and the impingement angle  is 80 • . The impingement location is 3 cm from the injectionpipe outlet and 0.8 cm off-center, closer to the right edge of the target. For the particle density and diameter we use the average values from the experiment. Table 2 shows the physical properties of the fluid and solid phases. The number of evolution steps is 5.

Computational domains, boundary conditions, and meshes
The flow computation domain is shown in Fig. 7. The origin of our coordinate frame is at the center of the target frontal surface. The boundary conditions are no-slip on the target and pipe outer surfaces, air jet injection velocity at the inflow boundary, zero-stress at the outflow boundary, and slip on all other boundaries. To reduce the computing time, we exclude the region downstream of the target. The mesh, shown in Fig. 8, has 4 millions of hexahedral elements, with a y + < 1 for all wall elements in the jet impingement area. The jet region has high refinement. The high-refinement zones we see in the interiors of the domain are mostly due to the block-structured nature of the mesh. We use the same mesh for the PCT computation.
We have 17 particle clouds, initially positioned at the injection-pipe outlet. Figure 9 shows the starting position for the cloud centers. The initial shape for all clouds is spherical with radius 1.87 mm. Each cloud has 1.0 × 10 5 particles.

Computational conditions
In the flow computations the time-step size is 1.0 × 10 −6 s. The number of nonlinear iterations per time step is 5, and the number of GMRES iterations per nonlinear iteration is 30. In the PCT computations the time-step size is 1.0 × 10 −5 s. In the first evolution step the flow computation is for 20,000 time steps. We extract the time-averaged quantities of Sect. 5 from the last 1000 steps. In the next 4 evolution steps the flow computation is for 10,000 time steps, and the time-averaged quantities are extracted again from the last 1000 steps. The PCT computation in all evolution steps is for 220 time steps. Figure 10 shows the flow field in the first evolution step. We see the turbulent nature of the impinging jet and a set of complex vortical structures. The jet enlarges immediately after the inlet, forming a conical volume with high turbulence production at the interface with the fluid at rest. Because the target is inclined, the flow is mostly upward after the impingement. We see a high-pressure region corresponding to the stagnation point, which is below the jet axis because the target is inclined. For the same reason, the elongation of the high-pressure region in the Y -direction is more in the upper part. Figure 11 shows the streamlines and (u F 3 ) 2 in the first evolution step. As expected, the turbulence is higher at the interface between the jet and the fluid at rest. Figure 12 shows the streamlines and cloud trajectories in the first evolution step. The trajectories are nearly symmetric with respect to the Y Z plane. The cloud centers travel basically straight, carried by the air jet, and strike the target. After the rebound, the clouds on each side cross the Y Z plane, while the clouds on the Y Z plane follow a straight  path. Figure 13 shows the time histories of the cloud ellipsoid semi-axis lengths in the first evolution step. We note that the time scale in the plots is the PCT computation time scale. We identify four phases.

Results
-Phase 1: t = 0 s to t 2 × 10 −3 s. All clouds are carried inside the jet streamtube. We see how the higher turbulence at the peripheries of the streamtube results in stretching of the clouds there, while the clouds near the tube center maintain their spherical shape until the impact. -Phase 2: t 2×10 −3 s to t 3.3×10 −3 s. The particles impact the target. After the impact the particles enter a high-turbulence region of the flow field. That increases the particle dispersion, with a rapid increase in all three semi-axis lengths. -Phase 3: t 3.3 × 10 −3 s to t 6.5 × 10 −3 s. The clouds, having exited the high-turbulence region, travel with little change in their shapes. -Phase 4: t 6.5 × 10 −3 s to t 8.5 × 10 −3 s. Some of the clouds impact the back wall of the chamber. That increases the particle dispersion again. -Phase 5: t 8.5 × 10 −3 s to the end of the computation.
The clouds, exiting again high-turbulence region, travel with little change in their shapes.

Erosion evolution
We set e thr = 95.6 μm, which is 1/5 of the spatiallymaximum erosion of 478 μm at the end of the test (see Fig. 6). Table 3 shows the scale-up factors and the spatiallymaximum eroded mass per unit area ((m e ) max ). Figure 14 shows the eroded mass per unit surface area at the end of each evolution-step PCT computation. Figure 15 shows the scaled-up erosion at the end of the last evolution step together with the corresponding test data. We can see from Fig. 14 that after an initial period of increase, the erosion area does not change much, and this is consistent with what we see in the experiment (see Fig. 5). Considering the complexity of the problem and the assumptions involved in the model, including the erosion-rate model, we think the discrepancy we see in Fig. 15 between the computed and experimental data is reasonable. Furthermore, we believe that some of the discrepancy is probably coming from not taking into account the dispersion caused by the particleparticle interactions.

Concluding remarks
We have presented an integrated method for computational analysis of particle-laden-airflow erosion. The analysis is valuable because erosion, e.g. on turbomachinery blades, over long periods of time, can degrade the aerodynamic performance. The analysis can help engineers have a better understanding of the erosion process, maintenance and protection of turbomachinery components. The analysis is challenging because it involves turbulent flows, large number of particles carried by the flow, turbulent dispersion of particles, and the target-geometry update due to the erosion has a very long time scale compared to the fluid-particle dynamics.
The main components of the integrated method are the VMS, a finite element PCT method, an erosion model based on two time scales, and the SEMMT. The turbulent-flow nature of the computational analysis is addressed with the VMS. The particle-cloud trajectories are calculated based on the time-averaged computed flow field and closure models defined for the turbulent dispersion of particles. One-way dependence is assumed between the flow and particle dynamics. Because the target-geometry update has a very long time scale, the update takes place in a sequence of "evolution steps" representing the impact of the erosion. A scale-up factor, calculated based on the update threshold criterion, relates the erosions and particle counts in the evolution steps to those in the PCT computation. The update threshold criterion is either a threshold erosion value that we expect to alter the aerodynamics of the target from its current operation pattern or a threshold operation period that we expect to be long enough to alter the aerodynamics of the target. As the target geometry evolves due to the erosion, the mesh is updated with the SEMMT, which not only protects the smaller elements from mesh deformation, but also protects the thin layers of elements near the target surface.
We presented a computation designed to match the sanderosion experiment we conducted with an aluminum-alloy target. We showed that, despite the problem complexities and model assumptions involved, we have a reasonably good agreement between the computed and experimental data. We also showed how updating the target geometry improves the quality of the computed data.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.

Appendix A: Physical properties of the particles and target
See Tables 4 and 5.