A strategy for the robust forecasting of gas turbine health subjected to fouling

Fouling represents a major problem for Gas Turbines (GTs) in both heavy-duty and aeropropulsion applications. Solid particles entering the engine can stick to the internal surfaces and form deposits. Components' lifetime and performance can dramatically vary as a consequence of this phenomenon. These effects impact the whole engine in terms of residual life, operating stability, and maintenance costs. In the High-Pressure Turbine (HPT), in particular, the high temperatures soft the particles and promote their adhesion, especially in the short term. Unfortunately, predicting the GT response to this detrimental issue is still an open problem for scientists. Furthermore, the stochastic variations of the components operating conditions increase the uncertainty of the forecasting results. In this work, a strategy to predict the effects of turbine fouling on the whole engine is proposed. A stationary Gas Path Analysis (GPA) has been performed for this scope to predict the GT health parameters. Their alteration as a consequence of fouling has been evaluated by scaling the turbine map. The scaling factor has been found by performing Computational Fluid Dynamic (CFD) simulations of a HPT nozzle with particle injection. Being its operating conditions strongly uncertain, a stochastic analysis has been conducted. The uncertainty sources considered are the circumferential hot core location and the turbulence level at the inlet. The study enables to build of confidence intervals on the GT health parameters predictions and represents a step forward towards a robust forecasting tool.


INTRODUCTION
Computer simulations are a powerful tool for the performance prediction of gas turbines. In many applications, a high fidelity computer simulation can substitute the real engine [1]. For example, in aeronautical applications, it is possible to simulate a complete mission in order to optimize the maintenance and the mission profile for the engine itself [2]. Instead, for heavy-duty applications, the load history (energy demand) can be analyzed [3]. In these simulations, the engine can be modelled at various levels of detail, from the full 3D description by using the computational fluid dynamics (CFD) to the thermodynamic description (0D) using the gas path analysis (GPA). One of the most useful applications of these tools regards the prediction of performance degradation of the engine during operations [3]. In this track, it is well recognized that the most critical cause of performance degradation is particle ingestion. Particles that enter the engine can deposit or erode its internal surfaces, increasing roughness and changing the aerodynamics geometries [4]. Even though this occurrence can affect either the cold and hot section, for aeronautical applications the latter results to be more susceptible to the phenomenon in the short term [5]. In particular, the most critical component is the high-pressure turbine vane, mainly due to the high temperature that occurs during operations.
In the literature, several studies have tried to analyze the effects of particle ingestion by using either GPA and CFD. For what concern the GPA, the only way to take into account degradation due to particles deposition and/or erosion is by imposing a degradation rate function [6]. Typically, this function is estimated by using data-driven techniques [6]. The main drawback of this approach is that the data-driven degradation rate lacks physical meaning. This may lead to unrealistic and un-physical results [7]. An example of this type of analysis is the work of Hanaki et al [6]. They estimated two correction factors for the compressor corrected mass flow rate and efficiency. These were used to calculate the modified compressor maps when degradation occurs. Another example is the work of Igie et al [8], where a reduction of fan capacity and efficiency were enforced to analyze the engine response during a complete mission. On the other hand, with the CFD a more realistic and physical prediction of performance degradation can be carried out. Here, the main issue is related to the use of deposition/erosion models. Even though many of these models are proposed in the literature [9], [10], [11], the physics underneath the phenomena still not completely understand. Works concerning particle deposition and erosion in HPT are mainly centred on the influence on film cooling and geometry variations. The first was deeply analyzed in the study of Borello et al [12]. They asserted that film cooling greatly affects deposition, mainly due to the effects of the jets on the particles' trajectories. The second theme was faced in the work of Casari et al [9]. They analyzed the particle deposition on the LS89 HPT vane, concluding that the deposit and the flow field are strictly coupled. Specifically, shifting of the shock and perturbation of pressure coefficients are the main consequences. More recently, Brandes et al [13] proposed an approach to estimate the severity of an engine mission due to erosion by combining computational and experimental results. They conclude that the taxi, take-off, and climb are the flight segments that impact the mission severity the most.
In this work, the authors tried to coupled CFD results with GPA to evaluate an aero-engine response to HPT fouling. The GT component treated as a reference is the HPT nozzle of the energy-efficient engine tested by Timko [14]. Furthermore, since several types of particles can enter the engine, the specific event regarding flight through the volcanic cloud was considered. The degradation rate was estimated by using a CFD analysis with particle ingestion on an HPT vane. Since the operating conditions at the exit of the combustor are known to be highly uncertain, a stochastic analysis was performed to have a probabilistic distribution of the output. This uncertain distribution was used to modify the turbine performance maps to fed the GPA.

