Next Article in Journal
Experimental Aerodynamic and Aeroelastic Investigation of a Highly-Loaded 1.5-Stage Transonic Compressor with Tandem Stator
Previous Article in Journal
Experimental Aerodynamic and Aeroelastic Investigation of a Transonic Compressor Rotor with Reduced Blade Count
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

Parallel CFD & Optimization Unit, Laboratory of Thermal Turbomachines, School of Mechanical Engineering, National Technical University of Athens, 15780 Athens, Greece
*
Author to whom correspondence should be addressed.
This paper is an extended version of our paper published in the Proceedings of the 14th European Turbomachinery Conference, Gdansk, Poland, 12–16 April 2021.
Int. J. Turbomach. Propuls. Power 2021, 6(2), 20; https://0-doi-org.brum.beds.ac.uk/10.3390/ijtpp6020020
Submission received: 28 May 2021 / Revised: 7 June 2021 / Accepted: 9 June 2021 / Published: 15 June 2021

Abstract

:
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 flow 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 flow constraints are included in the performed optimizations, keeping the cooling channels away from the airfoil sides and retaining the turbine inlet capacity and flow turning.

1. 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 gradient-based 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 shape optimizations for turbulent incompressible flows is developed in the OpenFOAM environment. There, the shape optimization of a 2D internally cooled turbine cascade and a 3D car engine cylinder cooling system are presented.
The internally cooled C3X turbine blade, reported in [3], is a well-known test case that has widely been used for validation purposes [4,5] and is used in this paper, as well. The optimization of the cooling system of the C3X blade can be found in [6,7], which both use evolutionary algorithms, and in [8] using the globally convergent method of moving asymptotes. In the latter, the objective function derivatives with respect to the design variables needed by the optimization algorithm are computed based on forward finite differences.
This paper is an extended version of [9] and deals with the continuous adjoint method for computing the gradient of the objective and constraint functions with respect to the design variables; the adjoint to the turbulence model PDE is also formulated and solved. Compared to the material presented in [9], the development of the continuous adjoint for CHT problems is presented in detail. The adjoint solver is used within a gradient-based algorithm to optimize the C3X turbine blade, aiming at minimizing the total pressure drop and the highest blade temperature, by controlling the blade shape and its cooling holes’ locations. This optimization is performed in two steps, in which constant heat transfer coefficients are imposed in the cooling holes. The latter are based on 3D CHT analysis runs instead of using the empirical expression of [3,9].
This paper is organized as follows. The governing equations and the adjoint development are presented in Section 2 and Section 3, respectively. The in-house CFD software PUMA [10], running on GPUs, is presented in Section 4 and used for the CHT analysis of the C3X turbine blade in Section 5. The adjoint-based computed sensitivity derivatives are compared with finite differences in Section 6. Section 7 presents the CHT optimization of the C3X turbine blade.

