Inverse design of internally cooled turbine blades based on the heat adjoint equation

Abstract

A method of solution of the inverse problem in heat conduction is presented. The method, based on an adjoint optimization procedure, is applied to the design of the pattern of circular cooling passages inside coated turbine blades. The general case of a non-homogeneous solid material is considered. The numerical solution of both the temperature field and the adjoint problem is based on a finite element method. The mathematical method is explained and the procedure is validated against theoretical and experimental data available in open literature.

Keywords:

• inverse problem
• heat conduction
• adjoint methods

1. Introduction

A wide range of engineering problems in thermal analysis and design has been formulated as inverse heat transfer problems 1. The inverse problem involves the estimation of the cause by utilizing the knowledge of the effect. Several studies focus on the cause–effect relation between heat flux and temperature. A classical example of this class of inverse problems deals with the estimation of an unknown boundary heat flux, by using temperature measurements taken below the boundary surface. Other studies estimate the steady or unsteady inlet temperature profile, or cover aerothermal effects by estimating steady/transient wall heat flux in laminar flow, or deduce inlet velocity estimation based on the temperature measurements. For conciseness, the interested reader can refer to 1–6 and the references therein. In fact, the abovementioned problems are out of the scope of this present work, which belongs to the class of the design problems.

The design problem can be formulated as an inverse problem in which some conditions are given at the boundary, while the shape of the body contour that realizes the imposed thermal features is unknown. Applications of this approach to heat conduction design problems have been proposed and successfully applied to turbine blade cooling in the past three decades. The heat conduction inverse problem of the geometry design of the internal cooling passages of a turbine blade has been solved for circular 7, super-elliptic 8 and generic geometries of the cooling passages in multi-holed turbine blade 9,10. The numerical procedures are mainly based on a direct solver driven by an optimization method. In one of the earliest applications 7, a boundary element method for heat conduction analysis and a gradient method, the steepest descent method (SDM) have been used to converge to the inverse problem solution. After that, the exponential growth of computational resources has allowed for an extensive use of more flexible and CPU consuming numerical approaches. Solvers based on finite element method (FEM) or finite volume method in two and three dimensions have been used to evaluate the thermal field, while, in the modelization of the physics, the convection and radiation effects have been included. Then the research evolved to conjugate heat transfer analyses by including the mutual interactions with the fluid flow and by modelling film cooling effects 11.

From the optimization counterpart, the improvements of gradient-based methods have led to the various formulations by adjoint methods for 2D/3D problems 10,12,13, while the use of genetic algorithms has also been introduced for the solution of the single objective and multi-objective optimization problems.

A new impulse to thermal design in aerospace propulsion is related to recent efforts of designing aeroengines which meet the Vision 2020 requirements on gasturbine emissions and efficiency from Advisory Council for Aeronautics Research in Europe. A straight way to enhance efficiency is obtained by increasing the maximum cycle temperature, i.e. the temperature into the combustion chamber. Moreover, other investigations introduce the interstage turbine burner concept to modify the thermodynamic cycle during flight, going towards variables cycle aeroengines. In both approaches, a closer control of the temperatures is required and the thermal design of the some engine components, as the burners and the turbine blades, becomes more aggressive. In this scenario, automated inverse problem solvers can help the designer to make choices based on a wider investigation of the design variables space.

In this work an approach to inverse problem solution is proposed. The method is based on the adjoint optimization and follows the footsteps of a technique of aerodynamic design 14,15. The mathematical treatment of the adjoint problem differs from previous adjoint approaches (see, e.g. 10,13). In fact, the proposed formulation does not need a sensitivity problem to be solved and it includes a penalization for imposing geometrical constraints. To allow for the treatment of coated blades, the heat conduction equation in a non-homogeneous material has been considered. As first step, circular cooling passages are treated and a way to impose non-intersecting contour constraints is developed. An FEM approach has been used both to compute the thermal field inside the turbine blade and to solve for the adjoint problem. The plan of this article is as follows: in the next section the mathematical model and the inverse problem solution are presented; then the adjoint problem is derived and the numerical technique is explained. Finally the geometric parametrization of the blade cooling passages is derived and the procedure is validated against an analytical solution and against a test case based on experimental data available in the open literature.

2. Governing equations and inverse problem formulation

Let us consider the heat conduction equation in a non-homogeneous material and Robin boundary conditions. This set of equations may be written in a compact integral conservative form as (1) (2) (3) where T is the temperature, T_g and T_cj are gas and coolant flow temperatures, h is the heat transfer coefficient and n is outward normal vector. For the thermal conductivity, we pose k = k(x) as a function of position vector x. In fact, in turbine blades with ceramic coating one can distinguish at least two regions with very different k-values. A sketch of a generic domain is shown in Figure 1. Γ_b is the external surface of the turbine blade, Γ_a is the interfacial surface between inner core blade (light grey) and coating material (dark grey), and Γ_c is the union of the Γ_j surfaces of the N inner cooling passages. In the classical thermal inverse problem, the boundary shape of the cooling passages Γ_j is regarded as unknown while some field variables are known at the external boundary Γ_b.

