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

Abstract

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.

1 Introduction

Turbomachinery blades are quite often subjected to particle-laden-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 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.

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 main components of the integrated method given in [3] are the Streamline-Upwind/Petrov-Galerkin (SUPG) [11] and Pressure-Stabilizing/Petrov-Galerkin (PSPG) [12] stabilizations, a finite element particle-cloud tracking (PCT) method [1,2,3] with one-way dependence, an erosion model based on two time scales, and the solid-extension mesh moving technique (SEMMT) [13,14,15,16,17].

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 stabilization 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, 22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43]). 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.

2 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 target-geometry 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.

3 Mathematical model

3.1 Navier–Stokes equations of incompressible flows

Let $\varOmega \subset {\mathbb {R}}^{nsd}$ be the spatial domain with boundary $\varGamma$, and (0, T) be the time domain. The Navier–Stokes equations of incompressible flows can be written on $\varOmega$ and $\forall t \in (0,T)$ as

$$\begin{aligned}&\rho \left( \frac{\partial {\varvec{u}}}{\partial t} + {\varvec{u}}\cdot \pmb {\nabla }{\varvec{u}} \right) - \pmb {\nabla }\cdot {\varvec{\sigma }} = {\varvec{0}}, \end{aligned}$$

(1)

$$\begin{aligned}&\pmb {\nabla }\cdot {\varvec{u}} = 0, \end{aligned}$$

(2)

where $\rho$ and ${\varvec{u}}$ are the density and velocity. The stress tensor ${\varvec{\sigma }} =- p{\varvec{I}} + 2\mu {\varvec{\varepsilon }}({\varvec{u}})$, where p is the pressure, ${\varvec{I}}$ the unit tensor, $\mu =\rho \nu$ is the viscosity, $\nu$ the kinematic viscosity, and the strain rate ${\varvec{\varepsilon }}({\varvec{u}})=(\pmb {\nabla }{\varvec{u}}+\pmb {\nabla }{\varvec{u}}^T)/2$. The essential and natural boundary conditions associated with Eq. (1) are represented as ${\varvec{u}}={\varvec{g}} \text{ on } \varGamma _g$ and ${\varvec{n\cdot \sigma }}={\varvec{h}} \text{ on } \varGamma _h$, where $\varGamma _g$ and $\varGamma _h$ are the subsets of the boundary $\varGamma$, ${\varvec{n}}$ is the unit outward normal vector, and ${\varvec{g}}$ and ${\varvec{h}}$ are given functions.

3.2 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:

$$\begin{aligned} {\varvec{x}}_c = \int _0^t {\varvec{v}}_c dt' + ({\varvec{x}}_c)_0. \end{aligned}$$

(3)

