## Abstract

Blade cooling technology will play a critical role in the next generation of propulsion and power generation gas turbines. Accurate prediction of blade metal temperature can avoid the use of excessive compressed bypass air and allow higher turbine inlet temperature, increasing fuel efficiency and decreasing emissions. Large eddy simulation (LES) has been established to predict heat transfer coefficients with good accuracy under various non-canonical flows, but is still limited to relatively simple geometries and low Reynolds numbers. It is envisioned that the projected increase in computational power combined with a drop in price-to-performance ratio will make system-level simulations using LES in complex blade geometries at engine conditions accessible to the design process in the coming one to two decades. In making this possible, two key challenges are addressed in this paper: working with complex intricate blade geometries and simulating high-Reynolds-number (*Re*) flows. It is proposed to use the immersed boundary method (IBM) combined with LES wall functions. A ribbed duct at *Re*=20 000 is simulated using the IBM, and a two-pass ribbed duct is simulated at *Re*=100 000 with and without rotation (rotation number *Ro*=0.2) using LES with wall functions. The results validate that the IBM is a viable alternative to body-conforming grids and that LES with wall functions reproduces experimental results at a much lower computational cost.

## 1. Background

Gas turbine technology plays a critical role in civilian as well as military propulsion and in power generation. One of the key challenges in the development of advanced gas turbines is effective thermal management in the high-pressure stage of the turbine. High combustion temperatures (approx. 1500–1700°C) mandate the use of active cooling technologies to cool the high-pressure stage nozzle vanes and blades to keep the metal temperature below the melting point [1]. In response to this need, the modern high-pressure gas turbine blade is designed not only for aerodynamics and structural strength but also for effective cooling. This is accomplished by circulating compressor bypass air at a relatively cooler temperature (about 600–700°C) through internal passages in the blade or through the skin (internal cooling) and ejecting into the hot gas mainstream to protect the blade surface from the hot mainstream gases (film cooling). Figure 1 shows a typical high-pressure aero-engine blade geometry consisting of internal passages fed by cooling air at the base through the rotor hub. It consists of four passages, one cooling the leading edge, another cooling the trailing edge, and with a two-pass passage at the centre of the blade. In practice, a variety of heat transfer augmentation techniques are used in the internal passages in the form of flow turbulators (ribs, pins, dimples and protrusions), and impingement cooling to promote turbulence and enhance heat transfer coefficients. In the example blade shown, the internal passages are equipped with angled ribs and coolant is ejected out on to the blade surface on the pressure side near stagnation, at the trailing edge and at the blade tip.

Accurate numerical predictions of blade temperature can have a direct impact on the life of the blade and fuel consumption. As a case in point, accurate prediction techniques can reduce the factor of safety required in designing for the quantity of compressor bypass flow, which is in the range between 10 and 15%. Reducing the bypass flow increases the amount of work extracted from the compressed air and increases stage efficiency. By the same token, accurate prediction techniques can allow higher turbine inlet temperatures for the same amount of bypass air flow and increase the thermal efficiency of the cycle. The realization of a hydrogen economy and using hydrogen to fuel land-based gas turbines with zero CO_{2} emissions depends on the ability of the high-pressure stage to withstand potential temperatures up to 1800–1900°C. Thus, the ability to predict the first-stage guide vane and blade temperatures accurately has a direct impact on the efficiency, fuel consumption and emissions. As an example, a single 10 h transatlantic flight of a Boeing 747 consumes about 109 000 kg of jet fuel and releases about 325 000 kg of CO_{2} into the environment. Even a 5% reduction in fuel consumption can reduce the amount of CO_{2} emitted by 16 250 kg per 10 h of flying time for a single large aircraft.

