Nomenclature

IGVs

inlet guide vanes

LE

leading edge

PS

pressure side

SS

suction side

TE

trailing edge

Symbols

$A$

area (m²)

c

chord (m)

${C_D}$

drag coefficient

${C_{LS}}$

slip-shear lift coefficient

${C_M}$

slip-rotation lift coefficient

${C_R}$

rotational coefficient

d

diameter of particle (m)

E

equivalent erosion rate (mg/g.cm²)

f

force by mass (N/kg)

g

gravity (m/s²)

H

pressure head (mm of water)

k

turbulent kinetic energy (m²/s²)

h

specific enthalpy (J/kg)

i

turbulence intensity level (%)

m

mass (kg)

$\vec n$

normal vector

p

pressure (Pa)

P

power (W)

r

radius (m)

Re

Reynolds number

$R{e_p}$

particle Reynolds number

$R{e_s}$

shear flow Reynolds number

$t$

time (s)

$\vec t$

tangential vector

$\vec{V}$

velocity (m/s)

Greek characters

$\beta$

impact angle (deg)

$\delta$

penetration (microns)

$\epsilon $

local erosion rate (mg/g)

$\varepsilon$

turbulence rate of dissipation (m²/s³)

$\phi$

particle shape factor

$\rho$

density (kg/m³)

$\mu$

dynamic viscosity (kg/m.s)

$\overrightarrow {\Omega} $

blade speed of rotation (rad/s)

$\overrightarrow {{\omega _f}} $

fluid rotation (s⁻¹)

$\overrightarrow {{\omega _{fr}}} $

relative fluid rotation (s⁻¹)

$\overrightarrow {{\omega _p}} $

angular velocity of particle (s⁻¹)

${\vartheta _s}$

step angle (deg)

τ

tip clearance (mm)

Subscripts

h

hub

f

fluid

i

inlet

o

outlet

p

particle

r

radial

rot

rotor

s

shroud, step

θ

tangential

t

total, tip

z

axial

1, 2

impact/rebound

2.0 Fan stage model

The studied turbomachinery is a variable speed axial fan stage. It includes 19 inlet guide vanes (IGVs) (Fig. 1(a)) and a rotor with 11 blades (Fig. 1(b)). The entire fan is made from an aluminium alloy. An optical coordinate measuring machine was used to obtain the 3-D geometry of the components. Table 1 provides the list of some design parameters of this axial fan stage.

Figure 1. Axial fan stage: (a) intake nose and IGVs; (b) rotor blades.

3.0 Flow field solution

The flow field and solid phase are solved independently since this particulate flow is assumed dilute and of one-way coupling. The 3D flow is modelled by using the simplified steady Reynolds Average Navier Stokes (RANS) equations and the Unsteady Reynolds Average Navier Stokes (URANS) equations, with the k-ω Shear Stress Transport (SST) [Reference Menter33] as the turbulence model.

(1a) \begin{align}\frac{{\partial \rho }}{{\partial t}} + \frac{\partial }{{\partial {x_j}}}\left( {\rho {W_j}} \right)\end{align}

(1b) \begin{align}\frac{\partial }{{\partial t}}\left( {\rho {W_i}} \right) + \frac{\partial }{{\partial {x_j}}}\left( {\rho {W_i}{W_j}} \right) = - \frac{{\partial p}}{{\partial {x_i}}} + \frac{\partial }{{\partial {x_j}}}\left[ {{\tau _{ij}} - \rho \overline {{w_i}{w_j}} } \right] + {S_{{M_i}}}\end{align}

(1c) \begin{align}\frac{\partial }{{\partial t}}\left( {\rho {h_t}} \right) - \frac{{\partial P}}{{\partial t}} + \frac{\partial }{{\partial {x_j}}}\left( {\rho {W_j}{h_t}} \right) = \frac{\partial }{{\partial {x_j}}}\left( {\lambda \frac{{\partial T}}{{\partial {x_j}}} - \rho \overline {{w_j}h} } \right) + \frac{\partial }{{\partial {x_j}}}\left[ {{W_i}\left( {{\tau _{ij}} - \rho \overline {{w_i}{w_j}} } \right)} \right] + {s_E}\end{align}

where $\overrightarrow {{S_M}} = - \rho\big( 2\vec \Omega \times \vec W + \vec \Omega \times ( {\vec \Omega \times \vec r} )\big)$ , $\vec \Omega $ rotational velocity, $\overrightarrow {W} $ relative flow velocity, ${\tau _{ij}}$ molecular stress tensor and $\rho \overline {{w_i}{w_j}} $ Reynolds stresses. ${h_t}=h+ \frac{W_i W_i}{2}+k$ with $h_t$ and $h$ are the total and static specific enthalpies, and the turbulent kinetic energy $k= \frac{1}{2} \overline{w_i^2}$ , and finally ${s_E}$ is the energy source term.

Table 1. Geometry parameters

The governing equations are integrated over each control volume defined by joining the centres of edges and the element centres surrounding each node. The volume integrals are discretised within each element sector and accumulated to the control volume to which the sector belongs, whereas the surface integrals are discretised at the central integration points of each surface segment. The flow properties are stored at mesh nodes and all approximations utilise a linear shape function $\varphi = \mathop \sum \limits_{i = 1}^{{N_{node}}} {N_i}{\varphi _i}$ . To prevent pressure field oscillations resulting from a non-staggered collocated grid arrangement, a coupled solver solves the flow equations as a single system. In the steady simulation, the time step control may use either a timescale based on the domain geometry and flow field or different time scales, or a fixed time scale $1/\Omega$ . The convergence is set at a root mean square (RMS) threshold of 10⁻⁶.

There are two options available for modelling rotating parts in CFX code: The frozen rotor interface is the most effective when there are significant circumferential variations in the flow properties compared to the component pitch. This interface connects the stationary and rotating frames in a way that ensures fixed relative positions, allowing for a steady solution to be obtained for the multiple frames of reference (MFRs) and accounting for interactions between them. On the other hand, the stage interface is more suitable when the circumferential variations of flow properties are of the order of the component pitch. This interface uses a quasi-steady algorithm that averages the discrete fluxes through the interface and applies them as boundary conditions to the adjacent zone. This interface is commonly used for predicting the aerodynamic performance.

3.1 Computational domain

