A numerical study of volcanic ash ingestion and erosion of the front components of a high bypass turbofan engine

Abstract

Airborne particles, such as dust and volcanic ash, pose a serious hazard to aircraft in flight due to their potential to cause erosion damage to engine components. It is crucial to anticipate and address the impact of erosion wear on engine performance and safety. This study aims to enhance our understanding of how volcanic ash particles behave when ingested through a high bypass turbofan engine (HBTFE) and assess the development of erosion wear in the front components. The effects of four different ash samples are assessed in various scenarios of encountering volcanic ash during cruise flight conditions. First, the flow solution is obtained for all front components, including the Pitot intake, spinner, fan, inlet guide vanes (IGVs), outlet guide vanes (OGVs), and connecting ducts. Based on the flow data, the particle motion equations are solved step by step using an in-house trajectory and erosion code. This latter adopts the Lagrangian approach, which incorporates a particle-eddy interaction model and includes probabilistic descriptions for the release positions of particles, sizes, and restitution factors. The finite element method (FEM) is used to track particles through the computational cells and determine impact positions and conditions. As a result, the Pitot intake design seems to prevent many ash particles from reaching the fan blade beyond 80% of the span. The fan blade leading edge (LE) exhibits extreme erosion on both sides. The blade’s pressure side (PS) displays erosion spreading practically on the entirety of the surface, especially near the trailing edge (TE). In contrast, the suction side (SS) has scattered erosion at lower rates. Furthermore, the rotor’s hub presents almost uniform erosion patterns, whereas the shroud depicts scattered erosion. This large fan appears to function as a separator, expelling a significant amount of ash particles through the secondary duct, thereby reducing the engine core’s susceptibility to erosion. Out of the four volcanic ash samples, those from the Kelud and Etna volcanoes appear to cause the highest hourly eroded mass, about twice as much as the samples from the Chaiten and Eyjafjallajokull volcanoes.

Nomenclature

AERD

average erosion rate density (mg/s.mm²)

FPR

fan pressure ratio

HBTFE

high bypass turbofan engine

HEM

hourly eroded mass (mg/hr)

IGV

inlet guide vane

OGV

outlet guide vane

LE

leading edge

PS

pressure side

SS

suction side

TE

trailing edge

Symbols

$A$

area (m²)

c

chord (m)

${C_D}$

drag coefficient

${C_{SL}}$

Saffman lift coefficient

c _p

specific heat at a constant pressure (J/kg.K)

d

diameter of particle (m)

E

erosion rate density (mg/s.mm²)

f

force by mass (N/kg)

g

gravity (m/s²)

h

specific enthalpy (J/kg)

i

turbulence intensity level (%)

k

turbulent kinetic energy (m²/s²)

m

mass (kg)

$\dot{m}$

mass flow rate (kg/s)

$\vec n$ , $\vec t$

normal and tangential unit vectors

$N,\;{N_d}$

rotational speed, design speed (rpm)

p

pressure (Pa)

r

radius (m)

$R{e_p}$

particle Reynolds number

$R{e_s}$

shear flow Reynolds number

T

temperature (K)

$t$

time (s)

$\vec V$

velocity vector (m/s)

Greek symbols

$\beta$

impact angle (deg)

$\epsilon$

local erosion rate (mg/g)

$\varepsilon$

turbulence rate of dissipation (m²/s³)

$\epsilon$

local erosion rate (mg/g)

$\phi$

particle shape factor

ρ

density (kg/m³)

μ

dynamic viscosity (kg/m.s)

$\overrightarrow \Omega $

blade speed of rotation (rad/s)

$\overrightarrow {{\omega _f}} $

fluid rotation (s⁻¹)

$\theta $

tangential component (rad)

Subscripts

abs

absolute

co

corrected

f

fluid

i

inlet

o

outlet

p

particle

r

radial

rel

relative

st

standard atmosphere

$\theta$

tangential

t

total

z

axial

1, 2

impact/rebound

1.0 Introduction

Gas turbine engines powering fixed or rotary-wing aircraft operate in various environments, ranging from arctic to desert regions and volcanic areas. The front components of these engines are especially vulnerable to erosion caused by airflows containing volcanic ash, dust, water droplets, and ice. As a result, the blades deteriorate, leading to a decrease in engine thrust and an increase in fuel consumption [1]. Particle ingestion happens at low and high altitudes. Low-altitude ingestion typically occurs during high-power settings at takeoff or thrust reverser use [2]. High-altitude ingestion occurs during events like dust storms or encounters with volcanic ash, which can lead to engine failure within minutes [3]. Since 1973, numerous volcanic ash encounter incidents have been reported [4]. Commercial aircraft that unexpectedly encountered volcanic ash in flight and at airports [5], have experienced in-flight power loss [6] and the melting of ash on hot sections [7–9]. The prolonged closure of European airspace after the eruption of the Eyjafjallajökull volcano in Iceland in 2010 led to significant economic losses [10]. Following this event, the European Commission established a concentration limit of 4mg/m³, which led to the creation of the Rolls-Royce ‘safe-to-fly’ chart [9, 11]. It has been suggested that ash concentration exceeding 2mg/m³ may pose safety risks [12]. Accordingly, a new approach proposed by Clarkson et al. [11, 13] considers the dose of dust ingestion event as the product of exposure duration and concentration.

Abrasive impacts are known to cause significant erosion near the LE and the tip of compressor blades [14]. Some experiments conducted by Batcho et al. [15] and Dunn et al. [16, 17] have noted that the loss in performance is attributed to factors such as increased tip clearance, surface roughness, blunting of the LE and shortening of the blade chord. Schmucker and Schaffler [18] found that even a 1% change in tip clearance can result in a 2% reduction in mass flow and efficiency and a 7.5% decrease in surge margin. Richardson and Smakula [19] examined a JT9D turbofan of extended service and identified that the changes in blade profile, surface roughness and tip wear are the causes of high-pressure compressor performance deterioration. Wulf et al. [20] observed that compressor blades developed a blunt leading edge and increased surface roughness after 6,000–8,000 hours of operation. Balan and Tabakoff [21, 22] investigated erosion damage in a NACA axial compressor and noted changes in the pressure rise coefficient and efficiency due to erosion of the LE and blade tip. Ghenaiet et al. [23] tested a high-speed axial fan stage that ingested sand particles at a concentration of 775mg/m³ during nine hours. They have reported a 7.1% decrease in efficiency and a 9.5% in pressure rise, while the surge margin dropped by 55% due to the blade blunting, chord reduction and increased tip clearance. Hergt et al. [24, 25] provided detailed information on the effects of LE blunting in an eroded transonic fan blade.

The limited data obtained from experimental campaigns highlight the importance of combining them with numerical tools to predict the particle trajectories and erosion in turbomachinery under various operating conditions. Tabakoff and his team [26, 27] were the first to conduct numerical studies on particle trajectory and erosion in axial turbomachines using experimentally determined particle restitution factors. Over the last decades, many studies have integrated 3D flows and physical effects into Eulerian-Lagrangian models. In this context, Ghenaiet carried out numerical studies of erosion in axial and radial configurations of fans [28, 29], compressors [30, 31] and turbines [32, 33]. In the first compression stage of a low bypass mixed turbofan, Ghenaiet [30] simulated particle trajectories and erosion considering a MIL-Spec sand with concentration ranging from 50 to 700mg/m³. As a result, erosion occurred at the LE and end-walls of the IGVs, while in the diffuser erosion spread out across the blade. The rotor erosion was mainly on the PS of the blade, with extreme erosion of the LE and TE added to blade tip distortion. De-Giorgi et al. [34] examined how volcanic ash affected erosion damage in the NASA Rotor 67 and compared the deterioration of total pressure rise between the baseline and eroded blade. Saxena et al. [35] detailed a CFD approach for tracking particles in a multi-stage axial compressor and illustrated radial centrifugation of large particles towards the casing. Vogel et al. [36] used a CFD package to track volcanic ash particles of different sizes ingested into a simplified design engine under three power conditions. Also, they analysed changes in particle size distribution within the core stream due to centrifugation by the fan.

