Simulation of Particle Trajectories in Gas Turbine Components and Assessment of Unsteady Effects Using an Efficient Eulerian-Lagrangian Technique

Abstract

In recent years, CFD has proven to be a very useful asset to help with predicting complex flows in a wide range of situations, including multiphase and gas-particle flows. On this track, numerical modelling of particle-laden flows in multistage turbomachinery has become an important step in helping to analyse the behaviour of a discrete phase in gas turbines. Furthermore, unsteady effects due, for example, to rotor–stator interaction may have an effect on trajectories and capture efficiencies of the discrete phase. Unfortunately, computational times for transient simulations can be exceedingly high, especially if a discrete-phase needs also to be simulated. For this reason, this work reports a new method for the efficient and accurate simulation of particle-laden flows in gas turbine engines components. The Harmonic Balance Method is exploited to gain orders of magnitude speedup exploiting the idea that once the flow field has been embedded in the spectral basis, it can be reconstructed at any desired time. In this way, not only can the computational time needed to reach convergence of the flow field be dramatically reduced, but there is also no need to keep simulating the flow field during particle tracking. On the contrary, the continuous phase field can be retrieved at any desired time through flow reconstruction. This technique is conceptually simple, but, to the authors’ knowledge, has never been applied so far in particle-laden flow simulations and represents a novelty in the field. First, the implementation of the method is described, and details are given on how phase-lagged boundary conditions can be applied to flow and particles to further speed up the calculation. Then, some relevant case studies are presented to highlight the performance of the method.

1. Introduction

The requirements for efficiency and safe operation of gas turbines have become of paramount importance in modern society. One of the main causes of performance degradation in gas turbines is solid particle ingestion, causing erosion and deposition on aerodynamic surfaces. Land-based units can ingest soot as a consequence of operation in highly polluted industrial environments, or sand particles during operation in desertic places. Filtration systems tend to limit the particles entering the machinery, but they are unable to completely prevent their ingestion [1]. Aircraft engines can ingest particles at takeoff and landing during sandstorms, or encounter ash clouds transported by volcanic eruptions to cruise altitude. Dunn et al. [2] showed that it is possible to consume the surge margin very quickly when the engine operates in a dust cloud. Clearly, investigation of fouling and erosion in gas turbines has had an important relevance for academic researchers and engine manufacturers. Many numerical and experimental investigations about these phenomena can be found in the literature. The interested reader is referred to [3,4,5] for comprehensive reviews.

In recent years, CFD has proven to be a very useful asset to help with predicting complex flows in a wide range of situations, including multiphase and gas-particle flows. Consequently, numerical modelling of particle-laden flows in turbomachinery has become an important step in analyzing the behaviour of a discrete phase in gas turbines. On the basis of this premise, the main driver of this work was the necessity of efficiently simulating particle-laden flows in gas turbine components. For turbomachinery problems, Eulerian-Lagrangian methods are the most commonly employed due to their accuracy and simplicity. On the other hand, such method requires a large number of particles to be tracked, so that its major drawback is the computational cost. This could translate into a massive deployment of resources for complex turbomachinery domains. Several authors used different methods to model solid particles in gas turbine passages. Ghenaiet [6] and Hamed et al. [7] studied erosion in an axial flow turbine using a steady simulation with a frozen—rotor interface to couple the different stages. It must be said that this method did not include unsteady effects and, in addition, the erosion pattern and the flow depended on the clocking between the turbine rows. Tabakoff et al. [8] studied erosion in a turbine stage. They also employed steady-state simulations, but they introduced random particles circumferential redistribution at the stator/rotor interface. This allowed taking into consideration the uniforming effect of reciprocal motion on time-averaged particle distribution at stator outlet. In a similar fashion, Mustafa [9] investigated droplet trajectories during online washing of a multistage compressor with the commercial code CFX-TASCflow. Yang and Boulanger [10] simulated the full annulus of an axial fan and compared the erosion rates obtained for steady and transient simulations. Zagnoli et al. [11] and Prenter et al. [12] used ANSYS FLUENT to simulate particle-laden flow in an axial turbine stage with steady and unsteady simulations. A different methodology was adopted by Aldi et al. [13] for analyzing particle ingestion in an axial compressor stage. In this case, the analysis was carried out with separate particle injections for the isolated rotor and stator rows. The results of the isolated rotor and stator characterisations were post-processed to include the interaction between the rotor and the discrete phase while analyzing deposition on the stator. This technique avoids the explicit treatment of particles crossing stator/rotor interfaces but includes only a first-order approximation of the effect of the upstream blade row. Recent advances in the simulation of particle-laden flows in multistage turbomachinery were carried out by Oliani et al. [14] by developing a unified and simple treatment for each type of turbomachinery interface. The baseline methodology is the same employed in this work.

Following the previous studies, the aim of this work is to propose a general methodology for the simulation of particle-laden flows in gas turbine components. We highlight some important features regarding particle trajectories and impacts on aerodynamic surfaces. Moreover, we investigate the impact of the flow unsteadiness caused by rotor–stator interaction on particle trajectories. This topic has still not been addressed in detail in the literature. Unfortunately, computational times to carry out fully transient simulations can be exceedingly high. This is all the more so when a discrete-phase needs also to be simulated. The reason for this is that, even after the convergence, the flow field needs to be continuously updated in order to include the unsteady effects on particle trajectories, often leading to unacceptable computational requirements. For these reasons, we propose an efficient method to overcome this issue. Specifically, we employ the Harmonic Balance Method (HBM) to tackle the problem. This method was devised relatively recently by Hall [15] and is based on casting a set of equations in the frequency domain in order to switch the time derivative with a spectral operator. This technique has gained popularity over the years due to its efficiency compared to fully transient calculations for the continuous phase. To gain orders of magnitude speed up, the idea that once the flow field has been embedded in the spectral basis, it can be reconstructed at any desired time, is exploited. In this way, not only can the computational time needed to reach convergence of the flow field be dramatically reduced, but there is also no need to keep simulating the flow field during particle tracking. On the contrary, the continuous phase field can be retrieved at any desired time through flow reconstruction. This technique is rather straightforward but has never been applied so far in particle-laden flow simulations and represents a novelty in the field.

The paper is organised as follows: In Section 2, the governing equations for the continuous and the discrete phase are described. Then, we briefly review the HBM and its implementation for turbomachinery problems. In Section 3, we introduce the HBM-reconstruction method for the efficient coupling between HBM and Lagrangian tracking of the discrete phase. The method is validated on an ad hoc test case in Section 4. Finally, in Section 5, we investigate the discrete phase behaviour in two representative gas turbine components. Particle trajectories and impacts on aerodynamic surfaces are discussed, and unsteady effects are assessed using the proposed HBM-reconstruction method.

2. Numerical Modelling

2.1. Governing Equations for the Continuous and the Discrete Phase

A new implicit density-based solver was implemented in OF by the authors to carry out the calculations of the continuous phase [16]. The unsteady three-dimensional compressible Reynolds-averaged Navier–Stokes (RANS) equations are written in integral form:

(1)${\int }_V\frac{\partial \mathit{Q}}{\partial t}dV+{\int }_{\partial V}\left({\mathit{F}}_{\mathit{c}}-{\mathit{F}}_{\mathit{v}}\right)dS=0$

where V and $\partial V$ are the control volume and the corresponding closed surface, respectively. The conservative variables $\mathit{Q}$, the convective flux ${\mathit{F}}_{\mathit{c}}$ and the diffusive flux ${\mathit{F}}_{\mathit{v}}$ arrays are equal to

(2)$\mathit{Q}=\left[\begin{array}{c}\rho \\ \rho \mathit{u}\\ \rho E\end{array}\right],{\mathit{F}}_{\mathit{c}}=\left[\begin{array}{c}\rho \mathit{u}·\mathit{n}\\ (\rho \mathit{u}\otimes \mathit{u})·\mathit{n}+p\mathit{n}\\ \rho \mathit{u}H·\mathit{n}\end{array}\right],{\mathit{F}}_{\mathit{v}}=\left[\begin{array}{c}0\\ \mathit{\tau }·\mathit{n}\\ (\mathit{\tau }·\mathit{u}+\mathit{q})·\mathit{n}\end{array}\right]$

where $\mathit{n}$ is the face normal vector, $\mathit{u}$ is the velocity, $\rho$ is the density, p is the static pressure, E is the total internal energy, H is the total enthalpy, $\mathit{\tau }$ is the viscous stress tensor and $\mathit{q}$ is the heat flux vector. To preserve the code structure and keep the freedom to choose among the many built-in OF turbulence models, turbulence equations are solved in a segregated manner. Finally, to relate the pressure and enthalpy to conservation variables, the system of equations is completed using the ideal gas law. Applying a finite-volume discretisation to Equation (1), a semi-discretised form is obtained:

(3)$V{\mathit{D}}_{\mathit{t}}\mathit{Q}=\mathit{R}\left(\mathit{Q}\right)$

where ${\mathit{D}}_{\mathit{t}}$ is the time derivative operator and $\mathit{R}\left(\mathit{Q}\right)$ is the residual term. The system of equation is linearised and solved using the GMRES solver [17] combined with a LU-SGS preconditioner [18]. In OF, the built-in solvers and the underlying code structure are based on SIMPLE formulation. For this reason, the new ICSFoam library [16] has been implemented to carry out the implicit solution of the system of equations using a density-based formulation. Now, it is possible to assemble an arbitrary number of blocks and equations, thanks to the generalisation of the structure for block-coupled matrices. As pointed out in Section 2.2, this is necessary for the solution of a fully implicit version of the HBM. Approximate Riemann solvers are used to calculate the inviscid fluxes at cell interfaces. Finally, the support for Multiple Reference Frame (MRF) solution in rotating frames has been included into the solvers. Second order accuracy in space is obtained with the MUSCL approach. The Van Leer limiter is applied to primitive variables for the reconstruction of the solution from cell average values to the faces of control volumes. A classic second-order accurate central difference formula is used for viscous fluxes.

For what concerns the discrete phase, we solve the dynamic equation of motion for each particle in the domain:

(4)$\begin{array}{c}\frac{d{\mathit{x}}_{\mathit{p}}}{dt}={\mathit{u}}_{\mathit{p}}\end{array}$

