Continuous Adjoint-Based Optimization of an Internally Cooled Turbine Blade—Mathematical Development and Application

: This paper presents an adjoint-based shape optimization framework and its demonstration in a conjugate heat transfer problem in a turbine blading. The gradient of the objective function is computed based on the continuous adjoint method, which also includes the adjoint to the turbulence model. Differences in the gradient resulting from making the frozen turbulence assumption are discussed. The developed software was used to optimize both the blade shape of the internally cooled linear C3X turbine blade and the position of cooling channels aiming at (a) minimum total pressure drop of the hot gas ﬂow and (b) minimum highest temperature within the blade. A two-step optimization procedure was used. A free-form parameterization tool, based on volumetric NURBS, controls the blade airfoil contour, while the cooling channels are free to move following changes in the coordinates of their centers. Geometric and ﬂow constraints are included in the performed optimizations, keeping the cooling channels away from the airfoil sides and retaining the turbine inlet capacity and ﬂow turning. optimization started with the best solution of the ﬁrst run, aiming at reﬁning the cooling channel locations for minimum highest temperature over the blade. The optimized conﬁguration has ∼ 20% less total pressure drop and ∼ 30 K lower highest blade temperature. The geometric constraints used in the second optimization run constrained the cooling channels from moving too close to the blade sides, increasing the highest temperature gradient magnitude inside the blade less than 0.3% of the baseline. The results show that the developed continuous adjoint CHT method, coupled with a gradient-based algorithm, provides a useful tool for the CHT optimization of internally cooled turbine blades. They also reveal possible losses in accuracy, regarding the gradient computation, from making the frozen turbulence assumption in the development of the adjoint method.


Introduction
The efficient design of gas turbines with improved performance and longer lifetime, for use in aerospace and marine propulsion, land power plants, and other industrial applications, is of great importance. Gains are expected from aerodynamically optimized compressors, reduced friction in mechanical parts, optimal performance of inlet ducts and nozzles, and/or higher turbine inlet temperatures. The latter is limited by the thermal resistance of the first turbine row blades which are currently able to operate at gas temperatures exceeding the limits set by the materials thanks to internal and/or external cooling. Relevant studies account for the interaction of blades and gas flow in gas turbine components. This can be done using conjugate heat transfer (CHT) analysis and optimization.
A CHT analysis involves the strongly or loosely coupled solutions of the governing equations in the adjacent fluid and solid domains. In a strongly coupled scheme, PDEs in different domains are solved simultaneously, while in a loosely coupled one, each discipline is solved separately and communicates with the other(s) by exchanging data along their interfaces [1]. For the shape optimization of cooled blades, both stochastic and gradientbased methods have been used. The latter make use of the adjoint method to compute the gradient of quantities of interest at a cost independent of the number of design parameters. In [2], the continuous adjoint method to compute the objective function gradient in CHT

Governing Equations
In a 3D CHT analysis problem, the flow equations in the fluid domain Ω F are coupled with the heat conduction one in the solid domain Ω S . In Ω F , the RANS (mean-flow; MF) equations for compressible fluid flows read where f inv k = [ρv k ρv k v 1 + pδ 1k ρv k v 2 + pδ 2k ρv k v 3 + pδ 3k ρv k h t ] T are the inviscid fluxes and f vis k = 0 τ 1k τ 2k τ 3k v τ k + q F k T the viscous ones. ρ, p and h t stand for the fluid's density, pressure, and total enthalpy, respectively. v k are the velocity components and δ km the Kronecker symbol. τ km = (µ + µ t ) ∂v k ∂x m + ∂v m ∂x k − 2 3 δ km ∂v ∂x − 2 3 δ km ρk is the stress tensor based on the Boussinesq assumption. The heat fluxes are q F k = κ F ∂T ∂x k , with κ F = C p µ Pr + µ t Pr t being the fluid's conductivity, where µ and µ t are the bulk and turbulent viscosity. C p is the specific heat under constant pressure, Pr = 0.72 and Pr t = 0.90. The dynamic viscosity depends on the fluid's temperature T F [11]. Turbulence is modeled either by the one-equation Spalart-Allmaras or the two-equation k−ω SST model. The transport equations of the Spalart-Allmaras variable (ν), the turbulent kinetic energy (k), and the turbulence frequency (ω) are 2 are the production and destruction ofν, and d stands for the distance from the nearest solid wall boundary. µ t can be computed via The constants κ, C b1 , C b2 , and C w1 along with the expressions of f v1 , f w , f t2 , and S are given in [12]. The last two quantities are functions of vorticity ζ = √ 2W km W km , . Constants β , α 1 , σ ω2 , the expression for the production of k P k = min(P k , 10β ρωk) , and those of σ k , σ ω , β, γ, F 1 , and F 2 can be found in [13].
In Equation (5), S stands for the strain rate magnitude. The distance field (d), which affects the production and destruction terms ofν and the solid wall condition for ω, is computed by solving the Eikonal equation, The last state equation, which stands for the heat conduction equation over Ω S , reads where q S k = κ S ∂T S ∂x k are heat fluxes and κ S is the thermal conductivity of the solid. The flow and heat conduction solvers communicate across the fluid-solid interface (FSI), where T F = T S and q F k n F k = −q S k n S k must be satisfied. n F k , n S k are the normal vectors to the FSI pointing toward Ω F and Ω S , respectively.