ENGINE DESCRIPTION
The aero-engine architecture chosen for the study is the twin-spool high-bypass unmixed flow turbofan, which is considered the best configuration for high subsonic commercial aircraft [15]. Its section view is reported in Fig. 1. In turbofan engines, a fraction of the total flow bypasses the core components (compressor, combustion chamber and turbine) before being ejected through a separate nozzle. Thus, the total thrust is the sum of two factors: the cold thrust, due to the fan, and the hot thrust, which results from the stream entering the engine core.
Since these engines aim to generate the desired propulsive thrust, the typical performance parameter chosen as representative is the Thrust Specific Fuel Consumption (SFC). It is mathematically defined as the engine fuel mass flow rate ratio to the amount of thrust developed. Thus, it represents the mass of fuel needed to provide the required thrust. Since the purpose of this study is to propose a general strategy, the reference turbofan engine chosen is the default one present in the commercial software GasTurb9. This choice is mainly guided by the difficulty to find data of real engines in the literature.

METHODOLOGY
The aim of this work is to propose a strategy for the robust forecasting of gas turbine performance degradation due to fouling. Specifically, the engine response during a flight across a volcanic cloud is considered. Ash clouds can carry particulate with concentrations as high as 250 mg/m 3 , which is considered a dangerous value for aircraft engines [16]. The flowchart of the proposed methodology is depicted in Fig. 2. As stated above, the component most influenced by this event is the first HPT vane. To compute the effects of volcanic ash deposition on this element, a 2D midspan CFD simulation was conducted. The results were then analyzed employing uncertainty quantification techniques to include the statistical variation of the operating conditions ( and ). The outcomes thus obtained were used to identify the rate of change of the HPT maps. In this respect, the scaling factor approach (rigid shifting) of the maps was used, even though several works consider it a strong hypothesis for degradation study purposes [17]. Since this paper aims to propose a strategy for robust forecasting, this assumption does not affect the generality of the approach. Finally, the engine response was evaluated by using the GPA with the modified maps. A detailed description of each step is reported in the following sections.

CFD simulations
The numerical flow field resolution was carried out using the sonicFoam solver of the open-source software OpenFoam-v1706. The solver has been added the support for particle tracking, allowing the exchange of momentum and heat with the carrier flow. Particle tracking was performed by using the Lagrangian-Eulerian coupling method. Turbulent fluctuations on particles were computed with the discrete random walk approach.

Geometry and computational domain
For this purpose, the HPT of the General Electric Energy Efficient Engine (EEE) was chosen as a reference [14]. Specifically, the midspan cross-section of the first vane was considered. Furthermore, since this study is a first attempt to tackle the robust forecasting in gas turbines, a simplified geometry was examined. Therefore, the cooling system was not taken into account. The computational domain was generated around the 2D geometry using Salome-v9.3. A 2D unstructured grid was used, with a 10-layers prism in the near-wall region. A mesh composed of 80,000 cells was chosen after a grid sensitivity analysis.

Continuous phase
The operating conditions considered for the study reproduces the design point described by the technical report written by Timko [14] and summarized in Tab. 1. The profile of the normalized temperature used for the study is the same reported in the work of Friso et al [18]. Since a high level of turbulence intensity is expected at the outlet of the combustion chamber, the Reynolds stress model proposed by Speziale et al [19] was used to resolve the turbulent flow fluctuations. This choice was guided considering the questionable accuracy that arises from the use of two-equations models when high turbulence occurs [20].