The computational domain depicted in Fig. 2 consists of sectors, including extended inlet domain, intake nose, two IGVs, one rotor blade and a step diffusing duct. The pitch ratio between the rotor blade and two IGVs is given by $\frac{{{s_{rot}}}}{{{s_{IGV}}}} = 0.8636$ , which is acceptable when using periodicity to reduce the computational effort in the particle trajectory computation.

Figure 2. Computational domain.

The fluid is air treated as a real gas with the reference properties set at 25°C. The total pressure and temperature values imposed at the inlet correspond to the standard atmospheric conditions, while a mass flow rate is imposed at the outlet of the diffusing duct. The turbulence intensity level [34] set at inlet is estimated by ${\rm{i}}\left( \% \right) = 100 \times 0.16\; Re_{{d_h}}^{ - 1/8}$ = 3.07 % for $R{e_{{d_h}}} = 5.44 \times 10^5 $ , but was increased to 5% as recommended in many turbomachinery simulations.

3.2 Mesh generation

Hexahedral mesh blocks are distributed to secure fine meshing near the the leading edge and trailing edge of the vane and blade, as well as the hub and shroud. There are 12 mesh lines around the walls to solve the boundary layers, while 15 meshlines are used in the tip gap for the leakage flow. Figure 3(a) displays the midspan meshes for two IGVs and one blade. The meshes of vanes, hub and extended nose are shown in Fig. 3(b) while those of blade, hub and step are depicted in Fig. 3(c).

Figure 3. Mesh details: (a) mid-span; (b) vanes and extended nose; (c) blade and step.

The first layer of nodes near walls corresponds to the distance $\Delta y $ targeting the parameter ${y^ + } = 2$ . The dimensionless parameter $y^+=\Delta y \frac{\rho}{\mu} U_t$ considers for the simplification the friction law of a flat plate $U_t=W_\infty \sqrt{c_f/2}$ , where ${c_f} = 0.027Re_c^{ - 1/7}$ . $W_{\infty}$ is the free velocity outside the boundary layer. The Reynolds number $R{e_c} = \rho {W_\infty }c/\mu $ is based on the mean chord of vane or rotor blade. When assesed at the nominal point the values of $R{e_c}$ are 172,531 and 630,157, for the vanes and blade, respectively. Figure 4(a) illustrates that the values of ${y^ + }$ in IGVs are low on both sides, except on the hub, not exceeding the value of 36.5. In the rotor blade, the values are very low except at the front corner of blade and the tip where the values are less than 13.4.

Figure 4. y⁺ values: (a) IGVs; (b) rotor blade.

The grid size dependency study at the nominal point demonstrates that the total-to-total isentropic efficiency stabilises (Fig. 5) beyond the fourth grid size of a total of 725,328 nodes. Moreover, the erosion assessed with sand MIL-E5007E at a concentration of 700 mg/m³ reveals slight changes (Fig. 5) in the hourly eroded mass from the vane and the blade above the fourth grid size.

Figure 5. Grid size dependency.

3.3 Flow solution

The high-resolution scheme is used in discretising the advection term and the k-ω based SST turbulence model. By setting the local time step factor to the value of $5$ a good convergence was reached at RMS values about 10⁻⁶. The aerodynamic performance parameters were calculated by adopting the mixing plane (stage interface), which utilises a quasi-steady algorithm where the discrete fluxes through the interface are circumferentially averaged and transferred to the adjacent zone. The mass averaging of different flow properties at the inlet and outlet of the fan stage allowed computing the total pressure head $\;{H_t} = \frac{{{{10}^3}\Delta {P_{to - i}}}}{{{\rho _w}g}}$ ( $mm{H_2}O)$ with $\Delta {P_{{t_{o - i}}}} = {P_{{t_o}}} - {P_{{t_i}}}\;$ and the total-to-total isentropic efficiency ${\eta_{is\;t - t}} = \dot m\frac{{\Delta \mathop {{P_{{t_{o - i}}}}}\limits^ \cdot }}{{{{\bar \rho }_a}{{\bf{P}}_{\rm{i}}}}}$ , with ${\rho _w}$ is the water density and $\overline {{\rho _a}} $ is the air average density. The input power ${{\rm{P}}_{\rm{i}}} = 11\smallint \overrightarrow {\Omega} .\vec r \times \left( { - p + \tau } \right)d\vec s\;$ is obtained from the torque integrated over the 11 blades. As shown in Fig. 6(a), this axial fan stage exhibits a wide operating range between 1 and 9 kg/s and a total pressure head of 56–533 mm of water column, while the rotational speed varies between 2,500 and 6,000 rpm. Figure 6(b) indicates that for each speed-line, there is a peak of efficiency, while the maximum peak efficiency of 78.42% (Fig. 6(b)) occurs at the nominal air mass flow rate of 3.716 kg/s and a rotational speed of 4,500 rpm.

Figure 6. Fan performance: (a) total pressure head; (b) total-to-total isentropic efficiency.

The flow solution allowed analysing the flow characteristics in this axial fan stage and provided the necessary data for the particle trajectory computations. Figure 7 depicts that the static pressure is decreasing between the vanes while it is increasing through the rotor blade passage. The static pressure continues to increase by diffusion behind the rotor step and through the duct, where the rotational kinetic energy converts into a static pressure. Figure 8 depicts that after a large stagnation region at the intake nose, the flow accelerates through the admission section and between the vanes, where the highest velocities are seen on the suction side. The wakes and vortices generated by the vanes are convected downstream and reach the leading edge of the blade, hence altering the flow incidence. When crossing the interface, the fluid gains a high prewhirl velocity component and the relative flow velocity increases significantly. Flow acceleration continues in the blade-to-blade passage to reach the highest velocities in the front part of blade. Figure 8(a) shows that at midspan the relative flow velocity equals 132.1 m/s and more near the blade tip (Fig. 8(b)), reaching up to 147.1 m/s. The flow is well guided in the mid-sections but becomes misaligned upper sections, as depicted in Figs. 8(b, c). The most complex flow features occur in the tip clearance region (Fig. 8(c)), where the flow is diverted transversally and mixes with the annulus wall boundary layer. Figure 9 illustrates the evolution of the relative Mach number. At midspan, the Mach number has a maximum of 0.383 at the front of suction side (Fig. 9(a)), while at the blade tip (Fig. 9(b)) has a maximum of 0.427. Figure 10 shows that the leakage flow caused by the pressure gradient creates a flow stream that tends to curl and forms a vortex structure known as the tip leakage vortex, which is convected by the mainstream. Additionally, a horseshoe vortex that emanates from the leading edge is enveloped by the leakage vortex and extends further along the suction side. These local flow details have a significant impact on the movement of the small solid particles.