Continuous Adjoint Formulation
In this work, two objective functions and two flow constraints are considered. The total pressure drop (F 1 ) and the highest temperature of the blade (F 2 ) are the two objective functions. The inlet capacity (F 3 ) and hot gas exit flow angle (F 4 ) are the two flow constraints. These quantities are where the highest blade temperature is approximated by the above differentiable expression where α takes on a large value.ṁ I ,p I t ,T I t are the massflow and the mass-averaged total pressure and total temperature at the hot gas inlet, respectively, andp O t ,v O 1 ,v O 2 stand for the mass-averaged total pressure and velocity components at the hot gas exit. Gradient-based optimization algorithms require the gradient of both the objective and constraint functions; thus, in the remainder of this section, the continuous adjoint method for a generic function F (standing for F 1 , F 2 , F 3 , or F 4 ) is formulated. In continuous adjoint, F is augmented (F aug ) by the field integrals of the product of the state equations' residuals with the adjoint variable fields over Ω F and Ω S . For the development of the adjoint method, only the Spalart-Allmaras turbulence model is used. Differentiating F aug with respect to the design variables b i gives where Ψ n ,ν a , and d a are the nth adjoint flow fields (n = 1,. . . ,5), the adjoint Spalart-Allmaras, and adjoint distance fields, respectively, all defined in Ω F . T a stands for the adjoint temperature and is defined in Ω S . By definition, correlates its total (δ/δb i ) and partial (∂/∂b i ) derivatives as well as grid sensitivities δx k δb i . Total derivatives of Φ with respect to b i and spatial derivatives permute according to [14] After a lengthy mathematical development, can be expressed as the sum of field and surface integrals that include variations in the state variables, and integrals depending on sensitivities of geometrical quantities with respect to b i . Setting the multipliers of variations in the state variables to zero gives rise to the field adjoint equations (FAE) along with the adjoint boundary conditions (ABC). In the next subsections, the four field integrals in Equation (9) are further developed in order to derive the FAE, the ABC, and the expressions of the sensitivity derivatives (SDs).

Development of the I MF Integral
Substituting the residuals of the MF equations in Equation (9), the first field integral of this equation takes the form Using the divergence theorem, its inviscid counterpart becomes Note that, in integrals over Ω F and Ω S or their boundaries (S F , S S ), superscripts denoting the fluid (F) or solid (S) domains are omitted in quantities such as q adj and T.
The adjoint stresses and adjoint heat flux are defined as Variations in µ are taken into account by differentiating the Sutherland law [11] with respect to the fluid's temperature. I MF µ t comprises field integrals over Ω F which include variations in µ t . Differentiating Equation (5) (first option), source terms to the adjoint Spalart-Allmaras equation arise, as follows

Development of the I SA Integral
Substituting the transport equation ofν, the second field integral of Equation (9) becomes where I SA C , I SA D and I SA S contain the convection, diffusion, and source terms of the Spalart-Allmaras model, namely The first two integrals can be rewritten using the divergence theorem as Integrals including variations inμ (whereμ = ρν) contribute to the field adjoint Spalart-A-llmaras equation (FAE SA ). Those Since the production and destruction ofν are functions of ρ, vorticity ζ, distance d, andμ, I SA S takes the form where Sν = Dν −Pν and the field integral including variations in distance contributes to the field adjoint distance equation (FAE D ). The partial derivatives of the production and destruction ofν with regard to ρ,μ, d, and ζ can be found in [15]. The last integral is rewritten as

