- Email: [email protected]
Aerodynamic optimization design for high pressure turbines based on the adjoint approach
Chinese Journal of Aeronautics, (2015), 28(3): 757–769
Chinese Society of Aeronautics and Astronautics & Beihang University
Chinese Journal of Aeronautics [email protected] www.sciencedirect.com
Aerodynamic optimization design for high pressure turbines based on the adjoint approach a,*
Chen Lei a b
, Chen Jiang
b
Heavy Engineering Research Association, Auckland 2104, New Zealand School of Energy and Power Engineering, Beihang University, Beijing 100191, China
Received 24 April 2014; revised 5 December 2014; accepted 2 March 2015 Available online 22 April 2015
KEYWORDS Adjoint method; Aerodynamic design; High pressure turbine; Optimization design; Objective function
Abstract A first study on the continuous adjoint formulation for aerodynamic optimization design of high pressure turbines based on S2 surface governed by the Euler equations with source terms is presented. The objective function is defined as an integral function along the boundaries, and the adjoint equations and the boundary conditions are derived by introducing the adjoint variable vectors. The gradient expression of the objective function then includes only the terms related to physical shape variations. The numerical solution of the adjoint equation is conducted by a finitedifference method with the Jameson spatial scheme employing the first and the third order dissipative fluxes. A gradient-based aerodynamic optimization system is established by integrating the blade stagger angles, the stacking lines and the passage perturbation parameterization with the quasi-Newton method of Broyden–Fletcher–Goldfarb–Shanno (BFGS). The application of the continuous adjoint method is validated through a single stage high pressure turbine optimization case. The adiabatic efficiency increases from 0.8875 to 0.8931, whilst the mass flow rate and the pressure ratio remain almost unchanged. The optimization design is shown to reduce the passage vortex loss as well as the mixing loss due to the cooling air injection. ª 2015 The Authors. Production and hosting by Elsevier Ltd. on behalf of CSAA & BUAA. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).
1. Introduction With the increasing need for high performance gas turbines to reduce the emission and the engine weight, the emerging trend * Corresponding author. Tel.: +64 92624751. E-mail addresses: [email protected] (L. Chen), chenjiang27@ buaa.edu.cn (J. Chen). Peer review under responsibility of Editorial Committee of CJA.
Production and hosting by Elsevier
is to use mathematical optimization techniques as an integral part of the aerodynamic design toolkit. Stochastic and gradient-based methods are the normal mathematical optimization algorithms. The stochastic method searches for a global optimal solution by monitoring the magnitude of the objective function, while the gradient-based method searches for a local optimal solution by monitoring the sensitivity of the objective function to the changes in the design variables. Both of the stochastic and the gradient-based methods tend to consume considerably computational resources for the cases with a large number of design variables. The adjoint method is another gradient-based approach especially for the optimization with numerous design variables.
http://dx.doi.org/10.1016/j.cja.2015.04.022 1000-9361 ª 2015 The Authors. Production and hosting by Elsevier Ltd. on behalf of CSAA & BUAA. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).
758 In this method, the adjoint system is conducted in a similar way as that in the optimal control problems. By introducing the adjoint variable vectors, the gradients of the objective function with respect to the design variables are calculated indirectly by solving the adjoint equations. Therefore, the sensitivity analysis is almost independent of the numbers of the design variables, and solving two sets of flow equations in one design cycle is nearly the total computing cost. When the flow equations and the adjoint equations are fully converged, the final gradients of the objective function with respect to the design variables can be obtained efficiently. Pironneau1 was the first to use the adjoint method in fluid mechanics, and then Jameson2 applied adjonit method in the aeronautical field. Combining the continuous adjoint method with CFD technology, Jameson developed the optimization design method which was applied to the transonic wing-body combinations.3 Moreover, the discrete adjoint method was also developed.4 ‘‘Continuous’’ and ‘‘discrete’’ methods symbolize two alternative approaches to deriving adjoint equations. As there is no clear quantitative comparison between the two approaches,5,6 they seem to achieve the same optimal goals. Both of the methods have performed well in the optimization for airfoils,7 wings,8 wing-body configurations9 and business jets.10 The adjoint method has been utilized in the area of turbomachinery. Li et al.11,12 used the continuous adjoint method based on Navier–Stokes and Euler equation respectively to conduct the aerodynamic optimization design for turbine blades, and the optimization system was validated by several numerical cases. Papadimitriou and Giannakoglou13 developed the continuous adjoint formulation to improve the aerodynamic performance of a 3D peripheral compressor blade cascade. Wang and He14 first proposed the adjoint nonreflective mixing-plane treatment method, and carried out the aerodynamic blading shape optimization design in a multi-stage turbomachinery environment. Luo et al.15 used the adjoint optimization method to reduce the secondary flow loss of turbine blades by redesigning the blade. Ji et al.16 combined the continuous adjoint method with thin shearlayer Navier–Stokes equations to construct an efficient sensitivity analysis optimization system for multi-stage turbomachinery blades, and the adjoint optimization code was validated through two compressor blades design cases. Zhang and Feng17 used the automatic differentiation tool to develop a discrete adjoint solver, and the optimization system was validated via a turbine cascade under the viscous flow environment. For the aerodynamic design of gas turbines, S2 surface design plays a crucial role in the entire design system.18 Blade design is based on the simulation results of S2 surface through flow calculation. This paper presents the results of the first study on the adjoint method applied to the S2 surface through flow calculation including the cooling air effect. The adjoint method is combined with the Euler equation with the source term to develop an efficient sensitivity analysis model for turbine blades (stagger angles and stacking lines) and the passage in a specified objective function. The validation of the optimization system is carried out via the case of aerodynamic optimization of a high pressure turbine.
L. Chen, J. Chen 2. Flow equations and solution methods The steady Euler equations of the curvilinear coordinates system are utilized to predict the aerodynamic performance of the gas turbine on S2 surface. The effects of the viscous losses, the leakages and the cooling air on the flow are concerned as the source terms in the right part of the Eq. (1). @ rU @ e @ e e þ Fþ G¼Q ð1Þ @t J @n @g e are the e G where U is the conservative flow variable vector; F; e convective flux vectors in the curvilinear coordinate systems; Q is the source term of the flow equations. @ni ; J ¼ detðKÞ Kij ¼ @xj e ¼ r Fnz þ Gnr þ H 1 nu F J r e ¼ r Fgz þ Ggr þ H 1 gu ; G J r 2
rm r m vz þ fz
e ¼ h þ h1 Q J 3
6 7 6 7 6 7 6 2 7 h ¼ 6 r m vr þ fr þ qðw þ xrÞ 7 6 7 6 r m v þ f qvðw þ 2xrÞ 7 u 4 5 u
r m H0 þ x2 r2 qv 2
3 0 6 r f @p þ p @ r f 7 6 J z @f @f J z 7 6 7 6r 7 @p @ r þ p f f 7 h1 ¼ 6 r @f J r 7 @f 6J 6 7 6 fu @p 7 @ fu 4 J @f þ p @f 5 J 0 where terms nz ; nr ; nu ; fz ; fr ; fu represent the partial derivatives @n @n @n @f @f @f ; ; ; ; ; ; F, G and H are the convective flux vectors in @z @r @u @z @r @u
the cylindrical coordinate system; m is the mass flow rate of the cooling air; vz ; vr and vu are the velocity components of the cooling air, and the injection of the cooling air is classified into nine types as shown in Fig. 1; H0 is the total enthalpy of cooling air; (fz ; fr ; fu ) are for accounting of viscous losses effects,
Fig. 1
Injection types of the cooling air.
Aerodynamic optimization design for high pressure turbines based on the adjoint approach and are calculated by the next relations: ðfz ; fr ; fu Þ ¼ u2 þvr/2 þw2 ðu; v; wÞ, / ¼ qT dr ; r is the entropy generation calculated by dt the loss coefficient, and there are four loss models (blade profile loss model, secondary flow loss model, trailing edge loss model and the angle of attack loss model) adopted in the code; u; v; w are the velocity in the axial, the radial and the tangential direction respectively; p is the static pressure and x is the angular velocity of rotor; h1 is for calculating the impact of is for the initial accelerablade thickness on the flow field; dx dt tion of the rotor. j is the number of the injection (1 6 j 6 9). j = 1: cooling air injection at the leading edge; j = 2: cooling air injection upstream the throat section; j = 3: cooling air injection at the trailing edge; j = 4: cooling air injection at the hub between the blade rows; j = 5: cooling air is injected upstream the blade row at the hub and then injected downstream the blade row; j = 6: cooling air is injected between the blade rows at the shroud; j = 7: cooling air is injected upstream the blade row at the shroud and then injected downstream the blade row; j = 8: cooling air is injected within the blade row at the hub; j = 9: cooling air is injected within the blade row at the shroud. Three types of boundary conditions are used to complete the numerical solution of the flow equation. Span-wise distributions of the total pressure, total temperature and velocity directions are set at the inlet, and the static pressure of the root is set at the outlet. It is the no flow through the boundary condition for the hub and the shroud. The Euler equation is numerically solved by using the highorder accurate Godunov scheme, and the explicit numerical scheme is employed first, followed by the implicit scheme. The implicit operator solves the linear systems of equations by scalar three diagonal solvers. Finally, the flow field is time-marched to a steady state in a pseudo time via the Euler time-marching approach. 3. Adjoint equations and solution methods Due to its lower memory requirements and easy implementation compared to the discrete adjoint method, the continuous adjoint approach is employed here. The detailed derivation of the adjoint equations from the Euler equations with source term (cooling air injection included) in a partial differential form is shown as below. 3.1. Adjoint equations The objective function can be expressed as an integral along the boundaries in the following general form, Z I¼ MdS ð2Þ S
where M is a function of the flow variables and the design variables, then Z Z dI ¼ MI dUdS þ MII ðddSÞ ð3Þ S
S @M . @U
where MI ¼ The subscripts I and II are used to distinguish the contributions due to the variations of @U and @S respectively.
759
The Euler flow Eq. (1) can be written as 1 @U þ RðU; SÞ ¼ 0 J @t
ð4Þ
For the fully converged flow equation, the variation of RðU; SÞ is dR ¼
@ e @ e e ¼0 dF þ dG dQ @n @g
ð5Þ
Then yields, 8 dR ¼ dRI þ dRII > > > e < fII ¼ AdU þ d F fII d F ¼ df FI þ d F e f f fII > d þ d ¼ BdU þ d G ¼ d G G G I II > > : e f f f f d Q ¼ d Q þ d Q ¼ Q dU þ d Q I
II
I
ð6Þ
II
where A and B are the Jacobian matrices, @e F G Q fI ¼ @ e f ¼ @e A¼f ; B¼G ; : FI ¼ @U Q I @U @U Integrating the Eq. (5) over the whole computational domain and introducing the adjoint variable vectors w into both sides of the equations Z Z @ e @ e e dV ¼ 0 ð7Þ wT dRdV ¼ wT dF þ dG dQ @n @g D D Regrouping Eqs. (6) and (7) Z wT dRdV ¼ wT ðdRI þ dRII ÞdV D ZD @ f @ f e I þ dRII dV ¼ 0 ð8Þ d FI þ d GI d Q ¼ wT @n @g D
Z
Performing integration by parts Z T @w f @wT f T f f w ðd FI nn þ d GI ng ÞdS d FI þ d GI dV @n @g S D Z Z f dV þ e wT d Q ZdV ¼0 I
Z
D
e ¼ wT where Z
D @ fII dF @n
ð9Þ
@ fII d Q f , and nn and ng are þ @g dG I
components of normal direction. Rearranging, Z Z T @w f @wT f fI ng ÞdS f dV wT ðdf d FI þ d GI þ wT d Q FI nn þ d G I @n @g S D Z @ f @ f f dV ¼ 0 ð10Þ d FII þ d GII d Q þ wT II @n @g D Subtracting Eq. (10) from Eq. (3) does not change the gradient in Eq. (3) Z Z Z fI ng ÞdS dI ¼ MI dUdS þ MII ðddSÞ wT ðdf FI nn þ d G S S S Z T @w f @wT f f dV d FI þ d GI þ w T d Q þ I @n @g D Z @ f @ f f ð11Þ wT d FII þ d G II d QII dV @n @g D Introducing the Jacobian matrices A and B of Eq. (6) into Eq. (11) Z Z dI ¼ MI dUdS þ MII ðddSÞ S S Z Z Z T e e CdV w ðAdUnn þ BdUng ÞdS þ XdV ð12Þ S
D
D
760
L. Chen, J. Chen
where 8 e I dU e ¼ @wT AdU þ @wT BdU þ wT Q
The subsonic outlet boundary condition is e @M @F wT nn ¼ 0 @q @q
where q ¼ q; u; v; w, and the adjoint variables w1 ; w2 ; w3 ; w4 are calculated by Eq. (19), while w5 is extrapolated from the interior fluid domain. At the wall boundary, Vn ¼ 0, and only one equation is required. The wall boundary is ð13Þ
D
where 8 T T ee > e I dU ee fII d Q f :C fII þ @ d G ¼ wT @ d F @n
@g
Eliminating the dependence of dI on dU achieves the adjoint equation and the boundary condition equation respectively @wT @wT eI ¼ 0 Aþ B þ wT Q @n @g
ð14Þ
MI wT ðAnn þ Bng Þ ¼ 0
ð15Þ
To solve the adjoint Eq. (14), a conservative form is firstly adopted and a pseudo time derivative term is added to conduct the Runge–Kutta time marching scheme. In addition, the characteristics of the adjoint equation are opposite to those of the Euler flow equation. Therefore, the adjoint equation shares the same Jacobian matrixes as those of the Euler flow equation with an additional negative sign, as shown in Eq. (16). @ 1 @w @w e T w AT BT Q w¼0 ð16Þ @t J @n @g 3.2. Adjoint boundary conditions The physical adjoint boundary conditions are defined by Eq. (15). Moreover, the propagation of the adjoint characteristics is in the opposite way to the physical information of the Euler equations. Thus, one physical boundary condition is required at the subsonic inlet and four are required at subsonic exit, and the other numerical boundary conditions (four conditions at the inlet and one condition at the outlet) are extrapolated from the interior fluid domains. The velocity vector k of the adjoint equation is also opposite to that of the flow equation ð17Þ
where Vn is the normal component of velocity, and a is the local speed of sound. Normally, it is assumed that the geometry of the inlet and the outlet is frozen during each optimization cycle; therefore, the subsonic inlet boundary condition is e @M @F wT nn ¼ 0 @p @p
w2 ng1 þ w3 ng2 þ w4 ng3 ¼
@M @p
ð20Þ
3.3. Objective function
II
k ¼ ½Vn ; Vn ; Vn ; ðVn þ aÞ; ðVn aÞT
ð19Þ
ð18Þ
and the inlet plane is assumed to be normal to the X axis. The adjoint variables w1 ; w2 ; w3 ; w4 are extrapolated from the interior fluid domain, while w5 is calculated through Eq. (18).
Entropy generation is defined as the objective function, which is a weighted sum of mass flow rate and it can be formulated as R 2 udA q~ s I ¼ þ # in 1 ð21Þ s0 m0 where s is the entropy generation given by Z Z p p 0 s¼ dA dA0 c c q q out in
ð22Þ
and s0 is the original entropy generation; m0 is the original n mass flow rate of the gas; ue ¼ unz þ vnr þ w ru ; # is a weight factor chosen to be 60 in the case. Then, the final expression of the objective function is Z @ f @ f f dV dI ¼ wT d FII þ d GII d Q II @n @g D Z Z @ fII dV þ wT d Q f dV fII þ @ d G dF ¼ wT II @n @g D D 0R 1 R Z fII nn dS @wT d F fII dVþ wT d F S D @n f dV A þ ¼ @ R wT d Q R @wT II T f f D w d GII ng dS D @g d GII dV S Z Z T @w f @wT f fII ng dS þ f dV d FII þ d GII þ wT d Q ¼ wT d G II @n @g D ZS gu g g d z w2 þ d r w3 þ d w4 pdS ¼ J J Jr S Z T nj @w @wT gj f dV d d Gj þ wT d Q Fj þ ð23Þ þ II @n J @g J D f is related to the variations of blade and passage where d Q II geometry, and its formulation is f ¼ ½d Q f1 ; d Q f2 ; d Q f3 ; d Q f4 ; d Q f5 dQ II and 8 1 f1 ¼ m dQ dJ > > > > > > f2 ¼ m vz d 1 þ fz d 1 @p d fz p @ d fz > dQ > r J J J J @f @f > > > 1 fr 1 1 > > 2 f3 ¼ m vr d þ d þ qw d þ > dQ > r J J Jr > < f f @ r 2qwx d 1J þ qx2 d Jr @p d Jr d p @f J @f > > > f > f > d Q4 ¼ m vu d 1J þ ru d 1J qvw d Jr1 > > > > > fu fu @p > 1 @ > > 2xqv d J @f d Jr p @f d Jr > > > : f d Q5 ¼ m H0 d 1J þ x2 qv d 1J
Aerodynamic optimization design for high pressure turbines based on the adjoint approach
761
3.4. Numerical solution
4. Optimization design
The finite difference method is employed. The time derivative term is time-marched by the second order Runge–Kutta time method, and the space derivative term is discretized by the Jameson’s spatial scheme employing the first and the third order artificial dissipative flux. The discretized equation is
Hicks-Henne hump functions19 are used as the parameterization method. In addition, the optimization algorithm adopts the quasi-Newton method of Broyden–Fletcher–Goldfarb–Sh anno (BFGS)20 together with a line search method.
T
4.1. Geometry parameterization
T
1 ðwnþ1 wn Þ Ai;j ðwiþ1;j wi1;j Þ Bi;j ðwi;jþ1 wi;j1 Þ ¼ þ J Dt 2Dn 2Dg T e e þQ w þd i;j
i;j
ð24Þ
where e d ¼ diþ12; j di12; j þ di; jþ12 di; j12 and diþ12; j ¼ tð2Þ Kiþ12;j ðwiþ1; j wi;j Þ tð4Þ Kiþ12;j ðwiþ2;j 3wiþ1; j þ 3wi; j wi1;j Þ
Three Hicks-Henne functions Eqs. (27)–(29) are used for the parameterization: one for passage perturbations, another one for blade stagger angle perturbations of different sections from hub to shroud, and the last one for the stacking line perturbations. K X dk fk ðxÞ
DyðxÞ ¼ di12;j ¼ tð2Þ Ki12;j ðwi;j wi1;j Þ tð4Þ Ki12;j ðwiþ1;j 3wi;j þ 3wi1;j wi2; j Þ
K X dk fk ðrÞ
DbðrÞ ¼ di; jþ12 ¼ tð2Þ Ki;jþ12 ðwi;jþ1 wi;j Þ tð4Þ Ki;jþ12 ðwi; jþ2 3wi;jþ1 þ 3wi;j wi; j1 Þ di; j12 ¼ tð2Þ Ki;j12 ðwi;j wi;j1 Þ tð4Þ Ki;j12 ðwi; jþ1 3wi;j þ 3wi;j1 wi; j2 Þ
tð4Þ ¼ max½0; ðeð4Þ tð2Þ Þ; eð4Þ normally equals 321 , and Kiþ12; j ¼ 23 ki 1 1 ki and ðk ; j þ k Þ; K 1 ¼ ðk þ k Þ, k ¼ 1 þ iþ1 i; j i; jþ1 i; j i i; jþ2 kj 2 2 23 k kj ¼ 1 þ kji kj , and k is the spectral radius of Jacobian
ð28Þ
k¼1
DstðrÞ ¼
tð2Þ is the first order dissipative flux, while tð4Þ is the third order dissipative flux. The formulation of ð2Þ tð2Þ is ti; j ¼ eð2Þ maxðri; j ; riþ1; j Þ; eð2Þ is the coefficient of the first order dissipative flux coefficient and unity is the typical jp j 2pi; j þpi1; j j . The formulation of tð4Þ is value for it and ri; j ¼ piþ1; iþ1; j þ2pi; j þpi1; j
ð27Þ
k¼1
K X dk fk ðrÞ
ð29Þ
k¼1
where
fk ðxÞ ¼ sincx ðpxeðkÞ Þ; eðkÞ ¼ lglg 0:5 ; xk
fk ðrÞ ¼ sincr ðprpðkÞ Þ;
dk is the kth design variable, and its initial value pðkÞ ¼ lglg0:5 rk is zero in every optimization cycle. x and xk are dimensionless axial coordinates, while r and rk are dimensionless radial coordinates; the indexes cx and cr are used to control the perturbation hump. The stagger angle b of the cascade and the stacking line of the blade are shown in Figs. 2 and 3 respectively. Once the values of the design variables dk are calculated, the perturbations (DyðxÞ; DbðrÞ and DstðrÞ) of the passage, the stagger angles and the stacking lines are achieved.
matrices. The discretization of the time term can be written as below
4.2. Optimization algorithm
dw þ Pi;j ¼ 0 dt
The quasi-Newton method of BFGS is employed to calculate the step values of the design variables when the flow equation and the adjoint equation are fully converged. The step values dI are defined as DS; ðDS ¼ e dS ; where e controls the step size and is determined by the line search methodÞ. To prevent drastic changes of the flow passage, stagger angles of the blade sections and the stacking lines, their
ð25Þ
where 3 AT i; j ðwiþ1; j wi1; j Þ þ 2Dn 7 6 7 6 Pi; j ¼ J6 BTi; j ðwi; jþ1 wi; j1 Þ þ Q eT w þ 7 5 4 i; j i; j 2Dg diþ12; j di12; j þ di; jþ12 di; j12 2
Then, the second order Runge–Kutta time method is 8 0 wi; j ¼ wni; j > > > > > > 1 0 0 > > < wi; j ¼ wi; j þ DtPi; j ð26Þ > 2 1 0 1 1 > w ¼ w þ Dt P þ P > i; j i; j i; j i; j 2 > > > > > : nþ1 wi; j ¼ w2i; j where n and n+1 represent the current and the new time steps respectively.
Fig. 2
Stagger angle of the cascade.
762
L. Chen, J. Chen
maximum perturbations are respectively set to be 0.2 mm, 0.2 and 0.15 mm in every design cycle. 5. Results and discussions 5.1. Validation of the S2 surface calculation
Table 1 Injection types
Velocity ratio of local gas velocity
3.28 1.31 0.57 0.57 0.41 0.41
700 700 700 700 700 700
0.3 0.3 0.3 0.3 0.3 0.3
Rotor
1 3 4 6
2.05 0.82 0.61 0.61
700 700 700 700
0.3 0.3 0.3 0.3
Table 2
Stacking line of the blade.
Total temperature (K)
1 3 4 6 8 9
Fig. 5
Fig. 3
Mass flow rate (kg/s)
Stator
5.1.1. Mesh independence study of S2 surface calculation The validation of the S2 surface calculation is conducted through a single stage high pressure turbine with cooling air injection. The blades for S2 surface calculation is shown in Fig. 4. The injection types of cooling air for the stator and the rotor are (1, 3, 4, 6, 8 and 9) and (1, 3, 4 and 6) respectively, as shown in Fig. 1. The parameters of the cooling air are shown in Table 1. The mesh independence study is conducted using two types of mesh as shown in Fig. 5. The comparison of the aerodynamic performance calculated from the two types of mesh above is shown in Table 2. The maximum deviation of the mass flow rate M calculated by Mesh 1 and Mesh 2 is 0.63%, and it is 0.42% and 0.09% for pressure ratio p and adiabatic efficiency g respectively. The distribution of the dimensionless velocity V calculated via the mesh above is presented in Fig. 6. It is obvious that the distribution of the velocity V in the radial direction simulated from the mesh is almost the same.
Parameters of the cooling air.
High pressure turbine mesh for S2 surface calculation.
Mesh independence study for high pressure turbines.
Mesh
M (kg/s)
p
g
1 2
89.21 88.65
2.878 2.890
0.8879 0.8871
The high gradient of the velocity area is concentrated at the leading edge of the rotor due to the rise in the diverging level of the passage. There are slight changes in the velocity contour of the stator. The zone with the lowest value of the velocity V is at the mid-span of the rotor, where the minor variations of the velocity only happen in the center of the contour lines. According to the discussion above, there is little variation in the results calculated from Mesh 1 and Mesh 2. Thus, Mesh 1 is adopted during the next optimization design step. 5.1.2. Comparisons of S2 surface and 3D simulation
Fig. 4
High pressure turbine blades for S2 surface calculation.
Three-dimensional steady numerical simulation of the high pressure turbine with cooling is carried out by the solver NUMECA. Reynolds averaged Navier–Stokes equations are
Aerodynamic optimization design for high pressure turbines based on the adjoint approach
Fig. 6 Comparison of the distributions of the dimensionless velocity on the S2 surface.
Fig. 7
Cooling air injection for 3D simulation.
763
Fig. 8 Distribution of static pressure on S2 surface from original design.
Fig. 9 Distribution of total pressure from hub to shroud at inlet and outlet for S2 calculation and 3D simulation.
Table 3 Aerodynamic performance from S2 surface and 3D simulation. Case
M (kg/s)
p
g
S2 calculation 3D simulation (S–A) 3D simulation (SST k x)
89.21 88.76 88.77
2.878 2.883 2.879
0.8879 0.8875 0.8878
numerically solved using the Spalart–Allmaras (S–A) turbulence model and SST (Shear Stress Transport) k x model to check the 3D simulation accuracy. The space discretization is based on a cell-centered finite volume approach, and the fourth order Runge–Kutta method is used for time marching. The types of the cooling air injection for the 3D numerical simulation are shown in Fig. 7, which should be the same as those for the S2 surface calculation.
The comparison of the simulation results of S2 surface calculation and 3D simulation has been made to validate the simulation accuracy of the S2 surface calculation. The aerodynamic performance parameters of S2 surface calculation and 3D simulation are shown in Table 3. As can be seen from the table above, the maximum deviation of the mass flow rate M calculated by S2 surface, 3D simulation S–A turbulence model and SST k x turbulence model is 0.51%, and is 0.17% and 0.05% for pressure ratio p and adiabatic efficiency g respectively. There is little change in the results for S–A and SST k x turbulence models. The comparison analysis of detailed information of the flow field from S2 surface calculation and the 3D simulation using S–A turbulence model is shown below. The static pressure distribution on the meridional plane is presented in Fig. 8. It is obvious that the overall distributions of the static pressure calculated by the two methods are the same. The static pressure
764
L. Chen, J. Chen
decreases when the gas passes the passage of the stator and rotor due to the expansion work from the gas on the turbine blade. The distribution lines of the static pressure of the turbine are uniform in the radial direction, and it is a clear echelon in the stream wise direction with non-uniform contour lines between the blade rows, where the gas starts expanding more greatly due to the rise in the diverging level of the shroud. The distributions of total pressure at the inlet and the outlet with S2 surface calculation and 3D simulation are presented in Fig. 9. There is a little discrepancy in the total pressure at the hub and shroud section for S2 surface calculation and 3D simulation. Meanwhile, the total pressure enjoys the same values at the sections from 20% to 80% span, with the averaged deviation being 0.26%. In addition, the total pressure curves of the two simulation methods have the same variation tendency at the outlet. There is a minor change in the total pressure between 10% and 40% span and the averaged deviation is 1.1%, which means that the results of S2 surface calculation have a good agreement with those of 3D simulation. 5.2. Validation of the aerodynamic optimization system The stagger angles of all the blades, the passage geometry and the stacking line of the stator are changed in the optimization, while the stacking line of the rotor is frozen due to the
Fig. 10 Variations of objective function and convergence history of adjoint equations during optimization.
Fig. 11
Distribution of adjoint variable vectors on S2 surface.
Aerodynamic optimization design for high pressure turbines based on the adjoint approach
765
Fig. 14 Comparison of passage geometry between the original design and the optimized design.
Table 4 Aerodynamic optimization.
performance
before
and
after
Case
M (kg/s)
p
g
Original design Optimized design
88.76 89.10
2.883 2.872
0.8875 0.8931
Fig. 12 Distribution of dimensionless static pressure of turbine on S2 surface.
Fig. 15 Comparison of blades at hub section, mid-section and shroud section before and after optimization.
Fig. 13 surface.
Distribution of dimensionless total temperature on S2
structural strength. Five design variables are chosen for hub, shroud, staggering angles of each blade and stacking line of the stator respectively. The variations of the objective function
during the optimization and the convergence history of adjoint equations during its iteration are presented in Fig. 10. The distribution of adjoint variable vectors w1 ; w2 ; w3 ; w4 ; w5 on the S2 surface is given in Fig. 11. The partial gradient information of flow variables with respect to the design variables can be denoted by the adjoint flow field. Due to the cooling air injection, the gradients of the adjoint variable vector are large at the trailing edge of the stator and the rotor (especially the area near the shroud), which indicates that the change of the design variables brings more variations of flow variables. The dimensionless static pressure distributions on the S2 surface before and after optimization are shown in Fig. 12. It is apparent that the distribution of the static pressure of the high pressure turbine is uniform in the radial direction, and
766
Fig. 16
L. Chen, J. Chen
Comparison of the contours of the static pressure and the limiting streamlines on the blades before and after optimization.
it is a clear echelon in the axial direction with only nonuniform contour lines between the blade rows due to the diverging passage. There are obvious changes in the passage after optimization. The radial coordinates of the shroud decrease greatly, while those of the hub rise between the blade rows, meaning that the converging level of the passage becomes stronger. According to the contour lines below, the static pressure contour lines from the optimized design bend to the high pressure side near the hub area of the blade, meaning an increase in the velocity of the main gas flow near the end-wall zone. It can restrain the accumulation of the fluid with low kinetic energy, resulting in carrying more of the fluid to the main gas flow and decreasing the secondary flow loss. The dimensionless total temperature distributions on the S2 surface before and after optimization are shown in Fig. 13. It can be seen that the total temperature of the stator has the following distributions: the high temperature zone concentrates on the mid-section and the area of the high temperature decreases obviously compared to that of the high temperature zone at the inlet, which is caused by the cooing air injection of different types. The high temperature area at the mid-section of the rotor decreases as well, with the reduction of the temperature of the first stator. The total temperature from the original design shows the distribution at the inlet: the temperature is higher at the mid-section and lower at the end-wall due to the non-uniform temperature boundary condition, and there is a large gradient of temperature at the shroud corner. With the optimized high pressure turbine, the more internal energy of the gas flow are converted into kinetic energy which makes the high temperature area at the inlet decrease because of the
changes of the diverging type in the passage, and the large gradient of the temperature is eliminated at the shroud corner to improve the total temperature distribution and reduce the gradient of the temperature at the end-wall. The original passage geometry and the optimized one are shown in Fig. 14, whereas the comparison between the original blades and the optimized ones of different sections is shown in Fig. 15. Normally, the design of the high pressure turbine passage is the ‘‘converging-diverging’’ type which can effectively restrain the development of the passage vortex and reduce the secondary flow loss including the loss caused by the cooling air injection. According to Fig. 14, there are minor changes in the hub geometry, while there are great variations of the shroud after optimization. The radial coordinate of the shroud falls apparently, and the converging level at the stator position increases with the reduction of the diverging level at the rotor position, which constrains the development of the boundary layer of the end-wall and the suction surface of the blade. The stagger angles and the stacking line of the stator change obviously after optimization, while there are slight variations of the stagger angles for the rotor as shown in Fig. 15. For the stator, the decrease of the stagger angles can make the outlet flow angle become bigger, so as to reduce the expansion level of the gas flow in the passage and lower the blade load whilst matching the mass flow rate. It is beneficial to constrain the development of the boundary layer of the blade suction surface. Moreover, the stacking line of the stator bends to the pressure surface direction to form the positive bowed blade, and the static pressure along the radial direction is thus redistributed, making the pressure lower at the
Aerodynamic optimization design for high pressure turbines based on the adjoint approach
Fig. 17
767
Contours of the Mach number and the streamlines on different spans of the stator before and after optimization.
mid-section and higher near the end-wall. The new type of static pressure distribution could bring more secondary flow near the end-wall to the main gas flow and thus reduce the secondary flow loss. 5.3. Three dimension numerical simulation verification As shown in Table 3, the calculation results of 3D simulation using S–A turbulence model and SST k x turbulence model are almost the same. Thus, the S–A turbulence model is adopted to conduct the verification in this section. The aerodynamic performance parameters of the original design and the optimized one are given in Table 4. As can be seen from Table 4, the adiabatic efficiency g has increased from 0.8875 to 0.8931 and the variation of the pressure ratio p is 0.38%, while the change of the mass flow rate M is 0.38%. The contours of the static pressure and the limiting streamlines on the blade before and after optimization are presented in Fig. 16. For the original design, there is apparent separation on the suction surface of both stator and rotor near the endwall, and the secondary flow of the rotor is stronger than that of the stator. A small separation occurs at the inlet of the
pressure surface of the rotor, and the flow of the stator pressure surface near the shroud and hub is impacted due to the cooling air injection. Moreover, the scale of the passage vortex (the upper and the lower vortex) at the suction surface of the stator is obviously smaller than that of the rotor because of the converging type design of the passage for stator. Meanwhile, the cooling air injection at the end-wall has a large effect on the flow of the suction surface, which contributes to the increase of the intensity of the secondary flow loss. According to the contour of the optimized design, the variations of the passage and the stacking lines allow the main gas flow to have more kinetic energy near the end-wall and thus to constrain the development of the boundary layer and the secondary flow. The decrease of the stagger angles of the stator reduces the turning angles of the gas flow in the passage, resulting in a higher static pressure of the peak at the suction surface. For the rotor, the reduction of the diverging level of the passage restrains the growth of the boundary layer at the end-wall and reduces the intensity of the passage vortex as well as the mixing loss caused by the cooling air injection, so as to improve the flow on the pressure surface of the rotor. The contours of the Mach number and the streamlines on the 5% span section, 50% span section and the 95% span
768
L. Chen, J. Chen
Fig. 18
Contours of the Mach number and the streamlines on different spans of the rotor before and after optimization.
section of the stator before and after optimization are presented in Fig. 17. According to the original design, the flow condition at the 5% and the 95% span section is greatly influenced by the passage vortex and the cooling air injection from the upstream, leading edge and the end-wall within the blade row. There is apparent separation at the upstream, leading edge and the blade surface,triggered by the cooling air injection. The mixing of the cooling air injected from the upstream end-wall with the main gas flow causes large separation, while the cooling air injection from the trailing edge has minor impact on the flow. The flow condition at the 50% span is better than that of the end-wall sections with only separation bubble at the leading edge. However, the Mach number of the 50% span is up to 1.3 running through the whole passage and the supersonic area is large. According to the optimized design, the mixing of the cooling air injected from the upstream and the end-wall within the blade row with the main gas flow is obviously reduced due to the variations of the passage, stagger angles and the stacking line. In addition, the separation bubble of the leading edge caused by the cooling air injection diminishes at the 95% span, while there are slight changes at the leading edge for the other sections. The maximum Mach number decreases and the supersonic area reduces apparently, lowering the shock wave loss. The contours of the Mach number and the streamlines on the 5% span, the 50% span and the 95% span of the rotor before and after optimization are shown in Fig. 18. It can be seen from the original design that the high intensity passage vortex makes a poor flow at the 5% span and the 95% span, and the mixing of the cooling air injected from the upstream and the leading edge with the main gas flow results in the
separation at the upstream and the leading edge, which contributes to the increase of the secondary flow loss. Moreover, the cooling air injection at the trailing edge has almost no influence on the flow field. The Mach number of the 50% span is up to 1.2 running through the whole passage. According to the optimized design, the mixing of the cooling air injection from the upstream and the leading edge with the main gas flow is reduced due to the new diverging passage, and the separation bubble at the leading edge near the end-wall diminishes. The intensity of the passage vortex is decreased as well and the Mach number at different sections changes slightly compared to that in the original design because of the minor variations of the stagger angles of the rotor. 6. Conclusions This paper presents an aerodynamic optimization design for high pressure turbines including the cooling air injection using the continuous adjoint method based on S2 surface Euler equation with source terms. From the analysis above, the following conclusion is made: (1) The aerodynamic optimization design system using the continuous adjoint method based on the S 2 surface is practicable and effective to optimize the high pressure turbine with cooling air injection for the case with a large number of design variables. In the optimization system, the flow equation and the adjoint equation are solved independently, and the objective function sensitivities are then obtained by the mesh perturbations efficiently.
Aerodynamic optimization design for high pressure turbines based on the adjoint approach (2) The cooling air injections from the leading edge, upstream position between the blade rows and endwall within the blade row have substantial impact on the main gas flow, causing huge separation near the leading edge and the blade surface, while the cooling air injection from the trailing edge has almost no influence on the flow field. (3) The optimization of the stagger angles, the stacking line and the passage geometry restrains the growth of the boundary layer near the end-wall and reduces the intensity of the passage vortex as well as the mixing loss caused by the cooling air injection. The optimized design redistributes the static pressure on the blade surface and improves the state of the attack angles, constraining the development of the boundary layer at the suction surface and the secondary flow loss drops.
Acknowledgments The authors are grateful to the anonymous reviewers for their critical and constructive review of the manuscript. This research work was funded by the Aeronautical Science Foundation of China – China (No. 2010ZB51023). References 1. Pironneau O. On optimum shapes in stokes flow. J Fluid Mech 1973;59(2):117–28. 2. Jameson A. Aerodynamic design via control theory. J Sci Comput 1988;3(3):233–60. 3. Jameson A. Aerodynamic shape optimization using the adjoint method. Brussels: Lectures at the Von Karman Institute; 2003. 4. Giles MB, Duta MC. Algorithm developments for discrete adjoint methods. AIAA J 2003;41(2):198–205. 5. Nadarajah SK, Jameson A. A comparison of the continuous and discrete adjoint approach to automatic aerodynamic optimization. Reston: AIAA; 2000. Report No.: AIAA-2000-0667. 6. Giles MB, Pierce NA. An introduction to the adjoint approach to design. Flow Turbul Combust 2000;65(3–4):393–415. 7. Jameson A. Re-engineering the design process through computation. J Aircr 1999;36(1):36–50. 8. Jameson A. Optimum aerodynamic design using CFD and control theory. Reston: AIAA; 1995. Report No.: AIAA-1995-1729. 9. Reuther J, Jameson A, Alonso JJ, Remlinger MJ, Saunders D. Constrained multipoint aerodynamic shape optimization using an adjoint formulation and parallel computers, Part 1. J Aircr 1999;36(1):51–60.
769
10. Reuther J, Jameson A, Alonso JJ, Remlinger MJ, Saunders D. Constrained multipoint aerodynamic shape optimization using an adjoint formulation and parallel computers, Part 2. J Aircr 1999;36(1):61–74. 11. Li YC, Yang DL, Feng ZP. Inverse problem in aerodynamic shape design of turbomachinery blades. New York: ASME; 2006. Report No.: ASME GT2006-91135. 12. Li YC, Feng ZP. Aerodynamic design of turbine blades by using adjoint-based method and N-S equation. New York: ASME; 2007. Report No.: ASME GT2007-27734. 13. Papadimitriou DI, Giannakoglou KC. Compressor blade optimization using a continuous adjoint formulation. New York: ASME; 2006. Report No.: ASME GT2006-904669. 14. Wang DX, He L. Adjoint aerodynamic design optimization for blades in multi-stage turbomachines: Part I- Methodology and verification. New York: ASME; 2008. Report No.: ASME GT2008-50208. 15. Luo JQ, Liu F, Ivan MB. Secondary flow reduction by blade redesign and endwall contouring using an adjoint optimization method. New York: ASME; 2010. Report No.: ASME GT201022061. 16. Ji LC, Li WW, Tian Y, Yi WL, Chen J. Multi-stage turbomachinery blades optimization design using adjoint method and thin shear-layer N-S equation. New York: ASME; 2012. Report No.: ASME GT2012-6853. 17. Zhang CL, Feng ZP. Aerodynamic shape design optimization for turbomachinery cascade based on discrete adjoint method; 2011. Report No.: ASME GT2011-45805. 18. Wu CH. A general theory of three-dimensional flow in subsonic and supersonic turbomachines of axial- radial, and mixed-flow types. Washington, D.C.: National Aeronautics and Space Administration; 1952. Report No. NACA-TN-2604. 19. Wu HY, Yang SC, Liu F. Comparison of three geometric representations of airfoils for aerodynamic optimization. Reston: AIAA; 2003. Report No.: AIAA-2003-4095. 20. Du L, Ning FF. Numerical optimization of com pressor blade profile based on the control theory. J Aerosp Power 2009;24(3):615–25. Chen Lei received the B.E. degree from Dalian University of Technology and graduated with the Ph.D. degree from Beihang University in 2007 and 2014 respectively, and then work at the Heavy Engineering Research Association in New Zealand. His main research interests are computational fluid dynamics, modeling of turbine and optimization of turbine based on adjoint method. Chen Jiang is a professor and Ph.D. supervisor at School of Energy and Power Engineering, Beihang University, China. His current research interests are computational fluid dynamics, aerodynamics, modeling of axial-centrifugal compressor, fan and turbine.
Chinese Society of Aeronautics and Astronautics & Beihang University
Chinese Journal of Aeronautics [email protected] www.sciencedirect.com
Aerodynamic optimization design for high pressure turbines based on the adjoint approach a,*
Chen Lei a b
, Chen Jiang
b
Heavy Engineering Research Association, Auckland 2104, New Zealand School of Energy and Power Engineering, Beihang University, Beijing 100191, China
Received 24 April 2014; revised 5 December 2014; accepted 2 March 2015 Available online 22 April 2015
KEYWORDS Adjoint method; Aerodynamic design; High pressure turbine; Optimization design; Objective function
Abstract A first study on the continuous adjoint formulation for aerodynamic optimization design of high pressure turbines based on S2 surface governed by the Euler equations with source terms is presented. The objective function is defined as an integral function along the boundaries, and the adjoint equations and the boundary conditions are derived by introducing the adjoint variable vectors. The gradient expression of the objective function then includes only the terms related to physical shape variations. The numerical solution of the adjoint equation is conducted by a finitedifference method with the Jameson spatial scheme employing the first and the third order dissipative fluxes. A gradient-based aerodynamic optimization system is established by integrating the blade stagger angles, the stacking lines and the passage perturbation parameterization with the quasi-Newton method of Broyden–Fletcher–Goldfarb–Shanno (BFGS). The application of the continuous adjoint method is validated through a single stage high pressure turbine optimization case. The adiabatic efficiency increases from 0.8875 to 0.8931, whilst the mass flow rate and the pressure ratio remain almost unchanged. The optimization design is shown to reduce the passage vortex loss as well as the mixing loss due to the cooling air injection. ª 2015 The Authors. Production and hosting by Elsevier Ltd. on behalf of CSAA & BUAA. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).
1. Introduction With the increasing need for high performance gas turbines to reduce the emission and the engine weight, the emerging trend * Corresponding author. Tel.: +64 92624751. E-mail addresses: [email protected] (L. Chen), chenjiang27@ buaa.edu.cn (J. Chen). Peer review under responsibility of Editorial Committee of CJA.
Production and hosting by Elsevier
is to use mathematical optimization techniques as an integral part of the aerodynamic design toolkit. Stochastic and gradient-based methods are the normal mathematical optimization algorithms. The stochastic method searches for a global optimal solution by monitoring the magnitude of the objective function, while the gradient-based method searches for a local optimal solution by monitoring the sensitivity of the objective function to the changes in the design variables. Both of the stochastic and the gradient-based methods tend to consume considerably computational resources for the cases with a large number of design variables. The adjoint method is another gradient-based approach especially for the optimization with numerous design variables.
http://dx.doi.org/10.1016/j.cja.2015.04.022 1000-9361 ª 2015 The Authors. Production and hosting by Elsevier Ltd. on behalf of CSAA & BUAA. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).
758 In this method, the adjoint system is conducted in a similar way as that in the optimal control problems. By introducing the adjoint variable vectors, the gradients of the objective function with respect to the design variables are calculated indirectly by solving the adjoint equations. Therefore, the sensitivity analysis is almost independent of the numbers of the design variables, and solving two sets of flow equations in one design cycle is nearly the total computing cost. When the flow equations and the adjoint equations are fully converged, the final gradients of the objective function with respect to the design variables can be obtained efficiently. Pironneau1 was the first to use the adjoint method in fluid mechanics, and then Jameson2 applied adjonit method in the aeronautical field. Combining the continuous adjoint method with CFD technology, Jameson developed the optimization design method which was applied to the transonic wing-body combinations.3 Moreover, the discrete adjoint method was also developed.4 ‘‘Continuous’’ and ‘‘discrete’’ methods symbolize two alternative approaches to deriving adjoint equations. As there is no clear quantitative comparison between the two approaches,5,6 they seem to achieve the same optimal goals. Both of the methods have performed well in the optimization for airfoils,7 wings,8 wing-body configurations9 and business jets.10 The adjoint method has been utilized in the area of turbomachinery. Li et al.11,12 used the continuous adjoint method based on Navier–Stokes and Euler equation respectively to conduct the aerodynamic optimization design for turbine blades, and the optimization system was validated by several numerical cases. Papadimitriou and Giannakoglou13 developed the continuous adjoint formulation to improve the aerodynamic performance of a 3D peripheral compressor blade cascade. Wang and He14 first proposed the adjoint nonreflective mixing-plane treatment method, and carried out the aerodynamic blading shape optimization design in a multi-stage turbomachinery environment. Luo et al.15 used the adjoint optimization method to reduce the secondary flow loss of turbine blades by redesigning the blade. Ji et al.16 combined the continuous adjoint method with thin shearlayer Navier–Stokes equations to construct an efficient sensitivity analysis optimization system for multi-stage turbomachinery blades, and the adjoint optimization code was validated through two compressor blades design cases. Zhang and Feng17 used the automatic differentiation tool to develop a discrete adjoint solver, and the optimization system was validated via a turbine cascade under the viscous flow environment. For the aerodynamic design of gas turbines, S2 surface design plays a crucial role in the entire design system.18 Blade design is based on the simulation results of S2 surface through flow calculation. This paper presents the results of the first study on the adjoint method applied to the S2 surface through flow calculation including the cooling air effect. The adjoint method is combined with the Euler equation with the source term to develop an efficient sensitivity analysis model for turbine blades (stagger angles and stacking lines) and the passage in a specified objective function. The validation of the optimization system is carried out via the case of aerodynamic optimization of a high pressure turbine.
L. Chen, J. Chen 2. Flow equations and solution methods The steady Euler equations of the curvilinear coordinates system are utilized to predict the aerodynamic performance of the gas turbine on S2 surface. The effects of the viscous losses, the leakages and the cooling air on the flow are concerned as the source terms in the right part of the Eq. (1). @ rU @ e @ e e þ Fþ G¼Q ð1Þ @t J @n @g e are the e G where U is the conservative flow variable vector; F; e convective flux vectors in the curvilinear coordinate systems; Q is the source term of the flow equations. @ni ; J ¼ detðKÞ Kij ¼ @xj e ¼ r Fnz þ Gnr þ H 1 nu F J r e ¼ r Fgz þ Ggr þ H 1 gu ; G J r 2
rm r m vz þ fz
e ¼ h þ h1 Q J 3
6 7 6 7 6 7 6 2 7 h ¼ 6 r m vr þ fr þ qðw þ xrÞ 7 6 7 6 r m v þ f qvðw þ 2xrÞ 7 u 4 5 u
r m H0 þ x2 r2 qv 2
3 0 6 r f @p þ p @ r f 7 6 J z @f @f J z 7 6 7 6r 7 @p @ r þ p f f 7 h1 ¼ 6 r @f J r 7 @f 6J 6 7 6 fu @p 7 @ fu 4 J @f þ p @f 5 J 0 where terms nz ; nr ; nu ; fz ; fr ; fu represent the partial derivatives @n @n @n @f @f @f ; ; ; ; ; ; F, G and H are the convective flux vectors in @z @r @u @z @r @u
the cylindrical coordinate system; m is the mass flow rate of the cooling air; vz ; vr and vu are the velocity components of the cooling air, and the injection of the cooling air is classified into nine types as shown in Fig. 1; H0 is the total enthalpy of cooling air; (fz ; fr ; fu ) are for accounting of viscous losses effects,
Fig. 1
Injection types of the cooling air.
Aerodynamic optimization design for high pressure turbines based on the adjoint approach and are calculated by the next relations: ðfz ; fr ; fu Þ ¼ u2 þvr/2 þw2 ðu; v; wÞ, / ¼ qT dr ; r is the entropy generation calculated by dt the loss coefficient, and there are four loss models (blade profile loss model, secondary flow loss model, trailing edge loss model and the angle of attack loss model) adopted in the code; u; v; w are the velocity in the axial, the radial and the tangential direction respectively; p is the static pressure and x is the angular velocity of rotor; h1 is for calculating the impact of is for the initial accelerablade thickness on the flow field; dx dt tion of the rotor. j is the number of the injection (1 6 j 6 9). j = 1: cooling air injection at the leading edge; j = 2: cooling air injection upstream the throat section; j = 3: cooling air injection at the trailing edge; j = 4: cooling air injection at the hub between the blade rows; j = 5: cooling air is injected upstream the blade row at the hub and then injected downstream the blade row; j = 6: cooling air is injected between the blade rows at the shroud; j = 7: cooling air is injected upstream the blade row at the shroud and then injected downstream the blade row; j = 8: cooling air is injected within the blade row at the hub; j = 9: cooling air is injected within the blade row at the shroud. Three types of boundary conditions are used to complete the numerical solution of the flow equation. Span-wise distributions of the total pressure, total temperature and velocity directions are set at the inlet, and the static pressure of the root is set at the outlet. It is the no flow through the boundary condition for the hub and the shroud. The Euler equation is numerically solved by using the highorder accurate Godunov scheme, and the explicit numerical scheme is employed first, followed by the implicit scheme. The implicit operator solves the linear systems of equations by scalar three diagonal solvers. Finally, the flow field is time-marched to a steady state in a pseudo time via the Euler time-marching approach. 3. Adjoint equations and solution methods Due to its lower memory requirements and easy implementation compared to the discrete adjoint method, the continuous adjoint approach is employed here. The detailed derivation of the adjoint equations from the Euler equations with source term (cooling air injection included) in a partial differential form is shown as below. 3.1. Adjoint equations The objective function can be expressed as an integral along the boundaries in the following general form, Z I¼ MdS ð2Þ S
where M is a function of the flow variables and the design variables, then Z Z dI ¼ MI dUdS þ MII ðddSÞ ð3Þ S
S @M . @U
where MI ¼ The subscripts I and II are used to distinguish the contributions due to the variations of @U and @S respectively.
759
The Euler flow Eq. (1) can be written as 1 @U þ RðU; SÞ ¼ 0 J @t
ð4Þ
For the fully converged flow equation, the variation of RðU; SÞ is dR ¼
@ e @ e e ¼0 dF þ dG dQ @n @g
ð5Þ
Then yields, 8 dR ¼ dRI þ dRII > > > e < fII ¼ AdU þ d F fII d F ¼ df FI þ d F e f f fII > d þ d ¼ BdU þ d G ¼ d G G G I II > > : e f f f f d Q ¼ d Q þ d Q ¼ Q dU þ d Q I
II
I
ð6Þ
II
where A and B are the Jacobian matrices, @e F G Q fI ¼ @ e f ¼ @e A¼f ; B¼G ; : FI ¼ @U Q I @U @U Integrating the Eq. (5) over the whole computational domain and introducing the adjoint variable vectors w into both sides of the equations Z Z @ e @ e e dV ¼ 0 ð7Þ wT dRdV ¼ wT dF þ dG dQ @n @g D D Regrouping Eqs. (6) and (7) Z wT dRdV ¼ wT ðdRI þ dRII ÞdV D ZD @ f @ f e I þ dRII dV ¼ 0 ð8Þ d FI þ d GI d Q ¼ wT @n @g D
Z
Performing integration by parts Z T @w f @wT f T f f w ðd FI nn þ d GI ng ÞdS d FI þ d GI dV @n @g S D Z Z f dV þ e wT d Q ZdV ¼0 I
Z
D
e ¼ wT where Z
D @ fII dF @n
ð9Þ
@ fII d Q f , and nn and ng are þ @g dG I
components of normal direction. Rearranging, Z Z T @w f @wT f fI ng ÞdS f dV wT ðdf d FI þ d GI þ wT d Q FI nn þ d G I @n @g S D Z @ f @ f f dV ¼ 0 ð10Þ d FII þ d GII d Q þ wT II @n @g D Subtracting Eq. (10) from Eq. (3) does not change the gradient in Eq. (3) Z Z Z fI ng ÞdS dI ¼ MI dUdS þ MII ðddSÞ wT ðdf FI nn þ d G S S S Z T @w f @wT f f dV d FI þ d GI þ w T d Q þ I @n @g D Z @ f @ f f ð11Þ wT d FII þ d G II d QII dV @n @g D Introducing the Jacobian matrices A and B of Eq. (6) into Eq. (11) Z Z dI ¼ MI dUdS þ MII ðddSÞ S S Z Z Z T e e CdV w ðAdUnn þ BdUng ÞdS þ XdV ð12Þ S
D
D
760
L. Chen, J. Chen
where 8 e I dU e ¼ @wT AdU þ @wT BdU þ wT Q
The subsonic outlet boundary condition is e @M @F wT nn ¼ 0 @q @q
where q ¼ q; u; v; w, and the adjoint variables w1 ; w2 ; w3 ; w4 are calculated by Eq. (19), while w5 is extrapolated from the interior fluid domain. At the wall boundary, Vn ¼ 0, and only one equation is required. The wall boundary is ð13Þ
D
where 8 T T ee > e I dU
@g
Eliminating the dependence of dI on dU achieves the adjoint equation and the boundary condition equation respectively @wT @wT eI ¼ 0 Aþ B þ wT Q @n @g
ð14Þ
MI wT ðAnn þ Bng Þ ¼ 0
ð15Þ
To solve the adjoint Eq. (14), a conservative form is firstly adopted and a pseudo time derivative term is added to conduct the Runge–Kutta time marching scheme. In addition, the characteristics of the adjoint equation are opposite to those of the Euler flow equation. Therefore, the adjoint equation shares the same Jacobian matrixes as those of the Euler flow equation with an additional negative sign, as shown in Eq. (16). @ 1 @w @w e T w AT BT Q w¼0 ð16Þ @t J @n @g 3.2. Adjoint boundary conditions The physical adjoint boundary conditions are defined by Eq. (15). Moreover, the propagation of the adjoint characteristics is in the opposite way to the physical information of the Euler equations. Thus, one physical boundary condition is required at the subsonic inlet and four are required at subsonic exit, and the other numerical boundary conditions (four conditions at the inlet and one condition at the outlet) are extrapolated from the interior fluid domains. The velocity vector k of the adjoint equation is also opposite to that of the flow equation ð17Þ
where Vn is the normal component of velocity, and a is the local speed of sound. Normally, it is assumed that the geometry of the inlet and the outlet is frozen during each optimization cycle; therefore, the subsonic inlet boundary condition is e @M @F wT nn ¼ 0 @p @p
w2 ng1 þ w3 ng2 þ w4 ng3 ¼
@M @p
ð20Þ
3.3. Objective function
II
k ¼ ½Vn ; Vn ; Vn ; ðVn þ aÞ; ðVn aÞT
ð19Þ
ð18Þ
and the inlet plane is assumed to be normal to the X axis. The adjoint variables w1 ; w2 ; w3 ; w4 are extrapolated from the interior fluid domain, while w5 is calculated through Eq. (18).
Entropy generation is defined as the objective function, which is a weighted sum of mass flow rate and it can be formulated as R 2 udA q~ s I ¼ þ # in 1 ð21Þ s0 m0 where s is the entropy generation given by Z Z p p 0 s¼ dA dA0 c c q q out in
ð22Þ
and s0 is the original entropy generation; m0 is the original n mass flow rate of the gas; ue ¼ unz þ vnr þ w ru ; # is a weight factor chosen to be 60 in the case. Then, the final expression of the objective function is Z @ f @ f f dV dI ¼ wT d FII þ d GII d Q II @n @g D Z Z @ fII dV þ wT d Q f dV fII þ @ d G dF ¼ wT II @n @g D D 0R 1 R Z fII nn dS @wT d F fII dVþ wT d F S D @n f dV A þ ¼ @ R wT d Q R @wT II T f f D w d GII ng dS D @g d GII dV S Z Z T @w f @wT f fII ng dS þ f dV d FII þ d GII þ wT d Q ¼ wT d G II @n @g D ZS gu g g d z w2 þ d r w3 þ d w4 pdS ¼ J J Jr S Z T nj @w @wT gj f dV d d Gj þ wT d Q Fj þ ð23Þ þ II @n J @g J D f is related to the variations of blade and passage where d Q II geometry, and its formulation is f ¼ ½d Q f1 ; d Q f2 ; d Q f3 ; d Q f4 ; d Q f5 dQ II and 8 1 f1 ¼ m dQ dJ > > > > > > f2 ¼ m vz d 1 þ fz d 1 @p d fz p @ d fz > dQ > r J J J J @f @f > > > 1 fr 1 1 > > 2 f3 ¼ m vr d þ d þ qw d þ > dQ > r J J Jr > < f f @ r 2qwx d 1J þ qx2 d Jr @p d Jr d p @f J @f > > > f > f > d Q4 ¼ m vu d 1J þ ru d 1J qvw d Jr1 > > > > > fu fu @p > 1 @ > > 2xqv d J @f d Jr p @f d Jr > > > : f d Q5 ¼ m H0 d 1J þ x2 qv d 1J
Aerodynamic optimization design for high pressure turbines based on the adjoint approach
761
3.4. Numerical solution
4. Optimization design
The finite difference method is employed. The time derivative term is time-marched by the second order Runge–Kutta time method, and the space derivative term is discretized by the Jameson’s spatial scheme employing the first and the third order artificial dissipative flux. The discretized equation is
Hicks-Henne hump functions19 are used as the parameterization method. In addition, the optimization algorithm adopts the quasi-Newton method of Broyden–Fletcher–Goldfarb–Sh anno (BFGS)20 together with a line search method.
T
4.1. Geometry parameterization
T
1 ðwnþ1 wn Þ Ai;j ðwiþ1;j wi1;j Þ Bi;j ðwi;jþ1 wi;j1 Þ ¼ þ J Dt 2Dn 2Dg T e e þQ w þd i;j
i;j
ð24Þ
where e d ¼ diþ12; j di12; j þ di; jþ12 di; j12 and diþ12; j ¼ tð2Þ Kiþ12;j ðwiþ1; j wi;j Þ tð4Þ Kiþ12;j ðwiþ2;j 3wiþ1; j þ 3wi; j wi1;j Þ
Three Hicks-Henne functions Eqs. (27)–(29) are used for the parameterization: one for passage perturbations, another one for blade stagger angle perturbations of different sections from hub to shroud, and the last one for the stacking line perturbations. K X dk fk ðxÞ
DyðxÞ ¼ di12;j ¼ tð2Þ Ki12;j ðwi;j wi1;j Þ tð4Þ Ki12;j ðwiþ1;j 3wi;j þ 3wi1;j wi2; j Þ
K X dk fk ðrÞ
DbðrÞ ¼ di; jþ12 ¼ tð2Þ Ki;jþ12 ðwi;jþ1 wi;j Þ tð4Þ Ki;jþ12 ðwi; jþ2 3wi;jþ1 þ 3wi;j wi; j1 Þ di; j12 ¼ tð2Þ Ki;j12 ðwi;j wi;j1 Þ tð4Þ Ki;j12 ðwi; jþ1 3wi;j þ 3wi;j1 wi; j2 Þ
tð4Þ ¼ max½0; ðeð4Þ tð2Þ Þ; eð4Þ normally equals 321 , and Kiþ12; j ¼ 23 ki 1 1 ki and ðk ; j þ k Þ; K 1 ¼ ðk þ k Þ, k ¼ 1 þ iþ1 i; j i; jþ1 i; j i i; jþ2 kj 2 2 23 k kj ¼ 1 þ kji kj , and k is the spectral radius of Jacobian
ð28Þ
k¼1
DstðrÞ ¼
tð2Þ is the first order dissipative flux, while tð4Þ is the third order dissipative flux. The formulation of ð2Þ tð2Þ is ti; j ¼ eð2Þ maxðri; j ; riþ1; j Þ; eð2Þ is the coefficient of the first order dissipative flux coefficient and unity is the typical jp j 2pi; j þpi1; j j . The formulation of tð4Þ is value for it and ri; j ¼ piþ1; iþ1; j þ2pi; j þpi1; j
ð27Þ
k¼1
K X dk fk ðrÞ
ð29Þ
k¼1
where
fk ðxÞ ¼ sincx ðpxeðkÞ Þ; eðkÞ ¼ lglg 0:5 ; xk
fk ðrÞ ¼ sincr ðprpðkÞ Þ;
dk is the kth design variable, and its initial value pðkÞ ¼ lglg0:5 rk is zero in every optimization cycle. x and xk are dimensionless axial coordinates, while r and rk are dimensionless radial coordinates; the indexes cx and cr are used to control the perturbation hump. The stagger angle b of the cascade and the stacking line of the blade are shown in Figs. 2 and 3 respectively. Once the values of the design variables dk are calculated, the perturbations (DyðxÞ; DbðrÞ and DstðrÞ) of the passage, the stagger angles and the stacking lines are achieved.
matrices. The discretization of the time term can be written as below
4.2. Optimization algorithm
dw þ Pi;j ¼ 0 dt
The quasi-Newton method of BFGS is employed to calculate the step values of the design variables when the flow equation and the adjoint equation are fully converged. The step values dI are defined as DS; ðDS ¼ e dS ; where e controls the step size and is determined by the line search methodÞ. To prevent drastic changes of the flow passage, stagger angles of the blade sections and the stacking lines, their
ð25Þ
where 3 AT i; j ðwiþ1; j wi1; j Þ þ 2Dn 7 6 7 6 Pi; j ¼ J6 BTi; j ðwi; jþ1 wi; j1 Þ þ Q eT w þ 7 5 4 i; j i; j 2Dg diþ12; j di12; j þ di; jþ12 di; j12 2
Then, the second order Runge–Kutta time method is 8 0 wi; j ¼ wni; j > > > > > > 1 0 0 > > < wi; j ¼ wi; j þ DtPi; j ð26Þ > 2 1 0 1 1 > w ¼ w þ Dt P þ P > i; j i; j i; j i; j 2 > > > > > : nþ1 wi; j ¼ w2i; j where n and n+1 represent the current and the new time steps respectively.
Fig. 2
Stagger angle of the cascade.
762
L. Chen, J. Chen
maximum perturbations are respectively set to be 0.2 mm, 0.2 and 0.15 mm in every design cycle. 5. Results and discussions 5.1. Validation of the S2 surface calculation
Table 1 Injection types
Velocity ratio of local gas velocity
3.28 1.31 0.57 0.57 0.41 0.41
700 700 700 700 700 700
0.3 0.3 0.3 0.3 0.3 0.3
Rotor
1 3 4 6
2.05 0.82 0.61 0.61
700 700 700 700
0.3 0.3 0.3 0.3
Table 2
Stacking line of the blade.
Total temperature (K)
1 3 4 6 8 9
Fig. 5
Fig. 3
Mass flow rate (kg/s)
Stator
5.1.1. Mesh independence study of S2 surface calculation The validation of the S2 surface calculation is conducted through a single stage high pressure turbine with cooling air injection. The blades for S2 surface calculation is shown in Fig. 4. The injection types of cooling air for the stator and the rotor are (1, 3, 4, 6, 8 and 9) and (1, 3, 4 and 6) respectively, as shown in Fig. 1. The parameters of the cooling air are shown in Table 1. The mesh independence study is conducted using two types of mesh as shown in Fig. 5. The comparison of the aerodynamic performance calculated from the two types of mesh above is shown in Table 2. The maximum deviation of the mass flow rate M calculated by Mesh 1 and Mesh 2 is 0.63%, and it is 0.42% and 0.09% for pressure ratio p and adiabatic efficiency g respectively. The distribution of the dimensionless velocity V calculated via the mesh above is presented in Fig. 6. It is obvious that the distribution of the velocity V in the radial direction simulated from the mesh is almost the same.
Parameters of the cooling air.
High pressure turbine mesh for S2 surface calculation.
Mesh independence study for high pressure turbines.
Mesh
M (kg/s)
p
g
1 2
89.21 88.65
2.878 2.890
0.8879 0.8871
The high gradient of the velocity area is concentrated at the leading edge of the rotor due to the rise in the diverging level of the passage. There are slight changes in the velocity contour of the stator. The zone with the lowest value of the velocity V is at the mid-span of the rotor, where the minor variations of the velocity only happen in the center of the contour lines. According to the discussion above, there is little variation in the results calculated from Mesh 1 and Mesh 2. Thus, Mesh 1 is adopted during the next optimization design step. 5.1.2. Comparisons of S2 surface and 3D simulation
Fig. 4
High pressure turbine blades for S2 surface calculation.
Three-dimensional steady numerical simulation of the high pressure turbine with cooling is carried out by the solver NUMECA. Reynolds averaged Navier–Stokes equations are
Aerodynamic optimization design for high pressure turbines based on the adjoint approach
Fig. 6 Comparison of the distributions of the dimensionless velocity on the S2 surface.
Fig. 7
Cooling air injection for 3D simulation.
763
Fig. 8 Distribution of static pressure on S2 surface from original design.
Fig. 9 Distribution of total pressure from hub to shroud at inlet and outlet for S2 calculation and 3D simulation.
Table 3 Aerodynamic performance from S2 surface and 3D simulation. Case
M (kg/s)
p
g
S2 calculation 3D simulation (S–A) 3D simulation (SST k x)
89.21 88.76 88.77
2.878 2.883 2.879
0.8879 0.8875 0.8878
numerically solved using the Spalart–Allmaras (S–A) turbulence model and SST (Shear Stress Transport) k x model to check the 3D simulation accuracy. The space discretization is based on a cell-centered finite volume approach, and the fourth order Runge–Kutta method is used for time marching. The types of the cooling air injection for the 3D numerical simulation are shown in Fig. 7, which should be the same as those for the S2 surface calculation.
The comparison of the simulation results of S2 surface calculation and 3D simulation has been made to validate the simulation accuracy of the S2 surface calculation. The aerodynamic performance parameters of S2 surface calculation and 3D simulation are shown in Table 3. As can be seen from the table above, the maximum deviation of the mass flow rate M calculated by S2 surface, 3D simulation S–A turbulence model and SST k x turbulence model is 0.51%, and is 0.17% and 0.05% for pressure ratio p and adiabatic efficiency g respectively. There is little change in the results for S–A and SST k x turbulence models. The comparison analysis of detailed information of the flow field from S2 surface calculation and the 3D simulation using S–A turbulence model is shown below. The static pressure distribution on the meridional plane is presented in Fig. 8. It is obvious that the overall distributions of the static pressure calculated by the two methods are the same. The static pressure
764
L. Chen, J. Chen
decreases when the gas passes the passage of the stator and rotor due to the expansion work from the gas on the turbine blade. The distribution lines of the static pressure of the turbine are uniform in the radial direction, and it is a clear echelon in the stream wise direction with non-uniform contour lines between the blade rows, where the gas starts expanding more greatly due to the rise in the diverging level of the shroud. The distributions of total pressure at the inlet and the outlet with S2 surface calculation and 3D simulation are presented in Fig. 9. There is a little discrepancy in the total pressure at the hub and shroud section for S2 surface calculation and 3D simulation. Meanwhile, the total pressure enjoys the same values at the sections from 20% to 80% span, with the averaged deviation being 0.26%. In addition, the total pressure curves of the two simulation methods have the same variation tendency at the outlet. There is a minor change in the total pressure between 10% and 40% span and the averaged deviation is 1.1%, which means that the results of S2 surface calculation have a good agreement with those of 3D simulation. 5.2. Validation of the aerodynamic optimization system The stagger angles of all the blades, the passage geometry and the stacking line of the stator are changed in the optimization, while the stacking line of the rotor is frozen due to the
Fig. 10 Variations of objective function and convergence history of adjoint equations during optimization.
Fig. 11
Distribution of adjoint variable vectors on S2 surface.
Aerodynamic optimization design for high pressure turbines based on the adjoint approach
765
Fig. 14 Comparison of passage geometry between the original design and the optimized design.
Table 4 Aerodynamic optimization.
performance
before
and
after
Case
M (kg/s)
p
g
Original design Optimized design
88.76 89.10
2.883 2.872
0.8875 0.8931
Fig. 12 Distribution of dimensionless static pressure of turbine on S2 surface.
Fig. 15 Comparison of blades at hub section, mid-section and shroud section before and after optimization.
Fig. 13 surface.
Distribution of dimensionless total temperature on S2
structural strength. Five design variables are chosen for hub, shroud, staggering angles of each blade and stacking line of the stator respectively. The variations of the objective function
during the optimization and the convergence history of adjoint equations during its iteration are presented in Fig. 10. The distribution of adjoint variable vectors w1 ; w2 ; w3 ; w4 ; w5 on the S2 surface is given in Fig. 11. The partial gradient information of flow variables with respect to the design variables can be denoted by the adjoint flow field. Due to the cooling air injection, the gradients of the adjoint variable vector are large at the trailing edge of the stator and the rotor (especially the area near the shroud), which indicates that the change of the design variables brings more variations of flow variables. The dimensionless static pressure distributions on the S2 surface before and after optimization are shown in Fig. 12. It is apparent that the distribution of the static pressure of the high pressure turbine is uniform in the radial direction, and
766
Fig. 16
L. Chen, J. Chen
Comparison of the contours of the static pressure and the limiting streamlines on the blades before and after optimization.
it is a clear echelon in the axial direction with only nonuniform contour lines between the blade rows due to the diverging passage. There are obvious changes in the passage after optimization. The radial coordinates of the shroud decrease greatly, while those of the hub rise between the blade rows, meaning that the converging level of the passage becomes stronger. According to the contour lines below, the static pressure contour lines from the optimized design bend to the high pressure side near the hub area of the blade, meaning an increase in the velocity of the main gas flow near the end-wall zone. It can restrain the accumulation of the fluid with low kinetic energy, resulting in carrying more of the fluid to the main gas flow and decreasing the secondary flow loss. The dimensionless total temperature distributions on the S2 surface before and after optimization are shown in Fig. 13. It can be seen that the total temperature of the stator has the following distributions: the high temperature zone concentrates on the mid-section and the area of the high temperature decreases obviously compared to that of the high temperature zone at the inlet, which is caused by the cooing air injection of different types. The high temperature area at the mid-section of the rotor decreases as well, with the reduction of the temperature of the first stator. The total temperature from the original design shows the distribution at the inlet: the temperature is higher at the mid-section and lower at the end-wall due to the non-uniform temperature boundary condition, and there is a large gradient of temperature at the shroud corner. With the optimized high pressure turbine, the more internal energy of the gas flow are converted into kinetic energy which makes the high temperature area at the inlet decrease because of the
changes of the diverging type in the passage, and the large gradient of the temperature is eliminated at the shroud corner to improve the total temperature distribution and reduce the gradient of the temperature at the end-wall. The original passage geometry and the optimized one are shown in Fig. 14, whereas the comparison between the original blades and the optimized ones of different sections is shown in Fig. 15. Normally, the design of the high pressure turbine passage is the ‘‘converging-diverging’’ type which can effectively restrain the development of the passage vortex and reduce the secondary flow loss including the loss caused by the cooling air injection. According to Fig. 14, there are minor changes in the hub geometry, while there are great variations of the shroud after optimization. The radial coordinate of the shroud falls apparently, and the converging level at the stator position increases with the reduction of the diverging level at the rotor position, which constrains the development of the boundary layer of the end-wall and the suction surface of the blade. The stagger angles and the stacking line of the stator change obviously after optimization, while there are slight variations of the stagger angles for the rotor as shown in Fig. 15. For the stator, the decrease of the stagger angles can make the outlet flow angle become bigger, so as to reduce the expansion level of the gas flow in the passage and lower the blade load whilst matching the mass flow rate. It is beneficial to constrain the development of the boundary layer of the blade suction surface. Moreover, the stacking line of the stator bends to the pressure surface direction to form the positive bowed blade, and the static pressure along the radial direction is thus redistributed, making the pressure lower at the
Aerodynamic optimization design for high pressure turbines based on the adjoint approach
Fig. 17
767
Contours of the Mach number and the streamlines on different spans of the stator before and after optimization.
mid-section and higher near the end-wall. The new type of static pressure distribution could bring more secondary flow near the end-wall to the main gas flow and thus reduce the secondary flow loss. 5.3. Three dimension numerical simulation verification As shown in Table 3, the calculation results of 3D simulation using S–A turbulence model and SST k x turbulence model are almost the same. Thus, the S–A turbulence model is adopted to conduct the verification in this section. The aerodynamic performance parameters of the original design and the optimized one are given in Table 4. As can be seen from Table 4, the adiabatic efficiency g has increased from 0.8875 to 0.8931 and the variation of the pressure ratio p is 0.38%, while the change of the mass flow rate M is 0.38%. The contours of the static pressure and the limiting streamlines on the blade before and after optimization are presented in Fig. 16. For the original design, there is apparent separation on the suction surface of both stator and rotor near the endwall, and the secondary flow of the rotor is stronger than that of the stator. A small separation occurs at the inlet of the
pressure surface of the rotor, and the flow of the stator pressure surface near the shroud and hub is impacted due to the cooling air injection. Moreover, the scale of the passage vortex (the upper and the lower vortex) at the suction surface of the stator is obviously smaller than that of the rotor because of the converging type design of the passage for stator. Meanwhile, the cooling air injection at the end-wall has a large effect on the flow of the suction surface, which contributes to the increase of the intensity of the secondary flow loss. According to the contour of the optimized design, the variations of the passage and the stacking lines allow the main gas flow to have more kinetic energy near the end-wall and thus to constrain the development of the boundary layer and the secondary flow. The decrease of the stagger angles of the stator reduces the turning angles of the gas flow in the passage, resulting in a higher static pressure of the peak at the suction surface. For the rotor, the reduction of the diverging level of the passage restrains the growth of the boundary layer at the end-wall and reduces the intensity of the passage vortex as well as the mixing loss caused by the cooling air injection, so as to improve the flow on the pressure surface of the rotor. The contours of the Mach number and the streamlines on the 5% span section, 50% span section and the 95% span
768
L. Chen, J. Chen
Fig. 18
Contours of the Mach number and the streamlines on different spans of the rotor before and after optimization.
section of the stator before and after optimization are presented in Fig. 17. According to the original design, the flow condition at the 5% and the 95% span section is greatly influenced by the passage vortex and the cooling air injection from the upstream, leading edge and the end-wall within the blade row. There is apparent separation at the upstream, leading edge and the blade surface,triggered by the cooling air injection. The mixing of the cooling air injected from the upstream end-wall with the main gas flow causes large separation, while the cooling air injection from the trailing edge has minor impact on the flow. The flow condition at the 50% span is better than that of the end-wall sections with only separation bubble at the leading edge. However, the Mach number of the 50% span is up to 1.3 running through the whole passage and the supersonic area is large. According to the optimized design, the mixing of the cooling air injected from the upstream and the end-wall within the blade row with the main gas flow is obviously reduced due to the variations of the passage, stagger angles and the stacking line. In addition, the separation bubble of the leading edge caused by the cooling air injection diminishes at the 95% span, while there are slight changes at the leading edge for the other sections. The maximum Mach number decreases and the supersonic area reduces apparently, lowering the shock wave loss. The contours of the Mach number and the streamlines on the 5% span, the 50% span and the 95% span of the rotor before and after optimization are shown in Fig. 18. It can be seen from the original design that the high intensity passage vortex makes a poor flow at the 5% span and the 95% span, and the mixing of the cooling air injected from the upstream and the leading edge with the main gas flow results in the
separation at the upstream and the leading edge, which contributes to the increase of the secondary flow loss. Moreover, the cooling air injection at the trailing edge has almost no influence on the flow field. The Mach number of the 50% span is up to 1.2 running through the whole passage. According to the optimized design, the mixing of the cooling air injection from the upstream and the leading edge with the main gas flow is reduced due to the new diverging passage, and the separation bubble at the leading edge near the end-wall diminishes. The intensity of the passage vortex is decreased as well and the Mach number at different sections changes slightly compared to that in the original design because of the minor variations of the stagger angles of the rotor. 6. Conclusions This paper presents an aerodynamic optimization design for high pressure turbines including the cooling air injection using the continuous adjoint method based on S2 surface Euler equation with source terms. From the analysis above, the following conclusion is made: (1) The aerodynamic optimization design system using the continuous adjoint method based on the S 2 surface is practicable and effective to optimize the high pressure turbine with cooling air injection for the case with a large number of design variables. In the optimization system, the flow equation and the adjoint equation are solved independently, and the objective function sensitivities are then obtained by the mesh perturbations efficiently.
Aerodynamic optimization design for high pressure turbines based on the adjoint approach (2) The cooling air injections from the leading edge, upstream position between the blade rows and endwall within the blade row have substantial impact on the main gas flow, causing huge separation near the leading edge and the blade surface, while the cooling air injection from the trailing edge has almost no influence on the flow field. (3) The optimization of the stagger angles, the stacking line and the passage geometry restrains the growth of the boundary layer near the end-wall and reduces the intensity of the passage vortex as well as the mixing loss caused by the cooling air injection. The optimized design redistributes the static pressure on the blade surface and improves the state of the attack angles, constraining the development of the boundary layer at the suction surface and the secondary flow loss drops.
Acknowledgments The authors are grateful to the anonymous reviewers for their critical and constructive review of the manuscript. This research work was funded by the Aeronautical Science Foundation of China – China (No. 2010ZB51023). References 1. Pironneau O. On optimum shapes in stokes flow. J Fluid Mech 1973;59(2):117–28. 2. Jameson A. Aerodynamic design via control theory. J Sci Comput 1988;3(3):233–60. 3. Jameson A. Aerodynamic shape optimization using the adjoint method. Brussels: Lectures at the Von Karman Institute; 2003. 4. Giles MB, Duta MC. Algorithm developments for discrete adjoint methods. AIAA J 2003;41(2):198–205. 5. Nadarajah SK, Jameson A. A comparison of the continuous and discrete adjoint approach to automatic aerodynamic optimization. Reston: AIAA; 2000. Report No.: AIAA-2000-0667. 6. Giles MB, Pierce NA. An introduction to the adjoint approach to design. Flow Turbul Combust 2000;65(3–4):393–415. 7. Jameson A. Re-engineering the design process through computation. J Aircr 1999;36(1):36–50. 8. Jameson A. Optimum aerodynamic design using CFD and control theory. Reston: AIAA; 1995. Report No.: AIAA-1995-1729. 9. Reuther J, Jameson A, Alonso JJ, Remlinger MJ, Saunders D. Constrained multipoint aerodynamic shape optimization using an adjoint formulation and parallel computers, Part 1. J Aircr 1999;36(1):51–60.
769
10. Reuther J, Jameson A, Alonso JJ, Remlinger MJ, Saunders D. Constrained multipoint aerodynamic shape optimization using an adjoint formulation and parallel computers, Part 2. J Aircr 1999;36(1):61–74. 11. Li YC, Yang DL, Feng ZP. Inverse problem in aerodynamic shape design of turbomachinery blades. New York: ASME; 2006. Report No.: ASME GT2006-91135. 12. Li YC, Feng ZP. Aerodynamic design of turbine blades by using adjoint-based method and N-S equation. New York: ASME; 2007. Report No.: ASME GT2007-27734. 13. Papadimitriou DI, Giannakoglou KC. Compressor blade optimization using a continuous adjoint formulation. New York: ASME; 2006. Report No.: ASME GT2006-904669. 14. Wang DX, He L. Adjoint aerodynamic design optimization for blades in multi-stage turbomachines: Part I- Methodology and verification. New York: ASME; 2008. Report No.: ASME GT2008-50208. 15. Luo JQ, Liu F, Ivan MB. Secondary flow reduction by blade redesign and endwall contouring using an adjoint optimization method. New York: ASME; 2010. Report No.: ASME GT201022061. 16. Ji LC, Li WW, Tian Y, Yi WL, Chen J. Multi-stage turbomachinery blades optimization design using adjoint method and thin shear-layer N-S equation. New York: ASME; 2012. Report No.: ASME GT2012-6853. 17. Zhang CL, Feng ZP. Aerodynamic shape design optimization for turbomachinery cascade based on discrete adjoint method; 2011. Report No.: ASME GT2011-45805. 18. Wu CH. A general theory of three-dimensional flow in subsonic and supersonic turbomachines of axial- radial, and mixed-flow types. Washington, D.C.: National Aeronautics and Space Administration; 1952. Report No. NACA-TN-2604. 19. Wu HY, Yang SC, Liu F. Comparison of three geometric representations of airfoils for aerodynamic optimization. Reston: AIAA; 2003. Report No.: AIAA-2003-4095. 20. Du L, Ning FF. Numerical optimization of com pressor blade profile based on the control theory. J Aerosp Power 2009;24(3):615–25. Chen Lei received the B.E. degree from Dalian University of Technology and graduated with the Ph.D. degree from Beihang University in 2007 and 2014 respectively, and then work at the Heavy Engineering Research Association in New Zealand. His main research interests are computational fluid dynamics, modeling of turbine and optimization of turbine based on adjoint method. Chen Jiang is a professor and Ph.D. supervisor at School of Energy and Power Engineering, Beihang University, China. His current research interests are computational fluid dynamics, aerodynamics, modeling of axial-centrifugal compressor, fan and turbine.