(5)$\begin{array}{c}{m}_p\frac{d{\mathit{u}}_{\mathit{p}}}{dt}={\mathit{f}}_{\mathit{I}}+{\mathit{f}}_{\mathit{NI}}\end{array}$

where ${\mathit{f}}_{\mathit{I}}$ and ${\mathit{f}}_{\mathit{NI}}$ are the inertial and non-inertial forces, respectively. Inertial forces can include different actions such as drag, gravity, Saffman lift and virtual mass. For turbomachinery applications in which a steady-state approach is adopted (use of MRF), ${\mathit{f}}_{\mathit{NI}}$ includes centrifugal and Coriolis forces, defined as

(6)$\begin{array}{c}{\mathit{f}}_{\mathit{centr}}=-m_p\mathbf{\Omega }\times \mathbf{\Omega }\times {\mathit{x}}_{\mathit{p}}\end{array}$

(7)$\begin{array}{c}{\mathit{f}}_{\mathit{Cor}}=-2m_p\mathbf{\Omega }\times {\mathit{u}}_{\mathit{p}}\end{array}$

where $\mathbf{\Omega }$ is the shaft angular speed in rotating domains and zero elsewhere. Depending on whether the continuous phase field is updated or not after each time-step of Equations (4) and (5), an unsteady or a steady Eulerian-Lagrangian algorithm is obtained, respectively. To track particles across the domain, we also need to know the current cell of the mesh where they reside. In OF, each cell is decomposed into tetrahedrons and the tracking is performed in barycentric coordinates. This significantly improves the accuracy and the robustness of the algorithm. At each step, the current cell is updated using a tracking routine; see, e.g., [19].

The treatment for particles crossing mixing plane interfaces in steady-state simulations is carried out according to the technique described in [20]. Particles are redistributed in a random manner in the tangential direction to simulate the reciprocal motion of stator and rotor. The velocity vector is rotated according to the redistribution angle to maintain the radial and tangential components. Since mixing planes usually divide one frame of reference from another, particle velocity is finally updated correspondingly. In addition, to ensure mass conservation in a full-annulus sense, the particle transfer function has been modified by adding or deleting particles basing on the coverage in the tangential direction of the two rows.

2.2. The Harmonic Balance Method

A detailed description of the implementation of a fully implicit HBM in OF can be found in [21]. In general, spectral techniques are based on a frequency-domain formulation of periodic unsteadiness. The main idea is that a nonlinear solution with $N_f$ harmonics is equivalent to $2N_f+1$ coupled steady flow problems. Therefore, the HBM leverages the signal sparsity in the Fourier basis to substitute an unsteady solution with a coupled set of steady-state calculations. In this work, we follow the mathematical formulation proposed in the time-spectral approach of Hall et al. [15].

Let us suppose that the vector $\mathit{Q}$ of conservative variables evolves with a known set of K frequencies ${\omega }_k$ in time. Now, we expand the vector ${\mathit{Q}}_{\mathit{l}}$ (representing $\mathit{Q}$ in the control volume l) and the residual term $\mathit{R}\left({\mathit{Q}}_{\mathit{l}}\right)$ in Fourier series

(8)${\mathit{Q}}_{\mathit{l}}\approx {\widehat{\mathit{Q}}}_{\mathit{l},\mathbf{0}}+\sum _{k=-K}^K{\widehat{\mathit{Q}}}_{\mathit{l},\mathit{k}}{e}^{-im{\omega }_{k}t},\mathit{R}\left({\mathit{Q}}_{\mathit{l}}\right)\approx {\widehat{\mathit{R}}}_{\mathit{l},\mathbf{0}}+\sum _{k=-K}^K{\widehat{\mathit{R}}}_{\mathit{l},\mathit{k}}{e}^{-im{\omega }_{k}t}$

where ${\widehat{\mathit{Q}}}_{\mathit{l},\mathbf{0}},{\widehat{\mathit{Q}}}_{\mathit{l},\mathit{k}}$ are the $N_T=2K+1$ Fourier coefficients of the time-average and the positive and negative frequencies terms, respectively. These coefficients are not known a priori. Substituting Equation (8) into Equation (1), we obtain, for the control volume l,

(9)$iV\mathit{A}\widehat{{\mathcal{Q}}_l}=\widehat{{\mathcal{R}}_l}$

where

(10)$\widehat{{\mathcal{Q}}_l}=\left\{\begin{array}{c}{\widehat{\mathit{Q}}}_{\mathit{l},\mathbf{0}}\\ {\widehat{\mathit{Q}}}_{\mathit{l},\mathbf{1}}\\ ⋮\\ {\widehat{\mathit{Q}}}_{\mathit{l},\mathit{K}}\\ {\widehat{\mathit{Q}}}_{\mathit{l},-\mathit{K}}\\ ⋮\\ {\widehat{\mathit{Q}}}_{\mathit{l},-\mathbf{1}}\end{array}\right\},\widehat{{\mathcal{R}}_l}=\left\{\begin{array}{c}{\widehat{\mathit{R}}}_{\mathit{l},\mathbf{0}}\\ {\widehat{\mathit{R}}}_{\mathit{l},\mathbf{1}}\\ ⋮\\ {\widehat{\mathit{R}}}_{\mathit{l},\mathit{K}}\\ {\widehat{\mathit{R}}}_{\mathit{l},-\mathit{K}}\\ ⋮\\ {\widehat{\mathit{R}}}_{\mathit{l},-\mathbf{1}}\end{array}\right\},\mathrm{and} \mathit{A}=\left[\begin{array}{ccccccc}0& \cdots & \cdots & \cdots & \cdots & \cdots & 0\\ ⋮& {\omega }_1& & & & & \\ ⋮& & \ddots & & & & \\ ⋮& & & {\omega }_K& & & \\ ⋮& & & & -{\omega }_K& & \\ ⋮& & & & & \ddots & \\ 0& & & & & & -{\omega }_1\end{array}\right]$

Equation (9) is a set of equations coupled through the nonlinear term $\widehat{{\mathcal{R}}_l}$. As reported by Hall et al. [15], the system of equations is more easily solved by modelling these terms implicitly, transforming them back into the time domain. This allows for working with time-domain solutions stored at $2N_f+1$ time levels, and the flux term can be computed as usual for each time level. The choice of the sampling points for the time levels is non trivial and influences the entries of the Fourier matrix $\mathit{E}$ that relates the Fourier coefficients vectors $\widehat{\mathcal{Q}}$ and $\widehat{\mathcal{R}}$ with the time-domain counterparts:

(11)$\widehat{{\mathcal{Q}}_l}=\mathit{E}{\mathcal{Q}}_l,\mathrm{and}\widehat{{\mathcal{R}}_l}=\mathit{E}{\mathcal{R}}_l$

where $\mathcal{Q}$ contains the variables $\mathit{Q}$ stored at $N_T$ subtime levels.

If the spectrum is composed of a single fundamental frequency $\omega$ and its harmonics (periodic flow), the Fourier matrix is always well-conditioned. Indeed, in such a scenario, a uniform sampling of the base period $T=2\pi /\omega$ leads to a condition number equal to 1 because of the orthogonality of the Fourier basis. On the contrary, if the frequency set is arbitrary, finding a proper set of time instants corresponding to a low condition number of the matrix $\mathit{E}$ is way more difficult.

Indeed, in many applications of industrial interest, almost-periodic signals are considered, where multiple discrete dominant frequencies exist. These frequencies need not be integral multiples of each other, a typical example being multi-rows turbomachinery where the number of blades and vanes are not multiples. When this is the case, we can still project the flow variables onto a set of non-harmonically related frequencies. It must be emphasised that now the frequency set does not form, in general, an orthogonal basis anymore. Additionally, there is no analytical expression for the matrix $\mathit{E}$. Instead, we define its inverse ${\mathit{E}}^{-\mathbf{1}}$ analytically by its components as $E_{n,k}^{-1}=e^{i{\omega }_k{t}_n}$, and then compute numerically $\mathit{E}$ [22,23]. The correct selection of the time levels is of paramount importance for the convergence of the simulation. Indeed, the higher is the condition number of the matrix $\mathit{E}$, the more errors arising from the HB term will be amplified during the iterative solution [22].

In the present paper, we follow the OptTP approach introduced by Nimmagadda et al. [24]. This method has negligible performance drop with respect to other methods, is robust, and is easier to implement. The basic idea is to compute the condition number for many possible periods and then uniformly sample the time period $T^{*}$ that minimizes the condition number of matrix $\mathit{E}$. Using Equation (11), and multiplying on the left Equation (9) by the inverse matrix ${\mathit{E}}^{-\mathbf{1}}$, we have

(12)$V\mathit{D}{\mathcal{Q}}_l={\mathcal{R}}_l$

where $\mathit{D}={\mathit{E}}^{-\mathbf{1}}i\mathit{A}\mathit{E}$. We can observe that Equation (12) is a system of steady-state equations. A comparison with Equation (3) shows that this procedure has allowed us to find an approximation of the time derivative operator ${\mathit{D}}_{\mathit{t}}{\mathit{Q}}_{\mathit{l},\mathit{j}}\approx {\left(\mathit{D}{\mathcal{Q}}_l\right)}_j$ in terms of the flow field snapshot j in each control volume l. The acceleration methods for steady-state simulations can be used (e.g., multigrid and local time-stepping), greatly reducing the computational cost of the solution.

Now, we must solve a set of equations coupling all the subtime levels through the HB term $V\mathit{D}\mathcal{Q}$. We can directly discretise this term implicitly since the operator $\mathit{D}$ is linear: $V\mathit{D}{\mathcal{Q}}^{n+1}=V\mathit{D}{\mathcal{Q}}^n+V\mathit{D}\Delta {\mathcal{Q}}^n$. Including this linearisation into Equation (12), one obtains:

(13)$\mathit{M}\Delta {\mathcal{Q}}^n={\mathcal{R}}^n-V\mathit{D}{\mathcal{Q}}^n$