Figure 7. Static pressure along the meridional plane.

Figure 8. Flow velocities: (a) at midspan; (b) near blade tip; (c) in tip gap.

Figure 9. Mach number: (a) at midspan; (b) near blade tip.

Figure 10. Tip vortical structure.

To compute the particle trajectories and assess the erosion development, while considering different blade positions, first the flow field was solved and then the data was transferred to the trajectory code.

The relative motion of the rotor blade is simulated using the transient rotor/stator interface. A total of four complete rotations of the rotor were simulated over a period of 53.33 ms when the fan operates at the nominal point (N = 4,500 rpm, m = 3.716 kg/s). The high-resolution advection scheme and the second order backward Euler transient scheme were used, while adopting a time step of 3.704x10⁻² ms equivalent to a 1 degree step angle. To initiate the unsteady computations the steady flow solution considering the frozen rotor interface were used. In the final round, five positions of the blade were chosen to analyse the flow patterns across the blade pitch. Thereafter, the particle trajectories and erosion were calculated, considering two types of sand particle seeded globally at the inlet of the fan. According to Fig. 11 the reference position ‘P0’ corresponds to the vane trailing edge facing the blade leading edge. The other pitchwise positions P1, P2, P3 and P4 are incremented by a step angle ${\vartheta _s} = $ 6.545 deg up and down of the reference position P0.

Figure 11. Blade positions.

Figures 12 present the evolution of static pressure with the blade defilement, indicating that the variation in pressure loading is more significant at the front part of the blade, and the diffusion is greatly influenced by the position of the blade relatively to the vane. The wakes convected downstream of IGVs and blades, along with the secondary flows at the hub and casing, are the major sources of entropy creation that changes with the blade position, mainly due to alterations of flow incidence and stagnation point, as well as the vortices produced by IGVs. Looking at the entropy patterns after the trailing edge of the blade for different positions, as depicted in Fig. 13, it is evident that the reference position P0, where the trailing edge of the vane is facing the leading edge of the blade results in minimum losses. On the other hand, position P1 experiences the highest losses, followed by position P4.

Figure 12. Evolution of static pressure.

Figure 13. Entropy downstream the blade TE.

4.0 Particle trajectory computation

There are two approaches commonly used in simulating particulate flows. The Eulerian approach treats the flow and dispersed phase as two interpenetrating continua, where the time averaged transport equations for the additional phase are solved. On the other hand, the Lagrangian approach focuses on tracking individual particles as they are driven by the airflow. For dilute gas-particle flows with particle volume fraction below 4x10⁻⁷ [Reference Sommerfeld35], a one-way coupling approach is typically used. A simple comparison indicates that the Lagrangian model is more cost-effective [Reference Abbas, Koussa and Lockwood36] and capable of providing a detailed physical description of the particle phase and interaction with walls.

The movement of particles is governed by the balance between the rate of momentum change and the external forces. Equation (1) involves the forces (N/kg) reduced by the mass of particle. The left side represents inertia, Coriolis, and centrifugal forces, while the right side includes external forces such as drag $\overrightarrow {\;{f_D}} $ , combined gravity and buoyancy force $\overrightarrow {\;{f_{GB}}} $ , slip shear layer induced (Saffman) force $\overrightarrow {{f_S}} $ , and Magnus force $\overrightarrow {{f_M}} $ . The pressure gradient force and added mass force are not considered due to the particle’s density being significantly higher than that of air.

(2) \begin{align}\left[ {\frac{{{d^2}{r_p}}}{{d{t^2}}} - {r_p}{{\left( {\frac{{d{\theta _p}}}{{dt}} + \omega } \right)}^2}} \right]\overrightarrow {{e_r}} + \left[ {{r_p}\frac{{{d^2}{\theta _p}}}{{d{t^2}}} + 2\frac{{d{r_p}}}{{dt}}\left( {\frac{{d{\theta _p}}}{{dt}} + \omega } \right)} \right]\overrightarrow {{e_\theta }} + \frac{{{d^2}{z_p}}}{{d{t^2}}}\overrightarrow {{e_z}} = \overrightarrow {\;\;{f_D}} + \overrightarrow {{f_{GB}}} + \;\overrightarrow {{f_S}} + \overrightarrow {{f_M}} \end{align}

The specific drag force is the dominant force. When developed in terms of drag factor and Reynolds number $R{e_p} = \frac{{{\rho _f}}}{{{\mu _f}}}{d_p}\left| {\overrightarrow {{V_f}} - \overrightarrow {{V_p}} } \right|$ this force becomes:

(3a) \begin{align}\overrightarrow {{f_D}} = \frac{3}{4}\frac{{{\mu _f}R{e_p}}}{{{\rho _{pd_p^2}}}}{C_D}\left( {\overrightarrow {{V_f}} - \overrightarrow {{V_p}} } \right)\end{align}

For very low $R{e_p} \lt 0.1$ the drag coefficient corresponds to the Stokes regime $\;{C_D} = 24/R{e_p}$ . For $R{e_p}$ between 0.01 and $2.6 \times 10^5$ , the general drag coefficient [Reference Haider and Levenspiel37] is given as follows:

(3b) \begin{align}{C_D} = \frac{{24}}{{R{e_p}}}\;\left( {1 + aRe_p^b} \right) + \frac{c}{{1 + \frac{d}{{R{e_p}}}}}\end{align}

The constants $\;a$ , b, c and $d$ are based on the shape factor $\phi = {A_S}/{A_P}$ ( ${A_S}$ area of the sphere of equivalent volume and ${A_P}$ is the actual area). The constants $\;a$ , $b$ , $c$ and $d$ are revaluated by Haider and Levenspiel [Reference Haider and Levenspiel37] based on the shape factor.

\begin{align*}a = exp\!\left( {2.3288 - 6.4581\phi + {{2.4486}\phi^2}} \right)\end{align*}