Discrete phase
Once the flow field was solved, solid particles were seeded from the inlet of the domain with a uniform distribution between 1 μm and 30 μm. These values are specific for volcanic ash clouds and can be found in the work of Talvalul et al [21]. In deposition studies, particle properties are of main importance. In this work, the values proposed by Ghosal et al [22] were considered. In Tab. 1 is reported the amount of mass of ash injected, which was derived considering a concentration of contaminant equals 250 mg/m 3 , representing a very dense ash cloud [16].
Particles trajectories were computed by integrating the force balance acting on each of them. Specifically, the Basset Boussinesq Oseen (BBO) equation with drag as the only force acting was used [23]. Furthermore, since the maximum particles volume fraction of 10 is small enough to consider negligible the effect of particles on turbulence, a one-way coupling model was considered between particles and fluid flow [24]. As already mentioned, an important parameter governing the adhesion is the surface temperature. For this reason, the heat transfer between the gas and the particles is also computed by means of the well-known Ranz-Marshall correlation.
Each particle is tracked through the domain until it impacts a wall or escapes from the domain. The particle-wall interaction was modelled by using the EBFOG [9]. It acts as computing a sticking probability at every impact, then a stochastic algorithm decides if the particle sticks or erode the impacted surface. Furthermore, the geometric variations due to this occurrence were counted through a mesh morphing technique. Faces of the vane grid where the impact take place are displaced outward normally to the wall surface. A subsequent Laplace smoothing was performed to preserve the mesh quality.

Uncertainty analysis
The stochastic analysis has been carried out considering two input uncertainties: the circumferential hot core location and the turbulence intensity at the inlet. Unfortunately, no data were available in the literature for the reference geometry considered to generate the statistical distributions of the inputs. This forced the authors to use the statistics already present in the literature for the same engine architecture chosen. Specifically, the statistics of the hot core location and the turbulence level measured by Montomoli et al [25] and reported in Tab. 2 was used.  Once the stochastic inputs were determined, the algorithm to perform the uncertainty quantification has to be established. Among the many techniques present in the literature, the probabilistic collocation method (PCM) was chosen. This was mainly guided by the necessity of cutting down the number of CFD simulations to obtain results as reliable as possible. The PCM is a technique that belongs to the family of polynomial chaos expansion (PCE) methods [26]. It consists of describing the stochastic output ( , ) by an orthogonal polynomials series in the following form: where is the stochastic part of the inputs, is the deterministic one, are the deterministic functions and Ψ are the multidimensional orthogonal basis polynomials. The expansion order ( ) represents also the minimum number of deterministic runs that have to be performed for the uncertainty evaluation and it is defined as follows: where is the number of the uncertain variables and the order of the polynomials Ψ . In this work, an order of = 2 was set as suggested by Montomoli et al [25]. In this track, the minimum number of CFD simulations required is = 6. To improve the robustness of the results, an oversampling was applied as suggested by Eldered [27]. Precisely, a total number of samples equals = 9 was chosen in order to fulfil the input space completely. Since Gaussian distribution was assumed for both the uncertain parameters, Hermite polynomials as the orthogonal basis are the best choice for modelling the uncertainty propagation as stated by Ghanem and Spanos [28]. Once the orthogonal basis was chosen, the sample points were determined by using the zeros of the orthogonal polynomials of order ( + 1) ℎ . These are called collocation points, and each of them is a specific CFD simulation.
To compute the polynomial coefficients, the Polynomial Chaos Projection approach was used. Therefore, the leastsquares method was applied. The final step of the stochastic analysis consists of the calculation of the statistics of the output. For this purpose, the method proposed by Hosder et al [29] was adopted.

Gas Path Analysis
In this section, the gas turbine thermodynamic model is presented. Its schematic is reported in Fig. 3. From a mathematical point of view, the model consists of many individual components: intake, fan, compressor (HPC), combustor, high-pressure turbine (HPT), low-pressure turbine (LPT), bypass, and nozzle, each of which requires input variables and generate the adequate output variables [30]. In the following paragraphs, after a brief description of how the components maps were integrated, all the equations used are discussed.