Here, subscript c refers to the cloud, ${\varvec{v}}_c$ is the velocity of the cloud, and $({\varvec{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 one-way dependence hypothesis according to Armenio and Fiorotto [67], reads as

$$\begin{aligned} \frac{d{\varvec{v}}_c}{dt} = \tau _R^{-1}\left( \langle {\varvec{u}}\rangle -{\varvec{v}}_c \right) +\left( 1-\frac{\rho }{\rho _p}\right) {\mathbf {a}}_\mathrm {GRAV}, \end{aligned}$$

(4)

where $\langle ~\rangle$ denotes ensemble average (defined later), $\rho _p$ is the particle material density, ${\mathbf {a}}_\mathrm {GRAV}$ is the gravitational acceleration, and $\tau _R$ is the particle relaxation time, which, for spherical particles, is

$$\begin{aligned} \tau _R^{-1}=\frac{3}{4d_p}C_D\frac{\rho }{\rho _p}\Vert \langle {\varvec{u}}\rangle -{\varvec{v}}_c\Vert . \end{aligned}$$

(5)

Here, $d_p$ is the particle diameter and $C_D$ is the drag coefficient based on the particle Reynolds number $Re_p = \frac{\Vert \langle {\varvec{u}}\rangle -{\varvec{v}}_c\Vert d_p}{\nu }$, first introduced in [68]. The Stokes number is defined as

$$\begin{aligned} Stk = \frac{\tau _R U}{L}, \end{aligned}$$

(6)

where U is the free-stream velocity and L is the flow length scale. The initial value of ${\varvec{v}}_c$ is given as ${\varvec{v}}_c(0) = \langle {\varvec{u}}\rangle |_{t=0}$.

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 $\theta$ it reads as

$$\begin{aligned} \langle \theta \rangle =\frac{\int _{\varOmega _c}\theta PDF({\varvec{x}},t)d\varOmega }{\int _{\varOmega _c} PDF({\varvec{x}},t)d\varOmega }. \end{aligned}$$

(7)

Here, $\varOmega _c$ is the cloud domain and $PDF({\varvec{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¹ as

$$\begin{aligned} PDF({\varvec{x}},t)=\frac{1}{\sqrt{(2\pi )^{nsd}}\sigma ^{nsd}}e^{-\frac{1}{2}\left( \frac{\Vert {\varvec{x-x_c}}\Vert }{\sigma }\right) ^2}, \end{aligned}$$

(8)

where $\sigma$ is the square root of the variance of particle position, and the cloud radius is $3\sigma$.

Here, we generalize the cloud shape to an ellipsoid with nsd independent semi-axis lengths. Then the particle-distribution PDF becomes

$$\begin{aligned} PDF({\varvec{x}},t)=\frac{1}{\sqrt{(2\pi )^{nsd} \det {\varvec{S}}}}e^{-\frac{1}{2}({\varvec{x}}-{\varvec{x}}_c){\varvec{S}}^{-1}({\varvec{x}}-{\varvec{x}}_c)}, \end{aligned}$$

(9)

where ${\varvec{S}}$ is the covariance matrix of particle position. We will define it in Sect. 5. The semi-axis lengths of the ellipsoid are $3\sigma _1$, $3\sigma _2$ and $3\sigma _3$, where $\sigma _1^2$, $\sigma _2^2$ and $\sigma _3^2$ are the eigenvalues of ${\varvec{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.

Combining Eqs. (4) and (5), we obtain

$$\begin{aligned} \frac{d{\varvec{v}}_c}{dt}&= C_D'\Vert \langle {\varvec{u}}\rangle -{\varvec{v}}_c\Vert \left( \langle {\varvec{u}}\rangle -{\varvec{v}}_c \right) +\left( 1-\frac{\rho }{\rho _p}\right) {\mathbf {a}}_\mathrm {GRAV}, \end{aligned}$$

(10)

where

$$\begin{aligned} C_D'=\frac{3}{4d_p}C_D\frac{\rho }{\rho _p}. \end{aligned}$$

(11)

4 Residual-based VMS

4.1 Formulation

The residual-based VMS formulation [17] of Eqs. (1) and (2) is given as

$$\begin{aligned}&\int _{\varOmega } {\varvec{w}}^h \cdot \rho \left( \frac{\partial {\varvec{u}}^h}{\partial t} + {\varvec{u}}^h \cdot \pmb {\nabla }{\varvec{u}}^h \right) d\varOmega + \int _{\varOmega } {\varvec{\varepsilon }} ({\varvec{w}}^h):{\varvec{\sigma }} d\varOmega \nonumber \\&\quad - \int _{\varGamma _h}{\varvec{w}}^h\cdot {\varvec{h}}^h d\varGamma + \int _\varOmega q^h\pmb {\nabla }\cdot {\varvec{u}}^h d\varOmega \nonumber \\&\quad + \sum _{e=1}^{n_{el}}\int _{\varOmega ^e} \frac{\tau _{SUPS}}{\rho } \left( \rho {\varvec{u}}^h\cdot \pmb {\nabla }{\varvec{w}}^h+\pmb {\nabla }q^h \right) \cdot {\varvec{r}}_M d\varOmega \nonumber \\&\quad + \sum _{e=1}^{n_{el}}\int _{\varOmega ^e}\nu _{LSIC}\pmb {\nabla }\cdot {\varvec{w}}^h \rho r_C d\varOmega \nonumber \\&\quad - \sum _{e=1}^{n_{el}}\int _{\varOmega ^e}\tau _{SUPS}{\varvec{w}}^h\cdot \left( {\varvec{r}}_M\cdot \pmb {\nabla }{\varvec{u}}^h\right) d\varOmega \nonumber \\&\quad - \sum _{e=1}^{n_{el}}\int _{\varOmega ^e}\frac{\tau _{SUPS}^2}{\rho }{\varvec{r}}_M \cdot \left( \pmb {\nabla }{\varvec{w}}^h \right) \cdot {\varvec{r}}_M d\varOmega = 0, \end{aligned}$$

(12)

where

$$\begin{aligned} {\varvec{r}}_M&=\rho \left( \frac{\partial {\varvec{u}}^h}{\partial t} + {\varvec{u^h\cdot \pmb {\nabla }u^h}}\right) -\pmb {\nabla }\cdot {\varvec{\sigma }}^h, \end{aligned}$$

(13)

$$\begin{aligned} r_C&=\pmb {\nabla }\cdot {\varvec{u}}^h \end{aligned}$$

(14)

are the residuals of the momentum equation and incompressibility constraint, and $\tau _{SUPS}$ and $\nu _{LSIC}$ are the stabilization parameters (see Sect. 4.2). The test functions associated with ${\varvec{u}}$ and p are ${\varvec{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.

4.2 Stabilization parameters

This section, included for completeness, is mostly from [3]. We first define the element length [69] in the advection-dominated limit:

$$\begin{aligned} h_{UGN} = 2\left( \sum _{a=1}^{n_{en}}|{\varvec{s}}\cdot {\varvec{\pmb {\nabla }}}N_a|\right) ^{-1}, \end{aligned}$$

(15)

where ${\varvec{s}}=\frac{{\varvec{u}}}{\Vert {\varvec{u}}\Vert }$, $n_{en}$ is the number of element nodes, and $N_a$ is the basis function associated with node a. In the diffusion-dominated limit, the element length [23] is defined as

$$\begin{aligned} h_{RGN} = 2\left( \sum _{a=1}^{n_{en}}|{\varvec{r}}\cdot {\varvec{\pmb {\nabla }}}N_a|\right) ^{-1}, \end{aligned}$$

(16)

where ${\varvec{r}}$ is the unit vector in the direction of the solution gradient:

$$\begin{aligned} {\varvec{r}} = \frac{\pmb {\nabla }\Vert {\varvec{u}}\Vert }{\Vert \pmb {\nabla }\Vert {\varvec{u}}\Vert \Vert }. \end{aligned}$$

(17)

The components of $\tau _{SUPS}$ corresponding to the advection-, transient-, and diffusion-dominated limits were defined in [23, 70] as

$$\begin{aligned} \tau _{SUGN1}&=\left( \sum _{a=1}^{n_{en}}|{\varvec{u}}\cdot {\varvec{\pmb {\nabla }}}N_a|\right) ^{-1} =\frac{h_{UGN}}{2\Vert {\varvec{u}}\Vert }, \end{aligned}$$

(18)

$$\begin{aligned} \tau _{SUGN2}&=\frac{\varDelta t}{2}, \end{aligned}$$

(19)

$$\begin{aligned} \tau _{SUGN3}&=\frac{h_{RGN}^2}{4\nu }. \end{aligned}$$

(20)

From these, $\tau _{SUPS}$ is defined as

$$\begin{aligned} \tau _{SUPS} = \left( \frac{1}{\tau _{SUGN1}^2} + \frac{1}{\tau _{SUGN2}^2} + \frac{1}{\tau _{SUGN3}^2}\right) ^{-\frac{1}{2}}, \end{aligned}$$

(21)

and $\nu _{LSIC}$ is defined [23] as

$$\begin{aligned} \nu _{LSIC} = \tau _{SUPS}\Vert {\varvec{u}}^h\Vert ^2. \end{aligned}$$

(22)

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 $\tau$ component, “$\tau _{SUGN4}$” [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.

5 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 $\varOmega _c = \bigcup _{e=1}^{n_{elc}}\varOmega _c^e$, where $\varOmega _c^e$ 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\sigma$. With that, the discretized version of ensemble averaging is written as

$$\begin{aligned} \langle \theta \rangle ^h=\frac{\sum _{e=1}^{n_{elc}}\int _{\varOmega _c^e}\theta PDF({\varvec{x}},t)d\varOmega }{\sum _{e=1}^{n_{elc}}\int _{\varOmega _c^e} PDF({\varvec{x}},t)d\varOmega }, \end{aligned}$$

(23)

where the element-level integration is performed by Gaussian quadrature.

5.1 Trajectory calculation

Spatially-discretized version of Eq. (10) is written as

$$\begin{aligned} \frac{d{\varvec{v_c}}^h}{dt} = {\varvec{a}}_c^h, \end{aligned}$$

(24)

where

$$\begin{aligned} {\varvec{a}}_c^h&=C_D'\Vert \langle {\varvec{u}}\rangle ^h -{\varvec{v}}_c^h\Vert \left( \langle {\varvec{u}}\rangle ^h -{\varvec{v}}_c^h \right) + \left( 1 - \frac{\rho }{\rho _p}\right) {\mathbf {a}}_\mathrm {GRAV}. \end{aligned}$$

(25)

Time discretization of Eq. (24) is performed with a predictor–multicorrector algorithm.

Predictor stage:

$$\begin{aligned} ({\varvec{v}}_c^h)_{n+1}^0 = ({\varvec{v}}_c^h)_n + ({\varvec{a}}_c^h)_n\varDelta t. \end{aligned}$$

(26)

Multicorrector stage:

$$\begin{aligned} ({\varvec{v}}_c^h)_{n+1}^{i+1} = ({\varvec{v}}_c^h)_n + \left( ({\varvec{a}}_c^h)_n + ({\varvec{a}}_c^h)_{n+1}^i\right) \frac{\varDelta t}{2}. \end{aligned}$$

(27)

Here the subscript n is the time level, and the superscript i is the counter for the multiple corrections. We stop the corrections when

$$\begin{aligned} \frac{({\varvec{v}}_c^h)_{n+1}^{i+1} - ({\varvec{v}}_c^h)_{n+1}^i}{({\varvec{v}}_c^h)_{n+1}^{i+1}} \le 2 \times 10^{-2}. \end{aligned}$$

(28)

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].

5.2 Parameters of the turbulent particle dispersion

For spherical clouds, in [3] one of the choices was

$$\begin{aligned} \sigma ^2&=2\langle \Vert {\varvec{u}}'\Vert ^2\rangle ^h \tau _L^2\left( 1-e^{-\frac{\tau _{SUPG}}{\tau _R}}\right) \left( \frac{t}{\tau _L} - \left( 1-e^{-\frac{t}{\tau _L}}\right) \right) \nonumber \\&\quad + \sigma _0^2, \end{aligned}$$

(29)

where $\tau _L = max(\tau _{SUPS},\tau _R)$. For ellipsoidal clouds, to define ${\varvec{S}}$, we first define ${\varvec{u}}^F={\varvec{u}}^h-\overline{{\varvec{u}}^h}$, where the overbar indicates time-averaging after a developed solution is reached. Then we write the components of ${\varvec{S}}$ as

$$\begin{aligned} S_{ij}({\varvec{x}},t)&=2\sqrt{\langle (u_i^F)^2\rangle \langle (u_j^F)^2\rangle }\tau _L^2\left( 1-e^{-\frac{\tau _G}{\tau _R}}\right) \nonumber \\&\quad \cdot \left( \frac{t}{\tau _L} - \left( 1-e^{-\frac{t}{\tau _L}} \right) \right) +(S_{ij})_0, \end{aligned}$$

(30)

where

$$\begin{aligned} \tau _L&= \max (\tau _G,\tau _R), \end{aligned}$$

(31)

$$\begin{aligned} \tau _G&=C_\mu ^{3/4}\sqrt{\frac{3}{2}}\max _i \left( \frac{\overline{(u_i^F)^2}}{\varepsilon ^D_{ii}}\right) \text{(no } \text{ sum } \text{ in } \varepsilon ^D_{ii}), \end{aligned}$$

(32)

with $C_\mu =0.09$ and $\varepsilon ^D_{ij}=2\nu \overline{\frac{\partial u_i^F}{\partial x_l}\frac{\partial u_j^F}{\partial x_l}}$ [72].

6 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.

6.1 Erosion thickness computation

The erosion thickness, calculated at the element level, is expressed as

$$\begin{aligned} e = \frac{m_e}{\rho _t}, \end{aligned}$$

(33)

where $m_e$ is the eroded mass per unit area, computed in the simulation associated with the evolution step, and $\rho _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 $\varDelta t$:

$$\begin{aligned} m_e = \sum \varDelta m_e, \end{aligned}$$

(34)

where, after a threshold particle impact count is reached,

$$\begin{aligned} \varDelta m_e = E \varDelta n_p m_p. \end{aligned}$$

(35)

Here E is the erosion rate, $\varDelta n_p$ is the particle impact count per unit area in $\varDelta t$, and $m_p$ is the particle mass. In the finite element PCT computation, $\varDelta n_p$ is calculated as

$$\begin{aligned} \varDelta n_p = C_{elem} v_{i,n,elem} \varDelta t, \end{aligned}$$

(36)

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.

Fig. 1

Impact angle in the erosion model

Fig. 2

Sketch of the erosion test bench: (1) reciprocating compressor; (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

Fig. 3

Dimensions of the injection pipe (mm)

Fig. 4

Dimensions of the erosion chamber (mm)

Fig. 5

Time evolution of the erosion zone

Table 1 Time evolution of the erosion rate

Fig. 6

Laser scan of the target at the end of the 40-min test. The lower picture shows the cross-section over the vertical plane passing through the injection-pipe axis

Table 2 Physical properties of the fluid and particle phases

Fig. 7

Flow computation domain. The impingement location is 3 cm from the injection-pipe outlet and 0.8 cm off-center, closer to the right edge of the target

Fig. 8

Flow computation mesh

6.2 Erosion rate calculation

For E used in Eq. (35), we adopt the expression given in [73]:

$$\begin{aligned} E=K(v_{i,elem})^n, \end{aligned}$$

(37)

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}$.

6.3 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\ge 1$ becomes

$$\begin{aligned} e_{SU}^{(i)}=F_{ACT}^{(i)} e^{(i)}= F_{ACT}^{(i)} \frac{E^{(i)} m_p}{\rho _t} n_p^{(i)}, \end{aligned}$$

(38)

with the scale-up factor $F_{ACT}^{(i)}$ determined from one of the two scale-up approaches described in Sect. 2. In the approach where we specify a threshold erosion value:

$$\begin{aligned} F_{ACT}^{(i)} = \frac{e_{thr}}{e_{max}^{(i)}}. \end{aligned}$$

(39)

In the approach where we specify a threshold operation period:

$$\begin{aligned} F_{ACT}^{(i)} = \frac{n_R^{(i)}}{n_{TOT}}. \end{aligned}$$

(40)

Here $n_R^{(i)}$ is the number of particles per unit area entering the PCT domain in the threshold operation period for the evolution step i, and $n_{TOT}$ is the total number of particles per unit area entering the PCT simulation domain.

Fig. 9

Starting positions for the cloud centers

Fig. 10

Flow field in the first evolution step. Isosurfaces of $\Vert {\varvec{u}}^h \times ({\varvec{\nabla }} \times {\varvec{u}}^h)\Vert$ (top) and $\Vert {\varvec{u}}^h\Vert$ (bottom-left) and $p^h$ (bottom-right) on the $X=0$ plane

7 Experiments

7.1 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.

7.2 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).

7.3 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).

7.4 Test description

The injection velocity is 15 m/s and the impingement angle is 80$^\circ$. The impingement location is 3 cm from the injection-pipe 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.

Fig. 11

Streamlines and $\overline{(u_3^F)^2}$ in the first evolution step, on the planes $X=0$ (top) and $Y=0$ (bottom)

Fig. 12

Streamlines and cloud trajectories in the first evolution step

Fig. 13

Time-histories of the cloud ellipsoid semi-axis lengths ($3\sigma _1$, $3\sigma _2$, $3\sigma _3$) in the first evolution step. We note that the time scale is the PCT computation time scale

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.

8 Computation

8.1 Problem setup

The computational setup reproduces the experimental setup. The injection velocity is 15 m/s and the impingement angle is $80^\circ$. The impingement location is 3 cm from the injection-pipe 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.

8.2 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 \times 10^{5}$ particles.

8.3 Computational conditions

In the flow computations the time-step size is $1.0 \times 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 \times 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.

Table 3 Scale-up factors and spatially-maximum eroded mass per unit surface area ($(m_e)^{(i)}_{max}$) ($\hbox {kg/m}^2$)

Fig. 14

Eroded mass per unit surface area (kg/m$^2$) at the end of each evolution-step PCT computation. The view is projected to the XY plane, and the black circle indicates the injection pipe

8.4 Results

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 $\overline{(u_3^F)^2}$ in the first evolution step. As expected, the turbulence is higher at the interface between the jet and the fluid at rest.

Fig. 15

Scaled-up erosion at the end of the last evolution step, cross-sectional view over the vertical plane passing through the injection-pipe axis

Figure 12 shows the streamlines and cloud trajectories in the first evolution step. The trajectories are nearly symmetric with respect to the YZ 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 YZ plane, while the clouds on the YZ 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.

• Phase 1: $t=0$ s to $t\simeq 2 \times 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\simeq 2 \times 10^{-3}$ s to $t\simeq 3.3 \times 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\simeq 3.3 \times 10^{-3}$ s to $t\simeq 6.5 \times 10^{-3}$ s. The clouds, having exited the high-turbulence region, travel with little change in their shapes.

• Phase 4: $t\simeq 6.5 \times 10^{-3}$ s to $t\simeq 8.5 \times 10^{-3}$ s. Some of the clouds impact the back wall of the chamber. That increases the particle dispersion again.

• Phase 5: $t\simeq 8.5 \times 10^{-3}$ s to the end of the computation. The clouds, exiting again high-turbulence region, travel with little change in their shapes.

8.5 Erosion evolution

We set $e_{thr}=95.6$ $\mu$m, which is 1/5 of the spatially-maximum erosion of 478 $\mu$m at the end of the test (see Fig. 6). Table 3 shows the scale-up factors and the spatially-maximum 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 particle–particle interactions.

9 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 sand-erosion 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.

Appendix A: Physical properties of the particles and target

See Tables 4 and 5.

Table 4 Physical properties of the red-brown corundum used in the experiments

Table 5 Physical properties of the target material Peraluman EN AW 5754