\begin{align*}b = 0.0964 + 0.5565\phi\end{align*}

\begin{align*}c = exp\!\left( {4.905 - 13.8944\phi + {{18.4222}\phi^2} - {{10.2599}\phi^3}} \right)\end{align*}

(3c) \begin{align}d = exp\!\left( {1.4681 + 12.2584\phi -20.7322 {\phi^2} + {{15.8855}\phi^3}} \right)\end{align}

The combined gravitational and buoyancy force:

(4) \begin{align}\overrightarrow {{f_{GB}}} = \left( {1 - \frac{{{\rho _f}}}{{{\rho _p}}}} \right)\vec g\end{align}

There are two lift forces: the slip shear layers induced force resulting from non-uniform pressure distribution known as the Saffman force, and the slip-rotation lift force known as Magnus force.

Saffman [Reference Saffman38] produced the first expression for the slip-shear lift force, which when developed in 3-D and reduced by particle mass becomes:

(5a) \begin{align}\overrightarrow {{f_{S}}} = \frac{3}{4}\frac{{{\rho _f}}}{{{\rho _p}}}\;{C_{LS}}\!\left( {\overrightarrow {{V_f}} - \overrightarrow {{V_p}} } \right) \times \overrightarrow {{\omega _f}} \end{align}

With the slip-shear lift coefficient [Reference Sommerfeld35]:

(5b) \begin{align}\;{C_{LS}} = \frac{{4.1126}}{{\sqrt {R{e_s}} }}f\!\left( {R{e_p},\;R{e_s}} \right)\end{align}

The correction function $f\!\left( {R{e_p},\;R{e_s}} \right)$ is proposed by Mei [Reference Mei39] for particle Reynolds number in the range $0.1 \le R{e_p} \le 100$ and $\beta = 0.5\frac{{\;R{e_s}}}{{\;R{e_p}}}$

\begin{align*}f\!\left( {R{e_p},\;R{e_s}} \right) = \left( {1 - 0.3314\sqrt \beta \;} \right){e^{ - R{e_p}/10}} + 0.3314\sqrt \beta \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;R{e_p} \le 40\end{align*}

(5c) \begin{align}f\!\left( {R{e_p},\;R{e_s}} \right) = 0.0524\sqrt {\beta R{e_p}} \,\,\;R{e_p} \gt 40\end{align}

The shear flow Reynolds number ( $R{e_s}$ ) is as follows

(5d) \begin{align}R{e_s} = \frac{{{\rho _{fd_p^2}}}}{{{\mu _f}}}\|\overrightarrow {{\omega _f}}\| {\rm{\;\;\;\;\;\;\;\;\;where}}\,\,\overrightarrow {{\omega _f}} = \nabla \times \overrightarrow {{V_f}} \end{align}

Crowe et al. [Reference Crowe, Sommerfeld and Tsuji40] developed an expression for the Magnus force, function of $R{e_p}$ , the Reynolds number of particle rotation ( $R{e_r}$ ) and the slip-rotation lift coefficient ( ${C_M}$ ). When reduced by particle mass, this force becomes as follows:

(6a) \begin{align}\overrightarrow {{f_M}} = \frac{3}{4}\frac{{{\rho _f}}}{{{\rho _p}}}\frac{{R{e_p}}}{{R{e_r}}}{C_M}\!\left[ {\overrightarrow {{\omega _{fr}}} \times \left( {\overrightarrow {{V_f}} - \overrightarrow {{V_p}} } \right)} \right]\end{align}

Where the Reynolds number of particle rotation:

(6b) \begin{align}R{e_r} = \frac{{{\rho _{fd_p^2}}}}{{{\mu _f}}}\|\overrightarrow {{\omega _{fr}}}\| {\rm{\;\;\;\;\;\;\;\;\;where}}\qquad \overrightarrow {{\omega _{fr}}}= 0.5\nabla \times \overrightarrow {{V_f}} - \overrightarrow {{\omega _p}} \end{align}

From the torque acting on a rotating particle the particle rotation speed, in terms of the rotational coefficient ${C_R}$ [Reference Sommerfeld35] and the flow gradient, is easily integrated through the time step as follows:

(7a) \begin{align}\frac{{d\overrightarrow {{\omega _p}} }}{{dt}} = \frac{{15}}{{16\pi }}\frac{{{\rho _f}}}{{{\rho _p}}}{C_R}\|0.5\nabla \times \overrightarrow {{V_f}} - \overrightarrow {{\omega _p}}\| \!\left( {0.5\nabla \times \overrightarrow {{V_f}} - \overrightarrow {{\omega _p}} } \right)\end{align}

According to Dennis et al. [Reference Dennis, Singh and Ingham41] the rotational coefficient is given by:

\begin{align*}{C_R} = \frac{{64\pi }}{{R{e_r}}}\qquad R{e_r} \le 32\end{align*}

(7b) \begin{align}{C_R} = \frac{{12.9}}{{Re_r^{0.5}}} + \frac{{128.4}}{{R{e_r}}}\qquad 32 \lt R{e_r} \le 1000\end{align}

4.1 Impact conditions

In the near-wall region (Fig. 14) the correct impact is set at a distance equal to half the particle diameter to prevent non-physical trajectories and impacts. Thus, an interpolated time step is used to get the exact impact. At a wall collision the kinetic energy implies erosion damage and the restitution coefficients are used to define the changes in the particle velocity amplitude and direction, knowing the local unit vectors $\vec n\;$ and $\vec t\;$ . The restitution coefficients for the aluminum 2024 are used for the vane, blade and casing, based on the experimental data obtained by Tabakoff et al. [Reference Tabakoff, Hamed and Murugan42], presented in terms of mean and standard deviation:

(8a) \begin{align}{e_v} = {V_{p2}}/{V_{p1}} = \mathop \sum \limits_{i = 0}^4 {a_i}\beta _i^i\;\qquad{\rm{and}}\qquad{\sigma _{\left( {{e_v}} \right)}} = \mathop \sum \nolimits_{i = 0}^4 {a_{\sigma i}}\beta _i^i\end{align}

(8b) \begin{align}{e_\beta } = {\beta _{p2}}/{\beta _{p1}} = \mathop \sum \limits_{i = 0}^4 {b_i}\beta _i^i\;\qquad{\rm{and}}\qquad{\sigma _{\left( {{e_\beta }} \right)}} = \mathop \sum \nolimits_{i = 0}^4 {b_{\sigma i}}\beta _i^i\end{align}