Components maps and algorithms
In their standard format, the fan and the compressor maps cannot be used directly in a performance calculation program. The solution to this problem is to introduce auxiliary coordinates ( -lines) as proposed by Kurzke [31]. These lines are parabolas parameterized by an index , that can span from 0 (choke line) to 1 (surge line). The parabolas intersect the speed lines and allow to access to the maps quantities. An example of how -lines looks like is reported in Fig. 3.
Intake. The intake module receives the altitude and the flight Mach number as inputs. The former is used to calculate the ambient conditions ( , ) at the inlet using the ISA relations [32], whereas the second is used to compute the flight conditions in terms of total temperature and pressure: Fan. The total quantities computed represents the inlet conditions for the fan. To find the output values of this element, the fan maps have to be used. Furthermore, the design operating conditions (supplied by GasTurb9) were assumed for this component. With the design pressure ratio ( ) and the design rotational speed ( ), the corrected mass flow rate ( ) and the efficiency ( ) can be obtained from the maps: After that, the estimation of the output quantities is straightforward: where ̇ is the air mass flow rate entering the engine, is the normalized pressure at the inlet and the normalized temperature at the inlet. In the turbofan, a portion of the total flow bypasses the core unit, forming the cold stream and then the cold thrust. A high-bypass configuration as the one considered, admits a bypass ratio (B) between 5 and 8. In order to proceed, this quantity has to be also provided by the user. The cold and the hot mass flow rate can be then computed as follows: Compressor. The high-pressure compressor (HPC) is the first component of the core engine. Its task consists of rising flow temperature and pressure before it enters the combustion chamber. Inputs of this component are the temperature, the pressure and the mass flow rate exiting from the fan. Here, the first interactive loop begins. The scope is to guess the value of beta ( ) that meets the continuity constraint using the quantities extracted from the component maps: As can be inferred, knowing the value of beta, all the output quantities can be computed. In the GPA constructed, the possibility of bleeding a portion of the HPC flow was also considered.
where is the thermodynamic power output of the HPC and the constant pressure specific heat of the air.
Combustion chamber. The air delivered by the compressor mixes with the fuel in the combustion chamber. Here, the second iterative process takes place. The value of the combustor outlet temperature ( 4 ) varies into a specified range in order to meet the continuity constraint at the HPT. Since the possibility of bleeding of the compressor flow, the computation of the temperature at the inlet of the HPT is computed as follows: Turbine. Given the similarity of the procedure adopted for either the HPT and the LPT, only the former is described. The inputs available for the turbine part are the mass flow rate, the thermodynamic power developed, the rotational shafts speeds, the temperature and the pressure. The first output calculated is the turbine outlet temperature by the use of the power compatibility: where is the thermodynamic power generated by the high-pressure components ( for the low-pressure ones). The next step involves the use of turbine maps. Specifically, the efficiency of the turbine can be extrapolated.
Finally, the pressure ratio and then the outlet pressure can be obtained: Nozzle and bypass. To calculate the output values for either the nozzle and the bypass, the expansion ratio is firstly compared with the critical one. If choked conditions occur, the thrust is being raised by the pressure contribution. Otherwise, only the velocity term contributes. =̇( 7 − 0 ) + 7 ( 7 − ) The whole algorithm implemented is shown in Fig. 5

RESULTS
In this section, the results obtained in the study are exposed. Firstly, the validation of the GPA against the GasTurb9 commercial software is treated. This is to assess the reliability of the computed engine response. After that, the outcomes of the stochastic analysis are briefly covered. A deeper discussion can be found in the paper of Friso et al [18]. Lastly, the engine response consequent to the turbine maps perturbation is evaluated. To compute the engine response, a fixed Fan operating point was considered as a constraint. This means a fixed mass flow rate entering the core engine, motivated by the time scale considered for the degradation (0.1 s of exposure).

Gas path validation
Before proceeding with the engine response, the validation of the GPA algorithm developed is reported. Since the commercial code GasTurb9 was used for this purpose, the default software inputs reported in Tab. 3 were chosen. The comparison was carried out in terms of temperature distribution along with the engine, SFC, and thrust. The Temperature results are reported in Fig. 6. As it can be seen, a perfect matching occurs for the two codes results. The same is for the other two quantities. The SFC and the thrust found by GasTurb9 is 19.2 g/(kN s) and 3.3 kN respectively, whereas the values found by the in-house code is 20.4 g/(kN s) and 3.7 kN. The authors ascribe the small differences to the hidden parameters in the commercial software, which cannot be integrated into the in-house code.