Figure 1. Domain geometry and nomenclature.

In the literature, whatever the method of solution used, various choices have been adopted in formulating the inverse problem. In 9, both the heat flux and the temperature are imposed in a discrete set of points on Γ_b. In 7, the temperature on Γ_b is imposed and the geometry that realizes a target heat flux distribution in Γ_b is sought for. In 10, mixed boundary conditions are imposed on both Γ_b and Γ_j, while looking for the geometry of the cooling passages that realizes the desired temperature in a discrete set on points on Γ_b. In this work, the approach is similar to the latter case, which is in our opinion the most suitable for a straight application, as well as for the experimental validation. We impose Robin boundary conditions both on the external blade surface Γ_b and on the internal cooling passages Γ_j, here supposed of circular shape, and we use the adjoint-based gradient method to find the hole locations and diameters which realize a desired temperature distribution T_b(x) along the boundary Γ_b. In this work, the temperature distribution T_b(x) is a known function related to a designer choice.

The numerical solution of system (1) is based on the FEM. The derivation of steady conduction equation is a standard exercise in classical textbook on FEM and will not be repeated here. As a guideline and using the terminology and notation of 16, pp. 87–92], let us assume for the temperature field the form (4) where [N] = {N₁, N_i, … , N_r} is the shape function matrix and {T } = {T₁, T_i, … , T_r}^T the nodal temperature vector. Equations (1) can be reduced to (5) where (6) (7) with (8) The system (3) may be linear or not, depending on the expressions of [K] and { f }. The parallel version of the SPOOLES package has been used to solve system (3).

The inverse problem is solved as an optimization problem, as outlined in the following section. The procedure iterates on a series of direct computations until all boundary condition and constraints are satisfied within an expected range of gradient and cost function residuals.

3. Variational formulation, adjoint equation and gradient

We define the cost function (9) where the 𝒫(Γ_c) is a penalization function added to enforce the geometric constraints to the optimization problem. The control variable is the set of relations defining the cooling passages geometry Γ_c. In two dimensions, each boundary has a parametric representation as follows: (10) where p_j(s), q_j(s) are shape functions and μ_j, ν_j are the control parameters that can be packed in a single vector α_i.

The optimal temperature field must satisfy the governing equation (1) and some geometric constraints. In order to solve such constrained extremum problem we introduce the Lagrangian function (11) where Λ(x) is a Lagrange multiplier. The Lagrangian will allow us to treat the problem as unconstrained. A stationary configuration is found when the variation of ℒ with respect to all its arguments, that are now considered independent functions, is 0. We compute δℒ as (12) All the contributions to δℒ must be 0 at the minimum. Hence, to find a stationary point we enforce In general this results in . To reach the minimum we take δΓ_j such that .

Note that the variations of ℒ with respect to the Lagrange multipliers Λ simply yield the heat conduction equation. The condition δℒ_T = 0 leads to the adjoint equation and its boundary conditions. Based on the second Green's identity the second integral can be formulated as (13) By perturbing the boundary conditions (1b,1c) we can write (14) Substituting in (11) and considering all this integral contributions must vanish, we have (15) The adjoint problem (13) is formally identical to system (1). By changing Λ ← T and with slight modifications to the set of boundary conditions, the numerical method used to solve the thermal problem can be applied to the adjoint problem.

To derive the variations we observe that, from Equation (8) they are formally equivalent to (and therefore replaced by) the variations of the Lagrangian function against the control parameters α_i (16) where (17) The final expression of the terms 𝒢_i(Γ_i, T, Λ) depends on the adopted parametrization of the cooling passages geometry. The formulation of 𝒢_i for circular cooling holes is derived in Section 3.1. If we update α_i with (18) by taking ρ > 0, then δℒ_α ≤ 0. By iterating such a procedure, the minimum is eventually reached. This method, namely the SDM, has a slow convergence. Better convergence rates can be obtained with the conjugate gradient method (CGM), in which the correction δα at the iteration n is given as (19) where βⁿ is defined by the modified Polak–Ribière formula (20)

3.1. Geometric parametrization and constraints

The objective of this section is to derive a formulation of 𝒢_i able to parametrize the geometric optimization variables α_j and to include the necessary set of constraints.

We consider M circular cooling passages. The boundary of each passage is defined by (21) where C = (a_j, b_j), R_j are the centre and radius of the j-th circle, respectively.