2. 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; M F ) equations for compressible fluid flows read
R n M F = f n k i n v x k f n k v i s x k = 0
where f k i n v   =   ρ v k   ρ v k v 1   +   p δ 1 k   ρ v k v 2   +   p δ 2 k   ρ v k v 3   +   p δ 3 k   ρ v k h t T are the inviscid fluxes and f k v i s   =   0   τ 1 k   τ 2 k   τ 3 k   v τ k + q k F 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 δ k m the Kronecker symbol. τ k m   =   μ + μ t v k x m + v m x k 2 3 δ k m v x 2 3 δ k m ρ k is the stress tensor based on the Boussinesq assumption. The heat fluxes are q k F   =   κ 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
R ν ˜   =   x k ρ v k ν ˜ ρ σ x k ν + ( 1 + C b 2 ) ν ˜ ν ˜ x k + ρ C b 2 ν ˜ σ x k ν ˜ x k P ν ˜ + D ν ˜   =   0
R k   =   x k ρ v k k x k μ + σ k μ t k x k P ˜ k + β ρ ω k   =   0
R ω   =   x k ρ v k ω x k μ + σ ω μ t ω x k γ ν t P k + β ρ ω 2 2 1     F 1 ρ σ ω 2 ω k x k ω x k   =   0
where P ν ˜ = ρ C b 1 1 f t 2 S ˜ ν ˜ and D ν ˜ = ρ C w 1 f w C b 1 κ 2 f t 2 ν ˜ d 2 are the production and destruction of ν ˜ , and d stands for the distance from the nearest solid wall boundary. μ t can be computed via
μ t = ρ ν ˜ f v 1   o r   μ t = ρ α 1 k m a x ( α 1 ω , S F 2 )
The constants κ , C b 1 , C b 2 , and C w 1 along with the expressions of f v 1 , f w , f t 2 , and S ˜ are given in [12]. The last two quantities are functions of vorticity ζ   =   2 W k m W k m , where W k m   =   1 2 v k x m     v m x k . Constants β , α 1 , σ ω 2 , the expression for the production of k P ˜ k = m i n ( 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,
R d F = d x k d x k 1 = 0
The last state equation, which stands for the heat conduction equation over Ω S , reads
R S   = q k S x k = 0
where q k S   = κ 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 k F n k F = q k S n k S must be satisfied. n k F , n k S are the normal vectors to the FSI pointing toward Ω F and Ω S , respectively.

3. 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
F 1 = p ¯ t I p ¯ t O , F 2 = s o l i d T S e α T S s o l i d e α T S , F 3 = m ˙ I T ¯ t I p ¯ t I , F 4 = t a n 1   v ¯ 2 O v ¯ 1 O
where the highest blade temperature is approximated by the above differentiable expression where α takes on a large value. m ˙ I , p ¯ t I , T ¯ t I are the massflow and the mass-averaged total pressure and total temperature at the hot gas inlet, respectively, and p ¯ t O , v ¯ 1 O , v ¯ 2 O 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 a u g ) 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 a u g with respect to the design variables b i gives
δ F a u g δ b i = δ F δ b i + Ω F   Ψ n δ R n M F δ b i d Ω I M F + Ω F   ν ˜ a δ R ν ˜ δ b i d Ω I S A + Ω F   d a δ R d δ b i d Ω I D + Ω S   T a δ R S δ b i d Ω I S
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, δ F a u g δ b i and δ F δ b i are the same. For any function Φ of the state variables, the expression
δ Φ δ b i   =   Φ b i + Φ x k δ x k δ b i
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]
δ δ b i Φ x   =   x δ Φ δ b i     Φ x k x δ x k δ b i
After a lengthy mathematical development, δ F a u g δ b i 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 ( S D s ).

3.1. Development of the I M F Integral

