Efficient Unsteady Simulations of Cooled Turbines using Harmonic Balance

17th International Symposium on Unsteady Aerodynamics Aeroacoustics and Aeroelasticity of Turbomachines ISUAAAT17 November 16-21, 2025, Melbourne, Australia ISUAAAT17-120 EFFICIENT UNSTEADY SIMULATIONS OF COOLED TURBINES USING HARMONIC BALANCE Michael M¨uller, Robin Sch¨ offler, Clemens Grunwitz, Christian Morsbach Institute of Propulsion Technology German Aerospace Center Cologne, Linder H¨ ohe, Germany ABSTRACT Efficient unsteady simulations of cooled turbines are essential for improving thermal load predictions and optimising turbine design. This study introduces a combined approach of the Harmonic Balance method, a differential Reynolds stress model, and a localised film cooling model, which captures the dominant unsteady effects in high-pressure turbine rotors. Traditional steady-state RANS simulations often underor over-predict adiabatic wall temperatures and heat transfer coefficients due to the neglect of unsteady interactions. Results show that unsteady effects significantly impact thermal loads, with frequency-domain harmonic balance simulations closely matching high-fidelity time-domain URANS results at a fraction of the computational expense. The localised film cooling model improves accuracy over a coarser volume source model, while the SSG/LRR-log(ω)turbulence model enhances turbulent mixing predictions. This can enable a much better capturing of relevant unsteady effects and turbulent mixing processes in both problem analysis and design, the impact of which would remain unknown using simpler methods. Keywords: Unsteady Aerodynamics, Cooled Turbines, Harmonic Balance, RANS Models 1 INTRODUCTION The drive for higher overall efficiency in civil aircraft propulsion has led to a significant rise of the turbine inlet temperature. Because the desired turbine inlet temperature is above the safe operating temperatures for the considered materials, active cooling is indispensable for meeting both efficiency and reliability targets [1]. But despite complex cooling used in modern jet engines, high thermal loads are still a critical factor governing reliability, maintenance intervals and, ultimately, the cost of ownership [2]. For rotors, this thermal load results from a complex, highly unsteady interaction between the hot main-stream, the cooling-air and the secondary flow structures. Ignoring these unsteady phenomena in a Computational Fluid Dynamics (CFD) simulation, as is customary in steady Reynolds-averaged Navier–Stokes (RANS) simulations, can lead to substantial underor over-prediction of wall temperatures and heat-transfer coefficients [3–5]. The present work aims to address the need for a computationally affordable yet sufficiently accurate method to capture the dominant unsteady effects in cooled turbine rotors. 1.1 Unsteady phenomena affecting thermal loads Several unsteady phenomena affect both the film cooling and the mixing and diffusion of the hot gas. When hot gas leaves the combustor it does not possess a uniform temperature field [6]. Instead, hot streaks are usually formed, which can drastically increase thermal loads. The transport and diffusion of those hot streaks is strongly affected by the positioning of the NGVs and the inlet turbulence [7]. At the same time, the blades generate a system of secondary flows that locally transport hot gas towards the blade and endwall surfaces and continuously disturb any protective cooling film [8, 9]. In addition, these vortices can also affect the hot streak. All of those phenomena are present in both stator and rotor rows. However, the flow is additionally unsteady in the rotor system. The hot streak, the wake and all vortices from the nozzle guide vanes (NGV) are deterministically unsteady phenomena, that also affect the thermal loads in the rotor. Some of the cooling films on the rotor will be excited to oscillate with each passing wake. This increases the spread of the coolant and can lead to a reduced time-averaged surface temperature in the affected regions. However, this oscillation also perturbs the boundary layer, 1 uhot gas Thot gas uwake Tcoolant Rotation u∞ urel,hot gas urel,wake Hot gas impacts pressure side Coolant impacts suction side FIGURE 1: Redistribution of thermal loads from the pressure side to the suction side through different incidence angles of coolant and hot gas causing an increase in the heat transfer coefficient (HTC) which can offset the reduction in average temperature. The most important unsteady effect for this paper is an uneven thermal loading on the pressure and suction sides of the rotors [10]. This effect is caused by a variation in the incidence in the rotor system [11] due to a speed difference between the cooling fluid and the hot gas and is shown schematically in Fig. 1. The hot gas leaving the combustion chamber is accelerated in the NGV passage, exiting at velocity uhot gas. However, as the NGV is film cooled, a cooler wake also leaves the NGV row. This wake has a similar flow direction, but is significantly slower. In the relative reference system of the rotor, subtraction of the rotation results in a different incidence. The hot gas from the passage moves towards the pressure side, while the cooling fluid from the NGV tends towards the suction side of the rotor. Consequently, the pressure side experiences a higher thermal load than the suction side [12, 13]. The presence of a hot streak in the NGV passage is not required for this effect. However, larger temperature differences increase the effect. 1.2 CFD for cooled turbines Overall, an unsteady system emerges, and reducing it to a steady state simulation carries significant error potential [4]. However, the relevance of unsteady simulations for highpressure turbines is evaluated inconsistently in the literature. On the one hand, it is argued that unsteady effects have no significant impact on total pressures and mass flows in an uncooled turbine with homogeneous inlet conditions [14]. On the other hand, it is reported that unsteady effects significantly impact the heat transfer coefficients, albeit not the temperatures, in an uncooled turbine [15]. Furthermore, the necessity of considering these effects in a fully cooled turbine with hot streaks is emphasized [3]. Notably, as the complexity of the flow conditions increases, the importance of considering unsteady interactions increases, too. In a real turbine, all unsteady flow phenomena occur simultaneously and should therefore be considered in the design process. In the literature, the unsteady consideration of turbines always results in higher thermal loads. However, the exact quantification of this increase varies greatly depending on the problem, boundary conditions, and methods in question. Thus, the heat transfer coefficient increases by 15% [15], 30% [3], or in an extreme case, by 230% [5] compared to the steady-state consideration. Similarly, the temperatures and the film cooling effectiveness can locally change by 40 K [3] or 50% [5], respectively. This wide range demonstrates the need for precise case-by-case consideration and the challenge of finding generally valid values for the underestimation of thermal loads by steady-state methods. Although URANS has obvious advantages, its widespread use in the design process is prevented by the required computational effort. One way to reduce this effort is to adjust the number of blades. Slightly changing the number of blades scales the problem so that periodicity is achieved in a smaller segment across all blade rows. This simplification usually results in only minor changes to global quantities, such as efficiency [16]. However, the ratio of hot gas streaks to stators and rotors, as well as their positioning, has a significant impact on heat loads in the rotor [17]. Consequently, scaling can have a negative effect on the prediction of thermal loads. Another approach to reducing computational effort for unsteady simulations is to use frequency domain methods, such as Harmonic Balance (HB) [18]. The most important phenomena in a turbomachine occur periodically, usually with frequencies known from the rotational speed and the number of blades. This can be exploited when the governing equations are transformed into the frequency domain. Rather than seeking a solution for each time step, the aim is to identify a small set of harmonic oscillations that can describe the deterministic disturbances of the periodic flow at all times. Consequently, each solution contains all the information necessary to describe the temporal evolution of the system. Phase shifts between individual passages can be easily represented in the Fourier coefficients. Therefore, usually only one passage per row needs to be calculated. This reduction in the computational domain, coupled with the fact that the solution no longer needs to be advanced in time, significantly reduces computation time [19]. 2 2 METHODOLOGY All simulations are performed with DLR’s turbomachinery CFD solver TRACE [20]. TRACE uses the finite volume approach to solve the RANS equations. The density based, compressible solver attains second order accuracy using Roe’s upwind scheme [21] and a MUSCL extrapolation. A van Albada flux limiter is used to smooth large gradients, e.g. in the vicinity of shocks [22]. Time integration for the URANS simulations is performed with a second-order Euler-backward scheme. Two different combinations of turbulence and turbulent heat flux models are used. Menter’s SST k-ω[23] model combined with a constant turbulent Prandtl number of 0.9 serves as the baseline. This is compared with the seven-equation differential Reynolds stress model SSG/LRR-ωproposed by Eisfeld et al. [24] with the algebraic turbulent heat flux model of Daly and Harlow [25] (DH). This model combination was chosen because it performed well in two previous studies [26,27]. For both turbulence models a logarithmic formulation of the ω-equation is used to improve their robustness in HB [28]. 2.1 Harmonic Balance The aim of the HB method implemented in TRACE is to solve the RANS equations for time-periodic solutions in the frequency domain. This transforms the problem of finding a vector of conservative variables qfor each point in the flow field for each time step to finding a set of KFourier coefficients ˆqkthat describe the complete temporal evolution of the flow field. This can now be solved as a steady-state problem in the frequency domain using an efficient steady-state solver scheme. The HB method as implemented in TRACE is a frequency domain formulation, as proposed by McMullen et al. [29]. The RANS equation in the frequency domain ikw ˆqk+d R(q)k=0,(1) where ˆqkare complex-valued Fourier coefficients for the conservative flow variables and with the complex-valued residual function d R(q)k, is solved with an alternating frequency-time domain approach. Therefore, d R(q)kis determined by calculating the residuals in the time domain at 4K+1 sampling points, which are then transformed into the frequency domain via a discrete Fourier transformation. Once the residuals are calculated, the K+1 equations given in (1) are solved by a pseudo-timestepping scheme. More details about the Harmonic Balance solver implemented in TRACE can be found in [19]. Two configurations with different numbers of harmonics are used in this paper to assess harmonic convergence. The HB-1-3 variant serves as a baseline, minimising the required computational resources, whereas a higher-resolution setup with four harmonics in the stator and seven in the rotor (HB-4-7) is also used. FIGURE 2: Example grids with resolved and modelled film cooling holes. Cells marked in red receive a source term. 2.2 Film cooling model There are two main advantages to modelling the film cooling. Firstly, the cell count is significantly reduced compared to a grid that resolves the inner structure of the blades and each film cooling hole. This reduction in grid cells can be seen in Fig. 2 for the local model used here. Secondly, variations in the cooling design can be easily implemented without the need to create a new grid. This is particularly useful when several designs need to be evaluated on the same blade geometry, for example during the optimisation of the cooling design. This paper compares two different cooling models. The localised film cooling model presented previously in [26] and a simpler volume source model presented in [30]. The local film cooling model will be denoted as FCM and the volume source model as VSM. Both models inject mass, momentum and energy into the domain, but they differ in where exactly they place those sources. The VSM injects cooling fluid into the complete O-grid around the blades and cannot resolve individual cooling holes. In contrast, the FCM only acts in grid cells that lie directly above unresolved hole exits. Nevertheless, using a homogeneous source term for each cooling hole deliberately keeps this model simple and robust, enabling its application over a broad range of operating conditions. More elaborate, carefully calibrated models may become unreliable when extrapolation is necessary. The governing equations for the two models are identical. However, the location and number of affected cells and the total volume differ. In the FCM, a source term is distributed volumeweighted among all cells used to model a particular cooling hole. In the VSM, the same distribution is applied to all cells in the Ogrid around the blade. Let Vtotal denote the sum of the volumes of the cells in question and let Vrepresent the volume of a particular cell. The homogeneous source vector ˆ Sis defined as  S= [S˙m, SM,SE](2) S˙m=˙m Vtotal V(3)  SM=S˙mu(4) SE=S˙mhc+1 2|u|2(5) 3 These terms appear on the right-hand side of the integral form of the conservative Navier–Stokes system ∂ ∂tZV  QdV+ZS  FdS=ZV  SdV,(6) where  Qis the vector of conservative variables and  Fthe combined Euler and viscous flux vector. The user supplies, for each hole of the FCM, the mass-flow rate ˙m, the total temperature (from which hcis obtained), and the entry/exit positions and shape of the cooling hole. The necessary information about the thermodynamic state of the coolant ( ˆu,hc) is calculated in conjunction with the static pressure just upstream of the hole exit from the CFD simulation, assuming an ideal gas and isentropic expansion. For the VSM the mass-flow rate is summed and the total temperature mass-flow averaged. 2.3 Calculation of metal temperatures The metal temperatures are calculated using a semiempirical software tool [31]. The external recovery temperature, heat transfer coefficient, and static pressure distributions derived from the CFD simulation are specified as boundary conditions. The cooling system, here consisting of two impingement baffles, a pin fin array in the trailing edge area, and film cooling holes, is parameterised by specifying the corresponding geometric dimensions. The total temperature and total pressure of the coolant are specified at the inlet of the cooling system of the NGV. Based on the flow boundary conditions and parameterisation of the cooling system, the resulting coolant mass flow rate is calculated and the internal heat transfer coefficients are determined using empirical correlations. Thus, internal and external boundary conditions are defined and metal temperatures are calculated by solving Fourier’s law. The NGV utilises the single-crystal nickelbased superalloy CMSX-4 coated with yttrium oxide-stabilised zirconia ceramics. The temperature-dependent thermal conductivity of the materials and the heating of the fluid within the internal coolant network are considered. 3 TEST CASE DESCRIPTION The high pressure turbine stage is derived from an ultrahigh-bypass-ratio geared turbofan engine concept [32], for which the entire turbine section and in particular the 3D geometry of the first stage of the high pressure turbine was performance optimized in a dedicated turbine design process [33]. A realistic film-cooling configuration is applied to both the NGV and the rotor. This turbine is scaled to a 1:2 configuration in order to keep the URANS simulations computationally affordable, which is not necessary for HB. Table 1 provides a summary of the averaged thermal design parameters, which describe the maximum take-off engine operation point, at which the highest thermal stresses occur. FIGURE 3: Total temperature distribution at the inlet TABLE 1: Boundary conditions Tt[K] pt[MPa] Tu [%] lTu [m] ˙m[kg s−1] Hot gas 2168.9 4.5437 20 0.0150 51.94 Coolant 971.6 4.760 5 0.0035 8.47 To represent the effect of a realistic hot-streak emerging from the combustor, the two-dimensional temperature distribution reported by Povey and Qureshi [6, 34] (Fig. 3) is imposed on the mass-flow averaged inlet total temperature while the total pressure is kept uniform over the inlet plane. The inflow turbulence characteristics are derived from engine-relevant combustor experiments. Following the overview of Fischer et al. [35], the high-estimate turbulence intensity of Schroll et al. [36] and the mean length-scale reported by Willert et al. [37] are adopted. At the outlet, a static pressure boundary condition enforces a prescribed mass-flow rate through the domain. 3.1 Grids All grids are block structured and do not resolve the internal geometry of the NGV or the film cooling holes. In the following discussion, grid levels are named according to the number of cells in the spanwise direction. The base grid has 146 cells in this direction. Two coarser and one finer grids are derived from the base grid by uniformly refining or coarsening the cells in all directions. The wall distance in the first cell is kept constant to ensure that the y+value is less than 1 for most of the surface. The highest y+values occur in a small area at the trailing edges, where they reach approximately 1.5. 4 (a) FCM steady SST (b) FCM HB SST (c) FCM HB SSG/LRR-DH FIGURE 4: Grid convergence of the area averaged film cooling efficiency FIGURE 5: Difference in adiabatic wall temperature between HB-1-3 and HB-4-7 on the rotor surface using SSG/LRR-DH 4 RESULTS AND DISCUSSION In a first step, the grid dependency is evaluated to determine if the base grid with 146 is sufficient to resolve the relevant phenomena. A previous study on flat plates with single cooling holes observed a significant grid dependency of the film cooling model [26]. Such a strong grid dependency would hinder the use of the model in the design process. Figure 4 shows the spanwise and, where applicable, temporally averaged film cooling effectiveness ηaw for different methods across the grid resolutions considered. As with a previous study using just the NGV [27], grid dependency on the rotor is significantly reduced compared to the individual film cooling holes on a flat plate. On the suction side, the effect of grid refinement is negligible. However, in contrast to the NGV, the averaged temperature at the leading edge remains virtually unaffected by grid refinement. Instead, the largest changes occur on the pressure side, primarily due to a slight increase in temperature in between the individual cooling films as the grid is refined. Overall, the grid dependency is of the same magnitude as observed on the NGV at the leading edge and is not significantly affected by using an unsteady method for the simulation of the turbine stage. Secondly, HB simulations can only be considered reliable if all relevant frequencies are resolved. However, the HB-4-7 configuration is significantly more expensive than a setup limited to three harmonics. Because omitted frequencies do not necessarily influence the time-averaged solution, a strong reduction of the number of resolved harmonics is attractive. In the NGV only weak unsteady activity is observed in all simulations, so a reduced configuration with a single harmonic is sufficient. Figure 5 shows the differences in the time-averaged temperature field on the rotor surface between HB-4-7 and HB-1-3. For reference, the temperature difference between coolant and hot gas is around 1000 K. The deviation between the two simulations are generally small on the suction side (Xrel >0). Isolated cooling films exhibit a temperature difference of up to 30 K, but most of the surface differs by less than 10 K. On the pressure side local differences reach up to 50 K, which is of the same order as the difference between URANS and the HB simulation. Nevertheless, the mean discrepancy between the two HB runs remains far below the difference between URANS and HB. Therefore, the cheaper HB-1-3 setup was chosen. 4.1 Comparison of rotor surface data Figure 6 shows the time-averaged adiabatic wall temperature on the rotor for different methods combined with the FCM. The panels are arranged from the least to the most detailed simulation. The steady SST case at the top exhibits the characteristic long, cold cooling films. This is a frequently encountered problem, where the common two equation turbulence models tend to under-predict the mixing and diffusion of the coolant [38, 39]. Compared with the HB results with the same turbulence model, the steady method cannot resolve the heat redistribution caused by unsteady effects. Consequently, a noticeably hotter region appears on the pressure side ( 1 ). The HB simulation reveals a clear lateral motion of the cooling films in the time resolved data on both suction and pressure sides, whereas the films near the leading edge ( 2 ) oscillate mainly in the axial direction. Com5 pared to the steady simulation, this reduces the cooling effectiveness of the cylindrical holes in this area. Near the trailing edge ( 3 ), in contrast, this lateral motion slightly enhances cooling. Switching the turbulence model leads to a marked increase of the predicted adiabatic wall temperature over the whole blade surface. The film cooling effectiveness is reduced due to the expected increase in mixing and diffusion, which leads to shorter and hotter cooling films. Additionally, the temperature of the uncooled blade regions close to the end-walls rises because more heat is transported from the hotter core of the hot streak to the periphery. This is again due to an increased turbulent heat flux. The differences between the URANS and HB results in the last two figures are small. The only easily visible difference is located at the trailing edge of the pressure side ( 3 ). Here URANS predicts a slightly weaker unsteady motion of the cooling films, which lead to a reduced film cooling effectiveness. To break down which method is responsible for which change in the adiabatic wall temperature, Fig. 7 shows the differences between the various simulations. The upper panel shows the difference between SST and a steady RANS method on the one hand, and HB with SSG/LRR-DH on the other. Red indicates that the HB solution predicts a higher adiabatic wall temperature. Generally, HB with SSG/LRR-DH predicts higher temperatures across the whole blade. Only isolated cooling films, where the cooling film adheres better to the surface, produce a local temperature reduction. The differences in the upper panel (steady SST versus HB with SSG/LRR-DH) can be separated into a turbulence model effect and an unsteady flow effect. The middle panel compares HB-SST with HB-SSG/LRR-DH, thereby isolating the turbulence model influence when both simulations employ the same unsteady method. Most of the differences in the upper image stem from the change of turbulence model, especially on the suction side and in the trailing edge region of the pressure side. A notable exception is the light-blue area around 1 , where the leading-edge cooling is more effective with SSG/LRR-DH because the film adheres better. Consequently, there is a reduction in temperature close to the leading edge. The effect of an unsteady method is visible in the bottom panel by comparing the steady calculation with the HB result. The redistribution of stator cooling fluid from the pressure to the suction side becomes evident in region 1 , where the unsteady solution exhibits a temperature increase that significantly outweighs the reduction seen in the middle panel, leading to an overall temperature rise. As expected, a slight temperature decrease is observed on the suction side. However, this does not compensate the increase when using SSG/LRR-DH, so the suction side also shows a net temperature rise in the top panel. Another unsteady effect is visible in the trailing edge region around 2 . A lateral oscillation of the cooling films helps to distribute the coolant and thus lowers temperatures. Moreover, unsteady effects cause a small shift of the passage vortex on the 2 1 3 2 1 3 2 1 3 1 2 3 FIGURE 6: Comparison of time-averaged adiabatic wall temperature on the rotor surface for the FCM suction side at 3 , deflecting some cooling films towards the tip and producing alternating temperature increases and decreases. However, these local variations result only in a minor change of the spatially averaged wall temperature. In addition to the temperature, the thermal loads are defined by the heat transfer from the flow to the structure. This can be 6 1 3 1 2 FIGURE 7: Differences of adiabatic wall temperature on the rotor surface quantified using Nusselt numbers, which are shown in Fig. 8 for three different simulations. The first two figures show the results of steady-state simulations using the two different turbulence models. The differences between SST and SSG/LRR-DH are comparatively small. Additionally, these steady-state simulations produce Nusselt numbers of a similar order of magnitude to those on the stator. However, significantly higher Nusselt numbers result if an unsteady method is used, as shown in the lower image. Notably, areas with a significantly higher Nusselt number stand out, particularly around the tip and the trailing edge. In time resolved data, a significant temporal variation of the temperature can be observed in those locations. This is due to the aforementioned movements in the cooling films. Lateral movement on both sides of the profile near the trailing edge, and an axial movement due to the unsteady incidence variation at the FIGURE 8: Nusselt numbers on the rotor surface leading edge. Consequently, the Nusselt number can locally increase threefold as this movement promotes the transport of heat through the disturbed boundary layer. In total, there is an increase in heat transfer in the spatial average, resulting in higher thermal loads when unsteady effects are considered. 4.2 Comparison of adiabatic wall temperatures To better evaluate the methods, it is useful to further reduce the surface data that has already been presented. Therefore, the adiabatic wall temperatures and the material temperatures are averaged over a section representing 20–80% of the span. This area represents the central part of the cooled blade. The temperature distribution obtained in this way can be directly compared with those of the 1D and 2D design tools. Additionally, this averaging allows for the complex effects at the end walls to be neglected. While these effects are relevant when assessing thermal loads, 7 FIGURE 9: Adiabatic wall temperatures and differences in adiabatic wall temperature in the center section of the rotor focusing on the middle area of the blade significantly simplifies the discussion at this point without losing relevance. Figure 9 shows the averaged adiabatic wall temperatures on the rotor surface in the top panel. The URANS simulation with the film cooling model and SSG/LRR-DH turbulence model is the highest quality simulation and serves as a reference for the steady and HB simulations. The temperature distribution determined by the VSM is qualitatively significantly different. Compared to estimates from the film cooling design tools [31], the FSM produces a significantly more plausible temperature distribution. This indicates that the VSM cannot correctly reproduce the mixing and diffusion processes that are relevant for film cooled turbines, as also discussed in the study on the NGV alone [26]. Even when excluding the VSM, there still are significant differences of almost 400 K between the remaining simulations. Those differences can be seen more clearly when switching to a representation of the temperature differences. Therefore, Fig. 9 depicts the difference between the average temperatures in the bottom panel. All simulations are compared to the highestquality simulation available, FCM URANS SSG/LRR-DH. Positive values indicate that the respective method overestimates the temperature compared to the reference. To better assess the methods, two literature values are used here, which provide an approximation of the temperature difference at which a doubling or halving of the service life can be expected in the event of failure due to creep or low-cycle fatigue [1]. These are plotted as two grey bands of 15 and 30 K in the representation of the temperature differences. Additionally, a halving of the service life at a temperature increase of 20 K is specified for failure due to corrosion [1]. However, this band is not shown to improve clarity. In general, no single failure mechanism dominates, with all contributing to turbine failure [2]. When discussing adiabatic wall temperatures, these values are for reference only as they have no physical significance for gas temperatures. However, the effect of the simulation methods on the estimated service life is also estimated later. The effect of the turbulence model is most evident when comparing the two steady simulations (blue and red) in the bottom panel of Fig. 9. The large differences at the trailing edge of up to 380 K have already been discussed. In general the temperature curves differ qualitatively and quantitatively between the two turbulence models, with on average higher deviations toward the trailing edges. Nevertheless, there are regions in which both turbulence models predict similar temperatures, especially around the leading edge on the pressure side. This means that the difference between SST and SSG/LRR-DH cannot easily be corrected by adding a constant offset or applying a simple function. Therefore, it seems unlikely that the SST model can be reliably used to accurately examine the thermal loads in cooled rotors. The effect of the unsteady methods can also be seen in the bottom panel of Fig. 9 as the red curve. This curve illustrates the difference between the steady RANS and the URANS method, with SSG/LRR-DH used in both cases. As previously discussed, redistributing hot gas from the suction side to the pressure side results in a temperature increase on the pressure side and slight decrease on the suction side when using an unsteady method. On average, this effect is smaller than that caused by changing the turbulence model. Nevertheless, when using a steady method, the temperature on the pressure side exceeds the 15 K failure limit due to creep by a factor of ten. In contrast, the HB simulation in orange can reproduce the URANS simulation temperature very well. The maximum difference is around 20 K, but significantly smaller differences occur over large areas. Therefore, HB achieves significantly better results than the steady simulation. 4.3 Comparison of metal temperatures The differences between the CFD methods can best be assessed by comparing the adiabatic wall temperature, which is an direct output of the simulations. However, this does not enable the effect on the estimated service life to be assessed. To this end, the material temperatures are compared with the aforementioned temperature threshold values of 15 to 30 K. The halving or doubling of the assumed service life would only occur if this temperature change takes place at the critical point in the material. This point’s location differs for each failure mechanism and is generally unknown. Accordingly, the changes in service life discussed here are estimates representing the worst-case scenario. 8 FIGURE 10: Metal temperature differences in the rotor center line Additionally, the observed temperature differences in the material between the steady-state and unsteady methods are similar to the same differences in conjugated heat transfer simulations presented by Carta et al. [40]. However, as these literature values can’t be used for a direct comparison due to a lack of information on the geometry, the cooling, and the boundary conditions, they are only suitable for a plausibility check. The impact of the unsteady effects on the service life can be estimated with the help of Fig. 10. Generally, due to the insulating effect of the ceramic protective layer and the internal cooling, the temperature differences in the metal are expected to be smaller than the differences between the methods in the adiabatic wall temperature. Firstly, HB reproduces the URANS results for material temperatures as well. For most of the profile, the temperature differences are well below 15 K. Only at the trailing edge of the pressure side are there larger differences. This is not due to significantly different wall temperatures, but rather to the higher heat transfer coefficients in the HB model, which are caused by stronger unsteady motion of the cooling films. Secondly, unsteady effects alone (red curve) result in a material temperature difference of around 100 K on the pressure side near the leading edge. On the suction side, however, there is a moderate reduction due to the redistribution of hot gas from the suction to the pressure side. Using a steady method results in an overestimation of the service life of a factor of 8 near the trailing edge and a factor of 32 near the leading edge due to creep (limit 15 K) on the pressure side. In contrast, the thermal loads are overestimated on the suction side. Consequently, the estimated service life doubles over large areas. This creates the potential for an optimised design to reduce coolant massflow and thus losses. Neglecting the unsteady effects results in the estimated material temperature of the steady SST simulation on the suction side being unexpectedly close to the reference results. However, this is likely a particular instance in which several errors cancel each other out. The temperature is underestimated by using SST and overestimated by neglecting the unsteady effects. These errors also influence the heat transfer, which is leads to these three errors largely cancelling out. The magnitude of these individual errors depends on the cooling design and the flow conditions. Therefore, it should not be concluded that steady RANS simulations with the SST turbulence model provide reliable results on the suction side of cooled rotors. Instead, emphasis should be placed on the fact that material temperature, and hence estimates of the service life, are highly sensitive to observed disparities between methods. 4.4 Computational cost In order to evaluate the methods, it is important to understand the computing resources required by each method. The HB method and the film cooling model are expected to significantly reduce the computing effort, whereas using SSG/LRR-DH will increase it. Figure 11 plots the cost of some simulations of the cooled turbine stage on AMD Epyc 7702 processors on a logarithmic scale. As these were not collected on a benchmarking system, they may be subject to greater fluctuations. Nevertheless, they can be used here to illustrate the orders of magnitude. The top part shows simulations of the entire stage from Fig. 9. URANS simulations were only performed on the base mesh (146). For completeness, the lower part shows the cost of simulations considering only the NGV of this turbine stage [27]. As the volume source model does not offer any advantage over the local film cooling model in terms of cost, it is not shown here for clarity. For SSG/LRR-DH, there is an increase in computational cost compared to SST: 20% for the NGV and 33% for the entire stage. While a significant difference between these two values is unexpected, due to the aforementioned fluctuations, this is within the margin of error. The increase of roughly 25% is significant, but small compared to the difference between modelled and resolved film cooling. In a previous study on this NGV [27] the resolved simulation required 175 times more computing resources than the simulation using the model, which illustrates the significant advantage of the film cooling model. Additionally, grid generation for the model is much simpler as the entire internal geometry can be disregarded. As can be seen in the top part of Fig. 11, the URANS simulation of the turbine stage took 87 times longer than the steady simulation. In contrast, the HB simulation consumed only 7.2 times as many computing resources as the steady simulation. This clearly demonstrates a significant advantage. Furthermore, HB can generally be used to simulate any blade number ratio of a stage with the same computing effort. For a classic URANS simulation, however, a full wheel simulation is usually necessary for real turbine blade counts. For this turbine, that would mean an additional factor of 25 for URANS. While these are only rough estimates, they nevertheless underscore the cost of unsteady sim9