where $M_{ii}=\left(-\frac{\partial {\mathit{R}}_{\mathit{i}}^{\mathit{n}}}{\partial {\mathit{Q}}_{\mathit{i}}^{\mathit{n}}}+V{D}_{i,i}\right)$ and $M_{i}j=V{D}_{i,j}\mathit{I}$. In the present work, the GMRES linear solver [17] with an LU-SGS preconditioner [18] are employed for the solution of the linearised system. Hence, the flow variables are updated for all the snapshots at the same time. The block-coupled solver structure described in [16] allows for a convenient assembly of the matrix $\mathit{M}$ for the linear system.

2.3. Turbomachinery Boundary Conditions

In multi-stage turbomachinery, flow unsteadiness can be expressed as a linear combination of the blade passing frequencies of the neighbouring rows [25]. The most common method that can be found in the HBM literature to deal with multi-stage turbomachinery is to resolve different time instants for each blade row. In this way, of course, the snapshots solved in each row do not, in general, match each other. To solve this problem, the data between adjacent zones are exchanged through spectral interpolation of the flow quantities from the donor side to the time levels of the receiving side [26], or via matching of the spatial Fourier coefficients [27]. The drawback of these techniques is that, if the interpolation is not combined with non-reflecting conditions [27] or oversampling [23,26], they can generate spurious wave reflections which may in turn corrupt the solution. In addition, if Lagrangian tracking must be coupled with the HBM, it is much easier to deal with synchronised time instants. Therefore, in this work, we follow the synchronised HBM strategy proposed by Crespo and Contreras [28]. This method ensures the flow continuity between stator and rotor interfaces by using the same time instants in all the blade rows.

Another important aspect concerns the exploitation of spatial periodicity in the azimuthal direction, as often carried out in turbomachinery simulations. It must be said, however, that, in unsteady calculations, it is necessary to have the same pitch for all different rows. For this reason, often half or even the full annulus of the machine needs to be simulated. The usage of phase-shifted boundary conditions [29] is a way to resolve the problem by exploiting the known spatio-temporal periodicity of the flow [30]. The method allows for reducing the computational domain to just one passage per row, independently of the blade number ratio. Indeed, there is a correspondence between the flow field inside a blade passage at a certain time t and the adjacent passage at a different time $t+\mathcal{T}$

(14)$\mathit{Q}(r,\theta +\Delta \theta ,z,t)=\mathit{Q}(r,\theta ,z,t+\mathcal{T})$

where $\mathcal{T}$ is called the time lag between the two passages, $\Delta \theta$ is the row pitch and $\theta$ is the tangential coordinate. Let us consider for simplicity a single stage, the results being generalizable to multiple stage. The time lag can be calculated as the phase of a wave traveling at a rotational speed ${\omega }_k=2\pi k{f}_{BP}$ in the azimuthal direction: $\mathcal{T}=k\sigma /{\omega }_k$, where $\sigma$ is the interblade phase angle (IBPA) and $f_{BP}$ is the blade passing frequency of the opposite row. To enforce phase-lagged conditions, one needs to update the Fourier coefficients on the boundaries at each time-step. Since HBM is intrinsically related to the Fourier decomposition of the solution, phase-shifted boundary conditions are well-suited and much simpler to implement in frequency-domain solvers. Equivalently, Equation (14) can be expressed in the frequency domain as

(15)$\sum _{k=-K}^K{\widehat{\mathit{Q}}}_{\mathit{k}}(r,\theta +\Delta \theta ,z)e^{i{\omega }_{k}t}=\sum _{k=-K}^K{\widehat{\mathit{Q}}}_{\mathit{k}}(r,\theta ,z)e^{i{\omega }_{k}t}{e}^{i{\omega }_k\mathcal{T}}$

With this in mind, we can reformulate the relation between the time instants on the two periodic boundaries as [23]

(16)$\mathcal{Q}(r,\theta +\Delta \theta ,z)={\mathit{E}}^{-\mathbf{1}}\mathit{SE}\mathcal{Q}(r,\theta ,z)$

where $\mathit{S}$ is diagonal with components $S_{kk}=e^{ik\sigma }$. The implementation of phase-lagged boundary conditions in the HBM solver is detailed in [21].

3. Coupling the HBM and Lagrangian Tracking

Previous studies showed that unsteady effects may have a more or less pronounced impact on the predictions of erosion and deposition in turbomachinery [10,12,31]. It is difficult to know a priori if the discrete phase will show an influence from rotor–stator interactions and to what extent for a particular case. Therefore, in general, one must know and take into account that steady simulations can experience a certain degree of error. The alternative is, of course, to perform unsteady Lagrangian tracking calculations, with the associated high computational cost. The term unsteady Lagrangian tracking denotes the calculations in which the continuous phase fields are updated after each Lagrangian step for the discrete phase. Therefore, after each particle $p_i$ in the domain is advanced of $\Delta {\mathit{x}}_{{\mathit{p}}_{\mathit{i}}}$ in space, the continuous phase fields are advanced from time t to $t+\Delta t$, where $\Delta t$ is the physical time-step of the calculation. The entire process is then repeated for as long as needed to track all the particles throughout the domain.

To gain orders of magnitude speed up the simple idea that once the flow field has been embedded in the spectral basis, it can be reconstructed at any desired time is exploited. In this way, not only can the computational time needed to reach convergence for the flow field be dramatically reduced, which is the advantage of using HBM over transient calculations, but there is also no need to keep simulating the flow field during particle tracking. On the contrary, the continuous phase field can be retrieved at any desired time through flow reconstruction. Therefore, the continuous phase quantities needed to calculate the forces acting on solid particles can be directly retrieved, without actually computing the flow field. This new approach will be called from now on the HBM-reconstruction method. An implicit assumption of this method is that particles do not affect the fluid flow (one-way coupling), but their trajectory is influenced by the continuous phase through fluid-dynamic forces. This is a typical assumption in turbomachinery flows since usually the particle’s volume fraction is very low ($<<$10%). This technique is rather straightforward but has never been applied so far in particle-laden flow simulations and represents a novelty in the field.

Some authors [10] have adopted a different approach by saving many timesteps of the converged flow field across one fundamental period of the current blade row. Then, after solid particles’ injection, the snapshot corresponding to the actual instant is loaded from the memory, instead of computing it. In this way, particle tracking can proceed as usual using the saved timesteps. However, this approach is not efficient and has several drawbacks. Firstly, it requires a lot of memory because the entire flow field needs to be stored for all the desired instants of tracking (which are typically tens or one hundred for one passing period). For larger simulations, even modern computers can not keep so many timesteps in the RAM due to memory limitations. Therefore, a continuous process of loading and unloading of the saved timesteps from the disk storage is needed. This, though much cheaper than recomputing the flow field at each timestep, results anyway in a very cumbersome and time-consuming process. Indeed, even with modern SSD, the memory bandwidth is limited compared to RAM, resulting in a pronounced access latency during the process. Just to give an idea, the memory bandwidth of RAM is typically around 10 GB s${}^{-1}$, depending on the computer, while SSD bandwidth is about 500 MB s${}^{-1}$. Things are even worse when a HDD storage is used, for which load/write speed is about 100 MB s${}^{-1}$. On the contrary, embedding the flow field with a ROM such as the HBM allows for keeping all the necessary data in the main memory and therefore there is no need to load data from the storage drive. Therefore, memory and computational efficiency are achieved. Of course, if the Fourier basis is chosen as the basis for the ROM, one could compress the data after performing a conventional transient calculation and reconstruct the field at any desired time. During the Lagrangian tracking step, this would be equivalent to the method presented above. Anyway, the HBM allows for obtaining accurate results in a reduced turnaround time also during the continuous phase calculation, and naturally embeds the results in the spectral basis without need for further post-processing. Hence, this strategy is preferable and has been used in this work.

Another major disadvantage of saving and loading the snapshots for the tracking is that the timestep of the calculation can not be modified. Indeed, since the snapshots of the flow field used for the tracking are saved and loaded from the storage, once the timestep has been selected, the only way to modify the snapshots is to recompute them. Instead, with the HBM, the flow field can be reconstructed at any time and one can therefore choose and modify the timestep used for Lagrangian tracking as preferred, with no restrictions. This is particularly useful if a suitable timestep for particle tracking can not be established a priori.

Of course, with the HBM being a reduced order model, accuracy losses can be expected if the number of harmonics retained in the spectrum is not sufficient. Indeed, errors arising from inaccurate reconstruction of the flow field can potentially accumulate along particle trajectory, leading to an increased computational error. A sensitivity analysis on the number of harmonics retained in the spectrum is carried out in Section 4 to further investigate this aspect.

To fully exploit the potential of the HBM, it is also desirable for the new method to allow particle tracking in reduced domains of a single passage per row. This can be achieved if the concept of phase-shift between periodic boundaries is extended to the discrete phase, and represents a further novelty in the methodology. The idea is based on the observation that, when a particle crosses a periodic boundary with a phase-shift boundary conditions imposed on it, another particle must reappear on the other side with the same velocity after a time lag:

(17)$\mathcal{T}={\sigma }_k/{\omega }_k\forall k$

where ${\sigma }_k$ is the IBPA associated with a travelling wave of rotational speed ${\omega }_k$. Remember that, even if multiple fundamental frequencies are present in the same row, Equation (17) must hold for each of them. In other words, the time-lag is a constant and can be calculated a priori for each row.

Please notice that, physically, the particle reappearing on the other side of the cyclic can not be the same particle that exits from the first side. Despite this fact, for computational purposes, the same particle is used on both halves of the cyclic since this is perfectly equivalent and easier to implement. Another way to see this is considering that, if a particle exits from a cyclic at a time t and is injected from the other half at time $t+\mathcal{T}$, then this is equivalent to continuing to track the same particle in the adjacent passage of the row at time t. Indeed, the flow fields of two adjacent passages are exactly the same after the time lag $\mathcal{T}$. This allows for taking advantage of single passage reduction also during the Lagrangian tracking step, thus avoiding reconstructing the flow field in adjacent passages. The underlying hypothesis in doing this is that particles’ injection inside the domain has the same periodicity of the flow. As a matter of fact, there is hardly any reason why this is not true in turbomachinery flows. In other words, this means that it is very unlikely that particles are injected with a temporal frequency different from the one of the flow field. In addition, particles are typically injected uniformly (in space and time) from the domain inlet, in which case this assumption is certainly true.