For the sake of conciseness, the set of control variables a_j, b_j, R_j can be collected into a single vector α_i of 3M elements defined as (22) The hole boundaries cannot intersect each other and must lie far from the coating region of the blade. These requirements are constraints of the optimization problem and are introduced using penalization. In a way similar to 7, considering M coolant passages, we write the penalization function as (23) where (24) is the minimum distance between the i-hole surface and the Γ_a contour (see Figure 1), and (25) is the minimum distance between two coolant passage contours Γ_i and Γ_j. The complete set of derivatives is (26) and (27) Finally, we derive the expression for the variation of the Lagrangian function with respect to the parametrization. Letting a variation α_i ← α_i + δα_i, neglecting higher order terms, the functional 𝒢_i can be reduced to (28) where the derivative of the penalization function 𝒫 are given by Equations (24) and (25). Note that in Equation (20) the vector α_i is composed by three set of control variables a_j, b_j, R_j, therefore we are now led to a formulation of 𝒢_i in which the three different contributions have been distinguished.

4. Numerical results

The proposed numerical procedure is a classical optimization approach based on the adjoint equation to compute the functional gradient. Basically it consists of the following steps:

1.	solve the direct problem (1),
2.	solve the adjoint problem (13),
3.	evaluate the functional gradient (14),
4.	compute the conjugate direction of search via (17) and march towards the extremum,
5.	update the solution and test the convergence criterion

Two numerical applications are proposed in the following sections. The first example is a test against an analytical solution. The second test is based on the experimental data of the Mark-II turbine blade 17.

4.1. Coated cylinder with internal heating

Let us consider the heat transfer problem on a composite hollow cylinder with coating and a single circular heating passage as shown in Figure 2. The cylinder radii are r_b = 58, r_a = 38, r₁ = 25 mm. The inner flow has a temperature T_c = 603.15 K, while the temperature of the external gas flow is T_g = 303.15 K. The thermal conductivities of the inner core and of coating are k₁ = 15 and k₂ = 0.2 W m⁻¹ K⁻¹, respectively.

Figure 2. Test 1: geometry.

Robin boundary conditions have been applied by specifying the heat convection rates for the inner flow passage (h_c = 400 W m⁻²) and for the exterior cylinder surface (h_g = 60 W m⁻²). The exact solution of the temperature field on the whole domain is obtained by classical analytical methods and it can be found in 18, p. 63, Example 3–11].

The inverse problem proposed here is inspired to a similar test case found in 7. From the knowledge of the temperature T_b(Γ_b) on the external cylinder surface, the correct geometry of the system is sought. Starting from a generic geometric configuration for the cylinder, as shown in Figure 3, we solve for the location and radius of the heating passage that realizes the target external wall temperature T_w = T_b(Γ_b) (= 337.230 K).

Figure 3. Test 1: geometry and temperature field.

The problem parametrization follows the guidelines given in Section 3.1. The control variables are i.e. the centre coordinates and the radius of the inner circular hole, respectively. The thermal field and the related adjoint problem are then solved at each optimization step and the domain geometry is updated following Equation (17) by the CGM strategy. We used parabolic triangular finite element. Adaptive mesh refinement with error control is obtained by using the bamg meshing and adaptation tool 19. Since the computational cost of the procedure still remains affordable, a very crude meshing strategy has been adopted. By a grid sensitivity analysis, an adequate level of mesh refinement is found in the initial computational domain. At each optimization step then:

1.	The domain geometry is derived from the parametrization used and the grid is generated.
2.	The thermal field is computed by the FEM solver.
3.	This solution is used to generate a new adaptivemesh by bamg meshing tool.
4.	The thermal field and gradient are computed on the second mesh.

At each step both thermal field and grid are computed twice. Although the resulting numerical procedure is not very efficient, it ensures a better grid resolution and a more accurate evaluation of the gradient functional.

Initial and final temperature field and geometry are shown in Figure 3. The system has shown a fast and monotonic convergence. In Figure 4, the control variables α_i and the functional residuals versus the optimization steps are plotted. As visible, after 20 iterations the functional residual is decreased of three orders of magnitude, even if a penalization has been applied. It must be noted, anyway, that there were very few control variables and that the influence of the penalty function was limited to the mitigation of asymmetric variations of the gradient components.

Figure 4. Test 1: convergence history of x_c = a₁, Y_c = b₁, R_c = r₁ (left) and of the objective function residuals (right) vs. optimization step.

A comparison between numerical and analytical results is given in Table 1. The close agreement between data confirms the high resolution of the computations: the highest relative error is below 0.3%. We also note the remarkable accuracy of the heat flux computation, which involves a derivative evaluation and therefore was supposed to be more prone to numerical error. Moreover, the heat flux is a dependent variable which plays a key role on the system since it is imposed as boundary condition.