Substituting the residuals of the MF equations in Equation (9), the first field integral of this equation takes the form
I M F = Ω F   Ψ n δ δ b i f n k i n v x k d Ω I i n v M F Ω F   Ψ n δ δ b i f n k v i s x k d Ω I v i s M F
Using the divergence theorem, its inviscid counterpart becomes
I i n v M F = S F   Ψ n n k δ f n k i n v δ b i d S I B Ω F   A n m k Ψ n x k δ U m δ b i d Ω F A E M F Ω F   Ψ n f n k i n v x x k δ x δ b i d Ω S D
where U   =   [ ρ   ρ v 1   ρ v 2   ρ v 3   ρ E ] (in a 3D flow) are the conservative flow variables, E the total energy per unit mass, and A n m k = f n k i n v U m the inviscid flux Jacobian matrix. An arrow (→) underneath any underbraced term indicates where this term contributes. For instance, field integrals including variations in U m contribute to the field adjoint mean-flow equations ( F A E M F ) . Integrals which contain sensitivities of geometric quantities contribute to the S D s . Finally, I B gathers all surface integrals. Using the divergence theorem
I v i s M F   = S F   Ψ n n k δ f n k v i s δ b i d S I B + Ω F   Ψ n x k δ f n k v i s δ b i d Ω I v i s , Ω M F + Ω F   Ψ n f n k v i s x x k δ x δ b i d Ω S D
Substituting the expression of the viscous fluxes, applying the divergence theorem, and rearranging the indices, I v i s , Ω M F becomes
I v i s , Ω M F = Ω F   τ k m μ + μ t Ψ m + 1 x k + v m Ψ 5 x k δ μ δ b i d Ω + Ω F   Ψ 5 x k T x k k μ δ μ δ b i d Ω Ω F   K k δ V k δ b i d Ω F A E M F + Ω F   τ k m μ + μ t Ψ m + 1 x k + v m Ψ 5 x k δ μ t δ b i d Ω + Ω F   Ψ 5 x k T x k k μ t δ μ t δ b i d Ω I μ t M F Ω F   τ k m a d j v k x + q m a d j T x x m δ x δ b i d Ω S D + S F   τ k m a d j n m δ v k δ b i + q k a d j n k δ T δ b i d S I B
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 a d j and T. Arrays V, K are
V = p   v 1   v 2   v 3   T F T K = 0   τ 1 m a d j x m Ψ 5 x m τ 1 m τ 2 m a d j x m Ψ 5 x m τ 2 m τ 3 m a d j x m Ψ 5 x m τ 3 m q k a d j , F x k T
The adjoint stresses and adjoint heat flux are defined as
τ k m a d j = μ + μ t Ψ k + 1 x m + Ψ m + 1 x k 2 3 δ k m Ψ + 1 x + Ψ 5 x k v m + Ψ 5 x m v k 2 3 δ k m Ψ 5 x v q k a d j , F = κ F Ψ 5 x k
Variations in μ are taken into account by differentiating the Sutherland law [11] with respect to the fluid’s temperature. I μ t M F 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
I μ t M F = Ω F   τ k m μ + μ t Ψ m + 1 x k + v m Ψ 5 x k + Ψ 5 x k T x k k μ t μ t μ ˜ δ μ ˜ δ b i d Ω F A E S A

3.2. Development of the I S A Integral