The implementation is as follows. For phase-lag cyclics, it is sufficient to transfer the particle from one side to the other and deactivate the particle for a time corresponding to the time lag. After a time $\mathcal{T}$ has passed, the particle is reactivated and the tracking proceeds as usual. It is also emphasised that, necessarily, the two halves of the phase-lag cyclic have IBPAs of the opposite sign (if half A is in advance with respect to half B, then half B is in delay of the same amount with respect to A). Although negative phase angles do not represent a problem, negative time lags are of no use for particle tracking. For this reason, it is always necessary to convert negative IBPAs into the corresponding positive time lags. For phase-lag AMI interfaces, things are slightly more complex. A strategy similar to the one described in [20] is used to find the correct position and cell on the other side of the interface. The particle is deactivated for a time $N_{Cp}\mathcal{T}$, where $N_{Cp}$ is the number of the copy on the receiving side and $\mathcal{T}$ is the time lag of the receiving blade row. Once again, $\mathcal{T}$ is always positive, depending on whether $N_{Cp}$ corresponds to a copy transformed in the forward or backward direction. The schematic procedure for particle transfer between phase-shifted interfaces is depicted in Figure 1.

4. Validation of the Method

The validation of the method is carried out on the rotor blade of the first stage of an $E^3$ axial turbine [32]. A fully transient simulation is compared with the results obtained with the HBM-reconstruction technique. The geometry is two-dimensional, and an artificial wake is generated and imposed at the domain inlet, according to:

(18)$\begin{array}{c}{p}_0\left(y\right)={\overline{p}}_0[1-\Delta p·e^{-0.693{\left(2\frac{y}{L}\right)}^2}]\end{array}$

(19)$\begin{array}{c}T\left(y\right)=\overline{T}[1+\Delta T·e^{-0.693{\left(2\frac{y}{L}\right)}^2}]\end{array}$

where ${\overline{p}}_0$ is the average total pressure and $\overline{T}$ is the average temperature at domain inlet, respectively equal to 240 $\mathrm{k}$$\mathrm{Pa}$ and 600 $\mathrm{K}$. The outlet total pressure is set equal to 144 $\mathrm{k}$$\mathrm{Pa}$. $\Delta p$ and $\Delta T$ are the total pressure and temperature variations equal to of 0.10969 and 0.015, respectively. L is the azimuthal wake width and has been selected as $20\%$ of the blade pitch. The average values of total pressure and temperature, as well as the flow inlet angle, are chosen according to a preliminary steady-state simulation of the entire stage with a mixing plane interface.

The standard k-$\epsilon$ turbulence model with its original constants and parameters is employed for the RANS closure of the turbulent viscosity ${\nu }_t$ [33]:

(20)$\begin{array}{c}\frac{\partial \rho k}{\partial t}+\nabla ·\left(\rho k\mathit{u}\right)=\nabla ·(\rho \frac{{\mu }_t}{{\sigma }_k}\nabla k)+P-\rho \epsilon \end{array}$

(21)$\begin{array}{c}\frac{\partial \rho \epsilon }{\partial t}+\nabla ·\left(\rho \epsilon \mathit{u}\right)=\nabla ·(\rho \frac{{\mu }_t}{{\sigma }_k}\nabla \epsilon )+\frac{C_1\epsilon }{k}P-C_2\rho \frac{{\epsilon }^2}{k}\end{array}$

(22)$\begin{array}{c}{\nu }_t=\frac{C_{\mu }{k}^2}{\epsilon }\end{array}$

where P is the turbulent kinetic energy production rate. The model constants employed are set according to the standard values: $C_{\mathsf{\mu }}=0.09$, $C_1=1.44$, $C_2=1.92$, ${\sigma }_k=1.00$, ${\sigma }_{\epsilon }=1.30$. The automatic wall treatment switches between the viscous and inertial sublayers values for the turbulent quantities by comparing the local $y^{+}$ with the estimated intersection of the viscous and inertial sublayers $y_{lam}^{+}$.

The mesh is composed of nearly $20,000$ elements. During the simulation, the wake profile described by Equations (18) and (19) is slid in the azimuthal direction at the domain inlet, to recreate an unsteady flow field similar to the one seen by the blade row in its relative frame during a simulation of the entire stage. The wavelength of the wake is twice the pitch of the rotor, to mimic an upstream vane row with a count ratio of 1:2 with respect to the blade row.

The fully transient calculation is performed on a domain including two blade passages (to take into account the wave length of the incoming wake) using periodic boundaries, while HBM simulations are carried out on a single passage domain exploiting phase-lag periodic boundary conditions. Figure 2a shows the instantaneous temperature field in the blade passages for the transient simulation. HBM solutions were obtained considering $4,5$ and 6 harmonics of the wake passing frequency. The number of sample points is obtained with the usual relation $2N_f+1$ and is respectively equal to: 9, 11, 13. Figure 2b shows instantaneous temperature profiles along the pitch in front of the rotor for the HBM solutions compared with the reference transient simulation. As can be seen, with six harmonics, the HBM results closely match the transient ones, while five harmonics are already sufficient to capture the shape of the profile well. Finally, a steady state calculation is also carried out using a mixed out state at the domain inlet with values of total pressure and temperature of ${\overline{p}}_0$ and $\overline{T}$.

For the unsteady calculations, particles are injected at the inlet following a probabilistic distribution in which the probability is high inside the incoming wake and very low outside. Particles are injected at the equilibrium velocity with the flow. In this way, a strongly non uniform incoming condition is generated for the particles. Indeed, the aim here is not so much to reproduce a realistic unsteady behaviour in the rotor as to create a test case in which unsteady effects are willingly enhanced to highlight the difference with a steady-state calculation. During the Lagrangian tracking steps of the fully transient and the HBM solutions, a timestep of $1\times {10}^{-6}$ is employed for the flow field, corresponding to a physical Courant number of 2. For the steady solution, instead, particles are injected uniformly from the domain inlet as if a mixing plane was placed between the stator and the rotor. In order to focus only on the correct implementation of the methodology, turbulent dispersion is not accounted for in the simulations of this case study. Particles are considered spherical and smooth, and the unique force acting on them is viscous drag. Three different diameters are injected to analyse the sensitivity of the results to particle size: $2,8$ and $16 \mathsf{\mu }\mathrm{m}$. Assuming a density of $2700 \mathrm{k}\mathrm{g}/{\mathrm{m}}^3$, the Stokes number of the particles based on the turbine axial chord is: $0.125,2.00$ and $8.00$, respectively. A perfectly elastic rebound is assumed when a particle impacts on the blade surface.

Figure 3 shows trajectories of $16 \mathsf{\mu }\mathrm{m}$ particles for transient and HBM-reconstruction methods at four equispaced instants in time. Transient results are shown on two blade passages while, of course, HBM results are available only on a single passage. The lower passage of the transient simulation is synchronised with the passage of the solution reconstructed using HBM. The qualitative comparison of the trajectories clearly shows that the new technique is able to reproduce the discrete phase movement across the domain. In particular, the non uniform injection and the variation in particles concentration due to wake stretching and distortion inside the passage are captured properly. Particles distribution is also very similar after they impact on the blade and exit the domain, meaning that impact angle and velocity are correctly reproduced.

Beside the instantaneous behaviour of the discrete phase, one is typically interested in obtaining reliable prediction of global parameters or time-averaged values quantifying impacts on the blade. Figure 4 reports the impact efficiency along the blade chord, defined as the ratio between the number of particles impacting on the blade and the total number of particles injected. For particles with Stokes $0.125$ and $2.00$, only the pressure side of the blade is shown, since almost no impacts were found on the suction side. Indeed, as also reported by Jacobs et al. [34], shock–particle interaction may create particle velocity components perpendicular to the main flow, producing transport perpendicular to the flow streamlines. This effect, combined with the curvature of the fluid flow streamlines is responsible for the very small number of impacts found on the suction side. Bigger particles instead tend to be projected towards the blade suction side after they rebound on the pressure side of the blade below (see Figure 3). A very good agreement is found for all diameters between HBM-reconstruction with six harmonics and transient cases (Figure 4a,c,e). Conversely, steady-state solutions are significantly different, showing that an unsteady treatment is necessary in this test case. A sensitivity analysis with respect to the number of harmonics is also carried out. A comparison of Figure 4b,d,f reveals that smaller particles are the most sensitive to the number of harmonics used for the simulation, being more influenced by the surrounding flow field due to the lower inertia. For these particles, six harmonics need to be retained in the spectrum to provide a good resolution of the impact efficiency. For $St=2.00$ and $St=8.00$, on the other hand, five harmonics are sufficient to obtain a good accuracy. From these observations, it can be inferred that, if the number of harmonics is adequate for the correct resolution of the flow field, a good approximation of the impact of unsteady effects on the discrete phase will also be obtained.

Inspection of computational times in Figure 5 reveals the potential of the new method for a remarkable speedup compared to conventional transient simulations. The CPU times are normalised with respect to the total simulation time for the steady case (notice the log scale on the y-axis). Total times are represented by the sum of the flow time (time to reach convergence of the flow field before injecting particles) and the Lagrangian time (time necessary to track all the particles until they exit the domain). It can be noticed that the time for the transient is two orders of magnitude greater than the steady simulation. Moreover, focusing just on the Lagrangian time (yellow bars), this is increased by three orders of magnitude with respect to steady-state. This observation clearly highlights the need for a reduction in the time required by this step of the calculation. This is accomplished by the HBM-reconstruction that reduces it by more than two orders of magnitude, making the bars diagram look very similar to the steady calculation, where the Lagrangian step represents only a minor contribution to the overall time. This is due to the combined effect of HBM and the reduction to a single passage with phase-lag BCs. This second contribution is expected to have a much bigger influence in real turbomachinery cases, where several passages typically have to be simulated. Overall, in this case, the total time for HBM-reconstruction is reduced by one order of magnitude compared to transient.

5. Results

