Digitisation of HPT blades

A fringe projection system is used (Betker et al., 2020) for the digitisation of the blade geometry. In order to automate the process and reduce systematic errors, the projector and camera of the measuring system are positioned with a serial robot (see Figure 2). The projector projects an alternating stripe pattern on the scraped and heavily deteriorated blades to measure the geometry from root to tip. The camera captures the deflection of the stripes. A measuring computer calculates the distance to the camera from these images and creates a point cloud from it. With a high dynamic range measurement, blades with different optical surface properties can be measured.

Figure 2.

Measurement system scanning a blade.

The blade is clamped using a workpiece carrier and a zero point clamping module, that is mounted on a rotation stage. On this stage, three projectors are mounted and project a stochastic pattern on the blade. This pattern, overlapping measurements, and the special geometry of the workpiece carrier are then used to merge the point clouds of the 13 single measurements to a combined point cloud.

The coordinate system of the measurement system is calibrated once with a normal. Because of the low absolute accuracy, as well as the assembly tolerances of the workpiece carrier, the position and orientation of the measured blade varies in the fixed coordinate system. An iterative closest point algorithm aligns the worn blade to a reference blade to determine the parameters of each blade on the same positions. After the alignment, blade slices of 1.2 mm around the planes that are used in the next steps are extracted and parameterisation (±600 µm) are extracted and converted into a 3D-model with surfaces. Additionally, a convex hull is used which allows the robust closing of holes from heavy damage and cooling air holes.

The alignment is performed by an Iterative Closest Point Algorithm. The alignment is carried out globally, since individual geometry parameters are not sufficient due to the geometry variations. Due to the strong wear in the blade tip, the alignment is only used in the range from the foundation to 5 mm below the blade tip. The iteration is converged if the change is less than 3 micrometers or 30 iterations are exceeded (Nübel and Denkena, 2022).

After the marcrosopic geometric properties of the investigated blades, the small-scale surface structures are optically measured with a con-focal laser scanning microscope (Keyence VK-X200). In the process, the measurement of surface roughness is performed at five positions on the suction side and at four positions on the pressure side at mid-span with a 20× magnification lens. For CFD simulations an equivalent sand-grain roughness ks is required to analyse surface roughness and compare different rough surfaces. In Nikuradse (1933) and Schlichting (1936) the equivalent sand-grain roughness was introduced, which describes the aerodynamic resistance of the surface roughness. In Bons (2010) a detailed overview of established correlations for converting a technical roughness to an equivalent sand-grain roughness is described. To assess the complexity of a three dimensional roughness topology, Hohenstein et al. (2013) and Gilge et al. (2019) propose using the shape and density parameter Λs of Sigal and Danberg (1990). Therefore, the approach of Bons (2005):

Parameterisation

The digitised blade is analysed in 20 different blade sections at constant radial positions based on the slices from the digitisation. Each blade section is parameterised by algorithms based on Heinze et al. (2014) and Ernst et al. (2016) and exemplified in Figure 3. The camber is approximated by the Delaunay triangulation (Aurenhammer, 1991), where the middle points of enveloping circles around the Delaunay triangles correspond to the camber line. This approach fails near regions with high curvature, e.g. leading- or trailing -edge. Therefore, different filters, based on the curvature of the camber, are applied in regions near the trailing- and the leading-edge in order to remove outliers. The leading-edge camber is approximated by a second order function attaching to the estimated camber based on the Delaunay triangulation and the trailing-edge by a first-order function, resulting in crossings with the blade´s macroscopic geometry. These crossings are the leading-edge and trailing-edge stagnation points which are used for rotating the blade around the origin and shifting the profile that the leading-edge stagnation point is placed at the coordinate system´s origin and the trailing-edge is horizontally aligned with the origin. This simplifies further calculations. The leading-edge radius rLE is approximated by fitting a circle within the given blade geometry using a least-square method solved by a Levenberg-Marquardt algorithm (Moré, 1978). This algorithm is used also for fitting an ellipse at the trailing-edge giving the parameter aTE along the camber line axis and bTE orthogonal to the camber. Furthermore, the curvature of the camber line has to be within certain defined tolerances. Stagger angle γ and chord length l are calculated by trigonometric functions. The geometric parameters considered in this work are shown in Table 1 including their absolute ranges and 95% confidence intervals. The absolute ranges are derived of the min-max value of each geometric parameter within the 36 digitised blades. The measurement accuracy of the digitisation and parameterisation is determined by measuring two blades in five different fixtures. The point clouds are aligned and the parameters extracted. The maximum deviation for each parameter is used as approximation of the measurement accuracy. Additionally, the range of surface roughness applied to the CFD simulations is presented with a relative error of the optical measurements based on a confidence level of 95% (see Gilge et al. (2019)).