According to the literature, the majority of research focused on compressors, turbines and occasionally industrial fans. However, limited information is available on how volcanic ash particles are ingested into HBTFEs and how the fan prevents larger particles from entering the engine core. Modern HBTFEs equipped with large transonic fans operate at high peripheral speeds and generate high mass flow rates, which portend significant erosion of front components due to abrasive particles carried by airflows.

The purpose of this research is to characterize the intricate flow and aerodynamic performance of a transonic fan stage, as well as to investigate how volcanic ash particles behave as they enter the front components of an HBTFE. Furthermore, the study aims to analyse the development of erosion of the surfaces that come into contact with ash particles, such as the Pitot intake, spinner, fan blade, OGVs discharging to the secondary nozzle, and IGVs leading to the engine core. The computations of particle trajectory and erosion used the validated in-house PARTRAJ code [23] which served in previous studies concerning turbomachinery erosion [28–33].

This paper is organised into eight sections, including an introduction. Section 2 outlines the HBTFE model, front components geometry, and CAD models. Section 3 discusses the computational domain, grid generation, boundary conditions, and solver settings, along with presenting the results of aerodynamic performance and flow structures. Section 4 explores the theoretical aspects of particle trajectory, Lagrangian tracking model including the random walk approach, and probabilistic models of particle release positions, sizes and restitution factors. The same section details the solution procedure, boundary conditions, and particle seeding distributions. Section 5 focuses on erosion modelling and characterisation. Section 6 addresses the code structure and validation. Section 7 presents and discusses the computed ash particles trajectories and describes their behaviour within the front components. Section 8 showcases erosion patterns, and presents in tables and graphs average erosion rates and material removal. Finally, Section 9 summarises the main findings with conclusions.

2.0 Fan stage configuration

Current commercial aircraft are powered by HBTFEs, which rely on a large fan as an essential component to achieve the overall engine performance. This fan is characterised by a low aspect ratio and highly loaded thin blades of transonic velocities. The studied components are from an old CF6-80A engine, with typical cruise flight operating parameters listed in Table 1, adapted from [37]. The rotor consists of 38 highly twisted transonic blades (Fig. 1(a)) made from titanium. A splitter located downstream of the fan rotor directs the majority of the air flow to discharge through the bypass duct via 80 OGVs, while the lower sections of the fan rotor are connected to the core of the engine through 78 IGVs, as illustrated in Fig. 1(b). The determination of the geometry of the fan blade and vanes used an optical scanning machine. The CAD models are shown in Fig. 2(a) for the rotor and in Fig. 2(b) for the complete fan stage. Detailed geometric data are listed in Table 2.

3.0 Flow solution

In dilute particulate flows, the gas phase is assumed to be a continuum phase while the particles are dispersed throughout. The flow field in this large transonic fan of low aspect ratio and highly twisted blades is 3-D, compressible and very complex. It involves flow features such as shock formation, boundary layer transition, separation, shock-boundary layer interactions, secondary flows, tip leakage, and vortices formation, as well as significant flow fluctuations. Moreover, the upper sections are prone to flow separation and often operate close to stall flutter in off-design conditions.

Table 1. Operating parameters

Figure 1. Fan stage components: (a) Fan blade (b) OGVs and IGVs.

Figure 2. CAD models: (a) Fan rotor (b) Fan stage components.

Table 2. Geometry parameters

Figure 3. Computational domain.

Although the flow complexity can be high, CFD simulations often rely on Reynolds-Averaged Navier-Stokes (RANS) approximations, since the Direct Numerical Simulation (DNS), capable of accurately representing all scales of turbulence, requires significant computating power. On the other hand, the Large Eddy Simulation (LES) remains prohibitively expensive and not commonly used in major engineering design problems due to computational resource limitations. Despite their inability to capture large separations, RANS models are widely used in industry due their simplicity and ability to achieve reliable results. The k- $\omega$ based SST (shear stress transport) turbulence model [38] combines the robustness of the standard k- $\omega$ model near walls with the capabilities of the k- $\varepsilon$ model, offering a good balance between grid resolution and flow complexity. Moreover, it can predict boundary layer thickness and separations based on wall function, reducing the need for large mesh size and making it a viable option for many industrial scenarios [39]. This model has been successfully verified in various instances [40], demonstrating good agreement between simulations and experimental results. Validation tests conducted on the NASA transonic Rotor 35 [41] demonstrated a high level of agreement between the predictions of the k- $\omega$ based SST turbulence model and LDV data collected at the design and 80% speeds. Another assessment focusing on ultra-high-lift turbine blades [42] concluded that the k- $\omega$ SST model accurately predicted aerodynamic performance. In the present research, the flow solution is based on the RANS equations with k- $\omega$ based SST as the turbulence model implemented in the CFX solver.

3.1 Computation domain

The front components of this turbofaninclude the Pitot intake, spinner, fan blade, two splitted downstream ducts. The first duct is connected to the engine core via IGVs, while the second one is connected via OGVs to the secondary nozzle, as shown in Fig. 3. Passage sectors are used through the lateral periodicity to reduce computing power. One blade is used for the fan rotor, and two vanes are used for both the IGVs and OGVs. Subsequently the pitch ratio between the rotor and IGVs $\frac{{{s_{rot}}}}{{{s_{IGV}}}} = 1.0263$ and that between the rotor and OGVs $\frac{{{s_{fan}}}}{{{s_{OGV}}}} = 1.0526$ , which are close to unity and justifying the true periodicity. Basically, there are two reference frames: stationary and moving. This leads to using the mixing-plane approach at the interface, where the flow data are averaged circumferentially. Although the mixing-plane approach involves simplifications, the solutions can still yield reasonable results for the time-averaged flow field [43].

The total pressure and temperature values set at the inlet correspond to the standard atmospheric conditions at the cruise altitude of 10,650m and flight Mach number of 0.8, thus resulting in total pressure of 36.335kPa and total temperature of 246.9K. The average static pressure at the outlet ducts was varied to cover the entire speedline along the choke limit. The computations were repeated for different speedlines. Upstream the Pitot intake the turbulence intensity was estimated by $i\left( \% \right) = 100 \times 0.16Re_{{d_h}}^{ - 1/8}$ [44], and thus $i\left( \% \right)$ = 2% based on the Reynolds number $R{e_{{d_h}}} = 1.6 \times {10^7}$ formed for the intake diameter.

3.2 Mesh generation

Hexahedral mesh blocks were distributed to secure fine meshing around the fan blade, spinner, vanes, hub and shroud, in addition to the Pitot intake. Near the walls, 12 meshlines guaranteed a minimum volume distortion and to solve for the boundary layers. In the tip clearance, 15 meshlines allowed to capture the leakage flow. The mesh obtained for the fan blade is shown in Fig. 4. The first layer of nodes near walls corresponds to a distance $\;\Delta y$ targeting ${y^ + } = 2$ . The dimensionless parameter $\;{y^ + } = \Delta y\;\frac{\rho }{\mu }{U_t}\;$ considers the friction law of a flat plate $\;{U_t} = {V_\infty }\sqrt {{c_f}/2}$ , where ${c_f} = 0.027Re_c^{ - 1/7}$ . ${V_\infty }$ is the free velocity outside the boundary layer. The Reynolds number $R{e_c} = \rho {V_\infty }c/\mu$ is based on the mean chord of blade or vane. At the cruise operating conditions the values of $R{e_{c\;}}$ of the rotor blade, IGVs and OGVs are equal to 2.5×10⁶, 5.3×10⁵ and 3.3×10⁵, respectively. The turbulence model k- $\omega$ based SST implements an automatic near wall treatment to switch between the low-Reynolds model for finer mesh and the wall-function for coarser mesh. Figure 5 demonstrates that the values of ${y^ + }$ on the surfaces of the fan blade, hub and spinner are below 4.5, with the highest values observed near the LE and blade tip. In the IGVs the values of ${y^ + }$ do not exceed 3.2, and that of OGVs is below 2.9, with a maximum occurring on the SS from the LE.