In this chapter, the behaviour of a discrete phase inside two gas turbine components is analysed: a transonic axial compressor stage and a 1.5 axial turbine stage. Despite the geometries being taken from experimental rigs rather than actual engines’ components, they can be considered as representative cases of turbomachinery applications of interest. Steady and unsteady results are presented, where the latter are obtained thanks to the HBM-reconstruction method. Mixing plane interfaces are used to connect different rows in the steady-state simulations. The HBM-reconstruction method is used to include transient effects and compared with steady results to assess their importance. The behaviour of the discrete phase is analysed both in terms of particle trajectories and global/local parameters referring to particle impacts on aerodynamic surfaces.

5.1. NASA Stage 37

The 3D geometry of the NASA stage 37 [35] is selected as the first case study. The stage is composed of 36 blades and 46 vanes. The rotor angular velocity is 18,000 rpm. One passage per row is modeled thanks to the single passage reduction technique described in Section 2.3. The computational domain for the stage is composed of a structured hexahedral mesh of 1.3 million elements generated using Ansys Turbogrid. The grid parameters were chosen according to a preliminary grid sensitivity study; see [13]. In the spanwise direction, about 100 grid points are used, and 45 grid points in the circumferential direction. The grid is realised by employing an O-grid around the rotor blade and the stator vane, with local refinements near the hub and shroud regions. In addition, 200 points are used around the rotor blade. The layer of cells on blade walls is such that the $y^{+}$ values range from 1 to 10. The tip clearance between the blade and the shroud is resolved with a mesh composed of six nodes across the gap span. At the inlet, an absolute total pressure of $101,325 \mathrm{Pa}$ and an absolute total temperature of $288 \mathrm{K}$ are imposed. A turbulence intensity of 3% and a turbulent length scale of $2\times {10}^{-3}$ $\mathrm{m}$ were also prescribed at the inlet. Two operating conditions are considered for the study, namely choked flow and the near stall. To achieve these conditions, the choking mass flow rate was initially computed by imposing a static pressure outlet boundary condition. The computed choking mass flow rate was $20.42 \mathrm{k}\mathrm{g}/{\mathrm{s}}^{-1}$, which is $2.3\%$ below the experimental value of $20.9 \mathrm{k}\mathrm{g}/{\mathrm{s}}^{-1}$. The static pressure was then increased until the mass flow rate achieved the desired value, equal to $p/p_{amb}=1.39$ and $p/p_{amb}=1.66$ for the choked and near stall conditions, respectively. For a detailed description and validation of this case in OF, see [16,21]. In each row, four harmonics of the passing frequencies are considered in the spectrum. Figure 6 shows the instantaneous pressure and entropy contours at midspan for the choked flow condition. The flow field has been obtained by reconstruction at the desired time. Four adjacent passages are shown, considering the phase-shift between them.

Spherical drag and Saffman–Mei [36] lift effect are considered to be the only forces acting on the particles. Here, a constant density of 2700 $\mathrm{k}\mathrm{g}$/${\mathrm{m}}^3$ for the discrete phase is assumed due to the nature of the air contaminants, made up of a large part of sand and very small particles of soil. Particles are injected at the same velocity as the air flow from the stage inlet, with randomly positioned injection points. Particle diameter $d_p$ varies in the range $0.5÷32 \mathsf{\mu }\mathrm{m}$, while the Stokes number, calculated at the rotor and stator inlet sections, are reported in Table 1.

Small particles (<2$÷10 \mathsf{\mu }\mathrm{m}$) are more prone to stick to the blades, causing fouling [37]. Larger ones tend to impact without adhesion, causing the erosion of aerodynamic surfaces. In this work, the attention is focused on particles’ impacts and deposition on aerodynamic surfaces. For each diameter, a sufficient amount of particles to obtain a statistically significant number of impacts on the blades surface must be injected. Since smaller particles tend to follow the fluid streamlines better, only a small number hits the walls. For this reason, a smaller number of particles were injected with growing diameter, as reported in the last row of Table 1.

Tian and Ahmadi [38] pointed out that different turbulence models have an impact on the particles velocity deposition in horizontal and vertical ducts. Their sensitivity analysis highlights the influence of turbulence models, wall mesh refinement and particle dimensions. They show that the most important parameter is the non dimensional particle relaxation time defined as

(23)${\tau }^{+}=\frac{\left({\rho }_p\right)/\rho d_p^2{u}_{\tau }^2}{18{\nu }^2}$

The turbulence model can be important in the resolution of particle trajectories near the wall. The turbulence model used for the NASA stage 37 calculation is the standard $k-\epsilon$ model and near-wall effects are modeled by means of scalable wall functions. In a previous work, Aldi et al. [13] demonstrated that the $k-\epsilon$ model is suitable to model particle dispersion near the walls for the NASA Stage 37 geometry. Their analysis shows that the $k-\epsilon$ turbulence model overpredicts deposition for particles in Brownian (${\tau }^{+}<{10}^{-2}$) regions, but in the inertial (${\tau }^{+}>10$) region, the predicted deposition velocity trend is similar to experimental data. As reported in Table 1, the ${\tau }^{+}$ for the compressor stage simulations is in the range corresponding to the inertial region in accordance with Tian and Ahmadi [38].

Particle rebound velocity and direction are computed with a specific particle–wall interaction model imposed on rotor blades, stator vanes and endwalls. The normal $e_n$ and tangential $e_t$ restitution coefficients are defined according to Forder et al. [39], as functions of the particle impact angle $\alpha$. The correlations for the restitution coefficients are based on impingement testing using AISI 4130 carbon steel and sand. The analysis of particle adhesion is instead performed using the experimental results from Poppe et al. [40]. Based on the trends of sticking probability, we can define representative functions to relate the normal impact velocity $u_n$ and the sticking probability $S_p$. According to the analysis carried out by Suman et al. [41], the trend can be represented by two equations. The first one refers to lower normal impact velocities (<$4\mathrm{m}/\mathrm{s}$)

(24)$S_p=-0.112u_n+0.99$

and the second one refers to normal impact velocities in the range (>$4 \mathrm{m}/\mathrm{s}$)

(25)$S_p=-6·{10}^{-5}{u}_n^2-6·{10}^{-4}{u}_n+0.545$

For the sticking decision, an auxiliary random number in the range $[0,1]$ is generated and compared with $S_p$. If greater, the algorithm rejects the sticking, and the restitution coefficients are used to calculate the new particle velocity. Otherwise, the particle sticks, and its velocity is set to zero. Two parameters can be defined in order to quantify the impact and adhesion characteristics. The impact efficiency (IE) has already been defined as the ratio between the number of impacting particles and the number of particles injected for a fixed diameter. The capture efficiency (CE) is defined as the ratio between the number of stuck particles and the number of particles injected for a fixed diameter. A time step of $1\times {10}^{-6} \mathrm{s}$ is used for the flow field reconstruction, corresponding to 100 instants per blade passing period. The computational time necessary to obtain the solution was about one order of magnitude higher than the steady-state case. This is considered as satisfactory, since fully unsteady calculations would have required the simulation of half annulus of the stage with an estimated computational time of two to three orders of magnitude higher than the steady-state solution.

The trends of the impact and capture efficiencies as a function of the Stokes number are shown in Figure 7 for the rotor and the stator. In the case of the rotor, only results for the choked flow condition are shown (Figure 7a) since the results were very similar in the near stall case. Results for both steady and unsteady calculations are shown. Steady results are discussed first (dashed lines). In general, it is possible to observe that the percentage of particles that hit the blade surface increases with the particle diameter for the rotor, while the stator has a maximum for $8 \mathsf{\mu }\mathrm{m}$ particles, more pronounced for the near stall condition. The number of particles that impact the blade and the vane is similar for small diameters. For large diameters, the blade has a greater number of impacts, and this difference increases with the particle diameter. As can be seen, smaller particles are more likely to stick when they make impact on the walls, but they undergo a lower number of impacts because they follow better the fluid streamlines thanks to the lower inertia. For larger particles, the opposite is true. The combination of these factors results in $8 \mathsf{\mu }\mathrm{m}$ particles having the highest capture efficiency (red line). The trends and the values are in line with those presented by Suman et al. [41] for the rotor and by Aldi et al. [13] for the entire stage. The most significant difference between the two operating conditions is represented by the vane impact efficiency. It can be observed that the near stall condition provides a higher number of impacts on the stator for the whole range of diameters, with $8 \mathsf{\mu }\mathrm{m}$ particles reaching a value almost double with respect to the choked condition. The reason for this behaviour can be explained in fluid dynamics terms, as described later in the text.

As expected, unsteady effects due to upstream running pressure disturbances caused by the vane are negligible, and the lines are almost overlapped. The situation for the near stall condition is equivalent and is not reported for the sake of brevity. Figure 7c refers instead to the stator. It can be observed that also for the vane the rotor–stator interaction effects are not extremely important. A slight increase in impact and capture efficiencies is observed in the choked case, but only of a few percentage points, with the maximum difference for $St=0.7$ ($8\mathsf{\mu }\mathrm{m}$ particles). For the near stall condition, results are almost unchanged, as can be noticed in Figure 7c. These findings are corroborated by the fact that the local impact pattern is also almost unchanged moving to the unsteady calculation, as testified by Figure 8, which compares the two calculations in terms of local impact efficiency on the pressure and suction sides of the blade and vane surfaces. Only $8 \mathsf{\mu }\mathrm{m}$ particles at a choked flow condition are shown since this size resulted in the largest discrepancy in the global IE. Nevertheless, only a small increase of impacts towards the trailing edge on the pressure side of the vane can be observed compared to the steady simulation. Therefore, it can be concluded that, for this geometry, transient effects do not produce important differences on particle behaviour and steady calculations are appropriate to capture the relevant features. This is an important remark that could not be argued a priori, and shows that, for a transonic axial compressor stage, the operating point is a more important parameter influencing deposition. Therefore, we investigate this aspect further.