Development of the I D Integral
Substituting the residual of the Eikonal equation in Equation (9) and applying the divergence theorem, I D becomes

Development of the I S Integral
In Ω S , using the heat conduction equation and the divergence theorem, the corresponding integral in Equation (9) takes the form as the adjoint heat fluxes in the solid domain,

Field Adjoint Equations (FAEs) and Adjoint Boundary Conditions (ABCs)
The multipliers of variations in the state variables give rise to FAE MF , FAE SA , FAE D , and FAE S . These read T a n k δq k δb i + q adj k n k The adjoint no-slip condition (Ψ 2 = Ψ 3 = Ψ 4 = 0) eliminates surface integrals with variations in τ km and p. Terms with variations in heat flux and T are eliminated by satisfying the adjoint adiabatic and FSI conditions. The former read q adj,F/S m n F/S m = 0 along the adiabatic walls of Ω F and Ω S . The adjoint FSI conditions are Ψ 5 = T a and q adj,F k n F k = −q adj,S k n S k . Finally, settingν a and d a to zero over the boundaries of Ω F eliminates the surface integrals which contain variations inν and d.

Expression for SDs
The remaining field and surface integrals comprise the formula for computing the gradient of F where Grid sensitivities along the parameterized surfaces (here, the FSI and the cooling channel boundaries) are computed directly by the differentiation of the parameterization tool. Then, the spring analogy method undertakes their propagation into Ω F and Ω S .

The PUMA Software
PUMA is a GPU-accelerated software which involves a RANS equations' solver, a heat conduction equation's solver, and their adjoint counterparts along with a set of shape parameterization and grid deformation tools. It runs on a GPU cluster using the shared on-board memory for data transfer between GPUs on the same computational node and the MPI protocol between GPUs on different nodes. Data communications overlap with computations.
The flow and their adjoint equations are discretized on unstructured grids following the vertex-centered finite-volume approach. A multi-stage Runge-Kutta scheme, applied in pseudo-time, with residual smoothing is used to solve the equations. A second-order accurate, central differences scheme with artificial dissipation discretizes the inviscid fluxes. Residual smoothing employs the first-order accurate Jacobian of the fluid's residuals, which is computed in double-, though stored in single-precision arithmetic. The mixed-precision arithmetic minimizes GPU memory transactions without jeopardizing solution accuracy. The heat conduction equation is also discretized using vertex-centered finite volumes and solved using a Krylov-based solver.
The flow and heat conduction solvers are loosely coupled by means of Aitken's dynamic relaxation formula. Figure 1 illustrates a schematic representation of the coupling between the flow and heat conduction solvers. Temperature fields computed by the heat conduction solver in Ω S after performing a number of internal iterations are scaled according to Aitken's formula and imposed as boundary conditions to the flow solver. The flow solver preforms a number of pseudo-time steps, and the computed heat fluxes on the FSI are communicated back to the heat conduction solver. The flow solver computes heat fluxes along the FSI by integrating the energy equation over the Ω F finite volumes next to it; for instance, for an FSI node i, where j are its neighbors and Φ inv energ , Φ visc energ the inviscid and viscous energy fluxes crossing the finite volume interfaces between i and j, respectively. The stability of the coupling scheme is enhanced by computing and communicating the Jacobians of the computed heat fluxes q F m n F m i with respect to T F i along the FSI to the heat conduction solver. Heat fluxes along with their Jacobians are computed and communicated at each pseudo-time iteration of the heat conduction solver. The CHT coupled adjoint equations are solved in a similar way.