Substituting the transport equation of ν ˜ , the second field integral of Equation (9) becomes
I S A = I C S A + I D S A + I S S A
where I C S A , I D S A and I S S A contain the convection, diffusion, and source terms of the Spalart–Allmaras model, namely
I C S A   = Ω F   ν ˜ a δ δ b i x k ρ v k ν ˜ d Ω I D S A   = Ω F   ν ˜ a δ δ b i ρ σ x k ν + ( 1 + C b 2 ) ν ˜ ν ˜ x k + ρ C b 2 ν ˜ σ x k ν ˜ x k d Ω I S S A = Ω F   ν ˜ a δ P ν ˜ δ b i d Ω + Ω F   ν ˜ a δ D ν ˜ δ b i d Ω
The first two integrals can be rewritten using the divergence theorem as
I C S A = Ω F   ν ˜ a x k v k δ μ ˜ δ b i d Ω F A E S A Ω F   ν ˜ a x k δ v k δ b i μ ˜ d Ω F A E M F + S F   μ ˜ ν ˜ a δ v k δ b i n k d S + S F   v k ν ˜ a δ μ ˜ δ b i n k d S I B Ω F   ν ˜ a x ρ v k ν ˜ x k δ x δ b i d Ω S D I D S A = Ω F   ν ˜ a ρ σ x k ν + ( 1 + C b 2 ) ν ˜ ν ˜ x k + ρ C b 2 ν ˜ σ x k ν ˜ x k δ ρ δ b i d Ω F A E M F + Ω F   ρ ν ˜ a σ C b 2 x k ν ˜ x k δ ν ˜ δ b i d Ω F A E M F   &   F A E S A + Ω F   ρ ν ˜ a σ C b 2 ν ˜ x ν ˜ x k x k δ x δ b i d Ω S D + Ω F   ρ ν ˜ a σ x ν + ( 1 + C b 2 ) ν ˜ ν ˜ x k x k δ x δ b i d Ω S D + Ω F   ρ C b 2 ν ˜ ν ˜ a σ x k δ δ b i ν ˜ x k d Ω Ω F   ρ ν ˜ a σ x k δ δ b i ν + ( 1 + C b 2 ) ν ˜ ν ˜ x k d Ω I g r a d ν ˜ S A
Integrals including variations in μ ˜ (where μ ˜ = ρ ν ˜ ) contribute to the field adjoint Spalart-A-llmaras equation ( F A E S A ) . Those including variations in ν ˜ contribute to both F A E S A and F A E M F as δ ν ˜ δ b i   =   ν ˜ ρ δ ρ δ b i + ν ˜ μ ˜ δ μ ˜ δ b i . The above field integrals which contribute neither to F A E M F , F A E S A nor to the S D s can further be analyzed using the divergence theorem, as
I g r a d ν ˜ S A =   Ω F   1 σ ν ˜ x k ρ ν ˜ a x k δ ν δ b i d Ω F A E M F + Ω F   1 + C b 2 σ ν ˜ x k ρ ν ˜ a x k δ ν ˜ δ b i d Ω F A E M F   &   F A E S A Ω F   x k 1 σ ν + ( 1 + C b 2 ) ν ˜ ρ ν ˜ a x k C b 2 σ ρ ν ˜ ν ˜ a x k δ ν ˜ δ b i d Ω F A E M F   &   F A E S A + S F   1 σ ν + ( 1 + C b 2 ) ν ˜ ρ ν ˜ a x k C b 2 σ ρ ν ˜ ν ˜ a x k δ ν ˜ δ b i n k d S I B S F   1 σ ρ ν ˜ a δ δ b i ν + ( 1 + C b 2 ) ν ˜ ν ˜ x k n k d S I B + S F   C b 2 σ ρ ν ˜ ν ˜ a δ δ b i ν ˜ x k n k d S I B + Ω F   1 σ ν + ( 1 + C b 2 ) ν ˜ ρ ν ˜ a x k C b 2 σ ρ ν ˜ ν ˜ a x k ν ˜ x x k δ x δ b i d Ω S D
Since the production and destruction of ν ˜ are functions of ρ , vorticity ζ , distance d, and μ ˜ , I S S A takes the form
I S S A = Ω F   ν ˜ a S ν ˜ ρ δ ρ δ b i d Ω F A E M F + Ω F   ν ˜ a S ν ˜ μ ˜ δ μ ˜ δ b i d Ω F A E S A + Ω F   ν ˜ a S ν ˜ d δ d δ b i d Ω F A E D + Ω F   ν ˜ a S ν ˜ ζ δ ζ δ b i d Ω I ζ S A
where S ν ˜     =     D ν ˜     P ν ˜ and the field integral including variations in distance contributes to the field adjoint distance equation ( F A E 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
I ζ S A =   Ω F   x m ν ˜ a S ν ˜ ζ 2 ζ W k m δ v k δ b i d Ω F A E M F Ω F   ν ˜ a S ν ˜ ζ 2 ζ W k m v k x l x m δ x l δ b i d Ω S D + S F   ν ˜ a S ν ˜ ζ 2 ζ W k m n m δ v k δ b i d S I B

3.3. Development of the I D Integral

Substituting the residual of the Eikonal equation in Equation (9) and applying the divergence theorem, I D becomes
I D = Ω F   x m 2 d x m d a δ d δ b i d Ω F A E D + S F   2 d x m n m d a δ d δ b i d S I B Ω F   2 d x m d x l d a x m δ x l δ b i d Ω S D

3.4. 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
I S   = S S   T a n k δ q k δ b i d S I B + Ω S   T a x k δ q k δ b i d Ω I q S + Ω S   T a q k x x k δ x δ b i d Ω S D
where, by defining q k a d j , S   =   κ S T a x k as the adjoint heat fluxes in the solid domain,
I q S = Ω S   T a x k T x k δ κ δ b i d Ω Ω S   q k a d j x k δ T δ b i d Ω F A E S + S S   q k a d j n k δ T δ b i d S I B Ω S   q k a d j T x x k δ x δ b i d Ω S D

3.5. Field Adjoint Equations (FAEs) and Adjoint Boundary Conditions (ABCs)

The multipliers of variations in the state variables give rise to F A E M F , F A E S A , F A E D , and F A E S . These read
R m Ψ = A n m k Ψ n x k K n V n U m + B μ F μ U m + B c F κ F U m + S m Ψ = 0 R ν ˜ a = v k ν ˜ a x k + 1 σ x k ν + ( 1 + C b 2 ) ν ˜ ρ ν ˜ a x k C b 2 σ x k ρ ν ˜ ν ˜ a x k 1 + C b 2 σ ν ˜ x k ρ ν ˜ a x k C b 2 σ ν ˜ a x k ν ˜ x k + S ν ˜ a = 0 R d a = x k 2 d x k d a + S d a = 0 R T a = q k a d j , S x k + B c S k S T S + F T S = 0
where B μ F , B c F , B c S stand for the multipliers of δ μ δ b i , δ κ F δ b i , and δ κ S δ b i , respectively. S Ψ are source terms added to F A E M F and result from the differentiation of the Spalart–Allmaras model. Similarly, S ν ˜ a and S d a are sources to F A E S A and F A E D , resulting from the differentiation of the mean-flow and turbulence model equations, respectively. It should be mentioned that F A E M F , F A E S A , and F A E S do not depend on the adjoint distance variable d a ; thus, F A E D is numerically solved after solving F A E M F or F A E S A and F A E S first and contributes to the computation of S D s . F T S is the contribution of F to F A E S , which is non-zero only for F 2 .
Gathering all surface integrals without variations in ν ˜ or d into a single integral gives
I B = S F   Ψ n n k δ f n k i n v δ b i Ψ n n k δ f n k v i s δ b i + τ k m a d j n m δ v k δ b i + q k a d j n k δ T δ b i d S S S   T a n k δ q k δ b i + q k a d j n k δ T δ b i d S
I B splits into integrals over (a) the inlet/outlet boundaries S I / O and (b) the solid wall boundaries S W of Ω F ; the latter include the FSI and adiabatic walls (hub and shroud). For the S I / O ,
I I / O B = S I / O   Ψ n A n m k n k δ U m δ b i + τ k m a d j n m Ψ 5 τ k m n m δ v k δ b i + q k a d j n k δ T δ b i d S
where variations in viscous stresses and heat flux are ignored. Eliminating the terms including variations in the flow variables gives rise to the adjoint inlet/outlet conditions, namely
Ψ n A n m k n k U m V j I / O + τ k m a d j n m Ψ 5 τ k m n m v k V j I / O + q k a d j n k T V j I / O + F V j I / O = 0
where V j I / O denotes any flow quantity extrapolated to S I / O from the interior of Ω F . For the inlet, V j I stands for the velocity magnitude, and for the outlet, V j O are the outgoing Riemann variables. For the solid walls of Ω F , taking into account the no-slip condition, Equation (15) becomes
I W B = S W F   Ψ m + 1 n m δ p δ b i d S S W F   Ψ m + 1 n k δ τ k m δ b i d S S W F   Ψ 5 δ q k n k δ b i d S + S W F   q k a d j n k δ T δ b i d S S W S   T a δ q k n k δ b i d S + S W S   q k a d j n k δ T δ b i d S + S W F   Ψ 5 q k δ n k δ b i d S + S W S   T a q k δ n k δ b i d S S D
The adjoint no-slip condition ( Ψ 2   =   Ψ 3   =   Ψ 4   =   0 ) eliminates surface integrals with variations in τ k m and p. Terms with variations in heat flux and T are eliminated by satisfying the adjoint adiabatic and FSI conditions. The former read q m a d j , F / S n m F / S = 0 along the adiabatic walls of Ω F and Ω S . The adjoint FSI conditions are Ψ 5 = T a and q k a d j , F n k F = q k a d j , S n k S . Finally, setting ν ˜ a and d a to zero over the boundaries of Ω F eliminates the surface integrals which contain variations in ν ˜ and d.

3.6. Expression for S D s

The remaining field and surface integrals comprise the formula for computing the gradient of F
δ F δ b i = I M F S D + I S A S D + I D S D + I S S D
where
I M F S D = Ω F   Ψ n f n k i n v x f n k v i s x τ k m a d j v m x q k a d j T x x k δ x δ b i d Ω + S W F   Ψ 5 q k δ n k δ b i d S I S A S D = Ω F   ν ˜ a x ρ v k ν ˜ x k δ x δ b i d Ω + Ω F   ρ ν ˜ a σ x ν + ( 1 + C b 2 ) ν ˜ ν ˜ x k + C b 2 ν ˜ x ν ˜ x k x k δ x δ b i d Ω + Ω F   1 σ ν + ( 1 + C b 2 ) ν ˜ ρ ν ˜ a x k C b 2 σ ρ ν ˜ ν ˜ a x k ν ˜ x x k δ x δ b i d Ω Ω F   ν ˜ a S ν ˜ ζ 2 ζ W k m v k x x m δ x δ b i d Ω I D S D = Ω F   2 d x m d x d a x m δ x δ b i d Ω I S S D = Ω S   T a q k x q k a d j T x x k δ x δ b i d Ω + S W S   T a q k δ n k δ b i d S
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 .

4. 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,
q m F n m F i = j N e i g h s ( i ) Φ e n e r g i n v , i j j N e i g h s ( i ) Φ e n e r g v i s c , i j
where j are its neighbors and Φ e n e r g i n v , Φ e n e r g v i s c 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 m F n m F i with respect to T i F 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.

5. 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.020716 T 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 t I   =   243700   Pa , T t I   =   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.
The pressure, temperature, and heat transfer coefficient distributions at the blade midspan were also computed by a 2D CHT analysis, in which the average temperature and coolant flow rate per cooling channel were replaced by heat transfer coefficient values (since the 2D analysis refrains from computing the coolant flow). The heat transfer coefficients per cooling channel result from h   =   κ Nu D / D , where κ is the coolant thermal conductivity, D the channel diameter, and Nu D the Nusselt number. For turbulent flows in smooth pipes, the latter is given by Nu D   =   0.022 Cr Pr 0.5 Re D 0.8 , where the Reynolds number Re D is 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.

6. 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.

7. 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 r e 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.

8. 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.

Funding

This research was funded by European Union’s Horizon 2020 research and innovation program under Grant Agreement No 769025 (MADELEINE).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Data sharing not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

CFDComputational fluid dynamics
CHTConjugate heat transfer
FAEField adjoint equation
FDsFinite differences
FSIFluid–solid interface
GPUGraphics processing unit
NURBSNon-uniform rational B-splines
RANSReynolds-averaged Navier–Stokes
SDSensitivity derivative

References

  1. Moretti, R.; Errera, M.; Couaillier, V.; Feyel, F. Stability convergence and optimization of interface treatments in weak and strong thermal fluid-structure interaction. Int. J. Therm. Sci. 2017, 126, 23–37. [Google Scholar] [CrossRef] [Green Version]
  2. Gkaragkounis, K.; Papoutsis-Kiachagias, E.; Giannakoglou, K. The continuous adjoint method for shape optmization in Conjugate Heat Transfer problems with turbulent incompressible flows. Appl. Therm. Eng. 2018, 140, 351–362. [Google Scholar] [CrossRef]
  3. Hylton, L.; Mihelc, M.; Turner, E.; Nealy, D.; York, R. Analytical and Experimental Evaluation of the Heat Transfer Distribution over the Surfaces of Turbine Vanes. In NASA-CR-168015; Detroit Diesel Allison: Indianapolis, IN, USA, 1983. [Google Scholar]
  4. Andrei, L.; Andreini, A.; Facchini, B.; Winchler, L. A decoupled CHT procedure; application and validation on a gas turbine vane with different cooling configurations. Energy Procedia 2014, 45, 1087–1096. [Google Scholar] [CrossRef] [Green Version]
  5. De Marinis, D.; de Tulio, M.; Napolitano, M.; Pascazio, G. A conjugate-heat-transfer immersed-boundary method for turbine cooling. Energy Procedia 2015, 82, 215–221. [Google Scholar] [CrossRef] [Green Version]
  6. Nowak, G.; Wroblewski, W. Optimization of blade cooling system with use of conjugate heat transfer approach. Int. J. Therm. Sci. 2011, 50, 1770–1781. [Google Scholar] [CrossRef]
  7. Mazaher, K.; Zeinalpout, M.; Bokaei, H. Turbine blade cooling passages optimization using reduced conjugate heat transfer methodology. Appl. Therm. Eng. 2016, 103, 1228–1236. [Google Scholar] [CrossRef]
  8. Wang, B.; Zhang, W.; Xie, G.; Xu, Y.; Xiao, M. Multiconfiguration Shape Optimization of Internal Cooling Systems of a Turbine Guide Vane Based on Thermomechanical and Conjugate Heat Transfer Analysis. J. Heat Transf. 2015, 136, 061004. [Google Scholar] [CrossRef]
  9. Trompoukis, X.; Tsiakas, K.; Asouti, V.; Kontou, M.; Giannakoglou, K. Continuous Adjoint Shape Optimization of Internally Cooled Turbine Blade. In Proceedings of the 14th European Turbomachinery Conference, Gdansk, Poland, 12–16 April 2021. [Google Scholar]
  10. Asouti, V.; Trompoukis, X.; Kampolis, I.; Giannakoglou, K. Unsteady CFD computations using vertex-centered finite volumes for unstructured grids on Graphics Processing Units. Int. J. Numer. Methods Fluids 2011, 67, 232–246. [Google Scholar] [CrossRef]
  11. Sutherland, W. The Viscosity of Gases and Molecular Force. Phil. Mag. 1893, 36, 507–531. [Google Scholar] [CrossRef] [Green Version]
  12. Spalart, P.; Allmaras, S. A one-equation turbulence model for aerodynamic flows. Rech. Aerosp. 1994, 1, 5–21. [Google Scholar]
  13. Menter, F.; Kuntz, M.; Langtry, R. Ten Years of Industrial Experience with SST Turbulence Model. Heat Mass Transf. 2003, 4, 625–632. [Google Scholar]
  14. Papadimitriou, D.; Giannakoglou, K. A continuous adjoint method with objective function derivatives based on boundary integrals, for inviscid and viscous flows. Comput. Fluids 2007, 36, 325–341. [Google Scholar] [CrossRef]
  15. Zymaris, A.; Papadimitriou, D.; Giannakoglou, K.; Othmer, C. Continuous adjoint approach to the Spalart-Allmaras turbulence model for incompressible flows. Comput. Fluids 2009, 38, 1528–1538. [Google Scholar] [CrossRef]
  16. Degand, C.; Farhat, C. A three-dimensional torsional spring analogy method for unstructured dynamic meshes. Comput. Struct. 2002, 80, 305–316. [Google Scholar] [CrossRef]
  17. Svanberg, K. A class of globally convergent optimization methods based on conservative convex separable approximations. SIAM J. Optim. 2002, 12, 555–573. [Google Scholar] [CrossRef] [Green Version]
Figure 1. Schematic representation of the coupling of the flow and heat conduction solvers. A similar coupling applies to the adjoint solver.
Figure 1. Schematic representation of the coupling of the flow and heat conduction solvers. A similar coupling applies to the adjoint solver.
Ijtpp 06 00020 g001
Figure 2. Computed T S fields (top left) at the solid blade’s hub, midspan, and shroud. Computed pressure (top right), temperature (bottom left), and heat transfer coefficient (bottom right) distributions across the blade midspan compared with measurements. Pressure is non-dimensionalized by the inlet total pressure of 243,700 Pa, temperature by 811 K , and the heat transfer coefficient by per-mode=symbol 1135.6   W / m 2 K .
Figure 2. Computed T S fields (top left) at the solid blade’s hub, midspan, and shroud. Computed pressure (top right), temperature (bottom left), and heat transfer coefficient (bottom right) distributions across the blade midspan compared with measurements. Pressure is non-dimensionalized by the inlet total pressure of 243,700 Pa, temperature by 811 K , and the heat transfer coefficient by per-mode=symbol 1135.6   W / m 2 K .
Ijtpp 06 00020 g002
Figure 3. Comparison of the heat transfer coefficients of the cooling channels computed by the two approaches.
Figure 3. Comparison of the heat transfer coefficients of the cooling channels computed by the two approaches.
Ijtpp 06 00020 g003
Figure 4. Temperature (left) and heat transfer coefficient (right) distributions along the blade contour obtained by 3D and 2D CHT runs. The 2D runs use the channels’ heat transfer coefficients computed by the two approaches presented in the text.
Figure 4. Temperature (left) and heat transfer coefficient (right) distributions along the blade contour obtained by 3D and 2D CHT runs. The 2D runs use the channels’ heat transfer coefficients computed by the two approaches presented in the text.
Ijtpp 06 00020 g004
Figure 5. Verification of the sensitivity derivatives computed by the adjoint method with the differentiated turbulence (DT) model and with the frozen turbulence (FT) assumption for F 1 (top left), F 2 (top right), F 3 (bottom left), and F 4 (bottom right). Finite differences are used as the reference values.
Figure 5. Verification of the sensitivity derivatives computed by the adjoint method with the differentiated turbulence (DT) model and with the frozen turbulence (FT) assumption for F 1 (top left), F 2 (top right), F 3 (bottom left), and F 4 (bottom right). Finite differences are used as the reference values.
Ijtpp 06 00020 g005
Figure 6. First optimization run.
Figure 6. First optimization run.
Ijtpp 06 00020 g006
Figure 7. Second optimization run.
Figure 7. Second optimization run.
Ijtpp 06 00020 g007
Figure 8. Computed Mach number fields for the baseline and optimized configurations.
Figure 8. Computed Mach number fields for the baseline and optimized configurations.
Ijtpp 06 00020 g008
Figure 9. Computed temperature (left) and temperature gradient magnitude (right) fields for the baseline and the optimized configurations.
Figure 9. Computed temperature (left) and temperature gradient magnitude (right) fields for the baseline and the optimized configurations.
Ijtpp 06 00020 g009
Figure 10. Baseline (red) and optimized (blue) blade’s shape and cooling system.
Figure 10. Baseline (red) and optimized (blue) blade’s shape and cooling system.
Ijtpp 06 00020 g010
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Trompoukis, X.; Tsiakas, K.; Asouti, V.; Kontou, M.; Giannakoglou, K. Continuous Adjoint-Based Optimization of an Internally Cooled Turbine Blade—Mathematical Development and Application. Int. J. Turbomach. Propuls. Power 2021, 6, 20. https://0-doi-org.brum.beds.ac.uk/10.3390/ijtpp6020020

AMA Style

Trompoukis X, Tsiakas K, Asouti V, Kontou M, Giannakoglou K. Continuous Adjoint-Based Optimization of an Internally Cooled Turbine Blade—Mathematical Development and Application. International Journal of Turbomachinery, Propulsion and Power. 2021; 6(2):20. https://0-doi-org.brum.beds.ac.uk/10.3390/ijtpp6020020

Chicago/Turabian Style

Trompoukis, Xenofon, Konstantinos Tsiakas, Varvara Asouti, Marina Kontou, and Kyriakos Giannakoglou. 2021. "Continuous Adjoint-Based Optimization of an Internally Cooled Turbine Blade—Mathematical Development and Application" International Journal of Turbomachinery, Propulsion and Power 6, no. 2: 20. https://0-doi-org.brum.beds.ac.uk/10.3390/ijtpp6020020

Article Metrics

Back to TopTop