The near stall condition results in general in higher impact efficiencies, especially for medium-sized particles. This behaviour can be explained by looking at Figure 9, which shows the trajectories inside the rotor and stator passages for the two operating conditions. A relative tangential velocity component is reported for the blade row, while absolute tangential velocity is shown for the vane row. The smallest particles are highly influenced by the flow field, and their velocity is very close to the local flow velocity (Figure 9a). Due to the different shock wave intensity, these particles experience a higher deceleration for the near stall condition and, in much the same way as the continuous phase, their incidence angle at the stator entry is increased (Figure 9a). However, these particles undergo a very small number of impacts on the blades, as reported in Figure 7 and, therefore, their impact efficiency is not significantly altered by the operating condition. On the contrary, larger particles are much less sensitive to the deceleration induced by the shock at the blade leading edge. Indeed, as illustrated in Figure 9c, trajectories are only marginally altered by the operating condition of the compressor. Nevertheless, the high impact efficiency amplifies these effects and results in higher local impacts on the rear part of the vane pressure side for the near stall condition. $8 \mathsf{\mu }\mathrm{m}$ particles are somewhat in between the two situations, since they are large enough to have a high impact efficiency on the blade, but their inertia still allows for experiencing the stronger deceleration caused by the stronger shock in the near stall condition. Indeed, looking at Figure 9b, it can be seen that a stronger deceleration is observed in the area enclosed inside the black ellipse, producing an increased incidence at the entrance of the stator. This is also evidenced by the wider particle “wake” observed in the rear part of the vane with respect to the choked condition. At this point, contrary to $0.5 \mathsf{\mu }\mathrm{m}$ particles, they are not able to follow the flow streamlines and the impact and capture efficiencies are significantly increased for the near stall condition. However, it is important to notice that the higher incidence angle increases the impacts on the pressure side, but in the leading edge area of the suction side, fewer impacts occur. For this reason, it can not be immediately concluded that the increased deposition for near stall condition will cause higher performance losses than the choked condition. Indeed, suction side and leading edge shape modifications are the most impacting on the losses inside the passage and could therefore be more detrimental to the global compressor performance [42,43]. Additional studies including blade surface roughness changes caused by deposition must be performed to assess the importance of these effects.

5.2. 1.5 Axial Turbine Stage

The second geometry analysed is the 1.5 axial turbine stage rig of the University of Bath [44]. The rig is composed of 32 upstream vanes, 48 rotor blades and 32 downstream vanes. A structured body-fitted grid of nearly 2.1 million elements was generated using ANSYS ICEM CFD. The mesh size is chosen according to a precedent study [44] showing that the cell number is appropriate for the simulation. The average $y^{+}$ value on the blades is about 14, allowing the usage of wall-functions.

The boundary conditions for the calculation are shown in Table 2. The flow coefficient is defined as $C_F=\frac{u_x}{{\mathbf{\Omega }}_{d}b}$, where $u_x$ is the axial velocity at the inlet section, ${\mathbf{\Omega }}_d$ is the rotor angular speed, and b is the turbine hub radius.

The timestep of the HBM-reconstruction calculation was chosen so that 100 instants per blade passing period were resolved. Figure 10a illustrates the time-averaged pressure field on stator and rotor walls and endwalls. Figure 10b shows the validation of the solution, in terms of the comparison between the pressure coefficient for numerical and experimental data. Illustrated is the time-averaged value of $C_p$ over two passages on a probe located on the hub downstream of the first vane.

The same structure of the previous section is followed here, by presenting a steady analysis to highlight the salient features of the discrete phase behaviour, followed by a comparison with the HBM-reconstruction technique. Spherical drag and Saffman–Mei lift force acting on the particles are considered during the calculation. A constant density of $2700 \mathrm{k}\mathrm{g}/{\mathrm{m}}^3$ is assumed for the discrete phase. Particles are injected from the domain inlet at the same velocity of the air flow. The variation of the particle diameter $d_p$ is in the range $0.5÷32 \mathsf{\mu }\mathrm{m}$, while the Stokes number, calculated at the three rows inlet sections, is reported in Table 3.

A mixing plane interface is placed between stators and the rotor for the steady analysis. The global impact and capture efficiencies for the first stator, the rotor and the second stator as a function of particle diameter are reported in Figure 11. It can be noticed that, for the first stator, these properties grow with particle size. This trend resembles the one seen for the rotor in the NASA stage 37 case (Figure 7a). Consequently, it can be deduced that, for the first upstream row of other axial turbomachinery, the trend will be roughly the same. In this case, however, the impact number is higher due to larger turning angle imposed to the flow by the vane row. The monotonic trend is interrupted after $d_p>8 \mathsf{\mu }\mathrm{m}$ because, for this size, all particles have already impinged on the vane surface. Indeed, it can be noticed that, for larger particles, the impact efficiency is higher than $100\%$, meaning that some particles rebound multiple times on the vane. On the contrary, the diagrams for the rotor and the second vane row present a maximum for middle-sized particles. For larger particles, the first stage acts as a “filter” with respect to the second stator. Indeed, despite the fact that the sticking probability is significantly lower, the larger inertia causes a high impact and capture efficiencies in the first stage, with the result that only a small fraction of large particles can reach the downstream vane row. This is testified by the fact that the impact and capture efficiencies of 8 and $32 \mathsf{\mu }\mathrm{m}$ particles for the second stator are orders of magnitude lower than for the other rows. This is not true for smaller particles since only a fraction of them actually make impact on the walls. To a lower extent, the first vane row acts as a filter for the rotor, explaining the decreasing number of impacts for large particles.

To better analyse these features, trajectories for a hundred particles are shown in a meridional plane for three diameters in Figure 12, coloured by absolute velocity magnitude. $1 \mathsf{\mu }\mathrm{m}$ particles undergo a higher acceleration in the first vane row, but their inertia is not sufficiently high and drag force is prevalent along their trajectory. When the Stokes number grows, particles are accelerated less by the flow, and they approach the following row with a lower incidence. Moreover, their inertia begins to be important and the average radial impact locations on the blade is shifted towards the shroud of the turbine. This is in agreement with the findings of Tabakoff et al. [8] who analysed the impact locations for particles of similar size on an axial turbine. Some of them are also able to move through the tip gap of the blade and exit on the suction side. It can be seen that $32 \mathsf{\mu }\mathrm{m}$ particles are centrifuged out so fast that almost no impacts occur on the pressure side of the blade; instead, they impact on the rotor case. As previously described, for larger particles, the almost totality of particles are not able to reach the downstream stator. However, it is important to note that, although Equations (24) and (25) predict a low sticking probability for large particles (due to the higher normal impact velocity), this value could be overestimated compared to other models typically used in the literature. This results in a higher number of particles being able to reach the downstream vane row in the end. It can be also noticed that the particles are not able to infiltrate inside the cavity, and their trajectory is not affected by the purge flow exiting from the rim gap, even for $1 \mathsf{\mu }\mathrm{m}$. The same simulation (not shown here) was also conducted for the high-purge flow condition, showing almost identical results.

Figure 11 shows the impact and capture efficiencies for the three rows. The only noteworthy differences between steady and transient simulations are due to smaller particles’ impacts. Up to $2 \mathsf{\mu }\mathrm{m}$, the impacts on the blade pressure side are increased and more uniformly distributed, while for 8 and $32 \mathsf{\mu }\mathrm{m}$, they are almost unmodified. Unsteady effects are negligible on the first vane row. This can be demonstrated also by comparing the global impact and capture efficiency for the upstream stator obtained with steady and HBM simulations (Figure 11a). As expected, the lines are superimposed to each other. For what concerns the rotor blade, Figure 11b shows that some differences are present for the transient solution. The most important are the total impact and capture efficiencies for $2 \mathsf{\mu }\mathrm{m}$ particles that are increased from $30.5\%$ to $40\%$ and from $24.5\%$ to $30.5\%$, respectively. Once again, rotor–stator interactions are more important for small particles due to the lower inertia, but it is the combination of this effect with the number of impacts that must be accounted for.

The impact efficiencies for $1 \mathsf{\mu }\mathrm{m}$ and $2 \mathsf{\mu }\mathrm{m}$ are shown in Figure 13. These are the particle sizes that resulted in the larger differences in the impact and capture efficiencies between the steady and transient calculations. Pressure and suction sides are reported for the three rows. The pressure side is the most struck by the discrete phase. It can be also noticed that a small percentage of $1 \mathsf{\mu }\mathrm{m}$ particles impact the suction side of the rotor near the tip. This is due to the tip leakage flow that transports some particles across the tip gap towards the suction side of the blade.

Unsteady effects could be potentially more important in the downstream stator, but the number of impacts is so small that no considerable percentage differences are observed, as illustrated in Figure 11c. As for the NASA stage 37, the very similar patterns found that comparing steady-state and transient approaches allows for arguing that steady-state solutions are a good approximation of the discrete phase physics, at least for the tested geometries and particle sizes. Unfortunately, this is not of much interest from the scientific standpoint, but allows for putting a higher trust in the results obtained with stationary simulations of discrete phases, which are still the most commonly employed in this research field. Please notice that under no circumstances should this remark be directly extended to the continuous phase, as already mentioned several times. Anyway, it is shown that unsteady effects on the discrete phase are small but not negligible for the rotor blade. Differences up to $10\%$ are present for $St\approx 0.1$ as shown in Figure 11b.

Unsteady analysis is also very useful to visualize the dynamic of particles approaching the rotating passage. To explore this feature, a small stream of particles exiting from the upstream vane is selected. Figure 14 shows the locations of 1 and $8 \mathsf{\mu }\mathrm{m}$ particles at four different time steps, coloured by absolute velocity magnitude. In the $1 \mathsf{\mu }\mathrm{m}$ case, particles are accelerated as they leave the vane passage and travel faster than the rotor blade, avoiding suction surface impacts. Larger particles instead experience a smaller acceleration and tend to concentrate in the wake region after impacting the vane walls. Differently from smaller particles, they travel slower and are chopped by the oncoming rotor blades. Having a smaller sticking probability, they rebound on the rotor blade suction surface as the blade keeps moving toward them.

6. Conclusions