CHT Analysis of the C3X Turbine
The geometry of the C3X turbine linear cascade can be found in [3]. The blade, made of stainless steel with density ρ S = 7900 kg/m 3 and heat capacity C S = 586.15 J/(kg K), is cooled by 10 straight channels. Its thermal conductivity is a function of T S as κ S = 6.8110 + 0.020716T S . Its span is 7.62 cm long, with an axial chord of 7.816 cm. According to Case 158, for the hot gas, p I t = 243 700 Pa, T I t = 808 K, and Tu I = 8.3% (turbulence intensity) at the inlet, whereas p O = 142.530 Pa at the outlet. For each cooling channel, experimental data for the average temperature and the coolant flow rate are available; to meet them, an iterative procedure for adjusting the total pressure and total temperature at the cooling channels' inlets is applied.
An unstructured 3D grid with ∼3.5 M nodes was generated with ∼2 M for the hot gas, ∼630 K for the coolant flow, and ∼950 K for the solid domain. The grid has structured layers close to the pressure and suction sides and the cooling channels to account for the large spatial gradients of T F and T S in these regions. The rest of the fluid and solid domain is filled with prismatic elements. Both the Spalart-Allmaras and k−ω SST turbulence models were used in the analysis phase. The T S fields computed at the hub, shroud, and midspan of the solid blade are presented in Figure 2 (top left), illustrating the effect of the cooling fluid flow at the blade temperature. Pressure, temperature, and heat transfer coefficient distributions along the blade midspan are compared with measurements in the same figure. A good agreement is observed for both models, though the temperature along the suction side and close to the leading edge is overestimated due to the absence of a transition model. determined by the channel's diameter, the measured flow rate, and the coolant's molecular viscosity. Cr is a correction for a fully developed thermal boundary layer to account for thermal entrance region effects and ranges from 1.03 to 1.12 for the 10 channels, [3]. An alternative computation is based on the post-processing of the heat flux and temperature fields at the blade midspan computed by the 3D CHT run. Figure 3 summarizes the heat transfer coefficients computed based on the two aforementioned approaches. The heat transfer coefficients based on the second approach, with both turbulence models, are slightly different, being ∼15% lower than those based on the first approach. To measure the effectiveness of the computed heat transfer coefficients, 2D CHT analyses were performed based on them with the Spalart-Allmaras turbulence model. Figure 4 shows the computed temperature and heat transfer coefficient distributions at the blade midspan resulting from two 2D CHT analyses employing the heat transfer coefficients from the two approaches and also includes the 3D CHT analysis results as well as measurements. It can be seen that, using the heat transfer coefficients computed by the second approach, the 2D CHT analysis produces almost the same results with the 3D run. On the other hand, with the first approach, lower T S is computed. Though it is not shown here, both 2D runs obtain practically the same pressure field at the blade midspan with the 3D run. This was expected since the pressure field is mostly affected by the blade shape.

Verification of the Adjoint Sensitivities
Prior to the optimization studies, the accuracy of the adjoint-based sensitivity derivatives is assessed. For this purpose, the sensitivity derivatives of the objective functions and flow constraints (F 1 , F 2 , F 3 and F 4 ) are computed by the developed continuous adjoint method and compared with finite differences (FDs). In this comparison, the 2D configuration with the heat transfer coefficients computed by the second approach and the Spalart-Allmaras turbulence model are used. The turbine blade shpae is parameterized using a 7 × 3 volumetric NURBS control grid. The control points at the leading and trailing edges are fixed while the rest are free to move in the axial and pitch-wise directions.
In total, 38 design variables are used. The cooling holes, though potentially controlled by the same morphing box, are fixed without being affected by changes in the design variables. Figure 5 presents the computed adjoint and FD-based sensitivity derivatives which are in very good agreement. The influence of the frozen turbulence assumption (which is frequently made in the literature) is investigated. As shown in Figure 5, avoiding including the turbulence model equation into the adjoint formulation does not affect the gradients of the total pressure drop (F 1 ), inlet capacity (F 3 ), and exit flow angle (F 4 ). However, this is not the case for the highest blade temperature (F 2 ), where the frozen turbulence assumption significantly reduces the accuracy of the computed derivatives. For instance, for some design variables, this leads to derivatives that even have a wrong sign.