Because of its paramount importance, the heat transfer and pressure drop characteristics of different surface augmentation geometries continue to be of great interest to the gas turbine industry. The flow and heat transfer characteristics depend on many different geometric parameters that characterize the shape of the cooling passage, the shape and configuration of roughness elements, and flow conditions. Computational modelling has played a very important role in identifying the underlying physics of cooling flows. Experiments have mostly measured the heat transfer coefficients without much regard to the underlying hydrodynamic mechanisms for a number of reasons, chief among them being the singular importance of heat transfer coefficients to the design process and the difficulty of making velocity measurements particularly in a rotating geometry, which is a challenge even to this day [2]. A vast number of numerical investigations have been performed to investigate the ability of numerical models, in particular turbulence models, to predict the heat transfer. Early research has mostly been conducted with the use of the Reynolds-averaged Navier–Stokes (RANS) method with a variety of turbulence models ranging from two-equation eddy-viscosity models such as the *k*–*ε*, *k*–*ω* SST, *v*^{2}–*f* model and algebraic stress models, to full Reynolds stress closure with varying degrees of success [3–12]. In the last two decades, with increasing computational power and capacity, these efforts have transitioned to more advanced time-dependent methods such as large eddy simulations (LES) [13–24] and hybrid methods using unsteady RANS (URANS)–LES or detached eddy simulations (DES) [25–28]. It is clear from these LES investigations that the predictions capture the underlying flow physics with good accuracy, and predictions of heat transfer coefficients are consistent and repeatable over a range of geometries and conditions. Similarly, DES investigations have also given encouraging results in predicting heat transfer in ribbed ducts.

Because of their time-dependent nature, these methods need large integration times and substantially finer grid resolutions than RANS. In spite of these challenges, they provide a much better ‘predictive’ capability than RANS and are frequently used in academe while also gaining traction in the gas turbine industrial community. It is envisioned that time-dependent methods will be the norm rather than the exception as the availability of high-performance computing (HPC) resources becomes more pervasive, and as the price-to-performance ratio drops further. In the late 1980s–early 1990s, a Cray-YMP costing about $10 million had a peak performance of about 4×330 MFlops (10^{6}), whereas today a peak performance of 1 Teraflop (10^{12}) can be purchased for about $2000 (NVIDIA-GPU card). If hardware costs follow the same exponential trajectory, which they are expected to, projecting outwards to 2030, $10 will buy approximately 1 Petaflop (10^{15}) of performance and a simulation is expected to run about one to two orders of magnitude faster. Under these conditions, it is envisioned that unsteady computations will transcend simple canonical internal cooling geometries to system-level analysis of heat transfer, just as they have transcended from simple channel geometries to more complex geometries in the last two decades.

Historically, computational hardware has advanced at a much greater pace than software. Thus, the real challenge will lie in our ability not only to exploit novel computer architectures to their full potential but also to look towards new methods and techniques for making system-level computations feasible. Considering the blade shown in figure 1 as our representative system with the goal of predicting the metal temperatures under representative engine conditions of external as well as internal cooling flow, a number of simulation and modelling challenges can be identified:

(1)

*Complex geometry*. Referring to the internal cooling geometry of the blade in figure 1, resolving the intricate geometrical features with a good-quality grid for stable time integration and suitable resolution of turbulence would be extremely time consuming using structured multiblock grids. Carefully designed unstructured grids with a judicious choice of elements can be used but they suffer from reduced performance compared with structured grids (between factors of 5 and 10) on modern hierarchical cache-based architectures. Additionally, even unstructured grid generation would be quite time consuming to incorporate all the relevant geometrical features. Thus, the current state of the art requires significant simplification of the geometry [29].(2)

*Large range of geometrical and turbulent length scales*. As a case in point, the overall dimensions of the blade in figure 1 is 61.5 mm in height and 20 mm wide at the hub, whereas the film cooling holes are 0.35 mm in size giving a length-scale range of two orders of magnitude. To resolve turbulence in the film cooling holes accurately, about 50–100 grid points would be needed in the radial direction of the hole. Extrapolating this requirement to the full blade geometry, accounting for the Reynolds number and internal cooling features, the total grid resolution is estimated to be of the order of a billion grid points to accurately resolve the internal flows in the blade. The computational complexity imposed by this requirement is not impossible to surmount with thousands of cores at our disposal. However, the orchestration of such a simulation from start to finish, including grid generation and long integration times, would be far from trivial. Such a calculation would be one of a kind in the present computing environment. The challenge then is to be able to perform multiple simulations of this type in a design environment.(3)