This work consisted of analyzing the behaviour of a discrete phase inside gas turbine components. A technique that combines HBM with Lagrangian tracking was presented, allowing for significant speedup during the continuous and discrete phase solution for unsteady calculations. The matter of particles crossing cyclic boundaries with a phase-shift relation was properly addressed. Moving to the results, firstly, a simple test-case has been used to assess the correct implementation of the methodology. Then, two representative applications of the Lagrangian tracking algorithm to gas turbine components were presented, namely a transonic axial compressor stage and a 1.5 axial turbine stage. A deposition model from the literature was used to determine which particles tend to stick on the blade surfaces. Impact statistics and particle trajectories were used to explain the observed trends in global impact and capture efficiencies. For the compressor stage, two operating conditions were investigated. Steady-state calculations with mixing plane interfaces and unsteady simulations with the HBM-reconstruction were performed and compared to each other for both case studies. Generally speaking, significantly different behaviours can be expected for particles with different Stokes number. For a transonic compressor stage, it was demonstrated that different shock intensities caused by different operating conditions have important effects on the vane impact efficiency according to particle size. For the axial turbine, it was demonstrated that impacts and deposition occur almost exclusively in the first stage that acts as a filter for the successive rows. Additionally, due to the significantly higher swirl component of the flow compared to the axial compressor, larger particles are centrifuged towards the casing of the turbine much faster. Interestingly, unsteady effects seem to affect the particle impacts only slightly. The most important differences are observed for the axial turbine case, in which the impact efficiency for $2 \mathsf{\mu }\mathrm{m}$ particles is increased by about $10\%$ on the rotor blade. However, these differences do not dramatically modify the conclusions drawn with the steady models. This is perhaps disappointing from a scientific standpoint but confirms that steady-state solutions are a good approximation of the discrete-phase physics at a reduced computational time, at least for the tested geometries and particle sizes. In any case, the HBM-reconstruction technique allows for obtaining more accurate results at an acceptable computational cost and therefore represents a valuable tool for the prediction of the discrete-phase behaviour.

Author Contributions

Software development and CFD simulations, S.O.; resources, review and supervision, N.C.; resources, review and supervision, M.P.; resources and review, M.C. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

Data presented in the paper can be obtained upon reasonable request by contacting the authors.

Conflicts of Interest

The authors declare no conflict of interest.

Nomenclature

The following abbreviations are used in this manuscript:

$\mathbf{\Omega }$	Domain angular velocity vector
$\mathit{F}$	Vector of fluxes
$\mathit{M}$	Matrix of the equations system
$\mathit{Q}$	Vector of conservative variables
$\mathit{R}\left(\mathit{Q}\right)$	Discretised fluxes residual term
$\mathit{u}$	Velocity
${\widehat{Q}}_m$	m-th Fourier coefficient of Q
${\mathit{D}}_{\mathit{t}}$	Time derivative operator
$\mathit{D}$	Harmonic balance operator
$\mathit{E}$	Discrete Fourier transform matrix
$\omega$	Angular frequency corresponding to a known flow periodicity
$\sigma$	Interblade phase angle
$\theta$	Azimuthal coordinate
$N_B$	Number of blades in a row
$N_T$	Number of snapshots for the HBM simulation
V	Control volume

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. Schematic procedure for particle transfer between phase-shifted interfaces. Particle trajectory is represented with a dashed blue line.

Figure 2. (a) Temperature field in the blade passages with unsteady wake imposed at the inlet; (b) temperature profiles in front of two blades obtained with the fully transient and HBM calculations. Reference pitch is twice the single blade pitch.

Figure 3. Comparison of [Forumla omitted. See PDF.] particle trajectories for transient simulation and HBM-reconstruction method at four equispaced instants.

Figure 4. (a,c,e) Impact efficiency on the blade surface for the steady, transient, and HBM-reconstruction methods; (b,d,f) impact efficiency for different number of harmonics in the HBM-reconstruction method. Results are shown for the various Stokes numbers: [Forumla omitted. See PDF.] (a,b), [Forumla omitted. See PDF.] (c,d) and [Forumla omitted. See PDF.] (e,f).

Figure 4. (a,c,e) Impact efficiency on the blade surface for the steady, transient, and HBM-reconstruction methods; (b,d,f) impact efficiency for different number of harmonics in the HBM-reconstruction method. Results are shown for the various Stokes numbers: [Forumla omitted. See PDF.] (a,b), [Forumla omitted. See PDF.] (c,d) and [Forumla omitted. See PDF.] (e,f).

Figure 5. Normalised CPU times (reference is steady-state calculation) for the [Forumla omitted. See PDF.] rotor case.

Figure 6. Contours of pressure and entropy at midspan for the NASA stage 37.

Figure 7. Comparison of efficiencies for the NASA stage 37 obtained with steady and HBM-reconstruction methods.

Figure 8. NASA stage 37 at choked flow condition: impact efficiency for [Forumla omitted. See PDF.] on the blade and vane surfaces obtained with steady and HBM-reconstruction methods.

Figure 9. NASA stage 37: particles trajectories for the steady calculations at choked and near stall conditions.

Figure 10. (a) Time -averaged pressure contours on the wall of the 1.5 axial turbine stage; (b) pitchwise distribution of time-averaged [Forumla omitted. See PDF.]: experimental and numerical results.

Figure 11. Comparison of efficiencies for the 1.5 axial turbine stage obtained with steady and HBM-reconstruction methods.

Figure 12. 1.5 axial turbine stage: particle trajectories for three different diameters.

Figure 13. Impact efficiencies comparison between steady and HBM simulations for [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.] in the three rows of the 1.5 axial turbine stage.

Figure 14. 1.5 axial turbine stage: particle positions at four equispaced time steps for two particle diameters. Particles are coloured by absolute velocity magnitude.

Table 1

Injection data for the NASA stage 37. $N_{inj}$ is the number of particles injected for each diameter.

${\mathit{d}}_{\mathit{p}}$ [$\mathsf{\mu }$m]	Rotor
	$0.5$	1	2	8	32
$St$	$2\times {10}^3$	$8\times {10}^{-3}$	$3\times {10}^{-2}$	0.5	8
${\tau }^{+}$	25	98	394	$6.3\times {10}^3$	$1.0\times {10}^5$
$N_{inj}$	$1\times {10}^6$	$1\times {10}^6$	$5\times {10}^5$	$1\times {10}^5$	$1\times {10}^5$
${\mathit{d}}_{\mathit{p}}$ [ $\mathsf{\mu }$ m]	Stator
	$\mathbf{0}.\mathbf{5}$	$\mathbf{1}$	$\mathbf{2}$	$\mathbf{8}$	$\mathbf{32}$
$St$	$3\times {10}^3$	$1\times {10}^{-2}$	$4\times {10}^{-2}$	0.7	11
${\tau }^{+}$	19	77	309	$4.9\times {10}^3$	$7.9\times {10}^4$
$N_{inj}$	$1\times {10}^6$	$1\times {10}^6$	$5\times {10}^5$	$1\times {10}^5$	$1\times {10}^5$

Table 2

Operating conditions for the 1.5 axial turbine stage.

Parameter	Value
Flow coefficient, $C_F$	0.41
Inlet total temperature	300 K
Outlet static pressure	99,300 Pa
Disc rotational speed	3000 rpm
Inlet turbulent intensity k	3%
Inlet turbulent length scale	$3\times {10}^{-3}$ $\mathrm{m}$

Table 3

Stokes numbers for the 1.5 axial turbine stage. $N_{inj}$ is the number of particles injected for each diameter.

${\mathit{d}}_{\mathit{p}}$ [$\mathsf{\mu }$m]	${\mathit{N}}_{\mathit{inj}}$	$\mathit{St}$ Stator1	$\mathit{St}$ Rotor	$\mathit{St}$ Stator2
$0.5$	$1\times {10}^6$	$1.7\times {10}^{-3}$	$7.2\times {10}^{-3}$	$1.8\times {10}^{-3}$
1	$1\times {10}^6$	$6.8\times {10}^{-3}$	$2.9\times {10}^{-2}$	$7.3\times {10}^{-3}$
2	$5\times {10}^5$	$2.7\times {10}^{-2}$	0.12	$2.9\times {10}^{-2}$
8	$1\times {10}^5$	0.44	1.85	0.47
32	$1\times {10}^5$	6.99	29.5	7.47

References