Figure 4. Rotor mesh: (a) Hub (b) Mid-span (c) Tip (d) Blade.

Figure 5. Contours of y⁺ on the surfaces of: (a) Fan blade (b) IGVs (c) OGVs.

The grid size independence verification carried out in cruise point shows (Fig. 6) that the fan pressure ratio (FPR), bypass ratio (BPR) and total-to-total isentropic efficiency tend to stabilise beyond the fourth mesh. Additionally, the evaluated hourly eroded mass per fan blade for the particle concentration of 2mg/m³ indicates (Fig. 6) slight variation beyond the fourth mesh size. As a result, separate grids of a total number of 1,611,925 nodes and 1,525,900 elements were generated and adopted in this study.

Figure 6. Grid size independence verification.

3.3 Flow results

The flow solutions based on the RANS equations and the stage interface allowed to compute the transonic fan stage performance maps and characterize the complex flow structures through the front components of this HBTFE. Moreover, the flow field provided the necessary data for calculating the trajectories of ash particles and erosion.

3.3.1 Performance maps

The mass averaging of flow properties at the inlet and outlet of the fan stage allowed calculating $FPR = \frac{{{P_{{t_o}}}}}{{{P_{{t_i}}}}}$ and the total-to-total isentropic efficiency ${\eta_{is\;t - t}} = {c_p}_i\;{T_{ti}}\left( {FP{R^{\frac{{{\rm{\gamma }} - 1}}{{\rm{\gamma }}}}} - 1} \right)/\left( {{h_{to}} - {h_{ti}}} \right)$ . These performance parameters are shown in Figs. 7(a) and 7(b), plotted against the corrected mass flow rate ${\dot m_{co}} = \dot m\sqrt {{T_{{t_i}}}/{T_{st}}} /\left( {{P_{{t_i}}}/{P_{st}}} \right)$ and speed ${\dot N_{co}} = N/\sqrt {{T_{{t_i}}}/{T_{st}}}$ , with ${P_{st\;}}$ and ${T_{st}}$ are the static values at standard sea level. The fan rotational speeds are 70, 80, 90, $\;100$ , $102.7$ and 110 ${\rm{\% }}{N_d}$ . From the computed performance (Fig. 7), the design point is determined for ${N_d} = 3,432{\rm{rpm}}$ at the corrected mass flow rate of 637.5kg/s corresponding to the true mass flow rate of 250kg/s. At that point the computed performance gives an FPR = 1.711 and a total-to-total isentropic efficiency $\eta_{is\;t - t} = 82.4$ 1%. Moreover, the cruise operating point corresponds to the rotational speed $\;N = 3,525{\rm{rpm}}$ and a corrected mass flow rate of 654.7kg/s equivalent to the true mass flow rate of 254.15kg/s. At the cruise point, the predicted values are FPR = 1.733, BPR = 4.287, and $\eta_{is\;t - t} = 79.96\% $ . These obtained values align with the data of the operating manual [37]. This latter states that in the flight cruise conditions, the true air mass flow rate is 257.6kg/s equivalent to a corrected mass flow rate of 663.5kg/s, while the FPR = 1.745 and BPR = 4.24.

Figure 7. Fan stage performance: (a) FPR (b) Total-to-total isentropic efficiency.

Figure 8. Vectors of flow velocity at span: (a) 10% (b) 50% (c) 75% (d) 98%.

Figure 9. Streamlines and blade tip flow structure.

Figure 10. Mach number at different blade spans: (a) 10% (b) 50% (c) 75% (d) 98%.

3.3.2 Flow field characterisation

In order to replicate the main flow patterns, the flow field was solved in the cruise flight conditions, Mach 0.8, altitude of 10,650m, and fan speed of $\;102.7{\rm{\% }}{N_d} = 3,525{\rm{rpm}}$ . This large fan is characterised by transonic flow velocities for more than 60% of the blade span. Figure 8 displays the flow velocity vectors at different locations, such as 10% span, mid-span, 75% span and near the tip of blade. Upstream and around the fan blade LE, the relative flow velocities are extremely high and increasing with the blade height. High flow velocities, reaching up to 465.9m/s, are observed on the SS of the fan blade, leading to the formation of shock structure. Towards the blade TE, the lower velocities are due to diffusion and wake formation. At the interface plane between the rotor and IGVs/OGVs, the flow exhibits drastic changes in flow direction, and the blade’s wake disturbs the angle-of-attack around the LE of the stationary vanes. This is clearer at the LE of IGVs since are closer to the fan blade. Across IGVs the flow accelerates, reaching high velocities on the front part of SS. In contrast, the OGVs, located farther, dissipate more the wake issued from the fan blade, hence the flow is well-aligned with the vanes’s LE.

The most intricate flow patterns are observed above the blade tip and through the tip gap (Fig. 9). In this typical region, the leakage flow is clearly depicted by the vectors of flow velocity crossing from the PS to SS, with their directions varying from the LE to TE of the blade. Figure 9 depicts the formation of the leading edge vortex (LEV) from the LE, and extending along the blade tip. The LEV rotates opposite to the tip leakage flow vortex (TLFV) and mixes out with the annulus boundary layer. At the exit from the rotor blade, a swirl velocity is established to ensure a radial equilibrium balanced by the axial velocity profile.

In this transonic fan the pressure loading is most pronounced at the front of the blade, and the diffusion process is greatly influenced by the radial position. The peak of static pressure occurs at the blade LE, while a noticeable depression is visible on the SS starting from mid-span onward. The formation of shock waves is responsible for the main increase in static pressure across the fan blade passage, with the highest aerodynamic loading at the upper sections from the mid-span.

The distributions of the relative Mach number are displayed in Fig. 10 at 10% span, mid-span, 75% span and near blade tip. As the flow moves from the inlet of the Pitot intake towards the fan rotor, the convergence imposed by the conical shape of the spinner caused an increase in Mach number before reaching the fan blade. Near the hub, the relative Mach number is practically subsonic and the flow velocity decreases along the blade passage. A supersonic bubble appears (Fig. 10(a)) at the SS close to the blade LE and extends towards 25% of the span. It is perceived that the transonic flow velocities dominate over more than one-third of the SS, while the subsonic velocities prevail on the remaining of the fan blade. As depicted in Fig. 10(b), a passage shock region develops from the LE over approximately the third of the SS and may reach the opposite fan blade. Downstream of the shock, the Mach number decreases until reaching a value around 0.6 at the exit of the fan blade. The distinctive characteristics of this transonic fan blade include an oblique shock sitting near the LE (Fig. 10(c)) and crossing the blade passage at 50% of the span, which interacts with the neighbouring blade. In the upper sections where the flow is entirely supersonic over the SS, a passage normal shock sits near the TE of the blade (Fig. 10(d)), with a sonic layer extending up to a third or two-thirds of the blade chord on the opposite PS. The strength of this shock is indicated by a jump in the Mach number from a value of 1.43 to 0.83. Near the TE of the fan blade, the boundary layer of the SS becomes considerably thicker compared to that from the PS, and eventually merges to form a thick wake diffusing downstream.

4.0 Particle trajectory