Figure 14. Impact conditions.

After a collision the model developed by Sommerfeld [Reference Sommerfeld43], extended to 3-D by Breuer et al. [Reference Breuer, Alletto and Langfeldt44], is used to determine the new angular velocity after an impact. ${e_n}$ and ${e_t}$ stand for the normal and tangential restitution coefficients. The non-sliding and sliding collision conditions depends on the relative particle velocity ${\vec V_{pr1}}$ and the limit velocity ${v_{lim}}$ , as well as the static and dynamic coefficients of friction [Reference Sommerfeld43].

(9a) \begin{align}{\vec \omega _{p2}} = {\vec \omega _{p1}} + \frac{{10}}{7}\frac{{\left( {1 + {e_t}} \right)}}{{{d_p}}}\left( {\vec n \times {{\vec V}_{pr1}}} \right)\;\;\;\;\;\;\;\;\;\left| {{{\vec V}_{pr1}}} \right| \le {v_{lim}}\end{align}

(9b) \begin{align}{\vec \omega _{p2}} = {\vec \omega _{p1}} - 5\frac{{\left( {1 + {e_n}} \right)}}{{{d_p}}}\left( {{{\vec V}_{p1}}.\vec n} \right)\;\frac{{{\mu _{dy}}}}{{\left| {{{\vec V}_{pr1}}} \right|}}\;\left( {\vec n \times {{\vec V}_{pr1}}} \right)\quad \left| {{{\vec V}_{pr1}}} \right| \gt {v_{lim}}\end{align}

4.3 Particle seeding

The trajectories of particles and erosion are computed when the fan operates at the nominal point. Sand particles are randomly released from various points at the inlet, distributed both radially and tangentially, following the specified particle size distribution and concentration. In compliance with the US standards and engine dust ingestion testing requirements, two sand particle concentrations were specified: a low concentration of 53 mg/m³ and a high concentration of 700 mg/m³, representing a highly polluted environment (MIL-E-8593E [Reference Shoemaker and Shurnate47]). Two distributions of sand particles (Fig. 15) were considered: AC-coarse sand (0–200 μm) and MIL-E5007E sand (0–1000 μm). The first distribution has a mean size of 37.7 microns and a standard deviation of 30.3 microns, while the second distribution has a mean of 237.1 microns and a standard deviation of 164.5 microns. The sizes of particles are randomly distributed conform to the particle size distributions (Fig. 15), while a subprogram iterates on particle diameters and release positions until the convergence in the total particle mass rate. Initial particle slip velocities are imposed based on experimental data [Reference Ghenaiet, Tan and Elder11] knowing the local flow velocities.

Figure 15. Size distributions: AC-coarse (0–200 μm) and MIL-E5007E (0-1000 μm).

5.0 Erosion assessment

Early research on erosion focused on the physical properties dependence of wear, such as the amplitude and direction of particle velocity during collisions, particle size and concentration, among other factors. Experimental results have indicated that erosion is affected by many factors such as the time duration, particle impact velocity and angle, and the restitution coefficients of the particle-target material. Finnie [Reference Finnie48] was the first to establish the basis for the sand erosion, which was later enhanced by Bitter [Reference Bitter49] to incorporate the effects of plastic deformation. The most effective erosion prediction models for aerospace materials originated from the experiments of Grant and Tabakoff [Reference Grant and Tabakoff50], which primarily involved aluminium-based alloys. The erosion rate, measured as the amount of material removed (in milligram) per mass of impacting particles (in gram), is influenced by two mechanisms: the first is dominant at low impact angles while the second is near normal impact angles.

(10a) \begin{align}\epsilon = {k_1}{\left[ {1 + ck\left( {{k_{12}}sin\left( {\frac{{90}}{{{\beta _0}}}{\beta _1}} \right)} \right)} \right]^2}V_{p1}^2{\cos ^2}{\beta _1}\!\left( {1 - R_\theta ^2} \right) + {k_3}{\left( {{V_{p1}}\sin {\beta _1}} \right)^4}\end{align}

The tangential restitution ratio ${R_\theta }$ :