Figure 3.

Blade section of deteriorated blade at mid-span.

Figure 4.

Heavily deteriorated HPT blade used in this project.

Reconstruction of HPT blade and sample generation

The blade reconstruction is performed using an original CAD blade geometry and altered based on the difference between reference and measured blades. This difference is applied for all blade sections along the span. In this work, the reference blade is defined as the least deteriorated one. Changing the geometry of the leading- and trailing edge requires smoothing between the new circles (for leading-edge) or ellipses (for trailing-edge) and the rest of the profile geometry for obtaining a continuous surface. Otherwise, the algorithm used is based on Ernst et al. (2016).

The ranges of the parameters are shown in Table 1 and a total of 250 blades with different geometries is reconstructed, each representing a combination of different geometry parameters and surface roughness. In order to carry out the numerical experiment effectively, a uniform Latin-hypercube sampling (McKay et al., 2000) is used. The tip gap variations are included during the meshing process. The roughness is implemented in the CFD setup on both, the suction and the pressure side. The digitised blades of the first rotor stage belong to a V2500-A5, while only original blade geometries of all rows of the V2500-A1 are available for usage in the computational models. To overcome this issue, the geometric variations are normalised with the blade height, assuming that both variants of the V2500 show geometrically similar wear.

CFD simulations of HPT with deteriorated first rotor

Based on the reconstructed blades, numerical set-ups are developed and simulated in a fully automated process. The CFD simulations are necessary to connect the effects of local deterioration, here first rotor blade, with the module performance. First, the meshes of the 2-stage HPT are constructed using NumecaAutogrid. The meshing parameters are set constant for all reconstructions to ensure mesh independence of the numerical solution. According to the DoE, the different tip gaps are included in the meshing. A total of 211 HPT simulations with a deteriorated first rotor blade are successfully meshed. For 39 combinations of geometric parameters, the meshing fails due to negative cell volumes or connectivity problems across the rows. In Figure 5 a mesh of an HPT is illustrated with the first blade highlighted in red. The original HPT geometry was already used in Goeing et al. (2020b) for evaluating combined compressor and turbine deterioration. Secondly, the simulations are performed using the non-commercial flow solver TRACE version 9.1 of the German Aerospace Center (DLR) (Nuernberger, 2004; Franke et al., 2005; Kugeler et al., 2008). The three-dimensional Favre-averaged Navier-Stokes equations are solved by a finite-volume method with structured multi-block meshes. Roe's second-order upwind scheme is used for the discretisation of the convective fluxes, while the diffusive fluxes are solved by a second order central difference scheme. In an automated pre-processing, the boundary conditions of the aerodynamic design point at take-off with a rotational speed of 13,972 rpm prescribed by the performance simulation are transferred to the CFD set-up. Next to the boundary conditions, the equivalent sand-grain roughness ks is specified at suction and pressure side of the first blade row. At the inlet, a total pressure of approximately 2,823 kPa and a total temperature of approximately 1,735 K is used. Additionally, a turbulence intensity of 4% and a turbulent length scale of 0.00208 m is set, which are typical for turbomachinery flows. At the outlet, a static pressure of 491.6 kPa is specified. Subsequently, periodic steady-state simulations with one pitch are conducted using the k−ω turbulence model of Wilcox (1988), assuming that the analysed blade represents the entire row. For the blade to vane interfaces mixing planes are used and the stagnation point anomaly fix by Kato and Launder (1993) is applied. The sidewalls are modelled through a wall function. In general, a non-dimensional distance of the first cell from the wall of y+<1 is obtained at the blade's surfaces, which allows resolving the viscous sublayer using the Low-Reynolds approach. Nevertheless, the first rotor blade is modelled by the wall function due to the increased wall shear stress caused by surface roughness. The logarithmic profile is described by:

Figure 5.

HPT simulation setup with studied first rotor blade highlighted in red.

Although an accurate initial flow solution is used, a CFD simulation of a mesh with approximately 8 M cells requires about 84 CPUh for a converged solution. A mass flow difference between inlet and outlet, a change in efficiency, and a change in pressure ratio of less than 0.1% satisfy the stopping criterion.

For the validation the original CAD blade geometry with a hydraulically smooth surface is used. The CFD simulations predict a polytropic efficiency of 88.68%, a pressure ratio of 5.2, and a massflow of 47.28 kg/s, while performance simulation gives 88.80%, 5.1 and 48.02 kg/s for same boundary conditions. With the same setup a mesh convergence study is conducted, resulting in a GCI of 0.136% for the mass flow, 2.59% for the polytropic efficiency, and 0.005% for the pressure ratio.

Performance simulation of the overall aircraft engine

In order to evaluate the influence of isolated and combined module variances on the whole turbofan engine a performance model of the V2500-A1 aircraft engine is developed to analyse the on- and off-design performance. The model is part of a model-based non-linear gas path analysis (NLGPA) to calibrate the performance model to real aircraft engines (here a factory-new and deteriorated V2500-A1 turbofan engine (see Table 2)).

This analysis investigates the ranges in which the miscellaneous compressors and turbines wear out in addition to the HPT. In conjunction with the CFD simulations of the HPT, the NLGPA and, thus, the virtual process allow a combination of the overall aircraft engine performance with the local effect of deterioration. This analysis is presented by the schematic flow chart illustrated in Figure 6 and is implemented in Matlab (Kurzke and Halliwell, 2018; Salomon et al., 2021).

Figure 6.

Global off-design calculation procedure for high-bypass turbofans. Used for non-linear GPA and to train the artificial intelligence.

The off–design calculation procedure consists of an iterative simulation process which requires various input data, such as environmental conditions. Moreover, further aircraft engine-specific boundary conditions are required, being crucial for the thermodynamic cycle at a reference point, also referred to as the engine's on–design operating point. This point represents the basic framework of the jet engine and defines the interaction between the miscellaneous turbomachines, secondary air system, geometry (e.g. nozzle area), and number of revolutions in order to fulfil the thermodynamic cycle. Further boundary condition, such as the steady-state performance maps of the compressors and turbines are required. Hereon, quantities related to the compressor and turbines are denoted with (.)C and (.)T sub-indexes, respectively. These characteristic diagrams describe the mass flow m˙, pressure ratio π, and efficiency η of the turbomachines by different rotational speeds. Furthermore, auxiliary coordinates, the so-called GL–lines, are placed through the diagram, which are necessary for the iteration of the algorithm. The on-design performance values and the boundary conditions for the V2500-A1 jet engine are based on Goeing et al. (2020b), as well as on the data of the IFAS-research aircraft engine.

After the aircraft engine–specific boundary conditions have been imported and the operating point to be reached has been defined, the analysis starts the iterations. A distinction is made between the outer loop, in which the rotational speed N2 is varied until the required aircraft engine thrust is reached, and the inner loop, in which a correct thermodynamic cycle, based on the input parameters, is matched through iterations. Matching in this context means iterating within the performance maps until: (1) the turbine power PT output matches the compressor power PC, as defined by the error function EP in Equation 5; (2) the mass flow is maintained, as defined by the error function Em˙ in Equation 6; and (3) the nozzle pressure pnozzle,t8 is equal to the downstream LPT pressure pLPT,t8 (including friction), as defined by the error function EN in Equation 7. The convergence limits are 1e-6.