Engine response
In this section, the engine response as a consequence of the turbine maps perturbations is analyzed. For this purpose, it has to be described how these maps were perturbed. Fig. 7 reports a schematic of the procedure adopted. First of all, the first two statistical moments for Γ , π , and were computed (see Tab. 4). Then, a representative probability distribution for the quantities was assumed by using the known statistics. At this point, the equally likely events were joined (green and blue dotted lines) to find the probability distribution of the shifting entity. As it can be seen, the perturbation is represented by a vector which components are the variations of the parameters due to fouling. As sketched in the figure, four shift entities were considered, i.e. the mean, and the 1, 2, and 3 events. These four altered maps were then being fed to the GPA. As stated above, the scaling factor approach has been adopted here. Different approaches would lead to a different value for the rate of shift of the maps.  The engine response to the perturbations is exposed in Tab. 5 in terms of deviations from the clean conditions of the SFC, the surge margin ( ), and the turbine entry temperature (TET). Specifically, SFC and TET are reported as deviation from the clean conditions, whereas are the actual values that correspond to the equilibrium points. As can be noted, all the parameters rise as the displacement increase. In particular, the increase of 1.4% in the SFC for a 3 event results in a rise in fuel mass flow rate and then total costs. The same for the TET, the increase of which lead to a reduction in the lifetime of the engine. Regarding the beta value, in off-design conditions, the equilibrium point tends to approach the surge line. For a 3 event, a drop of approximately 10 % of the distance between the original design point and the surge line occurs. The authors want to stress that these results correspond to exposure of 0.1 s to volcanic ash. This means that the change in the operating conditions can take place instantaneously for a flight across a volcanic cloud.
To better understand how the design point moves in the compressor map, Fig. 8 has been reported. Specifically, in Fig. 8a) the variations in terms of SFC, core rotation speed ( ), and compressor pressure ratio ( ) are depicted for a and 3 events. When volcanic ash deposits on the HPT vane, both HPC pressure ratio and rotation speed increase. This, with fixed Fan operating conditions, results in a shift upward of the HPC operating point (Fig. 8b)). The increase in speed leads to an increase in the SFC.
To summarize, when a flight through a volcanic cloud occurs, different engine responses must be expected. These are mainly due to its inherent uncertain operating conditions. Specifically, if fouling of the HPT vane takes place, the engine behaves following a sort of lognormal distribution. Restricting the analysis only up to 3 events, the operating point of the GT may end up working almost instantaneously in near surge conditions (a drop of roughly 10 % in the surge margin). Furthermore, an increase in the TET leads to an unexpected reduction in the engine lifetime. Finally, a significant rise in fuel consumption affects the total costs of the mission.

CONCLUSIONS
In this study, the stochastic engine response after a flight through a volcanic cloud has been analyzed. Specifically, an uncertainty analysis has been carried out with the results obtained by CFD simulations. The EBFOG deposition model has been utilized to account for the particles adhesion effects on the flow field. The results of the stochastic analysis have been used to feed a GPA to obtain the whole engine uncertain response.
From the analysis, it can be inferred that the hypothesis of deterministic operating conditions may result in a significant overestimation of engine degradation. In particular, considering the HPT vane as the only source of uncertainty, the following considerations can be done:  The drop in the distance between the surge line and the operating point goes from almost zero ( event) to 10 % (3 event). This may result in unstable operating conditions.  Since the engine lifetime is strictly related to the TET, its erroneous prediction may lead to an overestimation of the engine life. In this study, the maximum variation of the temperature (3 event) is equal to 0.263 %. In the work of Montomoli [25], an uncertainty of 0.6 % on the TET was found to be sufficient to reduce the HPT lifetime by 37 %. Thus, the value found in this study is not negligible, and may greatly affect performance prediction.
 As the SFC is strictly related to fuel consumption, then mission costs, its accurate prediction is of main importance. This study shows how this value may be highly uncertain, even only by considering the stochastic variations at the outlet of the combustor. An alteration of roughly 1.4 % was found in the most critical conditions (3 ), and this has direct consequences on the mission costs.
This study highlighted the importance of considering uncertainty in the prediction of gas turbine degradation. In its simplicity, this represents the first step towards a new strategy for the robust forecasting of GT performance degradation.