(10b) \begin{align}{R_\theta } = 1-0.0016{V_{p1}}\sin {\beta _1}\,\,{\rm{and}}\,\,ck = \left\{ {\begin{array}{*{20}{c}}{1\;\;\;\;if\;\;\;{\beta _1} \le 2{\beta _o}}\\[5pt] {0\;\;\;\;\;if\;\;\;{\beta _1} \gt 2{\beta _o}}\end{array}} \right.\end{align}

For aluminium alloy and quartz particle, the material constants are: $\;{\beta _0} = 20\;deg$ , $\;\;{k_1} = 3.67 \times {10^{ - 6}}$ , ${k_{12}} = 0.585$ , $\;{k_3} = 6 \times {10^{ - 12}}$ .

The local erosion rate evaluated from equation (10) is used to calculate the mass erosion at each impact point. The values of mass erosion per unit of time (s) are summed over the centred area ${A_{nod}}\;$ around the given node. To obtain the erosion rate density $\;E\left( {mg/g.c{m^2}} \right)$ , the unit of $c{m^2}$ is chosen instead of $m{m^2}$ to have better scales. Finally, the erosion rate density is distributed among the neighbouring nodes of the face elements.

(11a) \begin{align}{E_{nod}} = \frac{{\mathop \sum \nolimits_1^{ni} {\epsilon _i}{m_{pi}}}}{{{A_{nod}} \ }}\end{align}

At all initial positions of the rotor blade, the values of mass erosion of the rotor blade are summed along the pitchwise distance and averaged to get the erosion rate density as follows:

(11b) \begin{align}{\bar E_{nod}} = \frac{{\mathop \sum \nolimits_1^{{n_s}} \mathop \sum \nolimits_1^{ni} {\epsilon _i}{m_{pi}}}}{{{} \ \mathop \sum \nolimits_1^{ns} {A_{nod}}}}\end{align}

The depth of penetration for the given time Δt is estimated from the mass erosion at each node considering the material density and the centred area ${A_{nod}}$ around each node.

(12) \begin{align}{ \delta_{nod}} = \frac{\overline E_{nod}\Delta t}{\rho _{m}}\end{align}

The new coordinates are obtained using the depth of penetration and the normal unit vector. The deterioration of the blade, in terms of chord reduction and thickness were assessed for each blade section, hence resulting in average values for the entire blade.

The mean variation of blade chord:

(13) \begin{align}\Delta \bar c\!\left( \% \right) = 100\!\left[ {1-\frac{{\frac{1}{{{r_t} - {r_h}}}\mathop \smallint \nolimits_0^{{r_t} - {r_h}} \Delta {c}dr}}{{{c_m}}}} \right]\end{align}

The variation in the tip clearance:

(14) \begin{align}\Delta \tau \!\left( \% \right) = 100\!\left[ {1 - \frac{{\frac{1}{{{A_t}}}\mathop \sum \nolimits_1^n {\delta _{ti}}{A_{ti}} + \frac{1}{{{A_s}}}\mathop \sum \nolimits_1^n {\delta _{si}}{A_{si}}}}{\tau }} \right]\end{align}

6.0 Validation

The FORTRAN code PARTRAJ has been validated upon a high-speed axial fan stage (11,300 rpm) used in aircraft ventilation, made from aluminium as described in the author’s thesis [Reference Ghenaiet51]. The seeded particles consisted of MIL-E 5007E sand, (0 –1000 μm). The calculated trajectories showed impacts in IGVs concentrated at upper of pressure side due to the deflection imposed by the spherical fairing. In the rotor blade, the high frequency of impacts resulted in intense erosion spreading in a triangular region on the pressure side extending over the tip. Figure 16(a) and 16(c) depict the computed erosion patterns over the PS and SS of the rotor blade, which are in good agreement with the tested paint removal shown in Fig. 16(b) and 16(d). Furthermore, Table 2 presents steps of performance degradation, where sand particles were injected at a high concentration (775 mg/m³) while the fan operated near the nominal point of 0.8 kg/s. The tested performance degradation corresponded to the real eroded fan blade while the computed performance degradation corresponded to the predicted blade geometry deterioration. Initially, the fan stage pressure rise coefficient and efficiency were equal to 0.3952 and 0.802, respectively and the stall mass flow was equal to 0.616 kg/s. After 9 hours of sand ingestion, the measured performance decreased to the values of 0.3816 and 0.7487, respectively, while the computed performance based on the predicted eroded fan blade declined to 0.3402 and 0.7264, respectively. The same is observed for the stalling mass flow which become 0.710 kg/s as measured and 0.7384 kg/s as predicted. The predicted performance degradation aligned with the test results, however, at later periods the differences could be attributed to the eroded leading edge and the rounding of the tip corner that significantly influenced the stall limit and increased the leakage flow.

Figure 16. Comparison of blade erosion: (a) computed erosion rate density on the PS; (b) test paint removal from PS; (c) computed erosion rate density on the SS; (d) test paint removal from SS.

7.0 Results and discussions

Figure 17 shows a sample of trajectories related to released small sand particles (10 microns), distributed in both radial and circumferential directions. The small-size particles tend to follow the flow path (Fig. 17(a)) and are affected by the turbulence and secondary flows. Some of these particles deviate upwards near the trailing edge of the blade. Due to the effects of centrifugation and leakage flow, these particles move from the pressure side through the tip gap (Fig. 17(b)) and in the course induce erosion on the blade tip and shroud. In Fig. 18, it can be observed that impacts on the IGVs are caused by small particles spreading over the rear half of pressure side, while few impacts are seen on the suction side near the trailing edge, related to particles bouncing off the rotor blade. When the small particles pass through the rotor blade, they collide at high velocities near the leading edge of blade, but at lower velocities near the trailing edge. A distinct strip of impacts is observed on the suction side near the leading edge. Local erosion rates are the highest near the leading edge of the blade from both sides, but have lower values towards the trailing edge and hub.

Figure 17. Trajectories of 10 μm sand particles released: (a) at centreline, (b) near shroud.

Figure 18. Impacts coloured by velocities related to 10 μm sand particles.

Figure 19 illustrates paths taken by the mid-size particles of 50 microns which significantly deviate from the flow streamlines due to their higher inertia. They are shown to collide with the nose and bounce off to reach the outer duct and vanes. By crossing the interface plane, these particles reach the rotor blade at high incidence angles. When released close to the centreline (Fig. 19(a)), they impact the nose and bounce off twice before being deflected upwards to collide with the pressure side of the IGVs. Due to their high direction angle and imparted centrifugation, these particles hit the leading edge of the blade and repeatedly impact the pressure side at higher radii. When released near the shroud (Fig. 19(b)), the particles impact and bounce off the IGVs multiple times before reaching the rotor interface. By reaching the blade, they hit the blade tip and cross through the gap. Figure 20 provides examples of induced impacts and their corresponding velocities.

Table 2. Measured and computed performance degradation at the nominal point

Figure 19. Trajectories of 50 μm sand particles released: (a) at centreline, (b) near shroud.

Figure 20. Impacts coloured by velocities related to 50 μm sand particles.

Large particles, such as 200 microns (Fig. 21) of significant inertia deviate considerably from the flow streamlines. The particles arriving straight are shown to hit and bounce off the pressure side of IGVs, while those previously deflected by the nose shape hit the upper sections and the rear part of pressure side. Large particles released near the centreline (Fig. 21(a)), they first hit and bounce off the nose before being deflected upwards and impacting the pressure side of the IGVs. Due to high centrifugation, the particles are shown to hit the leading edge of the blade and repeatedly impact the pressure side at larger radii. Particles released near the shroud (Fig. 21(b)) repeatedly impact the vanes. By crossing the rotor interface, the particles are shown to hit the blade near the tip, and few are likely to pass through the tip gap. Figure 22 depicts impacts related to large-size particles spreading on the pressure side of the blade and cause high erosion rates. Owing to the relatively high rotational speed and flow velocities, as well as the pre-whirl imposed by IGVs, the large particles strike the blade leading edge at high velocities (Fig. 22), inducing impacts extending over the upper corner from the pressure side towards the tip and trailing edge of the blade. Many particles after hitting the leading edge of blade bounce back and reach the trailing edge of the vane.

Figure 21. Trajectories of 200 μm sand particles released: (a) at centreline, (b) near shroud.

Figure 22. Impacts coloured by velocities related to 200 μm sand particles.

Figure 23 displays the paths (coloured by particle size) followed by AC-coarse sand particles (0–200 microns), when released upstream the inlet duct. These particles are shown to strike the intake nose and deflect upwards, while others hit the intake lip and deviate downwards. Multiple impacts are seen on the pressure side of the IGVs. Through the rotor the large particles are likely to impact the pressure side of the blade, while impacts related to smaller particles are less. Small particles are seen to cross the tip gap as entrained by the leakage flow. Downstream the rotor, the large particles follow ballistic trajectories, while the small particles join the residual rotation in the diffusing flow. The sand particles, varying in sizes and release points, result in uneven impacts on the pressure side, as shown in Fig. 24. The leading edge and upper part of the blade pressure side experience high impact velocities and subsequently high local erosion rates. Similar trends are observed (Fig. 25) with MIL-E5007E sand (0–1000 μm), being highly deflected towards the shroud and the hub due to initial collisions with the intake nose and lip. Due to their larger sizes they cause more impacts with the leading edge and the pressure side, mainly at upper sections of the blade. Also, there are several particles shown (Fig. 25) to transversally cross via the tip gap. At the exit from the rotor blade, the particles highly deviate circumferentially and follow ballistic trajectories, while the smaller particles stick to the rotational flow downstream the rotor. As noticed from Fig. 26, the second type of sand particles caused high frequency of impacts and erosion rates, spreading on the pressure side, which is more evident in the upper corner region from the leading edge of the blade.

Figure 23. Trajectories of AC-coarse sand particles, released globally at the inlet.

Figure 24. Impacts coloured by velocities related to AC-coarse sand.

Figure 25. Trajectories of MIL-E5007E sand particles, released globally at the inlet.

Figure 26. Impacts coloured by velocities related to MIL-E5007E sand.

The next computations were conducted to analyse the erosion wear development for two types of sand particles: AC-coarse (0–200μm) and MIL-E5007E (0–1000 microns), globally seeded upstream of the intake at concentrations of 53 and 700 mg/m³. The initial position of rotor blade was varied according to Fig. 11 to investigate its effect on erosion patterns. When AC-coarse sand is seeded at a concentration of 53 mg/m³, the highest erosion rate density (mg/g.cm²) occurs at the leading edge of the blade and extends over the upper tip corner with a spotted erosion, as shown in Fig. 27(a–e) with blade positions. The development of erosion is related to the direct exposure to dense fluxes of particles entering the blade at a high incidence imposed by the large pre-whirl flow velocity, resulting in high frequency of impacts with high velocities. The blade tip exhibits evidence of erosion due to particles migrating from the pressure side to the suction side through the tip gap. There is no erosion observed on the blade suction side, except near the leading edge. The predicted erosion patterns in Fig. 27(a–e), show that the position of blade greatly influenced the erosion patterns and rate densities. All the depictions reveal a noticeable erosion of the leading edge and tip corner of the blade. When comparing between Fig. 27(b) and Fig. 27(c), it can be seen that as the blade moves from its initial position, the eroded area near the blade root extends from the leading edge. Similarly, when comparing between Fig. 27(b) and Fig. 27(d), the eroded area towards the trailing edge from the pressure side spreads more. In contrast, the erosion contours on the blade suction side remain mostly the same, except the slight variations near the root and tip of blade. In IGVs, the impacts are mainly on the pressure side, resulting in almost uniform erosion contours. The highest erosion rates within the four corners of the vanes are explained by the high frequency of impacts caused by the intake deflection. The erosion patterns of the vanes seem slightly affected by the variation of blade position; however, the particles bouncing from the leading edge of blade and shroud caused alteration to the erosion patterns. Erosion contours with low rate densities are seen upstream and around the blade at the downward sections characterized by fewer impacts since the vanes caused upwards deflection to the particles. Moreover, the erosion of the rotor’s shroud (Fig. 28(a–e)) spreads more from the lower part at blade entrance and the upper part of blade from the trailing edge. Despite the high centrifugation of particles while crossing the rotor blade passage, the erosion rate densities on the shroud are lower compared to the blade surfaces. When comparing between Fig. 28(b) and Fig. 28(d), it is evident that the erosion patterns on the hub and shroud are clearly affected due to variation of the blade position.

Figure 27. Erosion patterns caused by AC-coarse sand at the concentration of 53 mg/m³, with the blade positions.

Figure 28. Erosion patterns of rotor hub and shroud caused by AC-coarse sand at the concentration of 53 mg/m³, with the blade positions.

The erosion patterns differ when seeding sand particles (MIL-E5007E) at the concentration of 700 mg/m³. The erosion distribution looks more scattered with spotted regions of extreme rate densities compared to the first sand type. The patterns in Fig. 29(a–e) show that the highest erosion rate density (mg/g/cm²) occurs at the leading edge of the blade, as well as the upper corner from both sides in addition to the tip of blade. This is mainly due to the dense flux of particles arriving at high incidence added to the centrifugation of heavy particles. The blade tip experiences noticeable erosion typically by small particles crossing via the tip gap. In contrast, the blade suction side exhibits erosion only from the leading edge. Similarly, the predicted erosion patterns (Fig. 29(a–e)) clearly depend on the blade position. Near the blade root, the erosion patterns develop more, depending on the blade position, as observed when comparing Fig. 29(b) with Fig. 29(e). Also, an area of erosion towards the trailing edge from the pressure side tends to spread more, as revealed when comparing Fig. 29(a) with Fig. 29(e). On the blade’s suction side, however, the erosion contours remain almost similar, except the slight variations near the root and tip of blade where the spots of erosion become clearer in the case of Fig. 29(d) and Fig. 29(e) compared to Fig. 29(b) and Fig. 29(c). Impacts on IGVs occupy a large area of the pressure side, but the erosion rates are significantly lower compared to the rotor blade. Moreover, the erosion seems more scattered by comparing these patterns with those obtained with AC-coarse sand. Actually, the erosion patterns are more concentrated at the trailing edge from the hub corner, as the particles bounce back after impacting the blade, which is also revealed on the suction side near the trailing edge. On the rotor hub, the scattered erosion patterns (Fig. 30(a–e)) upstream and around the blade are related to the particles deflected after striking the IGVs’ shroud. The rotor’s shroud is characterised by areas of erosion spread, especially at the opposite of blade tip. Accordingly, the patterns of erosion of the hub and shroud are influenced by the blade position, which is more evident when comparing Fig. 30(a) with Fig. 30(e).

Figure 29. Erosion patterns caused by MIL-E5007E sand at the concentration of 700 mg/m³, with the blade positions.

Figure 30. Erosion patterns of the rotor hub and shroud caused by MIL-E5007E sand at the concentration of 700 mg/m³, with the blade positions.

It is evident that the erosion patterns and material removal are influenced by the concentration and sizes of particles, type of material being targeted, duration of exposure, and noticeably the position of the blade facing the trailing edge of the vane. The erosion process depicted in the above figures, resulted in blunting of the blade leading edge, gradual reduction of blade chord, increase in tip clearance, and alteration of the upper corner of blade that becomes rounded. These patterns of erosion have been confirmed by Ghenaiet et al. [Reference Ghenaiet, Tan and Elder11, Reference Ghenaiet, Tan and Elder18]. The results of particle trajectory and erosion, considering various particle size distributions, concentrations and rotor blade positions are presented in Tables 3 and 4 for the vane, and in Tables 5 and 6 for the rotor blade. It can be inferred from these findings that MIL-E5007E (0–1000 μm) sand led to higher erosion of the rotor blade compared to AC-coarse (0–200 μm) sand, while the opposite is true for the vane erosion. These tables demonstrate that the blade position significantly impacts the total eroded mass, overall erosion rate and blade geometry deterioration, all of which are expressed as percentages variations, after 10 hours of sand ingestion at concentrations of 53 and 700 mg/m³.

Table 3. IGV erosion caused by AC-coarse sand after 10 hours

Table 4. IGV Erosion caused by MIL-E5007E sand after 10 hours

Table 5. Rotor blade erosion caused by AC-coarse sand after 10 hours

Table 6. Rotor blade erosion caused by MIL-E5007E sand after 10 hours

The curves relating the mass removal (erosion) and percentages are shown in Fig. 31(a) for AC-coarse sand and in Fig. 31(b) for MIL-E5007E sand, at concentrations ranging from 53 to 700 mg/m³. These figures demonstrate that the eroded mass from the blade varies depending on the blade position along the pitch, while the erosion rates increase with the concentration. When seeding AC-coarse sand the maximum erosion is found at position ‘P4’, whereas the minimum occurs at position ‘P1’. On the other hand, with MIL-E5007E sand the highest erosion occurs at positions ‘P3’ and ‘P2’, while the lowest erosion is at the reference position ‘P0’. This suggests that the effect of blade position is also dependent upon the particle sand type. Figure 32(a) shows the overall erosion of the blade after 10 hours of sand ingestion and the percentage of mass reduction cumulated from the five blade positions. It is evidenced, that the erosion parameters of the blade increase with the particle concentration, where the highest values occur with MIL-E5007E sand. This is also observed for the vane in Fig. 32(b), but in contrast more erosion is obtained with AC-coarse sand. The deterioration of the blade’s geometry after 10 hours of sand ingestion at different concentrations is illustrated in Fig. 33(a). The changes in blade chord and tip clearance are noticeable and more significant at higher concentrations, with the most deterioration by MIL-E5007E sand. The noticeable variations in the computed values of the blade deterioration are attributed to the cumulated values obtained from the five blade positions, as well as the randomness in the release points, sizes and rebound factors. According to Figs. 32 and 33, after 10 hours of MIL-E5007 sand ingestion at extreme concentration, the eroded mass from the blade is about 0.29%, the tip chord reduced by a 0.45% and the tip clearance increased by a 0.23%. On the other side, AC-coarse sand caused 0.23% eroded mass, 0.4% decrease in tip chord and 0.16% increase in tip clearance. It seems that AC-coarse sand induces more erosion of the vane of about 0.0038% while that with MIL-E5007E is 0.0018%, which is significantly lower compared to the blade erosion.

Figure 31. Effect of blade position on blade erosion after 10 hours of sand ingestion: (a) AC-coarse; (b) MIL-E5007E.

Figure 32. Cumulated eroded mass after 10 hours of sand ingestion: (a) blade; (b) vane.

Ultimately, the prediction of blade material loss and geometry deterioration can be effectively used to generate an erosion fault model serving to simulate the fan stage performance degradation, based on the similar approach developed by Ghenaiet et al. [Reference Ghenaiet, Tan and Elder16, Reference Ghenaiet, Tan and Elder18].

8.0 Conclusion

The computations of the particulate flows and erosion helped in better understanding the erosion process in this category of axial fans.

The results indicate that the erosion is greatly influenced by the type, size, and concentration of particles. Moreover, the position of the rotor blades relative to the vanes has a significant impact on the development of erosion and the maximum amount of erosion wear.

• The MIL-E5007E sand caused more erosion on the blade compared to the AC-coarse sand. In contrast, AC-coarse sand caused more erosion of the IGVs.

• The highest frequency of impacts and erosion rates are observed at the leading edge of the blade due to dense fluxes of particles and pre-whirl imposed by IGVs.

• The pressure side of the blade is eroded at the upper sections and tip corners due to high centrifugation and leakage flow.

  

  Figure 33. Cumulated geometry deterioration after 10 hours of sand ingestion.

• The small particles are seen to cross through the tip gap and cause erosion of the blade tip.

• The blade suction side is only eroded along a narrow strip from the leading edge.

• The scatters of erosion observed at the shroud are caused by centrifugation and particles deflected upwards, while the hub revealed lesser erosion.

• The IGVs exhibit erosion spreading over the entire pressure side, but with low rate densities. The concentrated erosion near the trailing edge are related to the particles bouncing back after impacting the rotor blade.

• The present results clearly reveal that the blade position influences the erosion patterns and rates. Furthermore, the positions of rotor blade leading to the maximum and minimum erosion are well depicted.

• It is pertinent to consider the relative positions of the blades when computing the particle trajectory and erosion in order to get more realistic qualitative and quantitative results.

Given the limitations associated with experimental campaigns, the numerical tools become valuable to underline the process of erosion, assess the material loss and geometry deterioration and predict the performance degradation. By identifying the critical areas of erosion and their parameters, this help selecting adequate coating protection for this category of axial fan design.