A Methodology for a Coupled Structural–Computational Fluid Dynamics Analysis of Compressor Rotor Blades Subjected to Ice Impact With Uncertain Impactor Parameters

During operation, aircraft jet engines are exposed to certain weather scenarios and environmental conditions where ice is ingested into the engine core which may cause severe damage and loss of efficiency, see, for example, Refs. [1–3]. Ice ingestion into the compressor can be caused on the one hand by previous hail ball ingestion and its debris and on the other hand by detachment of accreted ice shreds from surfaces upstream of the compressor. The ice then travels downstream and enters the compressor. Here, especially the first stage rotor blades can suffer severe local mechanical damage by the ice impacting the leading edge. Resulting large plastic deformation of the rotor blades caused by high impact forces directly affects the aerodynamic performance of the rotor by changing stagger and blade inlet angles.

To enable a robust design of first stage compressor rotor blades which are vulnerable to ice impact scenarios, a consistent understanding of the caused structural damage and the related aerodynamic behavior of the system are essential but depend on several uncertain parameters: Major uncertainties exists for geometrical and physical properties of ice primarily caused by complex formation processes resulting in different crystalline structures. Many studies address the evolution of weather phenomena in which hail balls are formed by freezing water droplets, e.g., see Refs. [4–6] and thereby present variations of specific hail ball properties such as size and density, see Refs. [7,8]. These varying properties, in particular shape and mass, are crucial for its fundamental material behavior of the ice upon impact, which affects the occurring contact force on the impacted blade and its resulting structural response. Mechanical properties of ice on the other hand do not play a significant role for the resulting impact forces on high speed rotating compressor rotor blades as occurring inertial stresses are much higher than material strength so that a solid body shows hydrodynamic characteristics like a fluid under such high speed impact conditions as discussed in Ref. [9]. Several probabilistic models exist to predict property distribution of hail balls under severe hail storm conditions, for example, see Refs. [10,11]. A detailed correlation of these uncertain properties with their damaging potential and resulting effect on aerodynamic performance requires an efficient analysis methodology to consider all the relevant uncertain parameters of the impact. In Ref. [12], a probabilistic model was applied for the prediction of the maximum contact force due to ice impact.

Structural analysis of ice impact to assess the resulting mechanical damage is mostly computed by explicit finite element software, e.g., ls-dyna, which offers sophisticated modeling techniques to enable excessive deformation with defragmentation as shown in Ref. [13] and material models to capture fluid-like material behavior of the ice impactor during the impact, see, for example, Refs. [14–16]. Only a few recent studies focus on influence and change of flow and aerodynamic characteristics due to the deformed blade by foreign object impact events, in particular [17,18]. An optimization method for the multidisciplinary design process of compressor blades with regard to aerodynamics and ice impact worthiness has been presented by Schlaps et al. [19].

This present paper is structured as follows: At first, a probabilistic approach for those parameters governing ice impact such as impactor size and mass as well as impact location and rotational speed is presented. The following chapter describes the structural impact analysis method employed to predict mechanical damage and deformations of the impacted blade. Next, the computational fluid dynamics model to assess the effect of blade deformation onto the compressor flow is described. The paper ends with an uncertainty quantification study relating the uncertain input parameters to the mechanical and fluid dynamic effects based on an optimal Latin-Hypercube Monte Carlo analysis.

The hail ball diameter is approximated by a normal distribution based on observations over a 10-year period published by Changnon [7]. According to Fig. 1, the discrete frequencies H are taken from the literature to perform a least squares regression of the cumulative density function, which yields the mean μ and the variance $σ2$⁠. The density of the hail ball is assumed to be uniformly distributed. The corresponding limits are defined following Ref. [8]. Moreover, the impact height and the angular velocity are defined as uniformly distributed within the specified limits. All input parameters are assumed to be uncorrelated, which is ensured by applying the restricted pairing method according to Ref. [23].

For the probabilistic model, some simplifications were made that may affect the sensitivity of the system: e.g., uncorrelated input parameters with—except for the hail ball size—assumed probability density functions (PDFs). Under real world conditions, the input parameters would likely be correlated (e.g., impact height versus hail ball mass) and not uniformly distributed (e.g., impact height), which should be considered in future work.