Based on the on–design input parameters, the initial values of all necessary parameters are provided for the first iteration. The iteration parameters for the inner loop at any rotational speed N2 are the turbine entry temperature TTET, bypass ratio (BPR), the rotational speed N1, and the auxiliary coordinates GLT and GLC. Based on these, the operating point (constant thrust) of each individual aircraft engine is defined and a global cycle process calculation can be carried out. The next step is to validate whether Equations 5–7, are fulfilled and the error functions approaches zero. If this is not the case, then the iteration parameters are varied using the Newton–Raphson method until an operating point is found. The iteration parameter are varied throughout the ranges of all performance maps until a matching cycle according to Equations 5–7 is found. The next step is to check if this correct cycle also achieves the required thrust. If not, the rotational speed N2 is varied using the Newton–Raphson method. Here the convergence limit is 1e-5.

The global aircraft engine matching process is described with these iterations and is used in the DoE to physically calculate the relationship between state of jet engine and performance output. The state of the engine is described by the variation of the maps through the scaling factors (Δ) of π, m˙ and η.

In order to estimate the scaling factors for the LPC, IPC and HPC, as well as for the LPT, the performance model is calibrated to a factory-new and to the worn aircraft engine. The calibration is done via non-linear GPA, in which the performance map is iterated by scaling factors Δ via Newton-Raphson method until the performance model has the same performance output as the investigated aircraft engine. Here, the convergence limits are 2e-2 is assumed to be sufficient.

Sensitivity analysis

A fundamental knowledge of the effects of variations in the input variables of a model on the output variables is of paramount importance in a variety of engineering domains. Sensitivity analysis examines precisely these relationships. Saltelli et al. (2004) provide the following definition of sensitivity analysis: ‘The study of how uncertainty in the output of a model (numerical or otherwise) can be apportioned to different sources of uncertainty in the model input’. Therefore, an effective and wide-ranged group of sophisticated sensitivity instruments are variance-based methods (Saltelli et al., 2010). The Sobol indices (Sobol, 2001) are a widely utilised representative of this type of sensitivity analysis tools for model outputs. These indices not only allow to identify the direct effects of particular input factors on the variances of output variables, but also of all potential impacts due to interaction effects between the variances of the inputs on the outputs.

Metamodeling

In sensitivity analysis for the vast majority of cases, a relatively large sample size is mandatory. In aerodynamic engineering, however, many of the models employed are highly complex and simulations are consequently computationally demanding. As a result, generating sufficient samples is often too expensive or even unfeasible. A common approach to address this challenge is to develop and utilise metamodels that replicate the functionality of the original model, while being less complex and consequently less computationally intensive (see Shahsavani and Grimvall (2011) and Ben Abdessalem and El-Hami (2014)). In the context of aerodynamic engineering in general, and aircraft engine diagnostics in particular, such metamodels are increasingly based on artificial intelligence and deep learning approaches (see for example Khan and Yairi (2018) and Fentaye et al. (2019)). With these, larger sample sizes can be generated with less computational effort and thus, common sensitivity analysis tools can be applied.

In this work, artificial neural networks are trained to guarantee a sufficient sample size for the respective sensitivity analysis. A schematic layout of such a network is shown in Figure 7. It basically consists of an input layer, one or more hidden layers and an output layer. The number of neurons per layer, as well as the number of hidden layers, depends on the particular original model and cannot be generalised. The number of input neurons, as well as the number of output neurons, typically follows the number of input and output variables of the model to be mapped. Each neuron of a layer is weighted connected to each neuron of the predecessor and successor layer. For a detailed description of neural networks in general and fitting function networks in particular, as well as their training and handling, see Rafiq et al. (2001) or Samarasinghe (2006). Note, that the utilization of an artificial neural network, such as employed in this work, represents only one potential approach among various approaches, such as, e.g., polynomial chaos expansion, to perform sensitivity analysis under the condition of a small data set. Here, the methods were exclusively chosen for the purpose to present proof of concept further methods will be investigated in subsequent work of the authors.

Figure 7.

Architecture of a deep learning neural network.

HPT Sensitivities for single rotor deterioration

