Nomenclature
- CROR
-
contra rotating open rotor
- TSFC
-
thrust specific fuel consumption
- NLH
-
non-linear harmonics
- RANS
-
Reynolds Averaged Navier Stokes
- u-RANS
-
unsteady-Reynolds Averaged Navier Stokes
- LES
-
large eddy simulation
- SST
-
shear stress transport
- TVI
-
tip vortex interaction
- PC
-
probabilistic collocation
- CFD
-
computational fluid dynamics
- BPF
-
blade passing frequency
- SPL
-
sound pressure level
- OASPL
-
overall sound pressure level
- TR
-
torque ratio
- BL
-
blade loading
List of symbols
- ω1
-
rotational speed, front rotor
- ω2
-
rotational speed, rear rotor
- F x,1
-
axial force, front rotor
- F x,2
-
axial force, rear rotor
- M
-
mach number
- D e,1
-
external diameter, front rotor
- D e,2
-
external diameter, rear rotor
- D i,1
-
internal diameter, front rotor
- D i,2
-
internal diameter, rear rotor
1.0 Introduction
Current developments suggest that the aviation industry, responsible for 3.8% of human induced CO2 emissions in the EU [1], shall drastically reduce its’ environmental footprint. Specifically, the EU-funded Clean Sky 2 programme aims at reducing CO2, NOx and noise emissions per seat kilometre by 20–30% until 2035, through the implementation of new engine technologies. Whilst working towards the goals, industry has been looking toward novel aircraft engine concepts known as prop-fans or unducted fans, first researched in the 1980s. Depending on the position of the thrust producing fan row, these engines are classified as tractors and pushers [Reference Protopapadakis, Gkoutzamanis, Vouros, Kyprianidis and Kalfas2]. A tractor variant is well known as the contra rotating open rotor (CROR) and has received considerable research effort.
The CROR’s unducted architecture allows for a higher bypass ratio compared to conventional turbofans engine. Moreover, the second contra rotating blade row recovers remnant swirl aft of the first blade row. These two characteristics of a CROR lead to a step increase in propulsive efficiency and therefore a direct reduction in fuel consumption. A future CROR (entry into service 2030 and beyond) is expected to be 20% more efficient compared to a current state- of-the-art turbofan [3]. When compared to a geared turbofan and a direct drive turbofan as presented in the NASA Advanced Turbofan Study [Reference Hendricks and Tong4] an open-rotor engine can have 12.1% and 16% lower TSCF, respectively. However, one of the major challenges faced during the deployment of CROR is the increased noise emissions. The absence of the nacelle and the flow interactions between the two rotors lead to increased noise.
Therefore, it is an imperative need to mitigate the emitted noise to achieve their certification and decrease their impact on communities. Predicting the acoustic footprint of any aero engine is a challenging task requiring complex simulations and experiments, while the tools necessary are not always publicly available. Therefore, the first aim of this work is to develop a comprehensive methodology to study the aeroacoustics of a CROR with the tools available to the authors. The second, higher-level goal is to use the tools to assess the CROR’s sensitivity to variation in various operational parameters. Such an assessment is, up to now, not available in public literature and considering the criticality of this question, it was deemed appropriate to tackle it.
1.1 State of the art
After the reconsideration of the CROR concept in the recent years there has been a boom in literature regarding these engines. A large proportion of it is focused on the aerodynamic and aeroacoustic simulation of these engines and to experimental investigations. Colin et al. [Reference Colin, Blanc, Caruelle, Barrois and Djordjevic5–Reference Colin, Caruelle and Parry7] have published a series of three papers that are a go-to reference regarding the computational tools necessary for CRORs. They have compared various simulation methodologies, mesh types, turbulence models and have investigated the effect of those choices on the aeroacoustics. Research is also directed towards studying the installation effects of pylon-mounted configurations on the engine’s performance [Reference Fiore, Daroukh and Montagnac8].
Numerous works have delved into optimisation techniques aimed at enhancing the aeroacoustics and aerodynamics of propellers or other components of aeroengines [Reference Jaron9]. Recently, optimisation efforts, encompassing both aerodynamic and structural blade design, have been integrated with uncertainty quantification to yield a more robust design, better suited to diverse environmental conditions and less susceptible to manufacturing tolerances [Reference Hirsch, Wunsch, Szumbarski, Łaniewski-Wołłk and Pons-Prats10]. The significance of uncertainties in the optimisation process has also been underscored [Reference Mourousias, García-Gutiérrez, Malim, Domínguez Fernández, Marinus and Runacres11, Reference Mishra, Mukhopadhaya, Alonso and Iaccarino12].
In aerospace applications, there is a growing utilisation of computational methods that explicitly consider uncertainties. However, literature focused on the application of uncertainty quantification and optimisation in open rotors is limited whether for aerodynamics or aeroacoustics. Boisard et al. [Reference Boisard, Falissard, Chelius and Mauffrey13] have researched the impact of blade deformation due to operational loads and manufacturing tolerances on aerodynamics and acoustics. The authors constructed computational models for the blades based on strain gauges measurements. The research concluded that deviations up to 10% for noise and aerodynamics can be expected when the operational shape of the blade is considered in the simulations. Khurana et al. [Reference Khurana, El din and Yeo14] have also applied uncertainty quantification on the aeromechanics of the Onera 7A rotor, a four-bladed rotor of high solidity. They investigated the impact of uncertainties in blade stiffness on the performance of the rotor with computational methods. With uncertainties present, the output variables of interest lied within a 97% interval compared to the baseline case.
Therefore, it is evident that when considering the emerging uncertainties substantial performance differences are to be expected. Towards this, it is essential early in the design and performance analysis process to perform sensitivity analyses, to determine the impact of selected variables. Given the above, the first part of this work presents the steps of the developed methodology spanning the geometrical, aerodynamic and aeroacoustic characteristics of the selected engine configuration. Then, the results are presented and discussed.including a grid independence study, as well as the overall sensitivity analysis. Finally, the conclusions are provided.
2.0 Methodology
Noise is generated by pressure disturbances propagating in a medium. Therefore, to numerically derive the noise emissions from an engine, it is necessary to obtain the pressure field created by the rotation of the blades. This can be performed with large eddy simulation or Lattice Boltzmann solutions. However, these simulations are very expensive, especially when the sensitivity of the design with respect to operating conditions or manufacturing tolerances is to be studied. In that case, perturbations of the pressure field close to the blades are propagated to the far field using aeroacoustic solvers. Thus, computational fluid dynamics (CFD) and computational aeroacoustics are intertwined when investigating the aerodynamic and aeroacoustics of open rotors.
2.1 CROR geometry
The considered geometry is derived from a CROR design consisting of 14 front and 14 aft blades (14×14), with the nacelle geometry being representative of a mature open rotor configuration. The blades are sophisticated with high sweep and lean angles. Also, the aft rotor is clipped to avoid the tip vortex interaction (TVI) phenomenon and reduce noise emissions. An overview of the geometry can be seen in Fig. 1.
The original configuration had rotor diameters that would make numerical simulation difficult, therefore the geometry was scaled down by a factor of 6.35. To preserve the operational characteristics of the engine (tip speed, aeroacoustic footprint) the rotational speed of both rotors was increased by a factor of 6.35.
2.2 Aerodynamic analysis
2.2.1 Simulation type
Before proceeding to the aeroacoustic calculations, it is necessary to simulate the aerodynamic field around the engine. In published literature, three methods have been implemented: full annulus, passage, and non-linear harmonics. The full annulus method includes every single blade from every single blade row and is the most thorough method since it captures most unsteady phenomena along with incidence and installation effects. For example, Ref. [Reference Brandvik, Hall and Parry15] implemented the unsteady, full-annulus CFD simulation to quantify the effect of the engine’s angle-of-attack on the rotor aerodynamics. The passage approach takes advantage of the inherent periodicity in turbomachinery flows and therefore only a single passage of each bladerow is simulated. This leads to a step decrease in computational time but also ignores various unsteady phenomena due to the imposition of periodicity. Finally, non-linear harmonics is a method is a frequency-domain based methodology which has been successfully applied by Kritikos et al. [Reference Kritikos, Giordano, Kalfas and Tantot16] on CRORs.
For an isolated configuration, the differences between the full annulus and passage simulations are negligible, both in aerodynamic and aeroacoustic terms [Reference Colin17]. Considering that this work focuses on estimating the sensitivity of the noise emissions with respect the operating condition perturbations, a passage analysis was deemed adequate.
2.2.2 Solution methodology and parameters
Regarding the numerical solution scheme of the Navier-Stokes, u-RANS was deemed the most appropriate approach, having produced results in good agreement with experiments, as indicated by Refs. [Reference Smith, Filippone and Barakos18, Reference Soulat, Kernemp, Sanjose, Moreau and Fernando19]. The alternative solution method, LES, is more accurate but significantly costlier since the larger eddies are solved numerically rather than modelled. Not implementing LES leads to dropping of the turbulent noise source term, but considering the presented case is during take-off (low speed), the contribution of this noise source is expected to be low.
ANSYS CFX was chosen as the solver and every simulation was performed on the high-performance computing infrastructure of AUTH. Before commencing with the time-dependent solution is to initialise it with a steady state solution. In the steady state solution, the rotating domains do not move but the rotational effects are considered in the Navier-Stokes equations. For the transient simulations the following choices have been made, which are based on the theory from CFX manual [20]:
-
1. Advection scheme: High resolution, which applies non-linear methods to calculate the gradients more accurately at the nodes.
-
2. Transient scheme: Second order backward Euler for the discretisation of the transient term. Second order backward Euler is an implicit time-stepping solution scheme.
-
3. Turbulence numerics: High resolution since it increases accuracy without affecting the computational cost.
-
4. Energy equation: Energy equation is also solved, to assess the impact of temperature and compressibility in the flow field.
-
5. Convergence criteria: The convergence of the simulations is based on the MAX residuals and a convergence target of 10−3. MAX reflects the maximum residual for every iteration at any specific point of the computational domain. Furthermore, the torque and loading values were monitored and the simulation was considered converged once their values had a periodic behaviour.
One of the most influential parameters in the transient simulations is the choice of timestep. Timestep is defining the time discretisation during one revolution of the rotor. The timestep must be small enough to accurately capture the phenomena but not too small to require many iterations for one revolution to complete. The choice of second order backward method allows for bigger timesteps since it is numerically stable and of $O( {\Delta {t^2}} )$ accuracy. The solution is implicit therefore, the Courant number is not important for the numerical stability. The timestep calculation was based on the rotational speed of the faster rotor since for the same dt it results to the largest dθ. The target dθ = 0.22° leads to a timestep of dt = 8×10−6 (s). Normalising the timestep with the period of the slowest rotor yields a $d{t_{norm}} = \frac{{dt}}{{{T_{rotor,{\rm{\;}}fast}}}}=6.11 \times {10^{ - 4}}$ with a full revolution of the slowest rotor requiring 1,664 timesteps. The choice of dθ was based on the work of Boissard et al. [Reference Boisard, Falissard, Chelius and Mauffrey13] where they implemented a dθ = 0.18°.
To complete the simulation setup, a decision on the turbulence model was necessary. In this work, Menter shear stress transport (SST or k-ω SST) [1], was implemented. Moreover, the model was used with an automated wall function which decides on whether the low or high Reynolds variation of the turbulence model will be used, based on the solver y+. The choice of turbulence model was based on the works of Refs [Reference Colin, Blanc, Caruelle, Barrois and Djordjevic5, Reference Mieloszyk, Galiński and Piechna21], along with the inherent advantages for the SST model with a wall function for simulations where adverse pressure gradients with separation and reattachment are expected. Reference [Reference Mieloszyk, Galiński and Piechna21] tested the SST, Spalart-Allmaras (S-A), k-ω and Reynolds turbulence models against experimental data for a contra-rotating propeller. It was concluded that the S-A and SST models are the most efficient, with the SST model also leading to lower residuals compared to S-A. Colin et al. [Reference Colin, Blanc, Caruelle, Barrois and Djordjevic5] proved that Spalart-Allmaras, k-ω 2 equation, and SST models produce aerodynamic and aeroacoustic results with small discrepancies, when compared to one another, with a highly refined grid (30 million elements). Additionally, Ref. [Reference Tormen, Giannattasio, Zanon, Kühnelt and De Gennaro22] applied the SST model in ANSYS Fluent for unsteady simulation RANS of a CROR producing results in good agreement with experimental data. Overall, literature suggests the SST model is expected work well in open rotors, therefore making it a suitable choice for the presented case.
2.2.3 Computational domain and mesh
The computational domain for CRORs resembles those implemented for contra-rotating propellers, with two adjacent rotating domains to cover the blades. The computational domain is also expanded radially, axially upstream and downstream of the rotating domain, resulting in a large simulation domain. In Fig. 2 the computational domains implemented for the CFD simulation are presented. The two stationary domains are indicated as S1 and S2, while R1 and R2 are rotating with the rotational speed of the rotors. The stationary domains serve as the far field, so that the rotating domains and the subsequent calculations are not influenced by the freestream conditions.
With the multitude of interfaces generated between the computational domains, boundary conditions for the correct transition were defined. For the steady state simulation, the frozen rotor method of CFX [20] was implemented as it is computationally cheaper compared to other frame transition models. Frozen rotor is freezing the rotation of the rotor so both domains have a fixed relative position throughout the simulation. For the transient solution the relative position of the rotating domains changes and therefore the mesh is in that area is non-conforming. Therefore, a different interface model was utilised, the transient rotor stator (TRS) [20]. This is a common interface model for transient simulations which creates a new artificial interface, where the mesh is conforming, and the flow quantities are interpolated. For blade rows with a non-unity pitch ratio, TRS model can disturb the flow quantities in the interpolation interface. Since the blade rows of this case have a unity pitch ratio (14×14 configuration), this effect is expected to be minimised.
Meshing was one of the most challenging tasks during this research due to the unavailability of public tools that can effectively mesh CRORs. For this work, ANSA by BETA CAE systems was used since it allows for the creation of conformal grids between the rotating and stationary domains. Due to the different blade geometrical features, the meshes in the R1 and R2 domains are different and therefore non-conformal in this interface between the blades. Some numerical dissipation is expected to arise from the non-conformality in this region. Finally, all the domains were meshed with structured hexahedral elements. To test the mesh quality, quality criteria have been set in the CAE software. In particular, the minimum and maximum angles were set to 15° and 165° correspondingly. This range is set by empirical evidence and was set after consultation with industry professionals.
After careful consideration, it was decided that the blades should be configured with a tip gap. This decision was made because the first option would have necessitated three overlapping boundary conditions in the same area, which would have led to unrealistic or unphysical outcomes.
2.3 Aeroacoustic simulation
The next step in the computational setup is the aeroacoustic simulation. In this section, the aeroacoustic solution process along with the definition of the FW-H surfaces will be presented.
2.3.1 Aeroacoustic solver
To calculate the noise perceived by an observer located in the far-field solving the entirety of the domain to obtain the pressure fluctuations is computationally prohibitive. Therefore, the equation developed by Ffowcs Williams and Hawkings (FW-H) [Reference Williams and Hawkings23], as an extension to Lighthill’s acoustic analogy to include the effects of moving sources, is an appropriate tool. Utilising the FW-H equation is common practice in aeroacoustic simulations of open rotors. For this research work the aeroacoustic solver developed by Prof. Ghader Ghorbaniasl of Vrije Universiteit Brussels was implemented [Reference Ghorbaniasl and Hirsch24]. The code solves the FW-H equations using the Farassat 1A formulation. A moving medium formulation, developed by Ghorbaniasl and Lacor [Reference Ghorbaniasl and Lacor25] can be implemented to account for the effects of moving air on a stationary observer since the relative movement changes sound refraction and attenuation.
2.3.2 FW-H surfaces and noise sources
To solve the FW-H equations it is necessary to define the surfaces over which the numerical integration of the solution parameters will be performed. There are two types of such surfaces that can be elected: impermeable and permeable. Impermeable surfaces correspond to physical surfaces with no mass flux through them. When an impermeable surface is elected, the acoustic quadrupole term is neglected. On the contrary, permeable surfaces allow for mass flux through them and can be arbitrarily defined without necessarily corresponding to a physical surface. Choosing the type of FW-H surface is not bound by any standard process and is dependent on the scenario that is studied. For example, Stuermer and Yin [Reference Stuermer and Yin26] concluded that using the impermeable surface approach can give up to 10dB lower SPL across specific directivities due to the dropping of the quadrupole sources during a take-off scenario.
Since the solver handles both type of surfaces, the permeable approach was implemented, considering it allows for capturing more noise sources. The elected permeable surfaces enclose both rotors and allow for the direct capturing of all the noise sources (steady and unsteady) arising from the interaction of the front and aft rotors’ flow fields. The variables over the surfaces were considered in the stationary frame of reference. Finally, since the chosen surfaces were already part of the computational domain no extra effort for defining and meshing new surfaces were required. Finally, the choice of FW-H surfaces is coherent with the work of Deconinck et al. [Reference Deconinck, Capron, Hirsch and Ghorbaniasl27].
The selected surfaces can potentially capture every noise source of the engine. In this case, the limiting factor to the multitude of noise sources captured is the Navier Stokes solution method. By solving the u-RANS equations rather than applying LES broadband noise sources cannot be captured. Additionally, because no shock waves were identified during the aerodynamic analysis, no quadrupole sources are expected. Consequently, only monopole and dipole sources will be captured, and the solver configuration is set up accordingly.
2.3.3 Aeroacoustic simulation details
The intermediate step between the aerodynamic and aeroacoustic analysis, is the post processing of the aerodynamic results. The solution primitive variables are formatted accordingly before they can be inputted in the acoustic solver. In the following paragraphs, the methodology will be presented.
Firstly, the necessary solution variables and their sampling interval shall be defined. The necessary variables for the acoustic solver are pressure, density and the velocity components in the stationary frame of reference. These primitive values are extracted from the surfaces defined in Section 2.3.2. Regarding the processing to acquire the aeroacoustic results, a transient result file with the aforementioned values was exported every second iteration, corresponding to a Δt = 1.6 10−5s. Essentially, this timestep dictates the sampling period and frequency ${T_s}=1.6{\rm{\;}}{10^{ - 5}}{\rm{s\;}},{\rm{\;}}{f_s}=62.5{\rm{kHz}}$ accordingly.
According to Nyquist’s sampling theorem to capture a phenomenon of frequency ${f_N}$ the sampling frequency must be at least double of the phenomenon frequency or ${f_s}>2{f_N}$ . For open rotors the maximum frequency usually of practical interest is up to ${f_N}=5\left( {BPF + BPR} \right)$ . For a setup where both rotors have 14 blades, the front rotor rates at 4942RPM and the rear at 4,509RPM the maximum frequency is ${f_N}=11.026{\rm{kHz}}$ . Therefore, the proposed frequency is more than enough to capture the investigated phenomena.
2.4 Sensitivity analysis
Analysing the sensitivity of the open-rotor engine’s performance to changes in operating parameters is crucial to understand its behaviour beyond just the design point. This analysis helps assess how the engine’s performance fluctuates during operation, particularly during take-off. Furthermore, it offers insights into the model’s robustness by evaluating how well it adapts to variations in operating conditions.
2.4.1 Sensitivity analysis parameters
The decision for the parameters utilised for the sensitivity analysis was based on two factors. Firstly, the capacity of the computational model to account for these parameters was considered. Secondly, the relevance of the operational parameters on the performance of the engine were reviewed.
Due to computational restrictions three parameters were analysed: total pressure and temperature at the inlet and the rotational speed of the front rotor. Total pressure at the inlet was deemed more appropriate than the static pressure because it is already used as a boundary condition in the CFD. To avoid unphysical results due to lower pressure at the inlet compared to the outlet, only the dynamic fraction of the inlet total pressure was varied, while keeping the outlet pressure constant. Variation of the total temperature mainly affects the density of the air injected by the engine at the inlet and is therefore a crucial parameter on the performance of the CROR. The third parameter subject to variation and simultaneously impactful on the CROR’s operation, is the rotational speed of the two contra rotating blade rows. Variations can occur in the RPM since in a real-life scenario it is a controlled to maintain some other quantity constant, for example thrust. It is worth noting that this sensitivity study focuses on takeoff conditions, where noise emissions are particularly critical. Even a slight variation or divergence in operating conditions could result in noise production exceeding certification limits, posing significant challenges for manufacturers and operators.
Following the identification of the parameters, the variation range was defined. This was done rather arbitrarily and based on engineering judgement. Total pressure and temperature were varied by ±5% around their nominal values, while the front rotor’s RPM was varied by ±2.5%. The parameters nominal value and variation range are presented in Table 1 and have been non-dimensionalised with respect to their nominal value which is set to 1.
2.4.2 Sensitivity analysis methodology
For this research a local sensitivity analysis around the operating point of the engine during take-off is performed. To better understand how each operating parameter affects the aerodynamics and aeroacoustics of the engine the One - At - a -Time method was elected. With this method only one parameter is varied at a time and the main effects of each variation can be observed.
2.4.3 Evaluation parameters
To address the main goal of analysing the design’s sensitivity to small variations in operational parameters, it is essential to select appropriate output variables for evaluating the model’s performance. Hence, it is necessary to designate suitable aerodynamic and aeroacoustic output variables that accurately depict variations in the engine’s response.
Concerning the aerodynamics, torque ratio and blade loading will be the monitored values. This choice was based on the work by Kritikos et al. [Reference Kritikos, Giordano, Kalfas and Tantot16] to use these parameters as surrogates for the aerodynamic performance of the engine. Apart from being in direct correlation to the engine’s performance these values can serve as predictors for the noise emissions. The definition of the variables is presented in Equations (1) and (2).
For the aeroacoustics the overall sound pressure level (OASPL) was the metric of interest. OASPL deems itself appropriate since it is a single variable that can quantify the acoustic footprint of the engine. This is the reason why it is often used as the output variable for aeroacoustic sensitivity analyses or optimisation processes [Reference Wang, Diskin, Lopes, Nielsen, Lee-Rausch and Biedron28].
3.0 Results
3.1 Meshing and grid dependency study
In this section, the computational mesh used to resolve the flowflied is presented. The mesh consisted of 5,462,720 volumetric elements with their distribution across the domains being presented in Table 2. R1 and R2 have been more densely meshed since they enclose the blades. Consequently, a higher resolution is required to resolve all the phenomena. In addition, R2 has 200k more elements than R1 due to the existence of the wake of the front rotor (FR). The downstream movement of the FR wake and tip vortices create a more complex flow field in R2. To account for the added complexity of the flow, more volumetric elements were added.
The two rotor blades have been meshed with an O-grid topology (Fig. 3). This decision was made based on the work of Colin et al. [Reference Colin, Blanc, Caruelle, Barrois and Djordjevic5] who implemented the O-grid approach, after testing various mesh topologies, and deemed it the most appropriate. Schnell et al. [Reference Schnell, Yin, Voss and Nicke29] have also implemented O-grids for meshing open rotor blades with 2,9x106 elements. Regarding the placement of the outer far field, even when the outer boundary is placed nine radii away, the solution is still dependent to the distance according to Ref. [Reference Mieloszyk, Galiński and Piechna21]. In that work, it was indicated that five radii is a suitable distance to have minimal impact, reducing however the computational requirements. For this work, the far field is placed 3.3 radii of the blade tip to keep the computational domain at a manageable size.
With the mesh generated, its quality can be evaluated with the metrics described in Section 2.2.3. The quality evaluation indicated that only 250 elements (0.0045% of the total) are quality violating. In addition, no negative volume elements were observed. Given the complexity of the geometry and the manual creation of the grid, this small deviation from the set range is acceptable.
Before deciding on the total number of elements for the CFD cases, a grid dependency study was carried out. Torque ratio was the monitored value for the grid dependency since it is a non-dimensional variable. The mesh dependency was investigated through a steady state CFD, in a range of 3–12 m elements. It can be concluded that up to 12 m hexahedral elements the study is dependent of the mesh count as seen in Fig. 4. However, the deviations were small between the different scenarios (ΔTR<1%), as evident from the torque ratio comparison between the cases of Fig. 4, and since the sensitivities and not the absolute values were of importance, a mesh size of 5 m elements was elected. All in all, this choice represents a compromise between the accuracy and computational time.
3.2 Baseline case
The nominal or baseline case is representative for a CROR with 14 × 14 blade configuration during a take-off scenario. The most critical parameters are presented in Table 3.
Clipping refers to the reduction of the RR diameter to avoid the TVI from the FR. TVI is a phenomenon which additionally loads the tip of the rear blade and therefore leads to increased noise emissions. Clipping is defined in Equation (3).
3.3 Nominal case: aerodynamic analysis
For monitoring convergence, both RMS (root mean squared) and MAX residuals were assessed. The RMS residuals are an indicative metric for the overall quality of the simulation since they consider the average quality of the solution. On the other hand, the MAX provide information regarding the grid element with the highest residual value. Regarding the former, all values are below 1E-4 and therefore the overall mesh quality is sufficient. The mass residual, derived from the continuum equation shows the best convergence with values below 1E-6.
As expected, the MAX values are on average higher than the RMS, but still within acceptable limits. During the post processing, the location of the maximum residuals has been identified in the tip region of the FR and RR blades and can be attributed to a sharp mesh angle. Overall, all observed values have achieved periodicity, and the values are acceptable.
Regarding the two physical parameters monitored for convergence, torque ratio and blade loading, both achieved periodicity and thus convergence. The torque ratio converges after the 5,000th iteration while blade loading convergences earlier, at around 3,500 iterations.
3.3.1 Flowfield analysis
In this section fundamental flow properties and phenomena derived from the aerodynamic simulation of the investigated configuration will be discussed. First, the Mach number in the rotating frame of reference is presented at three different radial slices for both rotors in Figs. 5–7. The first observation is that the flow is always subsonic but compressible in the range of Mach 0.2–0.8. The flow is subsonic since the scenario is during a takeoff, where the freestream velocity is M = 0.2. Moreover, a non-uniform velocity profile at the inlet of the FR passage is evident, as a direct impact of the blocking that the boundary layer on the pressure side of the passage induces.
Additionally, the wake downstream of the rotors can be observed in the blade-to-blade Mach representations. Once the wake passes from the interrotor plane, it dissipates gradually. This dissipation is characterised more as a numerical phenomenon due to the difficulty to interpolate the velocity values across the mesh interface. Namely, the mesh lines do not align before and after the interface creating the wake dissipation. This phenomenon is observed in plenty CROR numerical studies, especially when coarser meshes is used, for example as in Ref. [Reference Protopapadakis, Gkoutzamanis, Vouros, Kyprianidis and Kalfas2]. The presence of numerical dissipation in the interface is expected to introduce inaccuracies in the results and could potentially affect both the aerodynamics and aeroacoustic sensitivities. Evidently, the wake of the front rotor influences the inflow of the aft rotor with an area of lower velocity and disturbed flow, like what has been noted by Ref. [Reference Stuermer and Yin26]. Overall, the influence of the front rotor’s wake is more pronounced near the tip area as indicated in Fig. 7 where it influences the entire inflow area of the aft rotor. This interaction leads to higher loading of the aft rotor and potentially higher noise emissions. A solution to decrease the wake induced loading of the aft rotor would be an increased spacing.
The wake effects can also be observed in the meridional views of Figs. 8, 9 where the Mach number contours are presented on the interstage plane and on a plane aft of the RR correspondingly. In Fig. 8 the effect of the FR blade is evident as a low-speed area that is shifted with the rotational direction. The deficit is most prominent in the blade area with the maximum LE curvature. Additionally, the effect of the tip vortex can be observed in the tip region as an area of lower speed, corresponding to the vortex core. Moving radially upwards the Mach number increases until the outer boundary of the surface, where the rotational domain meets the stationary. The overall increase of the velocity with the radius is explained by the higher tip speeds of the domain. The flowfield aft of the RR (Fig. 9) exhibits similar trends to the one of the FR, with a wake area of higher speed. However, significant differences can be seen above the tip area, with a highly disturbed flowfield. A larger area affected by the tip vortex of the downstream rotor can be observed with strong velocity gradients, indicating interaction with the vortex coming from the FR. The interaction of the tip vortices is also validated through the entropy of the flow in the examined regions, which is further explained in the following paragraph.
Analysis of the entropy generated by the powerplant can provide crucial insights to the loss mechanisms of the flow and the interactions that happen. It is also commonly used in open rotor literature to visualise the effects of the tip vortices, as suggested by Ref. [Reference Hall, Zachariadis, Brandvik and Sohoni30]. Figure 10 presents the entropy contours in the plane between the two rotors and in the plane aft of the rear rotor, where the Mach was previously analysed. In the interstage, the effects of the blade wake can be observed around the radial midspan location, where an area with increased entropy is present. This area coincides with the blade area that has the highest curvature. Moreover, the tip vortex from the front rotor can be seen in the tip region, as an inner circle with a higher static entropy which corresponds to the higher vorticity. Additionally, due to clipping of the aft rotor, the impingement of the front rotor’s tip vortex is mitigated as evident by the location of the elevated entropy area higher than the downstream rotor’s tip. Interactions with the front rotor tip vortex are still expected, although difficult to quantify which parts are associated with TVI and which with wake interactions. Due to the rotation of the domains, the vortical structures curl toward the rotational direction and this effect the impingement location, especially for the aft rotor, as also indicated by Ref. [Reference Hall, Zachariadis, Brandvik and Sohoni30].
Regarding the entropy generation downstream of R2, a similar pattern to R1 arises with two key differences. Firstly, the wake aft of the second bladerow has a larger area of higher entropy compared to the upstream rotor. This effect arises from the interaction with the wake of the upstream rotor. Secondly, the tip vortex of the second blade row has a higher core entropy indicating a stronger vortical structure. This is an expected phenomenon due to the distorted inflow at the input of the aft bladerow. As noted by Ref. [Reference Smith, Filippone and Barakos18], the aft tip vortex can be stronger compared to the fore and has significant contribution to the noise emissions of the powerplant.
3.4 Nominal case: aeroacoustic analysis
With the results of the aerodynamic analysis and a converged solution, the aeroacoustic analysis of the baseline case can be presented. This section aims at thoroughly presenting and analysing the aeroacoustic results obtained from the baseline case.
3.4.1 Microphone locations
Every noise source is directive, which means that the sound emitted towards different directions is different. Therefore, to capture the noise emitted by the CROR it is necessary to define microphones across different locations and directions, with each microphone representing an observer. According to the Noise Certification Workshop by ICAO [Reference Boettcher31] the measurements point necessary for the certification of any aero engine are both on the lateral and longitudinal axis. Although the investigation of the certification aspects is outside the scope of this work, it was decided to follow the reasoning presented by ICAO for the positioning of the microphones. Overall, microphones have been defined in three different distance clusters, at, 100 and 500m. Following this strategy of clustering the computational microphones in various distances, the effect of distance on SPL and any potential frequency-shifting can be observed. The positioning and exact coordinates can be seen in Table 4 and Fig. 11.
3.4.2 Aeroacoustic results
In this section, the results of the nominal case will be presented in detail only for the 500m microphone cluster. Noise generation mechanisms for open rotors are very similar to propellers, especially during take-off where the tip speed is not at the transonic region. Propellers have two main mechanisms generating tonal noise, thickness noise and loading noise [Reference Hubbard32]. The former is caused by the displacement of air by the blade, while the latter is caused by pressure fluctuations produced by accelerating forces acting on the surface of the blade due to lift and drag. Broadband noise can be generated from vortex shedding and generally the wake of the propellers. In this case, where the engine consists of two blade rows, interactions between them are expected to be prominent. The noise produced by the periodic interaction with a frequency equivalent to that of the count of the blades (BPF1 + BPF2) as they rotate in a counter direction is determined as unsteady loading noise [Reference Smith, Filippone and Barakos18].
At first the results for Microphone 1 are presented in Fig. 12(a). It is observed that the peak SPL value is found in the first interaction frequency BPF1 + BPF2 for this location. In frequencies lower and higher than the first interaction frequencies, the SPL drops rapidly indicating how strong this interaction is. The rotor alone tones BPF1, BPF2 are accompanied by a peak but their contribution is nowhere near that of BPF1 + BPF2. This is expected, as the unsteady loading generally has the biggest contribution(Reference Smith, Filippone and Barakos18). In higher frequencies, an SPL peak is observed in (BPF1 + 2BPF2), (2BPF1 + BPF2) region. After that peak a sudden drop in SPL is observed until the SPL rises again in 2(BPF1 + BPF2). For the higher frequencies, SPL exhibits more peaks, with 3(BPF1 + BPF2) being the most prominent.
Overall, for frequencies higher than 2(BPF1+ BPF2) the SPL magnitude starts to diminish since the strong aerodynamic interactions are concentrated in the lower frequency spectrum. Another reason for the behaviour of the SPL in the higher frequencies is the computational mesh used for CFD calculation [Reference Soulat, Kernemp, Sanjose, Moreau and Fernando19]. The aerodynamic mesh directly impacts the frequencies that can be captured by the aeroacoustic solver via a concept known as mesh frequency cut-off. According to the mesh frequency cut-off, there is a maximum solvable aeroacoustic frequency dictated by the CFD mesh. The mesh utilised is on the coarser side, to avoid overwhelming computational side and this impacts the maximum frequencies up to which SPLs can be accurately captured.
In general, it can be concluded that the SPL peaks coincide with the interaction frequencies as described in published literature such as by Guynn et al. [3]. The frequency spectrum analysis of the CROR indicates that the noise emitted is mostly due to the aerodynamic interactions between the two blade rows since the peaks coincide with the combinations of the fundamental BPFs.
What is different for every observer location is the SPL profile. For example, the microphone lying directly on the engine axis and behind the source, presented in Fig. 12(b), captures a radically different profile compared to Microphone 1. In this case, 2(BPF1 + BPF2) interaction is the main noise source.
It is interesting to discuss the results for the sideline microphone of Fig. 12(c). This microphone records three very distinct peaks at (BPF1 + BPF2), 2(BPF1 + BPF2), 3(BPF1 + BPF2) with continuously dropping SPL magnitude. Similar is the behaviour of the Mic. 4, in Fig. 12(d) which has only a z component. Finally, to investigate how the location of the microphone relative to the engine affects the noise captured, the results for Mic. 5 are displayed in Fig. 12(d). This microphone is upstream of the engine inlet contradictory to Mic. 1, which is downstream of it. The main difference for Mic. 5 stems from the 2(BPF1 + BPF2) interaction frequency, which has an SPL content as high as the BPF1 + BPF2. Moreover, the peak at BPF1 + BPF2 is 4dB lower than the Mic. 1.
To complete the aeroacoustic analysis, the overall sound pressure level (OASPL) for every microphone location is presented in Table 5. OASPL essentially weights the individual sound pressure levels across all computed frequencies, resulting to a single comparable quantity across all locations. For the 500m cluster, the highest OASPL is perceived by the observer at Mic. 2 location. All the microphones that are aft of the engine inlet, have similar OASPL values of around 180dB. The microphone front of the engine inlet (#5) has a lower OASPL value due to the lower magnitude of the SPL corresponding to the (BPF1 + BPF2) frequency. Moving closer to the source, the acoustic impact of the engine is stronger, and this can be concluded through the OASPL of the microphones in the 100m cluster, with the quietest location still louder than the loudest observer of the 500m. In the 100m the highest sound magnitude is at Mic. 7, which is like Mic. 2. In the final 10m cluster, Mic. 14 is the loudest.
To sum up, the microphones which lie directly on the engine axis and those who have only a z coordinate perceive the loudest sounds for every cluster of microphones. Microphones 1, 6, 11 are the farthest from the source for the corresponding clusters, and thus perceive the lowest OASPL magnitude from those aft of the engine inlet. The microphones fore of the engine inlet perceive lower SPL due to the directivity effects already described. The average OASPL for the 500m is 180.28dB, for the 100m 194.08dB and for the 10m 213.41dB.
To conclude, the SPL frequency diagrams for 5 microphone locations in the 500m region have been analysed. It has been concluded that the main interaction tones are well captured by the proposed methodologies. For most case the (BPF1 + BPF2) is the frequency leading to the highest SPL value. Also, the importance of the microphone location on the perceived SPL has been established. Finally, the three different cluster of microphones have been compared and it has been concluded that the close the microphone to the source the higher the SPL.
3.5 Sensitivity analysis
This section is dedicated to the in-depth analysis of the model’s response when uncertainties are present for various operating parameters. Briefly, the uncertainties have been modelled as distinct cases for each of the two extremes of the statistical distribution described in Section 2.4.1. Overall, six scenarios were simulated, two for each uncertain value. The scenarios will be named after the parameter that was varied and the deviation from the baseline values. For example, a scenario called “Total Pressure @ Inlet +5%” corresponds to a simulation where the total pressure at the inlet of the domain has been increased by 5% compared to the nominal value.
3.5.1 Aerodynamics
Following the steps of the modelling process, the impact of the considered uncertainties on the aerodynamics will be presented at first. To create an overview of the results, Table 6 is presented, where the values corresponding to the nominal case have been set equal to 1 and all the results have been non-dimensionalised with respect to the nominal case values.
Table 6 sums up the effects of the uncertainties on the monitored output variables, torque ratio and blade loading. Moreover, the values from each case are compared to the baseline case to check whether the uncertainties are damped or amplified. Amplification happens when the output of the model i.e., torque ratio and blade loading is varied more than the change of the model input. Namely for a 2.5% increase in the RPM of the front rotor, the TR is reduced by 7.7%. On the other hand, a 5% decrease in total inlet temperature is damped producing a 0.92% change in TR.
Important conclusions can be drawn from Table 6 regarding the sensitivity of the output variables and the physical relationships between the inputs and outputs. First, an increase of RPM leads to a reduction of TR and an increase of BL. Increasing the FR RPM by 2.5% will directly reduce the ratio of rotational speeds by 1/1.025. Moreover, the forces acting on the blades will be increased on the FR since the rotational speed is higher, in turn giving a lower forces ratio between the two rotors. The change of rotational speed and the higher axial forces on the first rotor, produce the reduction in TR. As already described, the blade loading is increased, due to the contribution of the higher loading of the first rotor. In this point, it must be reminded that blade loading value is the average for both rotors. For the decreased RPM simulation, the effects on TR and BL are the opposite. As a conclusion, the effect of RPM variation on the aerodynamic metrics is very prominent and is amplified through the computational model.
The second sensitivity studied is to the total inlet pressure. It is observed that an increase in total inlet pressure led to a slight reduction of TR and more considerable of BL. As the rotational speed ratio remains constant the only factor varied is the ratio between the axial forces between the aft and front rotor. The fraction of axial forces is decreased by 0.9% compared to the baseline case, but in general this is a slight reduction considering the 5% increase in total inlet pressure. As already mentioned, BL is more sensitive to uncertainties of the inlet pressure since it is reduced by 3.6%. Table 7 presents a comparison of the blade forces exerted on both rotors by the fluid, between the three different total pressure scenarios. It is evident that an increase in total inlet pressure leads to a reduction on forces acting on both blades. This behaviour is a direct consequence of the computational setup and boundary conditions. An increased total pressure at the inlet corresponds to a flow with higher potential energy and therefore the blades do not interact so heavily with the flow to meet the pressure requirements set at the outlet of the computational domain.
Finally, the effect of varying the total inlet temperature is analysed. From Table 6, it is evident that temperature has the most profound impact on blade loading. Increasing the total temperature at the inlet leads to a reduction in TR and blade loading. A 5% increase produces an 8.5% decrease in blade loading, while a 5% decrease leads to a value 10% higher than the baseline case. Therefore, the inlet uncertainty is amplified up to two times. TR is relatively insensitive to this variation, since the ratio of the forces exerted by the flow to the blades is almost constant and the rotational speeds remain unchanged. The cause of the amplification of the temperature change one can be attributed to the density distribution on the blade-to-blade view and the blade forces. When the temperature is increased, the density is reduced leading to lower blade forces and thus loading. This happens since the blades interact more lightly with the flow to achieve the pressure differential between the inlet and the outlet.
All in all, the aerodynamic performance is sensitive to the studied uncertainties since a variation of the total inlet temperature can have double the effect on the blade loading. Considering the criticality of blade loading design since it dictates both aerodynamic and structural performance of the engine an accurate variation range is necessary. On the other hand, torque ratio is not highly variable since it is a function of the RPM ratio, which can be set to constant (design choice) and the ratio of forces acting on the blades, which they result show to be almost constant.
Of course, the sensitivity of the aerodynamic model captured here is not necessarily representative of the real-life performance of the engine under such variations since the modelling assumptions have a high impact on the results. For example, the meshing and convergence influence the model’s accuracy and are two parameters which could be further improved. Overall, the qualitative and quantitative assessment of the design’s robustness has explored the relationships between the uncertain parameters and model outputs.
3.5.2 Aeroacoustics
This section aims to analyse the results from the aeroacoustic analysis and quantify the impact of the operational parameter’s uncertainties on the SPL and OASPL at 4 observer locations. Although the size of the engine is scaled down, as mentioned in the aerodynamics session, the rotating speed was adapted. Thus, the power of the noise generated is expected to be in normal scale, with the exception of the shift in frequency caused by the dynamic scale. Furthermore, as no absorption model is implemented, the amplitude should only be affected by the distance measured, while the distribution over the frequency should remain the same. The locations selected are clustered around 500m, as this is a reasonable distance for an actual observed being in the vicinity of an airport.
Firstly, the effect of varying the RPM of the front rotor is investigated. The expected outcome would be an increase of loading noise for increased RPMs. The reason is that loading is also expected to increase as shown in the previous paragraph. The same tendency is expected for the thickness noise, as with higher RPMs, a bigger air volume is expected to be displaced. The opposite behaviour would be expected for a reduction in RPMs.
From the SPL-frequency diagram in Fig. 13(a) for the corresponding case, it can be concluded that the SPL values in the interaction frequencies (BPF1 + BPF2), 2(BPF1 + BPF2) and 3(BPF1 + BPF2) are very close between the two cases for Microphone 1. Large discrepancies can be observed in the intermediate, non-dominant frequencies. This could be attributed to the CFD simulation. In the increased RPM case, the flow field is more turbulent and thus the existent mesh did not lead to a fully converged simulation, producing this tooth shape of the SPL function. Examination of Table 8 reveals the RPM variation does not have a large impact on the SPL perceived by an observer. Deviations of about 2.4dB are observed in location 2 for both cases in comparison to the nominal case with such differences considered minor.
It shall be noted that a reduction of the front rotor speed might reduce the noise generation for the front rotor but might increase the interaction noise in some cases. Such a case could be if, for example, the thickness of the boundary layer would increase, introducing more wake turbulence to the aft rotor. This would increase the unsteady loading noise of the aft rotor, and the overall noise level of the engine.
The presence of uncertainties on the total inlet pressure has demonstrated a non-negligible variation in the blade loading. The same does not stand for the aeroacoustic footprint of the CROR, as it would be expected. Reducing of blade loading should lead to lower loading noise. Table 9 indicates the OASPL levels for each microphone location with the largest divergence observed again in Microphone 2 (engine axis) at about 2.4dB when compared to the design point case. On every other location, the difference is below 1dB. Strikingly, although the increased pressure led to decreased blade loading and conversely, both uncertainties provided OASPL values that are larger compared to the nominal case. As presented in Fig. 14, both simulations yielded similar SPL distributions on the dominant interaction frequencies calculated the same SPL. Deviations are again observed in the intermediate frequencies, but they are not major enough to alter the OASPL to the nominal case. It can be concluded that the variation in the total inlet pressure and the subsequent variations in the aerodynamic performance of the rotors are damped by the model and thus not translated to any outstanding difference in the aeroacoustics.
In the aerodynamic analysis, the total inlet temperature was the variable with the most prominent change to the outputs of interest. The total temperature change is affecting mainly the density, and thus the loading as shown in the previous section. This effect was not observed for the aeroacoustic analysis, where the variations compared to the nominal case were similar to the ones mentioned afore. Investigating Fig. 13, the increased temperature produces a 5dB higher SPL peak for the (BPF1 + BPF2) and 3(BPF1 + BPF2). For the 2(BPF1 + BPF2), both simulations yield the same SPL value. Large deviations between the + and – simulations are observed in the region between (3BPF1 + 2BPF2) and 3(BPF1 + BPF2). This is trivial to explain since in general the two distributions are very close, but it is hypothesised to arise from the mesh cutoff frequency. To conclude, the variation of total inlet temperature follows the same pattern as the other parameters and does not appear to have a significant impact on the OASPL for the various observer locations.
To conclude, the variations in the aerodynamics presented afore do not significantly impact the aeroacoustic footprint of the engine. Variations of up to 2.4dB in OASPL were observed for all uncertain inputs in the microphone lying aft of the engine inlet. The sensitivities in the different parameters seem to affect the microphones in similar manners across the different cases. Moreover, the SPL-frequency comparison of the various cases indicated that the SPL peaks corresponding to the interaction frequencies remain almost unaffected when the uncertainties are present.
4.0 Conclusions
The purpose of this paper was to numerically analyse the sensitivity of the aerodynamic and aeroacoustic performance of a novel open rotor engine, to small variations of three operational parameters. Namely, the impact of variations in the rotational speed of the front rotor, in the total pressure and temperature at the inlet.
For the aerodynamic approach. unsteady Reynolds Averaged Navier Stokes simulation was implemented with a single-passage transient solution. This allowed to export the time-varying pressure field to the aeroacoustic solver, which solves the Ffowcs Williams – Hawking equations. In this way, the acoustic footprint of the engine can be calculated in the far field in an efficient manner.
The small variation analysis around the design point for take-off showcased significant sensitivity. Specifically, a 2.5% decrease in RPM could lead to 9% increase in the torque ratio, while a 5% decrease of inlet temperature would produce a 10% increase in blade loading. These aerodynamic alterations could have a profound impact on the operation and controllability of a CROR and should be considered during the design phase.
The aeroacoustic response of the engine, proved more resilient to the considered variations. Despite considerable variations of the aerodynamic properties, the overall sound pressure level differed by a maximum 2.4dB for the cases of reducing RPM, total inlet pressure and temperature, with the maximum sensitivity found in the microphone aft of the engine inlet. Regarding the SPL-frequency distribution, there was no frequency shift between the different cases and the SPL peaks coinciding with the interaction frequencies were of the same magnitude. Discrepancies in the intermediate frequencies were observed but they do not affect the final OASPL values due to their non-dominant nature.
Acknowledgments
Firstly, authors would like to express their gratitude towards the team of BETA CAE Systems in Thessaloniki, especially Antonis Karasavvidis, for providing the pre-processing software and continuous support. Finally, the authors would like to acknowledge the support provided by the IT Centre of the Aristotle University of Thessaloniki throughout the progress of this research work.