herein, $ΔFT∗$ denotes the relative change in thrust, $Δη$ is the change in isentropic efficiency, $ΔΠt∗$ is the relative change in total pressure ratio, $max(sc)$ is the maximum cup size, and f_i denotes the eigenfrequencies. The numerical models for structure and aerodynamics are described in Secs. 3 and 4 along with the definition of the mentioned result variables in Sec. 5. Subsequently, a sensitivity analysis is presented to quantify the contribution of a probabilistic input parameter x_i on the result variable y_i. For this purpose, the coefficient of importance (CoI, see Refs. [24,25] are employed.

A quarter-assembly of the whole compressor rotor is considered for structural analysis in order to reduce computational time effort, as depicted in Fig. 2. This model is based on the second stage of the transonic high speed compressor rig Rig250, which is representative for modern high pressure compression subsystems. More information about the rig may be found in Refs. [26] and [27].

The model contains seven individual blades, which are meshed with eight-node solid elements according to the approaches of Refs. [28] and [29]. The airfoil is modeled via three elements through the thickness. Each blade contains 58,538 elements and 71,190 nodes in total. As the model only reflects a quarter of the rotor, symmetric boundary conditions were applied to the corresponding surfaces in circumferential direction. The displacements of the assembly along the x-axis were fully constrained at the interfaces to the hub. As material of the rotor blisk a titanium alloy Ti-6Al-4V is used, modeled by the *MAT_JOHNSON_COOK material model based on Ref. [30], in combination with a Mie-Grüneisen equation of state [31]. This model expresses material stress as function of strain rate and plastic strain failure criterion. All parameters are derived from Refs. [32] and [33]. In this study, the explicit calculations were parallelized on the high performance computational cluster at the TU Dresden via 48 Intel(R) Xeon(R) CPU E5-2680v3. ls-dyna R12.0 MPP solver version was used with single precision for the explicit analysis and double precision for the implicit analyses.

In a first step, the rotor assembly is preloaded by assigning the rotational velocity ω in an implicit analysis step for each of the 50 eLHS-based case investigated impact scenarios Table 2 (see Appendix). Due to the occurring centrifugal force, an untwist deformation of the blade and prestresses arise, which need to be considered in a subsequent impact analysis. In a second explicit analysis step, the rotating assembly is impacted by a spherically shaped hail ball with diameter d_hb and density $ϱhb$ at impact height z_hb (Fig. 3).

The hail ball has no translational velocity and is modeled via 552 SPH particles for all scenarios which results in a varying interparticle distance. Material modeling of the ice material is carried out according to Ref. [14]. It should be noted that the variable density $ϱhb$ according to the different impact scenarios (Table 1) is represented through the particle mass definitions in the SPH model, whereas the material property definition is kept constant throughout the simulations. The hail ball impact event and its specific slicing process at the leading edge are illustrated in Fig. 4.

After impact, the rotor assembly is brought to rest in an implicit analysis step until the elastic strains are fully recovered. The resulting plastic deformations in the leading edge caused by highly localized impact forces (Fig. 5) can be expressed using the so-called cup size parameter s_c, visualized for selected scenarios in Fig. 6 referring to the leading edge positions along radius r.

In addition, modal analysis of the deformed single blade section is performed at the operational speed of $1153 rad/s$ to evaluate the sensitivity of the modeshape frequency. For the transfer of the deformed geometry to the CFD analysis, node coordinates and node displacements were extracted from the binary ls-dyna result files in an automated process using a user-defined script on the basis of the LASSO python library [34].

In order to be able to assess the flow changes caused by a single damaged rotor blade, the model for the CFD analysis contains the full 360 deg ring of the rotor, where one blade represents the damaged geometry and all other rotor blades are defined by the nominal blade shape. The damage shape is transferred from the FEA simulation to the CFD mesh using a morphing approach with radial basis functions. Additionally, a single stator passage is added upstream and downstream of the rotor domain to complete the CFD setup, as depicted in Fig. 7. This results in a block-structured mesh of 110 × 10⁶ nodes, which has been created using Rolls-Royce best practice guidelines and the results of previous mesh refinement studies.

In particular, the mesh of the rotor is sufficiently refined in radial direction in order to map the deformation shape from the structural analysis to the CFD mesh with high accuracy, see Fig. 8. The CFD simulation is carried out with the steady-state Rolly-Royce in-house solver Hydra, which solves the compressible Reynolds-averaged Navier–Stokes (RANS) equations. The turbulence is modeled by means of the turbulence model according to Ref. [35] with adaptive wall functions. Further details on the flow solver can be found in Ref. [36]. Mixing planes are used to couple the rotating and stationary domains. Boundary conditions for total pressure and temperature, flow angles, and turbulence intensity at stator inlet and radial equilibrium with prescribed reduced massflow at stator exit are employed for all simulations and were extracted from a CFD simulation of the full compressor domain.

In order to quantify the effect of the deformed rotor shape on the compressor flow, the deviations of the following flow quantities to the baseline are evaluated. In particular $ΔFT∗$ denotes the relative change in thrust, $Δη$ the change in isentropic efficiency, and $ΔΠt∗$ the relative change in total pressure ratio. All quantities are calculated using mass-averaged values of the stator inlet boundary upstream and the stator outlet boundary downstream of the rotor and extracted from each of the damaged setups. Exemplarily the result of the CFD solution for a severely damaged blade corresponding to configuration A is shown in Fig. 9, where the flow disturbance due to the large deformation is clearly visible, causing the observed drop in isentropic efficiency for this configuration.

R² denotes the coefficient of determination (COD) of a polynomial approximation with a certain polynomial order. In this context, $Ry,x2$ describes the COD with the entire polynomial basis and $Ry,x∖xi2$ the COD with the polynomial basis where all coefficients referring to x_i are removed. Therefore, the quality of the approximation should be considered when a CoI-based sensitivity analysis is performed. For this purpose, a k-fold cross-validation is performed, which yields $RCV2$⁠.

The application of eLHS allows to sequentially increase the number of model evaluations until certain quality criteria are met. Hereby, an increasing $RCV2$ is observed, while R² decreases until the total number of 50 model evaluations is reached. This effect is caused by over-fitting of the training data and is known from literature (e.g., see Refs. [37,38]). The final CoI sensitivity matrix is shown in Fig. 10, where the highly nonlinear system behavior is chosen to be represented by third-order polynomials (⁠$PO=3$⁠). Accordingly, the sample-to-coefficients ratio amounts to 3.8. Note that the shown sensitivity measures are based on the previously mentioned hypothesis of uncorrelated input parameters whose PDFs were—hail ball size excepted—estimated.

It becomes obvious that the most important parameter is the hail ball diameter d_hb. This can be demonstrated with sufficiently accurate polynomial approximation for aerodynamic result variables (e.g., $Δη$⁠) and the maximum cup size (⁠$max(sc)$⁠). In contrast, the approximation quality is rather low for almost all eigenfrequencies. However, the analysis leads to the conclusion that the impact position z_hb is a relevant parameter for eigenfrequencies whose importance occasionally exceeds that of the parameter d_hb (e.g., f₄). Given the low approximation quality in case of almost all eigenfrequencies, this should be understood as an indication rather than an evidence-based sensitivity measure. The respective dependencies, therefore, require a more detailed discussion regarding their causality, which will be discussed later in this section.

The influence of the hail ball diameter d_hb on selected result quantities is visualized in Fig. 11 by scatter plots along with the probability density functions f of the respective parameters. In the diagrams, the previously selected MCS runs are highlighted and labeled A, B, and C, respectively. Configurations A and C reflect an impact at about 60% span, while $dhb,A>dhb,C$⁠. In addition, large plastic deformations (⁠$max(sc)$⁠) result for A and B due to the relatively large hail ball diameter. In contrast to A, the impact in case B occurs near the blade tip at about 85% span.

According to Fig. 11(a), a critical hail ball size (⁠$dhb,crit≈15 mm$⁠) is required to create noticeable deformations, which is captured by the cup size s_c. Once $dhb>dhb,crit$⁠, a sharp increase in $max(sc)$ is observed, propagating through aerodynamic and structural quantities. This nonlinear and monotonically positive dependency is expressed by the rank correlation coefficient: $r̃=0.91$⁠. In case of $dhb<dhb,crit$⁠, the deformation that occurs is predominantly elastic. Mechanical and aerodynamic quantities therefore remain largely unaffected.

In terms of isentropic efficiency, the plastic deformations (⁠$dhb>dhb,crit$⁠) cause a slight increase (⁠$dhb≈20 mm$⁠) before a drastic decrease occurs (see Fig. 11(b)). The slight increase in efficiency could be due to mechanisms associated with shock control bumps known from literature (e.g., Ref. [39]). However, a clear explanation would require a more detailed aerodynamic analysis.

With respect to eigenfrequency f₄, the scatter plot (see Fig. 11(c)) shows a branching system, where each branch represents a specific causality between the inputs $x$ and f₄. Since third-order polynomials, as used for the CoI computation, are not capable of describing this system response, the approximation quality, expressed by $RCV2$⁠, is low. It is worth mentioning that the eigenfrequencies $f1,3,9$ do not exhibit such a distinct branching system. Therefore, the approximation quality is significantly higher.

To elaborate on this discussion, Fig. 12 shows the mode shapes of eigenfrequency f₄ for the selected MCS runs compared to the BL. This mode is characterized by remarkable amplitudes $||u||2$ in the blade tip region, especially at the leading edge. It can be inferred that an impact at about 60% span (A, C) causes a fundamental change in mode shape. Here, a hail ball diameter close to the critical value is sufficient to evoke this effect (⁠$dhb,C=16.6 mm$⁠). This change in mode shape is accompanied by an increase in f₄ of $81 Hz$ and $126 Hz$ for A and C, respectively. In contrast, f₄ decreases by $135 Hz$ due to ice impact near the blade tip (B). Although the mode shape is comparable to BL, moderately increased amplitudes are observed. To quantify the influence on mode shapes due to hail strikes more precisely, future work should incorporate the modal assurance criterion (MAC, see Ref. [40]).

In the case of a relevant ice impact, it is most likely that not only one impact occurs, but multiple ones. Consequently, several blades are deformed, which affects the eigenfrequencies accordingly and leads to probable frequency bands instead of deterministic solutions. In this context, Fig. 13 shows the frequency bands resulting from the present modal analysis. Here $Q0.5$ denotes the median of the frequency bands. The analysis is performed by varying the relative rotational velocity $ω/ωBL∈{0,0.5,1.0}$⁠. Note that the modal analysis is performed using single-sector models of the blisk. As a result, the following conclusions can be drawn:

Severe deformation due to impact would result in immediate overhaul and repair of the engine to avoid any risk. However, the broadening of the higher order frequency bands should be taken into account since the present study has shown that even small deformations at a sensitive location can significantly affect the eigenfrequency. The emergence of frequency bands due to geometric imperfections is known from literature, e.g., Ref. [41], where the influence of manufacturing scatter on the modal properties of a compressor blade was investigated.

In this study, a coupled structural-CFD analysis methodology has been applied on a compressor rotor blade row subjected to hail ball impact for the assessment of uncertain hail ball parameters and their effect on the system response. For this, in total 50 impact scenarios considering varying diameters and densities of the hail ball impactor, impact positions, and velocities of the rotor were derived from a probabilistic modeling approach. In a first step, structural ls-dyna analysis for all scenarios was carried out to obtain the impact-induced mechanical damage of the blades and respective vibration characteristics. Subsequently, the damaged blades are transferred to CFD analysis for each case to evaluate its effect on the aerodynamic properties.

It is shown that hail ball diameter may reach a critical size of $15 mm$ where considerable damage of the blade is caused by means of a deformation of the leading edge. In these cases, a sharp decline in isentropic efficiency is observed. On the other hand, the impact position plays a significant role in for the resulting vibration behavior, which may lead to a change in mode shape as well as a change in the eigenfrequency.

The system behavior is analyzed using the sensitivity measure coefficient of importance based on an MCS with 50 model evaluations. The analysis reveals the pronounced nonlinear properties and branching characteristics of the investigated system, which complicates the sensitivity analysis. However, the dominant influence of hail ball size on both aerodynamic and structural result quantities is demonstrated under the given conditions (e.g., uncorrelated input parameters, partially assumed PDFs). Moreover, the influence of the impact position on eigenfrequencies has been shown and discussed in terms of the underlying causality.

It should be noted, that the investigated ice impact due to hail balls is not primarily a representative loading scenario for compressor blades of the first stage, as undamaged hail balls impact the compressor in the rarest of cases. However, the developed and applied simulation approach provides the basis to focus on the more relevant impact scenario of ice sheds on the compressor blades due to icing. To achieve this, the current lack of information needs to be overcome in regard to the characteristics of the ice shed configurations in terms of shape, structure, and material properties. The presented approach is more relevant for the fan blades, which are more vulnerable to hail ball impact loading. Additionally, multiple impact scenarios where combined impact induced damages are considered should be subject of further investigations.

In addition, a comprehensive uncertainty quantification study could examine the statistics of the input parameters being considered to reveal their correlation structure and PDFs. For instance, the impact locations could be determined through the approach of Ref. [42] based on optical scans. This would contribute to a better understanding of the sensitivities under real world conditions.

In this context, other uncertain parameters such as the geometry of the compressor blades should also be included. For this purpose, Ref. [43] provides a parametric model that is widely used. In addition, specific geometric models (e.g., Ref. [44]) could be embedded to reflect the specific geometry of the damaged leading edge and consider their relevance for aerodynamic performance.

Finally, the nonlinear and branching characteristic of the system under investigation (eigenfrequencies) raised reliability issues of sensitivity measures. The limited capability to reflect this behavior through polynomial approximation could lead to ambiguous outcomes. The implementation of more sophisticated approaches for coefficient of importance computation (e.g., moving least squares, radial basis functions) could help to overcome this issue.

The authors wish to acknowledge the support of Rolls-Royce Deutschland (RRD) through the Lightweight Structures and Materials and Robust Design University Technology Centre (UTC) at the Technische Universität Dresden (TUD). The work presented in this paper was conducted within the framework of the PRESTIGE research project (20T1716A). Rolls-Royce Deutschland's permission to publish this work is greatly acknowledged. Additionally the authors would like to thank Sven Schrape and Jonas Sarrar for many fruitful discussions and Atilay Tamkan (RRD), Jonas Richter, and Johannes Gerritzen (TUD) for their help on postprocessing the structural FEA results. The authors thank the Center for Information Services and High Performance Computing [Zentrum für Informationsdienste und Hochleistungsrechnen (ZIH)] at the TU Dresden for providing its facilities for the calculations presented in this paper.

• Bundesministerium für Wirtschaft und Energie and Rolls-Royce Deutschland (Grant No. FKZ 20T1716A; Funder ID: 10.13039/501100006360).

Latin and Greek Symbols

CoI =

coefficient of importance

d =

hail ball diameter

f =

probability density function

f_i =

ith eigenfrequency

F_T =

thrust

H =

discrete frequency

$nsim$ =

number of simulations

$N(μ,σ)$ =

normal distribution

$Q0.5$ =

median of the frequency bands

r =

hail ball radius

$r̃$ =

rank correlation coefficient

$R2$ =

coefficient of determination

s_c =

cup size

$U(xmin,xmax)$ =

uniform distribution

x =

probabilistic input parameter

y =

result quantity

$y$ =

vector of result quantities

z =

impact height

η =

isentropic efficiency

ω =

angular speed

ϱ =

density

Π_t =

total pressure ratio

Subscripts, Superscripts, and Accents

$(..)hb$ =

hail ball

$(..)∗$ =

normed

$||..||2$ =

$ℓ2$ norm

CV =

cross validation

Abbreviations and Acronyms

BL =

baseline

Blisk =

blade integrated disk

CFD =

computational fluid dynamics

COD =

coefficient of determination

eLHS =

extended Latin hypercube sampling

FEA =

finite element analysis

HPC =

high performance computing

MAC =

modal assurance criterion

MCS =

Monte Carlo simulation

PDF =

probability density function

RANS =

Reynolds-averaged Navier–Stokes

SPH =

smooth particle hydrodynamics