Table 1. Composite hollow cylinder test. Comparison between numerical and theoretical value of main parameters.

	a1 (mm)	b1 (mm)	r1 (mm)	Q (W m^−2)	Tb (K)
Theoretical	0	0	25	2044.79	337.230
Numerical	0.303	0.012	24.933	2044.54	337.234
Error	5.22·10^−3(*)	2.07· 10^−4(*)	0.27%	0.012%	0.0012%
Note: *The absolute error has been scaled to the external cylinder radius rb.

4.2. Mark II turbine blade

Experimental data about the Mark II and CF3X internally cooled turbine blades are available in the open literature 17 and have been used as a reference test case in several work, on the conjugate heat transfer and optimal blade design. The Mark II stator vane geometry is depicted in Figure 5(b). The cooling system of this blade is composed by 10 circular passages (M = 10). In the original work a method of evaluating the heat transfer coefficient is outlined and the experimental results (wall pressure, blade surface temperature and heat flux coefficient) are given for a wide range of working conditions. Walking on the footsteps of the previous section, the experimental data are used here to formulate a design problem and to solve it by the inverse numerical procedure. The measured surface temperature for a selected working condition (namely, the run-15) is assumed as target temperature T_b(x). Then, starting from an arbitrary initial geometry, the correct coolant passage pattern is sought by solving the adjoint inverse problem.

First, we define the expression for the heat coefficients. For the coolant passages, the expression of the Nusselt number given in 17 is used (29) where Pr and Re are the Prandtl and Reynolds numbers based on the coolant flow rate, viscosity and temperature. Cr is an empirical correction parameter introduced in 17.

For the external blade surface a polynomial approximation of the heat transfer coefficients has been computed by fitting and smoothing the experimental data. The curve obtained is shown in Figure 6. Again, the control variables are the centre coordinates and radii of the circular cooling passages. These variables has been packed as Starting from the generic initial condition depicted in Figure 5(c), the inverse procedure has been applied until a satisfactory level of convergence is reached. The converged solution and temperature field are shown in Figure 5(d), while the experimental and the evaluated blade surface temperature are compared in Figure 6. The convergence of the objective functional and of the L²-norm residual of controls α_i is monitored in Figure 7.

Figure 5. Test Mark II: initial (a) and final (b) grid and cooling passage pattern, initial (c) and final (d) temperature field.

Figure 6. Test Mark II: heat transfer coefficient fitting (left) and comparison of the surface temperature obtained by the inverse problem and the target temperature distribution (right). Note: Symbols refer to the experimental data and solid lines to numerical results).

Figure 7. Test Mark II: convergence history of L²-norm of gradient residuals (left) and of the objective function residuals (right) vs. optimization step.

Although this test being quite severe (there are 30 control variables and strong penalty functions), the geometrical configuration obtained is close to the real one. For each α_i the relative error is less than 0.5%, a value slightly bigger than the error observed in the previous test T1.

It has also been observed that the penalization effect is twofold: it reduces the convergence rate dramatically and makes the solution less sensitive to the original objective functional. In the present case, it acts also as a smoothing filter for the target temperature peaks and experimental noise. Nevertheless, by using a parallel implementation of the FEM code, the full computation has run in about 1 h on a 8-core Intel i7 workstation.

5. Conclusions

An adjoint procedure for the solution of the inverse heat conduction problem has been proposed. Compared to similar approaches, the present method does not require the evaluation of an additional perturbation field and the functional gradient with respect to the control variables is computed from its parametrization. This simplification is supposed to increase the robustness of the whole numerical procedure.

The procedure has been validated against analytical and experimental solutions. The numerical tests carried out have shown an acceptable convergence rate, even if as much as 30 control variables were treated and a strong penalization is applied to avoid intersections between adjacent coolant flow passages or between a coolant passage and the external surface of the blade. The penalization is supposed to be responsible for the convergence worsening. In fact, it is well known that the number of control variables does not affect the performances of adjoint procedures. Nevertheless the computational cost is affordable on a medium level workstation.

The method has been specialized to the case of blade with circular passages holes. The extension to other geometries and to the three-dimensional case is straightforward but at the cost of 3D FEM computations and of a more complex parametric representation. From the numerical point of view a general purpose FEM engine can be used to solve both the thermal and the adjoint problem. Therefore, the method can be implemented in any free or commercial FEM engine that allows to impose user defined function as boundary conditions. Work on improvements on the grid adaptation strategy and on the parametrization of actual geometries of multi-holed turbine blades is in progress. In future work, the transient heat conduction/convection should be considered in the optimization procedure as a route to robust design and turbine blades life enhancement.