1. Wilcox, M.; Kurz, R.; Brun, K. Technology Review of Modern Gas Turbine Inlet Filtration Systems. Int. J. Rotating Mach.; 2012; 2012, 128134. [DOI: https://dx.doi.org/10.1155/2012/128134]

2. Dunn, M.G. Operation of gas turbine engines in an environment contaminated with volcanic ash. J. Turbomach.; 2012; 134, 051001. [DOI: https://dx.doi.org/10.1115/1.4006236]

3. Suman, A.; Casari, N.; Fabbri, E.; Pinelli, M.; Di Mare, L.; Montomoli, F. Gas Turbine Fouling Tests: Review, Critical Analysis, and Particle Impact Behavior Map. J. Eng. Gas Turbines Power; 2019; 141, 032601. [DOI: https://dx.doi.org/10.1115/1.4041282]

4. Suman, A.; Casari, N.; Fabbri, E.; Di Mare, L.; Montomoli, F.; Pinelli, M. Generalization of particle impact behavior in gas turbine via non-dimensional grouping. Prog. Energy Combust. Sci.; 2019; 74, pp. 103-151. [DOI: https://dx.doi.org/10.1016/j.pecs.2019.05.001]

5. Hamed, A.; Tabakoff, W.; Wenglarz, R. Erosion and Deposition in Turbomachinery. J. Propuls. Power; 2006; 22, pp. 350-360. [DOI: https://dx.doi.org/10.2514/1.18462]

6. Ghenaiet, A. Prediction of erosion in an axial turbine with initial position of blade. Proceedings of the 13th European Turbomachinery Conference, ETC2019-111; Lausanne, Switzerland, 8–12 April 2019.

7. Hamed, A.A.; Tabakoff, W.; Rivir, R.B.; Das, K.; Arora, P. Turbine Blade Surface Deterioration by Erosion. J. Turbomach.; 2004; 127, pp. 445-452. [DOI: https://dx.doi.org/10.1115/1.1860376]

8. Tabakoff, W.; Hamed, A.; Metwally, M. Effect of Particle Size Distribution on Particle Dynamics and Blade Erosion in Axial Flow Turbines. J. Eng. Gas Turbines Power; 1991; 113, pp. 607-615. [DOI: https://dx.doi.org/10.1115/1.2906284]

9. Mustafa, Z.; Pilidis, P.; Amaral Teixeira, J.; Ahmad, K. CFD Aerodynamic Investigation of Air-Water Trajectories on Rotor-Stator Blade of an Axial Compressor for Online Washing. Proceedings of the Turbo Expo: Power for Land, Sea, and Air; Barcelona, Spain, 8–11 May 2006; pp. 1385-1394.

10. Yang, H.; Boulanger, J. The Whole Annulus Computations of Particulate Flow and Erosion in an Axial Fan. J. Turbomach.; 2012; 135, 011040. [DOI: https://dx.doi.org/10.1115/1.4006564]

11. Zagnoli, D.; Prenter, R.; Ameri, A.; Bons, J.P. Numerical Study of Deposition in a Full Turbine Stage Using Steady and Unsteady Methods. Proceedings of the ASME Turbo Expo 2015: Turbomachinery Technical Conference and Exposition; Montreal, QC, Canada, 15–19 June 2015; Paper No: GT2015-43613

12. Prenter, R.; Ameri, A.; Bons, J.P. Computational Simulation of Deposition in a Cooled High-Pressure Turbine Stage with Hot Streaks. ASME J. Turbomach.; 2017; 139, 091005. [DOI: https://dx.doi.org/10.1115/1.4036008]

13. Aldi, N.; Morini, M.; Pinelli, M.; Spina, P.R.; Suman, A. An Innovative Method for the Evaluation of Particle Deposition Accounting for Rotor/Stator Interaction. J. Eng. Gas Turbines Power; 2016; 139, 052401. [DOI: https://dx.doi.org/10.1115/1.4034968]

14. Oliani, S.; Friso, R.; Casari, N.; Pinelli, M.; Suman, A.; Carnevale, M. Progresses in Particle-Laden Flows Simulations in Multistage Turbomachinery with OpenFOAM. J. Turbomach.; 2022; 144, 101007. [DOI: https://dx.doi.org/10.1115/1.4054076]

15. Hall, K.C.; Thomas, J.P.; Clark, W.S. Computation of Unsteady Nonlinear Flows in Cascades Using a Harmonic Balance Technique. AIAA J.; 2002; 40, pp. 879-886. [DOI: https://dx.doi.org/10.2514/2.1754]

16. Oliani, S.; Casari, N.; Carnevale, M. ICSFoam: An OpenFOAM library for implicit coupled simulations of high-speed flows. Comput. Phys. Commun.; 2023; 286, 108673. [DOI: https://dx.doi.org/10.1016/j.cpc.2023.108673]

17. Saad, Y.; Schultz, M.H. GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems. Siam J. Sci. Stat. Comput.; 1986; 7, pp. 856-869. [DOI: https://dx.doi.org/10.1137/0907058]

18. Yoon, S.; Jameson, A. Lower-upper Symmetric-Gauss-Seidel method for the Euler and Navier–Stokes equations. AIAA J.; 1988; 26, pp. 1025-1026. [DOI: https://dx.doi.org/10.2514/3.10007]

19. Macpherson, G.B.; Nordin, N.; Weller, H.G. Particle tracking in unstructured, arbitrary polyhedral meshes for use in CFD and molecular dynamics. Commun. Numer. Methods Eng.; 2009; 25, pp. 263-273. [DOI: https://dx.doi.org/10.1002/cnm.1128]

20. Oliani, S.; Casari, N.; Pinelli, M.; Suman, A.; Carnevale, M. Numerical study of a centrifugal pump using Harmonic Balance Method in OpenFOAM. Proceedings of the Taormina, Italy: The 34th International Conference on Efficiency, Cost, Optimization, Simulation and Environmental Impact of Energy Systems; Taormina, Italy, 27 June–2 July 2021.

21. Oliani, S.; Casari, N.; Carnevale, M. A New Framework for the Harmonic Balance Method in OpenFOAM. Machines; 2022; 10, 279. [DOI: https://dx.doi.org/10.3390/machines10040279]

22. Guédeney, T.; Gomar, A.; Gallard, F.; Sicot, F.; Dufour, G.; Puigt, G. Non-uniform time sampling for multiple-frequency harmonic balance computations. J. Comput. Phys.; 2013; 236, pp. 317-345. [DOI: https://dx.doi.org/10.1016/j.jcp.2012.11.010]

23. Gopinath, A.; van der Weide, E.; Alonso, J.; Jameson, A.; Ekici, K.; Hall, K. Three-Dimensional Unsteady Multi-stage Turbomachinery Simulations Using the Harmonic Balance Technique. Proceedings of the 45th AIAA Aerospace Sciences Meeting and Exhibit; Reno, NV, USA, 8–11 January 2007.

24. Nimmagadda, S.; Economon, T.D.; Alonso, J.J.; da Silva, C.R.I. Robust uniform time sampling approach for the harmonic balance method. Proceedings of the 46th AIAA Fluid Dynamics Conference; Washington, DC, USA, 13–17 June 2016.

25. Tyler, J.M.; Sofrin, T.G. Axial Flow Compressor Noise Studies. Proceedings of the SAE Technical Paper; SAE International: Warrendale, PA, USA, 1962.

26. Sicot, F.; Dufour, G.; Gourdain, N. A time-domain harmonic balance method for rotor/stator interactions. J. Turbomach.; 2012; 134, [DOI: https://dx.doi.org/10.1115/1.4003210]

27. Ekici, K.; Hall, K.C. Nonlinear Analysis of Unsteady Flows in Multistage Turbomachines Using Harmonic Balance. AIAA J.; 2007; 45, pp. 1047-1057. [DOI: https://dx.doi.org/10.2514/1.22888]

28. Crespo, J.; Contreras, J. On the Development of a Synchronized Harmonic Balance Method for Multiple Frequencies and its Application to LPT Flows. Proceedings of the ASME Turbo Expo 2020: Turbomachinery Technical Conference and Exposition; Virtual, 21–25 September 2020; Paper No: GT2020-14952

29. Erdos, J.I.; Alzner, E.; McNally, W. Numerical Solution of Periodic Transonic Flow through a Fan Stage. AIAA J.; 1977; 15, pp. 1559-1568. [DOI: https://dx.doi.org/10.2514/3.60823]

30. Gerolymos, G.A.; Michon, G.J.; Neubauer, J. Analysis and Application of Chorochronic Periodicity in Turbomachinery Rotor/Stator Interaction Computations. J. Propuls. Power; 2002; 18, pp. 1139-1152. [DOI: https://dx.doi.org/10.2514/2.6065]

31. Oliani, S.; Friso, R.; Casari, N.; Pinelli, M.; Suman, A.; Carnevale, M. A comparative analysis of particle-mixing plane interaction in multistage turbomachinery simulations. Proceedings of the European Turbomachinery Conference 14; Gdansk, Poland, 12–16 April 2021.

32. Thulin, R.D.; Howe, D.C.; Singer, I.D. High-Pressure Turbine Detailed Design Report; Technical Report CR-165608 NASA Report Lewis Research Center: Cleveland, OH, USA, 1982.

33. Launder, B.; Spalding, D. The numerical computation of turbulent flows. Comput. Methods Appl. Mech. Eng.; 1974; 3, pp. 269-289. [DOI: https://dx.doi.org/10.1016/0045-7825(74)90029-2]

34. Jacobs, G.; Don, W.S.; Dittmann, T. High-order resolution Eulerian–Lagrangian simulations of particle dispersion in the accelerated flow behind a moving shock. Theor. Comput. Fluid Dyn.; 2012; 26, pp. 37-50. [DOI: https://dx.doi.org/10.1007/s00162-010-0214-6]

35. Reid, L.; Moore, R.D. Design and Overall Performance of Four Highly Loaded, High Speed Inlet Stages for an Advanced High-Pressure-Ratio Core Compressor; Technical Report 1337 NASA Technical Paper Lewis Research Center: Cleveland, OH, USA, 1978.

36. Mei, R.; Klausner, J. Shear lift force on spherical bubbles. Int. J. Heat Fluid Flow; 1994; 15, pp. 62-65. [DOI: https://dx.doi.org/10.1016/0142-727X(94)90031-0]

37. Kurz, R.; Brun, K. Degradation in Gas Turbine Systems. J. Eng. Gas Turbines Power; 2000; 123, pp. 70-77. [DOI: https://dx.doi.org/10.1115/1.1340629]

38. Tian, L.; Ahmadi, G. Particle deposition in turbulent duct flows–comparisons of different model predictions. J. Aerosol Sci.; 2007; 38, pp. 377-397. [DOI: https://dx.doi.org/10.1016/j.jaerosci.2006.12.003]

39. Forder, A.; Thew, M.; Harrison, D. A numerical investigation of solid particle erosion experienced within oilfield control valves. Wear; 1998; 216, pp. 184-193. [DOI: https://dx.doi.org/10.1016/S0043-1648(97)00217-2]

40. Poppe, T.; Blum, J.; Henning, T. Analogous Experiments on the Stickiness of Micron-sized Preplanetary Dust. Astrophys. J.; 2000; 533, pp. 454-471. [DOI: https://dx.doi.org/10.1086/308626]

41. Suman, A.; Morini, M.; Kurz, R.; Aldi, N.; Brun, K.; Pinelli, M.; Ruggero Spina, P. Quantitative Computational Fluid Dynamic Analyses of Particle Deposition on a Transonic Axial Compressor Blade—Part II: Impact Kinematics and Particle Sticking Analysis. J. Turbomach.; 2014; 137, 021010. [DOI: https://dx.doi.org/10.1115/1.4028296]

42. Morini, M.; Pinelli, M.; Spina, P.R.; Venturini, M. Numerical Analysis of the Effects of Nonuniform Surface Roughness on Compressor Stage Performance. J. Eng. Gas Turbines Power; 2011; 133, 072402. [DOI: https://dx.doi.org/10.1115/1.4002350]

43. Kurz, R.; Brun, K. Fouling Mechanisms in Axial Compressors. J. Eng. Gas Turbines Power; 2012; 134, 032401. [DOI: https://dx.doi.org/10.1115/1.4004403]

44. Cosmo, G.D.; Scobie, J.A.; Lock, G.D.; Sangan, C.M.; Carnevale, M. Fluid Dynamics of Turbine Rim Seal Structures: A Physical Interpretation Using URANS. J. Eng. Gas Turbines Power; 2022; 145, 031009. [DOI: https://dx.doi.org/10.1115/1.4055752]