The performance sensitivity analysis of the overall aircraft engine requires significantly more than 211 performance samples of the HPT (see Figure 8). To expand the data, the described neural network as metamodel with three hidden layers, ten neurons per layer, eight input neurons, three output neurons, 15 epochs, and the Levenberg-Marquardt as back propagation algorithm is used. With these parameters a sufficient optimum was found, and for example, reducing the number of hidden layers results in a decrease in accuracy. Especially the epochs and neurons per layer showed an insensitivity for small changes in regards to prediction accuracy. For each neuron a sigmoid function is utilised as activation function. The final meta model reached a mean squared error of 1e-4 for all three outputs combined within a ten-times cross validation. Further improvement in accuracy would require more samples as well as further studies on different meta model approaches. The required extension of the geometric input is performed based on a Latin-hypercube sampling. Then, a variance-based sensitivity analysis is conducted on the input and output of the metamodel by means of Sobol indices (see Figure 9). This approach is validated by analysing the 211 original performance samples with the variance-based sensitivity approach according to Plischke (2010) and Plischke et al. (2013), so called EASI (Effective Algorithm for Computing Global Sensitivity Indices). This approach is purely based on existing data and does not require a specific input sampling. The sensitivity analysis is conducted on the three integral aerodynamic outputs in respect to the geometric inputs. First, a non-monotonic relationship is identified between mass flow and polytropic efficiency based on scatter plots. Classic correlation coefficients, e.g. Spearman's or Pearson's Rank correlation coefficient, are not sufficient in these cases. However, variance-based sensitivity studies are more suitable for detecting such non-monotonic relationships Saltelli et al. (2010). For the original CFD data the EASI indices of first-order indicate a direct sensitivity between stagger angle γ on pressure ratio π and mass flow m˙ approximately close to 1 (see Figure 8). This indicates negligible influence of the other geometric parameters for these output parameters. The EASI indices in Figure 8 show a sensitivity of around 0.1 for the influence of tip gap on pressure ratio and mass flow. However, due to the small sample size these results may include uncertainties and all other geometric parameters can be neglected. The polytropic efficiency has a sensitivity of 0.6 for γ and 0.2 for tip gap r. The influence of the tip gap is explained by the direct dependency of mass flow over the blade's surface resulting in higher pressure differences between pressure and suction side, and the leakage mass flow over the blades tip, which does not distribute to the turbines work conducted. The observation of the stagger angle agrees with the literature of Ernst et al. (2016) and Högner et al. (2016). The influence of the tip gap agrees with Yaras and Sjolander (1992).

Figure 8.

EASI indices for CFD data.

Figure 9.

Sobol indices after applying the HPT-metamodel.

The sample size is increased according to Saltelli et al. (2004) with N=n⋅(NumberofInputs+2) with seven input parameters and n=211 to N=18,432 samples for the Sobol sensitivity analysis by means of Latin-hypercube sampling. For all inputs and outputs under consideration, a convergence study of the mean as well as of the 2nd to 4th central moment was conducted. All examined variables converge according to the sample size N employed. Solely the output π exhibits slight leaps in the 3rd and 4th central moment. However, these have no influence on the results obtained, see Figure 10. The relationships of the original data are confirmed on the metamodel output, which is shown in Figure 9. Leading-edge radius, trailing-edge thickness, chord and surface roughness on both blade sides are negligible for all output parameters. Mass flow is completely determined by the stagger angle. The polytropic efficiency has a Sobol index of 0.6 for the stagger angle and 0.35 for tip gap, which is increased for the tip gap by 0.15 in comparison with the approximation of the original data set. For the pressure ratio the sensitivity regarding the tip gap is decreased further and the stagger angle is held constant versus in the original data.

Figure 10.

First and Total Kucherenko Indices of the scaling factors and overall performance output.

It is shown that the estimated variance-based sensitivity analysis of the original data is similar to the Sobol analysis with a metamodel based on the original data. However, for this study it is only possible to analyse the direct relationship between inputs and outputs, i.e., first-order indices. Higher-order effects, that can show the systems interactions of input variables and their combined influence on the output variables, require larger sample sizes. Subsequently, an analysis of total-order effects purely based on the metamodel can not be validated and is therefore not conducted within this study.

Performance analysis