*Complex flow physics*. The external flow is highly turbulent and compressible, with exit Reynolds numbers in the range 10^{6}–10^{7}. Resolving the turbulent boundary layers and predicting transition on the suction surface in the presence of free-stream turbulence and upstream wake effects is challenging. The presence of film cooling imposes additional challenges. The internal flow is mostly incompressible but could become compressible at high coolant flow rates. It is also characterized by complex turbulent flows with separation and reattachment, secondary strains giving rise to strong anisotropies, secondary flows which impact heat transfer, streamline curvature and a host of other non-canonical effects.(4)

*Large disparity in the characteristic fluid and solid transport time scales*. To obtain the metal temperature of the blade, it is necessary to perform a conjugate heat transfer analysis in the blade, which requires the coupling of the energy equation at the solid–fluid interface. While the flow time scale is of the order of milliseconds, the conduction time scale in the blade material is of the order of seconds to tens of seconds [30]. As the simulation time scale is limited by the fluid time scale, long integration times would typically be required for the solid temperature to reach a quasi-steady state, increasing the cost of already compute-intensive calculations.

In this paper, we address the first three challenges—eliminating the need for meshing the complex geometry and reducing the fluid grid requirements for resolving the time-dependent turbulent flow accurately. The method we propose for resolving the complex geometry of the blade is based on the immersed boundary method (IBM), whereby the full blade geometry is immersed into a volume background grid. In this framework, since the surface of the blade is not resolved exactly as in a body-conforming grid (BCG), very fine grids are required in the vicinity of the immersed boundary. The fineness of the grid is directly related to the Reynolds number of the flow. Because of this difficulty, IBM has mostly been applied to low-Reynolds-number flows. Time-dependent turbulent flows have also been calculated using LES but at low to moderate Reynolds numbers and in simple canonical geometries [31]. Because of the inability of the IBM to exactly resolve the immersed surface coupled with the fact that in wall-bounded flows resolving the small energy-containing eddies in the wall layer is extremely expensive versus the outer layer of the boundary layer, combining IBM with a wall model for LES (WMLES) has the potential of rendering the calculations more accurate and more tractable at high Reynolds numbers [32].

Towards realizing the goal of combining IBM with WMLES, preliminary work in this direction for ribbed internal cooling ducts is shown. First, results of IBM–LES for a ribbed duct at *Re*=20 000 are shown and compared to a BCG calculation. Secondly, the use of the WMLES is highlighted in a two-pass duct at *Re*=100 000 and the predictions are validated with experiments.

## 2. Large eddy simulation of fully developed flow in a ribbed duct using immersed boundary method

The IBM is implemented on a general curvilinear mesh on a non-staggered grid [33]. The major steps in the implementation are as follows: (i) identification of the location of the immersed surface and determination of solid and fluid nodes, (ii) solution of the governing equations at all nodes in the domain that are not in the immediate vicinity of the boundary and (iii) applying a special treatment to the nodes that are directly next to the surface (IB nodes) to reflect the presence of the immersed surface. For this purpose, each IB node is associated with a probe which is used for setting the boundary condition at the IB node. One of the key challenges in this method is the consistent identification of fluid, solid and IB nodes and the associated probe locations in complex geometrical configurations. As an example, the identification of the IB nodes for the turbine blade is shown in figure 1*c*. The identification uses a background grid of 32.8 million cells in a rectangular box around the blade and 6.2 million surface elements describing the blade features.