CHT Optimization of the C3X Turbine
The CHT optimization of the C3X turbine blade follows: The target was to minimize both the total pressure drop and the highest temperature of the blade while the inlet capacity and the hot gas exit flow angle were constrained to change up to 0.1% from their datum values. A first shape optimization was performed targeted at minimizing the total pressure drop under the above flow constraints. The 7 × 3 volumetric NURBS control grid of Figure 6a was used to control the blade shape. As in the verification study of the previous section, the control points at the leading and trailing edges remained constant, while the rest were free to move in the axial and pitch-wise directions. In this optimization run, in order to avoid the intersection of the cooling channels with the blade sides, the locations of their centers were controlled by the same control grid, while their radii were constant. In total, 38 design variables were used. To adapt the fluid domain grid to the changing boundaries, the spring analogy technique [16] was employed. The grid inside the blade was generated from scratch based on the new blade shape and cooling channel locations. Figure 6b shows the design variable bounds (used as extra geometric constraints). The method of moving asymptotes [17] was used, and its convergence is plotted in Figure 6b. The total pressure drop has decreased by ∼21% after 20 cycles by satisfying the flow constraints.  The best solution of the first optimization has almost the same highest blade temperature with the baseline, since this is mainly affected by the location of the cooling channels, which were only slightly changed by the control box of the first optimization. Thus, the outcome of the first optimization was further optimized for minimum highest temperature over the blade, aiming at refining the location of the channels. In this second optimization, the blade shape was not allowed to change, while the cooling channels were free to move leading to 20 design variables (2 per channel). Figure 7a shows the design variable bounds which were chosen so as to avoid channels overlapping with each other. To prevent the cooling channels from overlapping with the blade sides, the distance of each cooling channel from the airfoil sides was constrained to be higher than 30% of their radius. These 20 geometric constraints (2 per channel) ensured the generation of valid grids anew in each optimization cycle and kept the channels away from the blade sides, avoiding extremely high ∇T S and high thermal stresses. After 45 optimization cycles, the blade's highest T S was reduced by ∼30 K, still satisfying all geometric constraints. Figure 7b shows the convergence of the optimization algorithm, where the blade's temperature is normalized with T re f . For clarity, only the distance of the 10th cooling channel (the one closer to the leading edge) from the blade's sides is plotted in Figure 7b, where distances are normalized with the channels' radius. The best solution of the second optimization has a slightly increased (no more than 0.1%) total pressure drop compared to the best solution of the first run. This small increase was expected since reducing the highest T S on the blade results in a higher heat flux from the hot gas to the blade. At the same time, inlet capacity and hot gas exit angle differ no more than 0.1% from their datum values.    Figure 8 presents the Mach number fields computed for the baseline and the optimized configuration (i.e., the outcome of the second optimization run). It is clear that the reduction in total pressure drop is mostly achieved by reducing the flow velocity in the passage throat region. Figure 9 presents the computed T S and ∇T S fields inside the blade for the baseline and the optimized configuration. In the latter, the first two cooling channels moved toward the leading edge where T S is high, while the remaining channels moved toward the trailing edge where the highest T S exists. Thanks to the geometric constraints, the cooling channels did not move too close to the blade sides, and the highest ∇T S changed by less than ∼0.3% compared to the baseline. The optimized blade cooling system can also be seen in Figure 10.

Conclusions
The mathematical development of the continuous adjoint method for computing sensitivity derivatives in CHT problems involving turbulent compressible fluid flows was presented. The turbulence model equations are included in the adjoint formulation without making the frozen turbulence assumption, thus leading to an adjoint model that is consistent, in the continuous sense, to the one used in the CHT analysis. This was implemented in the in-house GPU accelerated software PUMA and used to optimize the internally cooled C3X turbine blade. Two optimization runs were performed. The first aims at reducing the total pressure drop while constraining the turbine inlet capacity and hot gas exit angle, allowing changes up to 0.1% compared to the baseline configuration. The second optimization started with the best solution of the first run, aiming at refining the cooling channel locations for minimum highest temperature over the blade. The optimized configuration has ∼20% less total pressure drop and ∼30 K lower highest blade temperature. The geometric constraints used in the second optimization run constrained the cooling channels from moving too close to the blade sides, increasing the highest temperature gradient magnitude inside the blade less than 0.3% of the baseline. The results show that the developed continuous adjoint CHT method, coupled with a gradient-based algorithm, provides a useful tool for the CHT optimization of internally cooled turbine blades. They also reveal possible losses in accuracy, regarding the gradient computation, from making the frozen turbulence assumption in the development of the adjoint method.
Author Contributions: Conceptualization, K.G.; software, X.T., K.T., V.A. and M.K.; validation, X.T. and M.K.; supervision, K.G. All authors contributed to writing and editing this paper. All authors have read and agreed to the published version of the manuscript.