In order to estimate the impact on the overall performance of the V2500-A1 aircraft engine, the ranges of the scaling factors of each turbomachine have to be estimated. These ranges are determined on the one hand by CFD of the HPT, and on the other hand by the NLGPA of a factory-new and operational stressed aircraft engine V2500 (see Figure 6 and Table 2). The resulting scaling factors are an identifier for the real range of degradation of the compressors and turbines in the turbofan. These parameters are validated by previous studies (Goeing et al., 2020a) and further literature (Maiwada et al., 2016; Cruz-Manzo et al., 2018) (see Table 3). Thus, the Δm˙ and Δπ of the HPC have the most significant impacts (−8%), followed by Δη of the HPT (−6%). Furthermore, the smallest impact was found for the LPT (−1% and −2%). Based on these scaling factors Δ, a uniform Latin-hypercube sampling is performed to simulate the performance of 200,000 different combinations of varying V2500-A1 turbofan engines at the same thrust level. Next, 70% of the results are used as a training set and 30% as a validation and testing set for the neural network. The mean square error of the metamodel is less than 1e-8.

This metamodel is used with 1e6 random parameter configurations within the parameter space to evaluate the isolated and combined sensitivities. The input variables are interdependent, due to the power balance and mass conservation (see Equations 5 and 6), so a variance-based sensitivity analysis according to Kucherenko which is able to analyse such interdependent input variables is performed.

In Figure 10 the isolated (first) and isolated plus interacting (total) sensitivities of measurable variables are shown. Therefore, the first- order indices are represented with blue for Tt5, green (SFC), light red (Tt3), yellow (N1), purple (N2) and brown (Pt5) bars. The stacked red bars represent the total indices. Further on, the 95% confidence interval is shown as a black error bar for the first-order and a grey error bar for the total order. Scaling factors Δ (e.g. IPC or Combustor) with an influence of less than 0.05 are neglected for the sake of clarity.

As shown in Figure 10, most of the peaks are located at the HPC and HPT. In particular, the EGT (0.44/0.38) and the SFC (0.36/0.32) are responded to the efficiency of the HPC/HPT. It is noticeable that the temperature Tt3 reacts on the efficiency of the HPC (0.56), as well as the mass flow of HPT (0.30), while the mass flow of the HPC has a significant influence on the HP-System rotational speed N2 (0.64). N2 also reacts to the efficiency of the HPC/HPT (0.14/0.13). Furthermore, there is a significant relationship between N1 speed and the pressure ratio and the mass flow rate of the LPC (0.57/0.42). In addition, the pressure Pt5 reacts sensitively to efficiency of the LPC and LPT (0.23/0.24).

In order to illustrate the virtual process and to demonstrate the significance of the interactions on the overall performance, the described methods are performed with the real HPT blade in Figure 4 as an example. The geometry of the worn blade, roughness parameters, and of their confidence limits are derived and supplied to the process. Subsequently, the metamodel of the CFD is used to obtain the module aerodynamics and scaling factors. The HPT efficiency is reduced by 2.4 ± 0.2%, the pressure ratio is reduced by 0.046 ± 0.003, and the mass flow is increased by 0.64 ± 0.08 kg/s compared to new blades. These substantial aerodynamic changes of the scrapped HPT blade confirm the heavy deterioration stated earlier.

In the next step, these results are integrated in five aircraft engine configurations Θ by means of the performance metamodel. The configurations represent the following health conditions

Figure 11.

Impact on Tt5 of miscellaneous aircraft engine health conditions Θ using the example of the real HPT blade in Figure 4.

For this purpose, the Tt5 difference between the real and the new turbine blade over the various configurations Θ is shown. The red bar represents the influence of the HPT, the blue bar the influence of all other deteriorated modules on the Tt5 temperature. As expected, Tt5 increases the most due to a deteriorated HPT and HPC. It is particularly noticeable that the impact of a the HPT blade on the Tt5 is not constant for each configuration. Especially, for the aircraft engines Θ3 and Θ4, the difference decreases from 50 K to 41 K. This effect can be explained by the interdependence of the various modules, e.g. the throttling of the various compressors due to degradation. It agrees with the high global sensitivity previously described and already shown in Figure 10. The strength of these interaction effects is indicated by the red bars in Figure 10. The interaction effects are particularly significant for modules whose isolated influence is low (see IPC and Θ3). The relationship and level between HPT deterioration and Tt5 is comparable to studies on real aircraft engines (see Zaita et al. (1997)).