To test the capability of IBM to simulate the turbulent flow in a ribbed duct, a square duct with square ribs with rib height to hydraulic diameter ratio *e*/*D*_{h}=0.1 and pitch to rib height ratio *P*/*e*=10 is simulated under the assumption of fully developed flow at *Re*=20 000 (*Re*=*UD*_{h}/*ν*; *U* is the mean velocity; *ν* the kinematic viscosity). The dynamic Smagorinsky subgrid stress model is used with a second-order central difference discretization [19]. The IB geometry is defined by 2.8 million surface elements and is placed in a non-uniform background grid of 128×136×136 in the *x*-, *y*- and *z*-directions, respectively, as shown in figure 2. Approximately 128×128×128 cells are available for the fluid domain calculation, which is similar to that used in the BCG. The triangulated surface is immersed in the background grid to identify the fluid IB nodes adjacent to the surface. The calculation was run for about nine flow-through time units with mean (RMS) and turbulent statistics calculated for the last 3.5 units.

Figure 3 compares the mean flow field and the turbulent kinetic energy (TKE) in the centre of the duct (*z*=0.5). The shear layer formed on the rib separates at the leading edge of the rib as part of it attaches back on the rib forming a small recirculating zone on top of the rib with a large primary recirculation region behind the rib. This region is characterized by unsteady vortex shedding from the separated shear layer and is the primary source of turbulence augmentation in the duct. A weak sympathetic secondary recirculation zone is also formed between the primary recirculation zone and the rib. The recirculating zone at the front of the rib is a result of instantaneous unsteady junction eddies that form in this region. High levels of turbulent fluctuations are present in the separated shear layer (TKE between 14 and 18%) as a result of unsteady vortex shedding. These high levels persist in the RMS components of velocity at *x*=1.0 and *z*=0.5. In addition to these features, localized secondary flows are formed which result in lateral flow impingement near the ribs [20]. The IBM captures all the major features of the flow with some slight quantitative differences in the predicted turbulent quantities.

## 3. Wall model for large eddy simulation in a two-pass duct

The geometry, based on the experiments of Iacovides *et al*. [34], consists of a square cross-sectioned two-pass internal cooling duct with 180° U-bend as shown in figure 4. Square-sectioned ribs are employed in a staggered arrangement along the inner and outer walls of the duct. The rib pitch to height ratio is *P*/*e*=10 and height to hydraulic diameter ratio is *e*/*D*_{h}=0.1. Two cases are examined at *Re*=10^{5}: stationary case and positive rotation at a rotation number *Ro*=0.2, where the rotation number *Ro*=*ω*_{z}*D*_{h}/*U* (*ω*_{z} is the angular velocity about the *z*-axis).

The calculations are performed on a boundary conforming grid with each rib pitch meshed with a grid of 378 000 cells using 14 blocks around the rib giving a total grid size of 7.78 million computational cells in the domain. The first layer of cells is spaced to give a of approximately 20 at walls. The dynamic Smagorinsky model is used together with a two-layer wall model [35]. In the wall model, a simplified one-dimensional boundary layer equation is solved in the layer between the first grid point and the wall (inner layer) on a virtual embedded mesh by using the instantaneous velocity boundary condition from the outer LES calculation. The coupling between the inner and outer layers is completed by using the wall shear stress obtained from the inner layer calculation to set the boundary condition for the outer layer LES. Details of the procedure are given in Patil and Tafti [35]. While the wall model plays a dominant role in conventional shear-dominated boundary layer flows as experienced in the U-bend, it plays a relatively modest role in the separated shear-layer-driven flow dynamics in the ribbed sections. The calculations are run for approximately seven flow-through time units, with averaging performed over the last five units to obtain the mean and turbulent statistics, at which point the turbulent statistics were found to reach a stationary state.

### (a) Stationary duct

Figure 5*a* shows the mean velocity vectors at the centre of the duct. On entering the duct, the flow reaches a fully developed state by the third or fourth rib. As the bend is approached, the velocity profile is distorted and the flow accelerates around the inner corner and separates as it traverses the tight inner radius. Downstream of the bend, the inner recirculating flow region forces the bulk of the flow to the outside and takes substantially longer to recover back to a fully developed state. Iso-surfaces of coherent vorticity representative of the turbulent intensity in the duct show that the flow coming out of the bend is highly turbulent compared with the flow upstream of the bend. This leads to higher heat transfer augmentation in the second leg.

Figure 5*c*–*e* shows comparisons of mean and turbulent quantities (*U*_{RMS} and ) at select locations (figure 4) in the flow. The distortion of the streamwise velocity is captured very well by the simulations as the flow enters and exits the bend. There is a slight acceleration on the inner side as the flow enters the bend, while the flow coming out of the bend is concentrated on the outer side with negative velocities on the inside characterizing the separated zone. The distribution of *U*_{RMS} entering and exiting the bend is also captured with good accuracy except some under-prediction of the peak values measured on the rib going into the bend and on the rib coming out of the bend. The turbulent shear stress at the first and second rib locations coming out of the bend also agree well with the experiments. The low values of shear stress on the outside at location B are consistent with the small gradients in the streamwise velocity.

### (b) Rotating duct

Coriolis forces resulting from blade rotation profoundly influence the heat transfer coefficient. The rotation number *Ro* is typically less than 0.3 for aero-engines but can be greater than unity in large power generation turbines. In a rotating blade, Coriolis forces act to stabilize and attenuate turbulence along the leading side of a duct (aerofoil suction surface) and to destabilize and augment turbulence on the trailing side (pressure surface). Rotation also induces cross-stream pressure gradients, resulting in the development of secondary flow cells in the cross section of the duct, which develop strong upwash and downwash regions. The direct effect of turbulent structure causes the heat transfer on the trailing wall to increase significantly, while heat transfer on the leading wall is decreased because of the attenuation of turbulence. Secondary flows resulting from the strong movement of fluid from the trailing to the leading wall increase heat transfer along the smooth side wall.

The effect of rotation on the mean and instantaneous turbulence in the duct is shown in figure 6*a*,*b*. The orientation of the rotational axis along the positive *z*-direction destabilizes the flow at the inner wall, which acts as the leading wall in the first pass and the trailing wall in the second pass. Comparing the mean velocity vectors in the bend region to that without rotation, the inner bend separated shear layer is pushed towards the outer wall creating a much larger recirculation region in the downstream leg. The instantaneous coherent vorticity clearly shows the augmentation and attenuation of turbulence at the inner and outer wall, respectively, in the first leg. Rotation also results in higher levels of turbulence than without rotation immediately downstream of the bend. Figure 6*c*–*e* compares streamwise velocity, *U*_{RMS}, and profiles with experiments. The mean velocity entering the bend shows some differences with the measurements, whereas very good agreement is obtained coming out of the bend. The peak value of *U*_{RMS} is slightly under-predicted at location C coming out of the bend but is very well predicted further downstream at location D. The predictions of turbulent shear stress just downstream of the exit to the bend C and location D further downstream also agree very well with the experiments. At the location D, the augmentation and attenuation of turbulence at the trailing (inner) and leading (outer) wall, respectively, is clearly discernible in the distribution of turbulent shear stress.

## 4. Concluding remarks

Effective active cooling of turbine nozzle vanes and blades in propulsion and power generation gas turbines has a direct impact on energy efficiency and pollutant emissions. A key component in their design is the accurate prediction of metal temperatures under engine conditions, which in turn requires that the highly turbulent external and internal flow fields are predicted accurately. This raises a number of challenges in modelling and simulation, chief among them being the resolution of complex geometries over a multitude of length scales, anisotropic turbulent flows with transition, complex secondary flows induced by non-canonical flow phenomena, and conjugate heat transfer and the associated disparity in time scales of heat flow between solid and gas. It is envisioned that as the price-to-performance ratio of computational hardware drops further and high-performance computers become more pervasive, time-dependent simulations based on large eddy simulations and hybrid URANS–LES, which are superior in their predictive capability than steady-state methods, will provide the avenue for full blade simulations under engine conditions. Towards this goal, it is proposed to develop the IBM to resolve the intricate internal blade geometry together with appropriate turbulence wall modelling capability for LES (WMLES). Preliminary results using IBM–LES at *Re*=20 000 for a ribbed duct, and WMLES at *Re*=100 000 show the feasibility of the proposed approach. Future work will focus on predicting the heat transfer coefficients and combining the IBM with WMLES, predicting conjugate heat transfer and applying the method to more complex blade geometries.

## Funding statement

The authors would like to acknowledge the support of Advanced Research Computing (ARC) at Virginia Tech.

## Footnotes

One contribution of 13 to a Theme Issue ‘Aerodynamics, computers and the environment’.

- © 2014 The Author(s) Published by the Royal Society. All rights reserved.