There are two common approaches in particle-laden turbulent flows simulations, which differ in the treatment of the dispersed phase. The Eulerian approach considers both the continuous and dispersed phases as interpenetrating and solves the time-averaged transport equations to predict the field distributions of momentum and volume fraction for each phase [45]. In the present computations, as the volume fraction at the highest concentration is about 1.5×10⁻⁹ below the limit of 4×10⁻⁷ [46], the dispersed phase is assumed to be dilute, and the Lagrangian particle tracking model has to be considered. This model is used to track large numbers of particles using a stochastic approach to achieve a significant solution [47]. Compared to the Eulerian approach, the Lagrangian approach is more cost-effective [48] and better at handling particle-wall interactions [49]. After obtaining the flow field, the flow data are exported to the trajectory code to compute the particle trajectory and impacts in step by step.

4.1 Particle motion modelling

Particle motion is governed by the balance between the rate of momentum change and the superposition of external forces. The left side of the vector equation of particle motion (1) takes into account inertia, centrifugal force, and Coriolis force, while the right side sums the external forces reduced by the mass of particle, including the drag force $\overrightarrow {\;{f_D}} $ , combined gravity and buoyancy force $\overrightarrow {\;{f_{GB}}} \;$ , pressure gradient force $\overrightarrow {{f_P}} $ , and Saffman force $\overrightarrow {{f_S}} $ .

(1) \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_P}} + \overrightarrow {{f_S}} \end{align}

The dominant external force is the drag force, expressed in terms of the drag coefficient and Reynolds number $R{e_p} = \frac{{{\rho _f}}}{{{\mu _f}}}{d_p}\left| {\overrightarrow {{V_f}} - \overrightarrow {{V_p}} } \right|$ , which when reduced by particle mass becomes:

(2a) \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}

At very low $R{e_p} \lt 0.1,$ the drag coefficient corresponds to the Stokes regime $\;{C_D} = 24/R{e_p}$ . For the wider range 0.01 $ \le {\rm{R}}{{\rm{e}}_{\rm{p}}} \le $ 2.6×10⁵, the four-parameter drag coefficient [50] is considered as follows:

(2b) \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 revaluated by Haider and Levenspiel [50] based on the shape factor.

(2c) \begin{align}a &= \exp\left( {2.3288 - 6.4581\phi + {{2.4486}}\phi}^2 \right)\nonumber\\[5pt]b &= 0.0964 + 0.5565\phi\nonumber\\[5pt] c &= \exp\left( {4.905 - 13.8944\phi + {{18.4222}}{\phi}^2 - {{10.2599}}\phi}^3 \right)\nonumber \\[5pt] d &= \exp\left( {1.4681 + 12.2584\phi - {{20.7322}\phi}^2 + {{15.8855}}\phi}^3 \right)\end{align}

Since ash particles are of irregular shapes, the constants $a$ , $b$ , $c$ and $d$ are based on the shape factor $\phi = {A_S}/{A_P}$ ( ${A_S}$ area of a sphere of equivalent volume, and $\;{A_P}\;$ actual area of the particle). De-Giorgi et al. [51] by analysing the Etna volcano ash particles with a scanning electron microscopy found that the shape factor is equal to 0.68, in the same range of 0.6–0.8 recommended by Riley et al. [52] based on three samples of volcanic ash. In this work the value of 0.7 is opted for.

The combined gravitational force and buoyancy force per unit of particle mass is as follows:

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

The pressure gradient force relates the effect of the local fluid pressure gradient around a particle, which when reduced by particle mass becomes:

(4) \begin{align}\overrightarrow {{f_P}} = - \frac{1}{{{\rho _p}}}\vec \nabla p\end{align}

The slip shear layers induced force is due to non-uniform relative flow velocities. Saffman [53] established the first expression, which when developed in 3-D and reduced by particle mass becomes:

(5a) \begin{align}\overrightarrow {{f_{SL}}} = \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 lift coefficient [46] given by:

(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 [54] for particle Reynolds number in the range of $0.1 \le R{e_p} \le 100$ and $\beta = 0.5\frac{{\;R{e_s}}}{{\;R{e_p}}}$ .

(5c) \begin{align}\begin{array}{l@{\quad}l}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\\[5pt]f\left( {\;R{e_p},\;R{e_s}} \right) = \;0.0524\sqrt {\beta R{e_p}} & R{e_p} \gt 40\end{array}\end{align}

The shear flow Reynolds number ( $R{e_s}$ ) is given by:

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

4.2 Turbulence effect

The stochastic random walk implemented within the Lagrangian particle tracking model is based on the eddy interaction model of Gosman and Ioannides [55]. This latter assumes that a particle and an eddy interact as long as the interaction time $\Delta {t_{in}} = {\rm{min}}\left( {\Delta {t_e},\,\,\Delta {t_r}} \right)$ [56] is shorter than the eddy lifetime $\Delta {t_e} = 0.37\frac{k}{\epsilon }\;$ and the particle displacement is smaller than the eddy’s length ${l_e} = 0.3{{{k^{3/2}}} \over \epsilon }$ (Shirolkar et al. [49]). The transit time $\Delta {t_r}$ is the duration taken by a particle to cross through an eddy, which has a linearised form [55] function of ${l_e}$ and relaxation time ${\tau _p}$ as follows:

(6) \begin{align}\Delta {t_r} = - {\tau _p}Ln\left[ {1 - {{{l_e}} \over {{\tau _p}\|\overrightarrow {{V_f}} - \overrightarrow {{V_p}}\| }}} \right]\;\,{\rm{and}}\,\,{\tau _p} = {4 \over 3}{{{\rho _p}} \over {{\rho _f}}}{{{d_p}} \over {{C_{D \| \overrightarrow {{V_f}} - \overrightarrow {{V_p}} \| }}}}\end{align}

When a particle is trapped inside an eddy the turbulent components $u',v',w' = \xi\sqrt {\frac{2}{3}k} $ are added to the mean velocity components of the fluid through a Gaussian random number $\xi $ . The turbulence effect is updated when the particle either passes into a new mesh cell or exceeds the eddy lifetime which is a function of the specific dissipation rate.

Figure 11. Impact conditions.

4.3 Boundary conditions

The true impact on a wall is determined by half the diameter of a particle (Fig. 11) to avoid non-physical trajectories and impacts. Interpolation is used to determine a more precise time step, and the calculation is repeated. At an impact point, the local surface unit vectors $\vec n\;$ and $\vec t\;$ (Fig. 11) are determined and used to calculate the normal and tangential components of impact velocity. Tabakoff et al. [57] first studied the post-impact properties of solid particles using high-speed photography to analyse the rebound correlations for annealed 2024 aluminium alloy. Then followed the work of Grant et al. [58] carried out in an erosion wind tunnel. By using the transit laser anemometer, Tabakoff et al. [59] determined the restitution factors for typical turbomachinery materials such as 2024 Aluminium, TiA1-4, AM355 and Rene 41. The particle impact is influenced by its hardness, irregular shape and surface roughness, leading to the rebound characteristics which are dispersed around the averaged values, considering the values of standard deviation [59]:

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

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

At the symmetry boundary condition, such as the axis of rotation, the velocity components are updated by pure reflection. When it comes to the periodic boundary condition, the particle coordinates and velocity components are transferred to the opposite side using the pitch angle. A FEM-based search algorithm is used to determine the natural coordinates and update for the cell number. At the interface between the stationary domain and rotating domain, the coordinates of a particle are transferred locally to the next frame and the FEM algorithm searches for the next cell number. The particle velocity is updated using the velocity composition.

(8) \begin{align}{\overrightarrow{V_{p_{{}_{abs}}}}}= {\overrightarrow {V_{p_{{}_{rel}}}}} + \overrightarrow {\Omega} \times \overrightarrow {r\;} \end{align}

4.4 Solution procedure

The Runge-Kutta-Fehlberg seventh-order technique is used to integrate the set of non-linear differential Equations (1). Initial particle slip velocities are set based on the available experimental data [23]. The integration time step is determined from the computational cell and the local flow properties, while keeping the leading truncation error within a specified tolerance. If a particle interacts with an eddy, the integration time step is adjusted to $\;\Delta {t_{step}} = {\rm{min}}\left( {\Delta {t_{cell}},\;\Delta {t_{in}}} \right)$ . The particle tracking algorithm is using the FEM to interpolate for the local flow properties and to find the exact impact locations and conditions. To locate a particle in a mesh cell, the physical coordinates are transformed into the local coordinates using a non-linear system of equations solved via the Newton-Raphson method. If a local value exceeds unity, the particle is considered outside and the cell number is updated. Upon updating for the new cell, the coefficients of the Jacobian matrix are modified to interpolate for the flow properties at the next step of computation. The trajectory code can handle particle release positions that follow a prescribed concentration profile of particles randomly distributed within their size distribution, or grouped based on similar physical characteristics. Injection modes, such as point, line or global injection could be selected to simulate particle ingestion. However, the global injection seems more appropriate in the present study of turbomachinery components erosion.

4.5 Particle seeding

Particle trajectory and erosion calculations were performed considering four samples of volcanic ash collected from various volcanic regions around the globe, and close to the eruption source. The first sample was collected from the Eyjafjallajokull volcano in Iceland (North Europe, 2010), and particle sizes determined through sieving (Fig. 12(a)) to mimic steady-state flight conditions [60]. The second sample (Fig. 12(b) was collected from the Chaiten volcano in Chile (South-America, 2008) [61]. The third sample (Fig. 12(c)) was collected from the Etna volcano in Italy (South Mediterranean, 2002) [51]. Lastly, the fourth sample (Fig. 12(d)) corresponds to the Kelud volcano in Indonesia (South Asia, 2014) [62]. The cumulative frequency curves of particle sizes of the four samples are shown in Fig. 13. The mean size and standard deviation (Table 3) of each ash sample are calculated from the size density derived from the cumulative frequency curve. The density of volcanic ash (Table 3) is influenced by the quartz concentration, as reported by Lau et al. [63]. The main component of volcanic ash is quartz, ranging from 47 to 73.2% [63]. The Chaiten sample, from a rhyolitic eruption, has the highest SiO₂ percentage of 73.9%, while the samples from Eyjafjallajokull, Kelud, and Etna have percentages of 57.8, 54.07, and 47.1% respectively, which are mainly consisting of basaltic material.

Figure 12. SEM images of different volcanic ash samples: (a) first sample from Eyjafjallajokull [60] (b) second sample from Chaiten [61] (c) third sample from Etna [51] (d) fourth sample from Kelud [62].

Figure 13. Cumulative frequency of ash particle sizes.

The maximum concentration in the vicinity of an aircraft is at least equal to 4mg/m³ [9, 11], although the level of confidence in the estimated concentrations is low and uncertainties of an order of magnitude are possible [64]. In this study, the number of ash particles and seeding positions, radially and tangentially, in the Pitot intake are randomly set according to the concentration of particles, size distribution (Fig. 13), intake geometry and flow conditions. The mass rate of volcanic ash particles is obtained by multiplying the air volume flow rate and particle concentration of successively 1, 2 and 4mg/m³. A subprogram continuously iterates on the particle sizes and release positions until the convergence of the total particle mass rate.

5.0 Erosion model

Early studies on erosion have concluded that this phenomenon is affected by various factors, such as the characteristics of erodent particles, size, concentration, impingement velocity and angle, and the properties of the target material. Finnie [65] established the basis for the sand erosion, which was later improved by Bitter [66] to account for the effect of plastic deformation. Grant and Tabakoff [67] developed successful erosion models for aerospace materials, beginning with aluminium alloys and quartz particles. Their correlation of erosion rate incorporates two mechanisms predominant at low impact angle and high impact angle, and expressed as material removed (milligram) per mass of impacting particles (in grams). Tabakoff et al. [68] developed further erosion models for aluminium, stainless steel, and titanium impacted by ash particles, which are expressed as follows:

(9) \begin{align}\epsilon = {k_1}{\left[ {1 + ck.{k_{12}}{\rm{sin}}\left( {90/{\beta _0}} \right){\beta _1}} \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 } = 1 - 0.0016{V_{p1}}\sin {\beta _1}$ , and the constant $ck = \left\{ {\begin{array}{*{20}{c}}{1\;\;if\;\;{\beta _1} \le 3{\beta _o}}\\{0\;\;\;if\;\;\;{\beta _1} \gt 3{\beta _o}}\end{array}} \right.$

Table 3. Volcanic ash particle size and density

The Pitot intake and spinner are made from aluminium alloy, whereas the fan blades, IGVs and OGVs are made from titanium. For the 2024 aluminium alloy and ash as erodent particles, the angle of maximum erosion $\;\;{\beta _0} = 20\;{\rm{deg}}\;$ and the material constants are: $\;{k_1} = 0.156988 \times {10^{ - 5}}$ , ${k_{12}} = $ 0.3193, ${k_3} = 2 \times {10^{ - 12}}$ . For the titanium alloy 6Al-4V and ash as erodent particles, the angle of maximum erosion $\;\;{\beta _0} = 20\;{\rm{deg}}\;$ and the material constants are: $\;{k_1} = 0.1564951 \times {10^{ - 5}}$ , ${k_{12}} = $ 0.173636, and ${k_3} = 3 \times {10^{ - 12}}$ .

At each impact point, the eroded mass is calculated using the erosion rate evaluated from Equation (9). An area around the concerned node, denoted centred area $\;{A_{node}}$ , is formed by taking a quarter of area of the four faces that share this node. The values of eroded mass per unit of time (second) are summed and divided by the area $\;{A_{node}}$ to calculate the equivalent erosion rate $E$ (mg/s.mm²) according to Equation (10) and distributed among the nodes of impacted faces.

(10) \begin{align}E = \frac{{\mathop \sum \nolimits_1^{ni} {\epsilon _i}{m_{pi}}}}{{{A_{node}}}}\end{align}

6.0 Code structure and validation

The organisation and functionalities of the FORTRAN code PARTRAJ, depicted in Fig. 14, are summarised as follows: COMMON BLOC contains all the arrays and matrices allocations, and the declarations and common variables. INPART reads the control and initial parameters. The geometry, grids, and flow data are imported from the flow solver. PGEODATA reads the geometry and topology, while FLODATA reads the flow data. ARDIST distributes particles at the release plane according to the selected injection mode, size distribution, and concentration. BOUNDARY sets all boundary conditions. INTVEL interpolates for the flow data. EXFORCE calculates all activated external forces. RKUTTA integrates the particle motion equations. JACOBI computes the shape function and the Jacobian matrix. TRANSF converts the global coordinates to the local coordinates. LOCATE searches for particle position and updates for the cell number. REBOUND calculates the rebound properties. STATISTIC generates random numbers and calculates the mean and variance for the given distributions. SPLINT interpolates functions. EROSION calculates erosion rate and cumulates eroded mass using CUMEROS , and calls DISTEROS to calculate erosion rate density and its distribution. DEGGEOM calculates depths of penetration and deteriorated geometry. OUTPUT generates files for the plots of trajectories, erosion and deteriorated geometry.

Figure 14. PARTRAJ flowchart.

This code has been validated based on a high-speed axial fan stage of 11,300rpm, made of aluminium and used for aircraft equipment ventilation, as detailed in the author’s thesis [69]. Silica sand MIL-E 5007E of particle sizes ranging from 0 to 1,000 microns was used to test and simulate particle trajectories and erosion. Figure 15(a) shows sand particle paths deviating from the streamlines due to inertia and radial spread of particles. Large particles follow ballistic trajectories after hitting walls, while the small ones closely follow the streamlines. In IGVs, the impacts mainly occur at the aft upper corner of the PS, while the remaining particles traverse completely without collisions. Figure 15(b) shows frequent impacts around the blade’s LE caused by sand particles at high velocity and incidence. High frequency of impacts result in intense erosion spreading in a triangular region on the PS, and over the tip of the blade. Particles migrating from the blade’s PS to the SS, through the tip clearance, cause clear erosion of the tip. Figures 15(c) and 15(d) illustrate the computed erosion rate densities on the PS and SS of the blade, which are consistent with the paint stripping test results. Figures 16(a) and 16(b) reveal the measured performance degradation based on the tested eroded blade profiles compared to the predicted performance degradation based on the computed eroded blade profiles. Initially, at the design point defined for the mass flow rate of 0.8kg/s, the fan stage had an efficiency of 0.802 and a pressure rise coefficient of 0.395. After nine hours of sand ingestion, the measured performance droped to the successive values of 0.748 and 0.381, while the predicted performance declined to the values of 0.740 and 0.360, respectively.

Figure 15. Sample of trajectories (a) sample of impacts (b) computed erosion rate density on the blade PS (c) and SS (d) tested erosion patterns showing the PS (e) and SS (f).

Figure 16. Compared performance degradation: (a) Pressure rise coefficient (b) Efficiency.

Figure 17. Samples of trajectories of particle size 1μm coloured by velocity: (a) Meridional view (b) Top view of fan blade, OGVs and IGVs.

7.0 Particle trajectories results

The results of the particle trajectory calculations include information about the particles coordinates, velocities, collision conditions, and impact frequency. The behaviour of ash particles trajectories can be seen in different components, including the Pitot intake, spinner, fan blade, splitter, IGVs, OGVs and outlet ducts. It is evident that most of released particles do not reach beyond 80% of the fan blade span due to the shape of the Pitot intake lip and inner contour. Many ash particles are observed to deflect downward as they collide with the intake lip and strike the spinner, rotor’s hub, and lower sections of the fan blade. The paths of small ash particles closely align with the streamlines, as revealed by the trajectories of particles circulating around the intake lip. Small particles released centrally accelerate along the intake duct and impact the spinner while others move away alongside the spinner, and by entering the rotor they become strongly affected by the centrifugal force. On the other hand, large ash particles of higher inertia than the drag, acquire an outward radial velocity after striking the spinner, while others after collisions bounce off and curve upwards to collide with the upper sections of the fan blade and shroud. Generally, the deviation of particle paths from the streamlines increases with the particle size. Heavy ash particles are more likely to strike the blade surfaces and outer walls at high velocities, bouncing, and bending trajectories.

Figure 17(a) displays the paths of 1-micron ash particles, with their velocities indicated by colours. These particles generally follow the flow streamlines and are influenced by the turbulence and secondary flows. These tiny particles circulate and collide repeatedly around the Pitot intake lip and travel alongside the inner contour. Evidence from Fig. 17(a) reveals that these particles adhere to the streamlines and closely follow the flow details, to reach velocities up to 466m/s. As observed from the top view (Fig. 17(b)), small ash particles travelling around the blade are entrained by the leakage flow over the blade tip, thereby causing multiple impacts and erosion. Figure 18(a) depicts the trajectories of 10 microns ash particles. Many of these particles are seen to circulate around the intake lip a short distance before inducing multiple collisions and being diverted slightly. After recovering kinetic energy from the flow, they move directly towards the fan blade. The same particles when released near the center of intake, by reaching the front of the fan they circulate away from the spinner and only induce a few collisions. Owing to the high peripheral blade speed and flow velocity combined with the blade twist, many high-velocity particles impact the blade’s LE, spanning from the root upward, as depicted in Fig. 18(b). Added to that, there are multiple impacts spreading over the PS of the fan blade. Figure 18(b) shows particles reaching the top of the blade at high velocities up to 457m/s, which when they hit the blade PS at low-impact angles become practically aligned with the streamlines. Other particles at the inner sections of the fan blade, are seen to impact and bounce off at high rebound angles but at relatively lower velocities. Due to the shape of the intake lip and inner contour pushing ash particles away, the majority of them reach the fan blade along approximately 80% of the span, while the rest exhibits scattered impacts. Upon leaving the rotor blade, the particles change their directions and pass to the IGVs and OGVs. The tiny ash particles such as 1 micron are observed to travel through the stationary vanes and follow closely the streamlines while exchanging kinetic energy from the fluid. The large particles of 10 microns seem to deviate notably (Fig. 18(b)) from the streamlines and impact more frequently the PS of the stationary vanes. As also observed, the small ash particles issued from the lower sections of the fan blade are redirected by the secondary flow to reach the LE of IGVs at high incidence and travel farther while inducing several impacts.

Figure 18. Samples of trajectories of particle size 10μm coloured by velocity: (a) Meridional view (b) Top view of fan blade, OGVs and IGVs.

Figure 19. Samples of trajectories of particle size 100μm coloured by velocity: (a) Meridional view (b) Top view of the fan blade, OGVs and IGVs.

Figure 19 illustrates the trajectories of large ash particles (i.e. 100 microns), which are shown to collide repeatedly around the intake lip and diverte downward in sort of ballistic trajectories and impact the lower sections of the fan blade. As shown in Fig. 19(a), the particles released centrally travel along the intake, strike the spinner, bounce back, and then turn to follow paths in a parabolic manner. Moreover, owing to the influence of flow structure and imparted centrifugation, large particles are forced to move upwards and reach the upper parts of the fan blade. The high flow velocities in this transonic fan, leading to high drag force, combined with the significant centrifugal force, cause the large particles to impact the LE of the fan blade from both sides and rebound to reach the opposite PS at the aft of the surface. As observed in Fig. 19(b), the particle velocity increases significantly while moving across the fan blade, reaching up to 434m/s. These large particles are seen to impact the PS of the blade near the optimum impact angle and rebound off from the surface, and after multiple collisions they exit the fan blade as coalesced particles. By crossing the rotor interface, their directions change significantly. Through IGVs or OGVs, they seem to travel almost straight, as depicted in Fig. 19(b), and impact the stationary components at high velocities and angles, primarily at the LE from the SS and more frequently on the PS of the vanes. It should be noted that the majority of ash particles, after undergoing significant centrifugation by the fan blade, experience higher frequency of impacts on the OGVs' surfaces and the shroud compared to the small particles. Furthermore, owing to intense centrifugation imparted by the fan blade, a large amount of ash particles bypass outward, preventing them from traversing the IGVs and entering the engine core.

To demonstrate the behaviours of ash particles of various sizes randomly released, trajectory simulations were performed for the ash samples Eyjafjallajokull (0.4–80 microns) and Kelud (0.5–300 microns), which are characterized by relatively small and large particle sizes, respectively. These particles, randomly injected upstream of the Pitot intake combine all the discussed aspects associated with small, medium and large ash particles. Their paths along the Pitot intake and through the fan blade, OGVs and IGVs, and outlet ducts are depicted in Fig. 20, and coloured by the particle diameter. As observed in Fig. 20(a), the small particles belonging to the first sample tend to bypass the spinner without collisions and move away from the rotor’s hub, and then deviate upward to reach the upper sections of the fan blade. On the other hand, the large particles released centrally are shown to strike the spinner, rebound, and deflect upward to reach the upper sections of the fan blade. Even particles moving near the lower sections of fan blade are redirected upwards and projected against the shroud, due to the high centrifugation. The trajectories of particles belonging to the fourth ash sample, of medium and large sizes (Fig. 20(b)), appear more diverted while crossing the intake. After multiple collisions with the intake lip, these particles follow ballistic paths towards the lower sections of the fan blade. When ash particles collide with the spinner, they rebound backward and then curve upward in parabolic trajectories, striking the upper regions of the fan blade and shroud. For the same reasons, larger particles (in red colour) repeatedly collide with the shroud and are more likely to infiltrate through the tip gap. Upon leaving the fan blade, ash particles enter the OGVs/IGVs at high incidence angles, while others deflect upward and move alongside the shroud. Heavy ash particles are predominantly diverted upward by the fan rotor, pass through the OGVs and are finally expelled through the bypass nozzle. On the other hand, the light particles are less diverted and may reach the engine core. Consequently, this large fan operates as a separator for the large-size particles.

Figure 20. Trajectories coloured by particle diameter: (a) Sample 1 Eyjafjallajokull (b) Sample 4 Kelud.

8.0 Erosion results and discussions

The local erosion rates and patterns are obtained in cruise flight conditions, considering ash particle concentrations 1–4 ${\rm{mg}}/{{\rm{m}}^3}$ and the four samples collected from the major volcanic regions worldwide. Due to high number of impacts on each face of a mesh element, the local erosion rates $\left( {{\rm{mg}}/{\rm{g}}} \right)$ determined from Equation (9) were used to evaluate the erosion rate density ( ${\rm{mg}}/{\rm{s}}.{\rm{m}}{{\rm{m}}^2})$ based on Equation (10). To facilitate the discussion, erosion patterns of the four ash samples are presented only for the concentration of 2mg/m³. Contours of erosion rate density are displayed with the same scales to compare between patterns of different samples. As noticed, higher erosion rates are obtained with the Kelud’s ash sample, followed by that of the Etna volcano.

Figure 21. Erosion rate density on the fan blade caused by volcanic ash, at concentration of 2mg/m³: (a) Sample 1 Eyjafjallajokull (b) Sample 2 Chaiten (c) Sample 3 Etna (d) Sample 4 Kelud.

Figure 22. Erosion rate density on the fan hub and shroud caused by volcanic ash, at concentration of 2mg/m³: (a) Sample 1 Eyjafjallajokull (b) Sample 2 Chaiten (c) Sample 3 Etna (d) Sample 4 Kelud.

Figure 23. Erosion rate density in OGVs and IGVs caused by volcanic ash, at concentration of 2mg/m³: (a) Sample 1 Eyjafjallajokull (b) Sample 2 Chaiten (c) Sample 3 Etna (d) Sample 4 Kelud.

Figure 24. Erosion rate density in Pitot intake, spinner and fan blades caused by volcanic ash, at concentration of 2mg/m³: (a) Sample 1 Eyjafjallajokull (b) Sample 2 Chaiten (c) Sample 3 Etna (d) Sample 4 Kelud.

The patterns of erosion, depicted in Figs. 21–24, are globally similar, except the slight variations in erosion rates, spread, and some specific details between the four ash samples. The first and second samples from the Eyjafjallajokull and Chaiten volcanoes caused dense erosion spreading more on the fan blade’s PS, but with lesser erosion intensity on the SS. On the other hand, the third and fourth samples from the Etna and Kelud volcanoes resulted in scattered erosion on both sides of the fan blade due to dispersed high-velocity impacts, leading to higher local erosion rates compared to the first and second ash samples. As depicted in Fig. 21(a)–(d), there are noticeable areas of erosion spreading on the PS of the fan blade, typically from the hub to 80% of the span, where the frequency of impacts and velocities are very high. The remaining of the top sections of the fan blade from the PS are characterized by fewer impacts and subsequently lower erosion rates. The maximum of local erosion rates on the PS of the fan blade are equal to 1.352×10⁻⁵, 1.974×10⁻⁵, 1.106×10⁻⁴ and 1.791×10⁻⁴ ${\rm{mg}}/{\rm{s}}.{\rm{m}}{{\rm{m}}^2}$ for the first, second, third and fourth ash samples, respectively. It is evident, that the fourth and third ash samples caused the highest erosion densities, due to higher proportion of large size particles. The LE and front of the fan blade were directly impacted by high-velocity particles, resulting in a band of extreme erosion at the LE that extends to both sides. Figures. 21(a)–(d) exhibit pronounced erosion extending upward from the hub and covering a large portion of the PS until the TE, where the third and fourth samples led to extensive erosion rates. The scattered areas of dense erosion observed at the top sections of the fan blade are caused by the centrifuged heavy particles added to the high-velocity particles reflected upward by the spinner. Moreover, the SS of the fan blade reveals impacts spreading from the LE up to two-thirds of the surface and resulting in scattered erosion. The front of the blade’s SS presents concentrated erosion, but of lower rates compared to the PS. As observed, the spread of erosion and rates are the highest with the third and fourth ash samples, characterised by a broader range of particle sizes, higher density and quartz content. According to Figs. 22(a)–(b), the rotor’s hub is severely eroded from its entrance alongside the PS whereas the opposite side has scattered erosion due to irregular impacts and particles deflected away around the blade root. Furthermore, the rotor’s shroud presents scattered erosion due to irregular impacts caused by the centrifuged and reflected heavy particles from the spinner and blade root sections, and subsequently, the cumulated erosion appears more scattered compared to the hub.

Erosion of the IGVs, as illustrated in Figs. 23(a)–(d), extends across both the PS and SS of the vanes, and are characterized by erosion rate densities reaching up to 5.108×10⁻⁵, 9.601×10⁻⁵, 1.426×10⁻⁴ and 1.351×10⁻⁴ ${\rm{mg}}/{\rm{s}}.{\rm{m}}{{\rm{m}}^2}$ , between the first, second, third and fourth ash samples, respectively. However, the third and fourth ash samples caused higher erosion rates at the aft of the vane’s PS until the TE, while the SS is eroded across the entire surface until the TE, but at lower rates. The IGVs’ hub is eroded slightly since ash particles are diverted upwards due to preceding centrifugation while the shroud is more eroded near the splitter. The front junction between IGVs and shroud appears the region more prone to erosion wear. On the other side, the OGVs display patterns of erosion spreading across the entirety of the PS, with the highest rates at the rear of the PS, while the lowest rates are observed near the root. Moreover, the main eroded areas of the SS are concentrated on the LE and the front part of OGVs, while the upper sections seem unaffected by erosion and the same for the aft of this surface. The OGVs’ shroud displays erosion patterns, particularly, at the front part, generated by previously centrifuged particles and those following ballistic trajectories. The hub, on the other hand, is eroded predominantly from the splitter entrance, and onward the erosion becomes weak.

Figures 24(a–d) display the front and back views of erosion patterns formed in the Pitot intake, spinner, fan blades, hub and shroud. The highest erosion rates in the Pitot intake is observed on the intake’s lip, unlike the inner contour of much lower erosion. The Pitot intake experiences erosion rate density of up to 2.811×10⁻⁸, 2.805×10⁻⁸, 7.469×10⁻⁸, and 7.865×10⁻⁸ ${\rm{mg}}/{\rm{s}}.{\rm{m}}{{\rm{m}}^2}$ , for the first, second, third and fourth ash samples, respectively, but of lower rates compared to the spinner and the fan blades. The spinner displays erosion spreading more near the fan hub due to the high frequency of impacts, and subsequently, the erosion rate densities are high and equal to 1.248×10⁻⁶, 1.191×10⁻⁶, 3.291×10⁻⁶ and 9.005×10⁻⁶ ${\rm{mg}}/{\rm{s}}.{\rm{m}}{{\rm{m}}^2}$ between the first, second, third, and fourth ash samples, respectively. As noticed, the fourth ash sample caused more erosion of the spinner. The back view shows that the fan blades suffer more erosion mainly from the blades’ PS, with extreme erosion seen at the rear part and the TE. Accordingly, the fourth ash sample induced the highest erosion rate density of up to 1.791×10⁻⁴ ${\rm{mg}}/{\rm{s}}.{\rm{m}}{{\rm{m}}^{2{\rm{\;}}}}$ . On the other side, the front view of the components displays erosion extending across the intake lip. Also, reveals scattered erosion on the blades’ SS, with the highest rates observed near the LE. One may conclude that the most of erosion scatters, of higher rates, are induced by the fourth ash sample.

The average erosion rate density (AERD), in milligrams per second and square millimetre (mg/s.mm²), is another metric for assessing the overall erosion and characterising how erosion develops in each component. This parameter is determined using the local erosion rates, mass of impacting particles, and component’s area. To evaluate the AERD, the erosion rate density $E$ based on Equation (10) is integrated over the whole surface of the concerned component and then divided by its area. In addition, the hourly eroded mass (HEM), expressed in milligrams per hour (mg/hr), provides a comprehensive assessment of the material removal from a component. The HEM is calculated by multiplying the value of AERD by the component’s area, for duration of 3,600 seconds.

Table 4. Erosion parameters caused by four volcanic ash samples

Table 4 summarises the obtained results in terms of ash particle concentration, particle mass rate, AERD, and HEM. At the highest ash particle concentration, the fan blades’ AERD, reached the maximum values of 0.161×10⁻⁶, 0.162×10⁻⁶, 0.266×10⁻⁶ and 0.276×10⁻⁶ ${\rm{mg}}/{\rm{s}}.{\rm{m}}{{\rm{m}}^2}$ between the first, second, third and fourth ash samples, respectively. It is evident that the fourth and third samples led to the highest values of AERD. Moreover, the HEM (mg/hr), for the totality of 38 blades, is equal to 5.871×10³, 5.910×10³, 9.705×10³, and 10.067×10³mg/hr, where the highest values occur for the third and fourth ash samples. The fan’s hub exhibits AERD values of 0.175×10⁻⁶, 0.178×10⁻⁶, 0.335×10⁻⁶, and 0.352×10⁻⁶ ${\rm{mg}}/{\rm{s}}.{\rm{m}}{{\rm{m}}^2}$ , for the first, second, third and fourth ash samples, respectively. Accordingly the HEM values are 0.514×10³, 0.523×10³, 0.983×10³ and 1.033×10³mg/hr, where the highest values are caused by the third and fourth samples. The shroud reveals AERD values of 0.064×10⁻⁶, 0.583×10⁻⁶, 0.525×10⁻⁶, and 1.430×10⁻⁶ ${\rm{mg}}/{\rm{s}}.{\rm{m}}{{\rm{m}}^2}$ , between the first, second, third and fourth ash samples, respectively, leading to HEM of 0.4952×10³, 4.487×10³, 4.041×10³ and 11.013×10³mg/hr, with the highest values corresponding to the fourth ash sample. Overall, the fourth ash sample induced the most severe erosion of the fan blades, followed by the third sample. The rotor’s shroud, despite its lower erosion rates, has a higher HEM than the hub because of larger surface.

Figure 25. HEM (g/hr) of the complete fan rotor: (a) blades (b) hub (c) shroud.

Figure 26. HEM (g/hr) of complete: (a) IGVs (b) OGVs.

Figure 27. HEM (g/hr) of complete: (a) spinner (b) Pitot intake.

The entire IGVs are characterized by AERD values equal to 0.404×10⁻⁶, 0.386×10⁻⁶, 0.506×10⁻⁶, and 0.514×10⁻⁶ ${\rm{mg}}/{\rm{s}}.{\rm{m}}{{\rm{m}}^2}$ , respectively, and the subsequent HEM are equal to 2.337×10³, 2.233×10³, 2.923×10³, and 2.968×10³mg/hr, respectively. Moreover, the IGVs’ shroud, presents HEM values of 0.522×10³, 0.551×10³, 0.983×10³, and 1.008×10³mg/hr, while the erosion of the hub is negligible. The fourth and third ash samples exhibit higher values compared with the two other samples. As for the OGVs, the values of AERD are respectively equal to 0.095×10⁻⁶, 0.098×10⁻⁶, 0.136×10⁻⁶, and 0.141×10⁻⁶ ${\rm{mg}}/{\rm{s}}.{\rm{m}}{{\rm{m}}^2}$ , and their respective HEM values are 2.277×10³, 2.341×10³, 3.245×10³, and 3.375×10³mg/hr. Furthermore, the OGVs’ shroud depicts values of HEM equal to 0.164×10³, 0.165×10³, 0.565×10³, and 0.817×10³mg/hr, respectively, while the hub has insignificant erosion.

The spinner is characterized by values of AERD equal to 0.341×10⁻⁶, 0.375×10⁻⁶, 0.973×10⁻⁶, and 1.099×10⁻⁶ ${\rm{mg}}/{\rm{s}}.{\rm{m}}{{\rm{m}}^2}$ , and the subsequent HEM values are 0.736×10³, 0.811×10³, 2.10×10³, and 2.373×10³mg/hr between first, second, third and fourth ash samples, respectively. Despite the low erosion rate density in the Pitot intake, its large geometry resulted in high values of HEM equal to 0.115×10³, 0.113×10³, 0.161×10³, and 0.163×10³mg/hr between the four ash samples, respectively, where the fourth sample produced more erosion.

Finally, Figs. 25–27 compare between the values of HEM of the Pitot intake, spinner, entire fan blades, and entire IGVs and OGVs, depicting essentially a marked increase with ash concentration. The highest values of HEM correspond to the third and fourth ash samples. As outlined in Table 5, in these front components of an HBTFE, the fourth volcanic ash sample (Kelud), followed by the third sample (Etna), caused the most severe erosion due to their relatively broader size range (0.5–300 microns) and (0.5–140 microns), combined with higher density and quartz content.

Table 5. Compared HEM at high ash concentrations

The continual removal of material from the surfaces of this transonic fan stage, particularly the rotor blades, causes a rapid decline in aerodynamic performance and lifespan. By simulating various operating situations, ash particle types, and concentration, a correlation can be formed between particle type and concentration, operating conditions, blade geometry deterioration, performance degradation, and lifespan.

9.0 Conclusion

Large HBTFEs are well-known for their high blade speed and transonic flow velocities, making them susceptible to erosion damage. This research is intended to replicate scenarios of encounters with volcanic ash and assess how erosion affects the front components of an engine.

The findings effectively demonstrate the behaviour of ash particles and the severity of erosion in each component. The Pitot intake design seems to cause the large-size ash particles to deflect downwards, strike the spinner, bounce back, and then travel in ballistic trajectories to reach the upper sections of the fan blade and shroud. As a result, a large quantity of ash particles impact the fan blade from the hub to around 80% of its span. The fan blade’s PS is almost entirely eroded, with the highest erosion rates along the LE and aft of the surface. In contrast, the SS of the fan blade displays erosion scattered over two-thirds of the surface. The rotor’s hub presents almost uniform erosion patterns, while the shroud depicts scattered erosion of lower rates. The spinning motion of the blades appears to fling the heavy ash particles into the secondary duct, preventing them from entering the engine core. The smaller ash particles issued from the lower part of the fan blade pass via the IGVs to the engine core, eroding the entire PS and SS of the vanes. Furthermore, erosion is shown to spread across the whole PS of the OGVs, with higher erosion rates concentrated in the second half of the surface, as opposed to the SS which erodes in the first half.

The highest levels of HEM correspond to the volcanic ash samples from the Kelud and Etna volcanoes, most likely due to larger particle sizes, higher density, and quartz content. The presence of volcanic ash in the atmosphere poses a major concern to flight routes in southern Europe, Asia, America, and northern Europe. Therefore, forecasting the effects of these abrasive particles on engine performance and mechanical integrity is critical for assuring flight safety and developing solutions to minimise erosion and extend engine lifespan.