# S2 Stream Surface Approximation Approach for Quasi Three-Dimensional Turbine Time-Marching Design Method

Authors: Zhuyu Jiang; Zhaolin Fan; Ming Qiu; Cheng Gong
Journal: International Journal of Aerospace Engineering
 (2023)

---

## Abstract

A fast and reliable turbine design method is necessary in aeroengine development practice. As Q3D time-marching method iteration process can be time consuming and unstable, a novel S2 stream surface approximation approach coupled with flow tangential condition is proposed to reduce the computation requirement. The surface approximation uses a quadratic function of axial coordinatex to formulate distribution of circumferential coordinate θ. The flow tangential condition inherently represents inviscid blade force effect, and the formulation of inviscid blade force is not needed. A corresponding personal Q3D turbine computer code was developed, which was suitable for both design and analysis applications. The Q3D design method solves Euler’s equation through third-order Godunov’s scheme with TVD property in finite volume method. Semi-implicit Crank-Nicolson’s temporal scheme is implemented. Profile, secondary, and tip clearance energy loss models are added to predict viscous losses. Blade row exit swirl is required as the convergence target for Q3D design. A two-stage turbine and a single-stage high-pressure turbine were designed and analyzed by the Q3D method. Viscous 3D CFD was utilized to check design performances. The results showed that the Q3D method could finish a design case within 1 min. At design point, Q3D mass flow rate error was no more than 1.25%, expansion ratio error was no less than -0.022, and isentropic efficiency error was no more than 0.37 percentage points. The Q3D design method is fast and accurate. The stream surface approximation approach is suitable for Q3D design and analysis.


---

## Body

## 1. Introduction

Turbine is one of the key components in gas turbine aeroengine. In aeroengine development practice, a fast and reliable turbine aerodynamic design method is necessary. At present, the turbine aerodynamic design process still follows the system of one-dimensional (1D), two-dimensional (2D) or quasi three-dimensional (Q3D), and fully three-dimensional (3D) computation, of which Q3D design and analysis work is an important phase. The Q3D design theory was first proposed by Wu in 1952. As Wu suggested, flow patterns in turbomachinery could be divided into two types, of which one was the flow along S1 stream surface and the other along S2 stream surface [1]. Wu’s idea inspired many researchers to develop numerical methods for Q3D or so-called through flow computation on S2 stream surface. The main purpose of Q3D design work is to obtain radial distribution of aerodynamic parameters, from which blade shape profiling are supported. The Q3D design method is still the core technology in turbine design process. Two classical calculation methods are streamline curvature (SLC) method proposed by Novak and matrix stream function (MSF) method developed by Marsh [2, 3]. These two methods, especially SLC method, were prevalently used and improved for many years. However, their biggest drawback is that they are not suitable for transonic and supersonic flows, and chocking mass flow rate cannot be predicted [4]. To cope with this problem, the Q3D time-marching method based on computational fluid dynamics (CFD), which solved the Euler or Navier–Stokes (N–S) equations, was established. Through time-marching method, shock wave and chocked flow can be predicted automatically. Spurr [5] firstly developed Q3D time-marching method by solving the 2D Euler equation in stagnant cylindrical coordinate, and the result showed the potential of shock wave prediction. Numbers of researchers were devoted to the development of Q3D time-marching method [6, 7]. In general, the applications of Q3D time-marching method consist of two parts. They are design and analysis applications.For analysis application, the Q3D time-marching computation procedure is in accordance with 3D CFD technique. At first, turbomachinery configuration is needed, from which computational domain and grid is generated. In most cases, blade camber surface shapes and flow path geometry are required [8, 9]. Usually, it is thought that fluid moving direction is closely tangential to camber surface. However, because of incidence and deviation angles at blade leading and trailing edges, the real S2 stream surface deviates from camber surface more or less [10, 11]. 2D twisted grid can be constructed along stream surface in blade zone, with meridian plane gird in nonblade zone [9, 12]. The derived flow equation in blade zone is supposed to satisfy the flow tangential condition on stream surface. Alternatively, from another perspective, the entire grid can be regarded as a meridian plane one, and the inclination of fluid direction in blade zone is imposed by an additional inviscid blade force term fb [13]. In other cases, both blade suction and pressure side surfaces are demanded, and a single grid cell is set circumferentially in the flow path filling the whole flow passage [11, 14]. Generally, the requirement of a single cell in circumferential direction is related to the corresponding flow solver, which is modified from a 3D solver.For design application, since the profiles of blades are unknown, the formulation of S2 stream surface along with the correspondingfb becomes the key point. Controllable velocity distribution is usually used in turbomachinery design [15, 16]. One way to produce stream surface and fb is to use distribution of swirl vθr or circumferential velocity vθ during blade zone [8, 17]. Stream surface geometry evolves during iteration to meet the desired vθ from some initial value. For example, Damle et al. [18] specified vθr as design target, and stream surface were updated through flow tangential condition formula. Then, source vector about fb was renewed. Yao and Hirsch [19] implemented the Euler equation in Cartesian system in his design work, where relative speed wθ was given. The calculation was conducted in meridian plane, and the speed component wx or wy was renewed in each step. After velocity vector w⟶ was regenerated, fb was then recalculated. Both Simon [20] and Hirsch and Sturmayr [21] used target vθ and its computed value to formulate the partial derivative of fb with respect to a pseudo time parameter. As iteration convergence progressed, the difference between computed vθ and the target became less, and fb became convergent. Taddei and Bertini and Taddei and Larocca [22, 23] set spanwise distribution of specific work as the design target, which was proportional to the difference between rotor inlet and exit vθr. They used meridian convective equation to replace the original energy equation in Euler’s system to update fb through vθr. S2 surfaces waved around their fastened leading edges until they reached the steady shape that was consistent with the specified vθr distribution.Apart from the above examples, direct calculation is also used in optimum design. For instance, Ivanov et al. [24] used an analysis method in turbine optimization process with output power as the design target. However, detailed optimization methodology was not stated in his paper. Pasquale et al. [25] incorporated axisymmetric ZFLOW code with metamodel assisted genetic algorithms. They established a pseudo-objective function that must be minimized to reduce deviation from target tangential velocity, mixed-out total enthalpy, and entropy production. Blade row exit flow angle was taken as a design parameter, and inlet flow angle was obtained from upstream condition. After the best stream surface was found, geometric surface of blade was then calculated via incidence and deviation angle model.Basically, the majority of inverse Q3D time-marching design methods contain the iteration ofvθ and fb throughout blade zone. However, the iteration procedure can be time consuming and unstable. As Jin et al. [17] said, the convergence problem often arose when iteratively solving fb. The contribution of the present work is to develop a fast and reliable Q3D turbine design method. To achieve this, a novel S2 stream surface approximation approach is proposed. In this way, no more formulation of fb is needed, and the iteration of vθ is simplified. The design method determines the stream surface circumferential coordinate θ in blade zone first in each iteration step by means of quadratic function of axial coordinate x, and then calculation is conducted. Radial distribution of blade row exit vθr is required as convergence target. A two-stage turbine and a single-stage high-pressure turbine were designed and analyzed by the Q3D method, and their 3D aerodynamic configurations were built. 3D viscous CFD software NUMECA was used to simulate turbine performances for comparison. Finally, the accuracy of Q3D design method is found to be reliable, and the approximation approach of S2 stream surface is proven to be reasonable.


## 2. Quasi Three-Dimensional Design Method

### 2.1. Governing Equations

The Euler equations in orthogonal curvilinear systemξ,η,ζ in the relative frame are used here, which are derived from cylindrical system x,r,θ. The coordinate system is illustrated in Figure 1.Figure 1
Coordinate system.The equations are written in conservative form [12]
(1)∂∂trU¯J+∂∂ξrJF¯∂ξ∂x+G¯∂ξ∂r+H¯∂ξr∂θ+∂∂ηrJF¯∂η∂x+G¯∂η∂r+H¯∂ηr∂θ+∂∂ζrJF¯∂ζ∂x+G¯∂ζ∂r+H¯∂ζr∂θ=h¯J,where
(2)U¯=ρρvxρvrρwθe,F¯=ρvxρvx2+pρvxvrρvxwθe+pvx,G¯=ρvrρvxvrρvr2+pρvrwθe+pvr,H¯=ρwθρvxwθρvrwθρwθ2+pe+pwθ,h¯=00p+ρwθ+ωr2−ρvrwθ+2ωrω2r2ρvr.J is the Jacobian determinant, which is in the form of
(3)J=∂ξ,η,ζ∂x,r,θ.vx, vr, and wθ are velocity components in relative cylindrical coordinates. ω is the angular speed of turbine spool. Then, the absolute circumferential speed vθ can be expressed as wθ+ωr. ρ denotes density, p is the pressure, and e is the total energy. The equation can be simplified under the Q3D hypothesis. In the axial adjacent blade row gap, i.e., nonblade zone, the flow is assumed to be axisymmetric, S2 stream surface is meridian plane, and
(4)ζ=θ.Another three relations are deduced(5)∂ξ∂θ=0,∂η∂θ=0,∂∂ζ=0.Consequently, Formula (1) can be simplified as
(6)∂∂trU¯J+∂∂ξrJF¯∂ξ∂x+G¯∂ξ∂r+∂∂ηrJF¯∂η∂x+G¯∂η∂r=h¯J.While in blade zone, fluid flow direction is tangential to S2 stream surface, flow surface is twisted, and flow tangential condition must be satisfied(7)vx∂ζ∂x+vr∂ζ∂r+wθ∂ζr∂θ=0.Thus, Formula (1) becomes
(8)∂∂trU¯J+∂∂ξrJF¯∂ξ∂x+G¯∂ξ∂r+H¯∂ξr∂θ+∂∂ηrJF¯∂η∂x+G¯∂η∂r+H¯∂ηr∂θ=h¯J+h¯1,where
(9)h¯1=0−r∂ζJ∂x∂p∂ζ−p∂∂ζr∂ζJ∂x−r∂ζJ∂r∂p∂ζ−p∂∂ζr∂ζJ∂r−∂ζJ∂θ∂p∂ζ−p∂∂ζ∂ζJ∂θ0.−r∂ζ/J∂x∂p/∂ζ, −r∂ζ/J∂r∂p/∂ζ, and −∂ζ/J∂θ∂p/∂ζ can be equivalently seen as components of inviscid blade force vector f⟶b in the x, r, and θ directions. Formulas (6) and (7) can be written in a universal form
(10)∂U∂t+∂F∂ξ+∂G∂η=H,where
(11)U=rU¯J,F=rJF¯∂ξ∂z+G¯∂ξ∂r+H¯∂ξr∂θ,G=rJF¯∂η∂z+G¯∂η∂r+H¯∂ηr∂θ,H=h¯J+h¯1.To avoid erroneous pressure fluctuations generated by the conservative wave model [26], Equation (8) is next transformed to nonconservative form [27]
(12)∂U~∂t+∂F~∂ξ+∂G~∂η=H,~where
(13)U~=rJρvxvrwθp,∂F~∂ξ=∂U∂U~−1∂F∂U~∂U~∂ξ,∂G~∂η=∂U∂U~−1∂G∂U~∂U~∂η,H~=∂U∂U~−1H.Since different flow conditions are assumed in blade and nonblade zones, special flow data transfer is needed between blade and nonblade zones. Li et al.’s [27] method is used to conduct the data transfer at blade leading and trailing edges, which mainly utilizes the conservation of mass, energy, and entropy, along with momentum equation. For the conservation of energy, it means the total enthalpy is invariant.


### 2.2. Numerical Schemes

Finite volume method is applied to solve the differential equations. Semi-implicit second order Crank-Nicolson’s scheme is used in the time orientation in Equation (12), and the temporal discrete difference equation is
(14)U~n+1−U~n∆t=12∂U~∂tn+1+12∂U~∂tn.Substituting Equation (12) into the above formula, the result is
(15)I+Δt∂An2∂ξ+Δt∂Bn2∂ηΔU~n=−Δt∂F~n∂ξ+∂Gn∂η−H~n,where
(16)ΔU~n=U~n+1−U~n,A=∂F~∂U~,B=∂G~∂U~.Then, approximate factorization (AF) method is implemented to reduce the computation requirement, and Equation (15) becomes
(17)I+∆t∂An2∂ξI+Δt∂Bn2∂ηΔU~n=−Δt∂F~n∂ξ+∂Gn∂η−H~n.Since eigenvalues of Jacobian matricesA and B can reflect wave propagation direction, A and B can be split with positive and negative eigenvalues
(18)A=A+−A−,(19)B=B+−B−,where
(20)A+=SξΛξ+Sξ−1,A−=SξΛξ−Sξ−1,B+=SηΛη+Sη−1,B−=SηΛη−Sη−1,Λξ+=12Λξ+Λξ,Λξ−=12Λξ−Λξ,Λη+=12Λη+Λη,Λη−=12Λη−Λη,A=SξΛξSξ−1,B=SηΛηSη−1.Λξ and Λη are diagonal matrices whose diagonal elements are the eigenvalues of A and B, while Λξ and Λη stand for diagonal matrices whose diagonal elements are the absolute values of those of Λξ and Λη. Here, forward and backward differences are defined as
(21)∆f=fi+1−fi,∇f=fi−fi−1.Substitute (18) and (19) into Equation (17), and utilize upwind scheme to discretize A+, A−, B+, and B− with local matrix freezing treatment at the same time. For the right-hand side terms, flux difference splitting (FDS) method is used. Finally, Equation (17) turns to
(22)SξI+Δt2Δξ∇Λξ+n−∇Λξ−nSξ−1SηI+Δt2Δη∇Λη+n−∇Λη−nSη−1ΔU~n=−ΔtF~i+1/2n−F~i−1/2nΔξ+G~i+1/2n−G~i−1/2nΔη−H~n.Although the approximate local linear treatment of matricesA+, A−, B+, and B− reduces the time accuracy, it is not a problem for steady computation in the present research. F~i+1/2, F~i−1/2, G~i+1/2, and G~i−1/2 are calculated by the primitive parameter fluxes U~i+1/2 and U~i−1/2. The Riemann exact solution is used to compute U~i+1/2 and U~i−1/2 at the grid interface. To realize this, the left and right values of the interface need to be reconstructed. Here, a piecewise parabolic distribution of arbitrary primitive parameter u is commenced
(23)ui+1/2n−=uin+1+φ4Δu^jn+1−φ4∇u^jn,ui−1/2n−=uin−1−φ4Δu^jn−1+φ4∇u^jn,where Δu^jn are ∇u^jn restricted increments to meet the monotonicity condition or total variation diminishing (TVD) property
(24)Δu^jn=Δujnχ∇ujn∆ujn,∇u^jn=∇ujnχ∆ujn∇ujn.The limiterχ is written as
(25)χr=0,r≤0,r,0<r<1,1,r≥1.Whenϕ is 1/3, the scheme has third order accuracy in space orientation which is a higher order Godunov’s scheme. Finally, since Equation (22) is the 5×5 block tridiagonal matrix, it can be solved by the chasing method in ξ and η orientations alternately.


### 2.3. Aerodynamic Design Method

To carry out Q3D turbine aerodynamic design, several necessary boundary conditions and design parameters are needed, including turbine inlet total pressurepin∗ and total temperature Tin∗, outlet static pressure pe, and spool rotational speed N. Geometry information of meridian flow path shape and number of blades Nb in each row are also needed. In addition, radial distribution of vθr at blade row exit should be prescribed as well. For the time-marching method, the numerical procedure of inverse design is essentially the same as direct calculation, even though the input conditions are different. Based on the CFD methodology, the computational domain and grid are in need at first. Therefore, the formulation of the physical S2 stream surface in blade zone is a key point in Q3D time-marching design work. A prevalent design process is to establish an original stream surface according to some initial conditions, and then automatically adjust surface shape during the iteration loop until vθ or vθr goal inside blade row is fulfilled. At the same time, inviscid blade is also recalculated repeatedly in the design loop, which forces flow to make a turn and influence vθ in return.A novel S2 stream surface approximation method is proposed here to decrease the iteration requirement. Since rotor specific work is only related to thevθr difference between inlet and exit, the radial distribution of vθr at blade exit is regarded as the design goal instead of the whole internal blade zone. The S2 stream surface in blade zone is formulated through distribution of θ, and θ is estimated via function of x coordinate. Noera and Satta [28] used the quadratic polynomial function of x to calculate θ of blade camber line in his SLC analysis work. Deviation angle is predicted by the Ainley and Mathieson model. Considering that blade camber line is close to streamline, Noera and Satta’s research proved the rationality of quadratic distribution of θ of streamline to some degree. To investigate the geometric characteristics of turbine blade camber line a step further, a series of camber lines of reference blades with normalized coordinates are calculated herein, as shown in Figure 2. The blade profiles come from NASA uncooled core turbine mean blade section [29], NASA high temperature engine turbine mean blade section [30], and E/TU-1 and T106 cascade section [31]. In Figure 2, camber line axial and circumferential coordinates are normalized by axial chord length Cx. The leading edges of lines are all set at the origin point in the coordinate system.Figure 2
Turbine blade camber lines.As depicted in Figure2, it can be speculated that the turbine blade camber line might be approximated by the quadratic function of x. Therefore, the distribution function of θ in blade zone is defined here similar to Noera and Satta’s, which is in the form of
(26)θx=a2x2+a1x+a0,xLE≤x≤xTE,where a0, a1, and a2 are constant coefficients. xLE and xTE denote the axial coordinates of the blade leading and trailing edges, respectively. To calculate these three constants, three additional conditions are added. First, the blade leading and trailing edge flow angles, βLE and βTE, are given, and the function should satisfy
(27)dθxLEdx=2a2xLE+a1,dθxTEdx=2a2xTE+a1.In addition, thex coordinate x0, where θ is 0 is also given. Then,
(28)θx0=a2x02+a1x0+a0=0.a0, a1, and a2 can then be obtained. Another parameter in need for Q3D design is the blade blockage factor, which reflects the influence of metal blade blocking on flow. Blade geometry information is typically used to define blockage factor, and both circumferential and normal blockage factors are often used [4, 10, 14]. Here, to accommodate the flow governing equations in the design method, ∂ζ/∂θ−1 is taken as blockage [10]. The physical meaning of ∂ζ/∂θ−1 is illustrated as follows. In orthogonal curvilinear system, the projected area Aηζ of flow unit volume on the ηζ plane is [12].
(29)Aηζ=∆VΔξ=r/J∂x/∂ξ2+∂r/∂ξ2+r∂θ/∂ξ2=1∂ζ/∂θr∂x/∂ξ∂x/∂η∂r/∂ξ∂r/∂η∂x/∂ξ2+∂r/∂ξ2+r∂θ/∂ξ2,where ΔV is the unit volume and Δξ is the edge length along the ξ axis of unit volume. Aηζ reflects the through flow area in turbine blade row passage and mass flow capability. All partial derivatives in the above formula can be calculated by numerical difference once the computational grid is given except for ∂ζ/∂θ−1. Therefore ∂ζ/∂θ−1 is regarded as the blade blockage factor, which needs to be specified manually. In nonblade zone, ∂ζ/∂θ−1 is always equal to 1. In blade zone, ∂ζ/∂θ−1 is less than 1. Blockage is distributed through an arbitrary smooth function in blade zone. Normally, the minimum blockage is located at the half Cx position. For inviscid blade force fb, many other design studies add a step of formulation and update of fb. For the present research, since the flow tangential condition is incorporated into the flow equation within blade zone, Formula (8) inherently contains the effect of fb, which no longer requires an additional special formulation anymore.After the flow field is computed, blade geometry parameters can be calculated. Incidence anglei is assigned 0, thus inlet blade angle βk1 equals inlet flow angle β1. Deviation angle δ is estimated by the empirical model in reference [32]. This model comes from the generalization of linear turbine cascade experiment results. δ result depends on blade row exit Mach number Mar2 and flow angle β2. The model can be expressed in the form of a chart as shown in Figure 3. The definitions of i and δ are
(30)i=β1−βk1,δ=β2−βk2.Figure 3
Deviation angle model.Other geometrical information of blade leading and trailing edge radius, stagger angle, thickness distribution, rotor blade tip clearance, etc., is also needed to complete blade profiling. The values of these parameters have a broader tolerance to maintain turbine cascade flow performance thanβk1 and βk2. The latter two essentially dominate mass flow rate and expansion level of cascade flow. Consequently, extra geometrical parameters are decided manually as long as the cascade passage has favorable convergence.The major Q3D design process is focused on the iteration of blade row exitvθr. When applying the Q3D design method, the radial distribution of the target swirl vθr∗ should be specified as well as the initial βLE and βTE for each blade row. The design process is illustrated in Figure 4. There are mainly two levels of loops. ε1 and ε2 denote convergence tolerance of the two loops in Figure 4. The inner loop is iteration of βLE, aiming to eliminate incidence at blade leading edge. In each step, βLE is renewed by β1:
(31)βLE′=β1,where βLE′ is the updated blade leading edge flow angle. In this way, βLE finally converges to β1. The outer loop is the iteration of blade row exit vθr. βTE reflects the degree of flow turning in blade row and has direct impact on vθr. The larger βTE is, the higher vθr is. As a consequent, βTE is corrected proportionally to the difference between vθr∗ and vθr in every step:
(32)βTE′=βTE+Kvθr∗−vθr,where βTE′ is the updated blade trailing edge flow angle, and K is a constant that is assigned 0.1 °·s/m2. To avoid an endless loop, the maximum number of iterations Imax is set. Since βLE and βTE are related to the stream surface shape, S2 surface continues to transform until convergence.Figure 4
Q3D design process.Off-design performance prediction is also important in the turbine design phase. Therefore, the Q3D method is also applied to off-design analysis with some modifications. The analysis process based on the geometry information produced by the design process is shown in Figure5. In the analysis process, the outer loop is replaced by an iteration of βTE. After Q3D design is finished, βk1 and βk2 are obtained. They are used as the initial βLE and βTE. In Figure 5, βTE′ stands for updated βTE, which is calculated by deviation angle model and βk2(33)βTE′=βk2+δ.Figure 5
Q3D analysis process.In the context of maintaining radial stage work distribution, the present design process only needs blade row exitvθr to satisfy the design target, simpler than those demanding the iteration of parameters throughout the whole blade row passage. Simultaneously, modeling and iteration of fb is no more needed, which may cause convergence problem as Jin et al. [17] pointed out. In this way, iteration requirement can be much reduced which is beneficial to the computation efficiency and robustness.


### 2.4. Loss Models

Viscous dissipative effect is modeled via empirical loss models when solving the Euler equation [33]. Generally, shock loss is thought to be accurate in the time-marching method. Extra losses, including blade profile loss, secondary flow loss, and blade tip clearance leakage loss, are supplemented. Energy loss Z models are used here. First, profile loss Zp is calculated in the following way [32]:
(34)Zp=0.0030.09kp/cosβ2+0.46cosβ1−cosβ2+0.085+0.017,where
(35)kp=−cosβ1sinβ2−β1,β1−β2≥90°,−cosβ1sinβ2−β1,β1−β2<90°,Secondary lossZs is calculated through [32]
(36)Zs=−0.0474cosβ2cosβ1tanβ1−tanβ2H/Cks+0.0118where
(37)ks=1,HC≥2,1−0.252−HC,HC<2.H stands for blade height, and C is the blade chord length. Blade tip clearance leakage loss is estimated through the Balje and Binsley model [34]
(38)Zk=0.0696tanh13dkCHCcosβ∝tanβ2−tanβ1,where dk is the rotor blade tip clearance, and β∝ is the vector mean flow angle in cascade. Next, the overall loss Zall is
(39)Zall=Zp+Zs+Zk.Empirical loss model is vital for performance prediction [35]. For viscous loss prediction in Q3D computation, it is not enough to just use the aforementioned models because these models are summarized from linear cascade flow performance. Therefore, loss parameter should be redistributed radially to reflect 3D effect. Secondary loss is mainly caused by cross flow and boundary layer separation near flow path endwall. Therefore, secondary loss at hub and tip should be higher. Tip clearance leakage loss is generated from rotor blade tip clearance leakage flow and vortex, thus loss near tip must be higher. Sullerey and Kumar [36] proposed that secondary loss could be redistributed radially in a quadratic function style. Under this assumption, secondary loss at midspan was 0, and loss increased from midspan to the endwall. Yang et al. [10] and De Jesus et al. [37] also adopted this synthesis in their Q3D computation. Based on the above work, the same quadratic distribution of Zs is used here. Furthermore, rotor blade tip clearance loss should also be redistributed. According to reference [38], redistribution of Zk is adopted from 0.7 to 1.0 span. The cubic function is used to redistribute Zk. Leakage loss at 0.7 span is 0, and increases to maximum value at 1.0 span. Finally, the radial distribution formulae of Zs and Zk are shown as follows:
(40)ZsH¯=12Z¯sH¯−0.52,0≤H¯≤1,ZkH¯=40000Z¯k81H¯−0.73,0.7≤H¯≤1,where Z¯s and Z¯k are the average secondary and tip clearance leakage losses calculated by the models. H¯ denotes relative spanwise location. Z¯s and Z¯k are separately equal to the radial definite integral of the corresponding distribution function:
(41)Z¯s=∫01ZsH¯dH¯.Z¯k=∫0.71ZkH¯dH¯.Inside blade row, loss coefficient also needs to be distributed axially or streamwise. Usually, loss coefficient predicted by the loss model is converted to an entropy increment, and then entropy gradient forms the function of viscous blade forcefv [8, 10, 13]. The entropy field in blade zone, however, needs an assumption. Persico and Rebay [13] suggested that smooth function distributions could be applicable, including linear, sinusoidal, and exponential functions. Taddei and Larocca [39], on the other hand, pointed out that the dominant contribution to entropy production came from total pressure loss when streamline curvature was neglected, and he formulated fv as a function of total pressure gradient. However, the distribution of total pressure in blade zone still requires assumption. Therefore, it is apparent that some assumptions are inevitable when modeling viscous loss. A simpler treatment is used here, which assumes that the loss coefficient is linearly distributed axially from 0 at blade leading edge to Zall at trailing edge. Further investigation to support this method is not available at present. However, loss result of the computation example could later prove some rationality. After energy loss Z is distributed at each grid cell, velocity loss coefficient φ is computed [38]:
(42)ϕ=1−Z.φ is equal to the real velocity versus isentropic velocity, and φ is used to correct the velocity field instead of adding a source term fv in flow equation.


### 2.5. Performance Validation Method

After fully 3D aerodynamic configuration of turbine is designed by Q3D method, its performance is validated through 3D viscous CFD tool. Here, NUMECA software is used. The 3D computational grid is constructed by AutoGrid software. The grid is in a single channel, and its topological structure is O-H type, as shown in Figure6. Grid y+ near wall is approximately 1, the maximum cell aspect ratio is approximately 3000, and the maximum expansion ratio is around 2.5. Orthogonal angle of the grid cell is greater than 30°. The flow solver is FINE/Turbo software, which solves N–S equation at steady condition. The temporal scheme is fourth-order Runge–Kutta’s approach, and turbulence model is the Spalart-Allmaras model. The flow data transfer at stator/rotor interface is solved by the mixed-plane approach. The same boundary condition used in Q3D computation is applied in 3D simulation.Figure 6
3D computational grid in blade row.


## 2.1. Governing Equations

The Euler equations in orthogonal curvilinear systemξ,η,ζ in the relative frame are used here, which are derived from cylindrical system x,r,θ. The coordinate system is illustrated in Figure 1.Figure 1
Coordinate system.The equations are written in conservative form [12]
(1)∂∂trU¯J+∂∂ξrJF¯∂ξ∂x+G¯∂ξ∂r+H¯∂ξr∂θ+∂∂ηrJF¯∂η∂x+G¯∂η∂r+H¯∂ηr∂θ+∂∂ζrJF¯∂ζ∂x+G¯∂ζ∂r+H¯∂ζr∂θ=h¯J,where
(2)U¯=ρρvxρvrρwθe,F¯=ρvxρvx2+pρvxvrρvxwθe+pvx,G¯=ρvrρvxvrρvr2+pρvrwθe+pvr,H¯=ρwθρvxwθρvrwθρwθ2+pe+pwθ,h¯=00p+ρwθ+ωr2−ρvrwθ+2ωrω2r2ρvr.J is the Jacobian determinant, which is in the form of
(3)J=∂ξ,η,ζ∂x,r,θ.vx, vr, and wθ are velocity components in relative cylindrical coordinates. ω is the angular speed of turbine spool. Then, the absolute circumferential speed vθ can be expressed as wθ+ωr. ρ denotes density, p is the pressure, and e is the total energy. The equation can be simplified under the Q3D hypothesis. In the axial adjacent blade row gap, i.e., nonblade zone, the flow is assumed to be axisymmetric, S2 stream surface is meridian plane, and
(4)ζ=θ.Another three relations are deduced(5)∂ξ∂θ=0,∂η∂θ=0,∂∂ζ=0.Consequently, Formula (1) can be simplified as
(6)∂∂trU¯J+∂∂ξrJF¯∂ξ∂x+G¯∂ξ∂r+∂∂ηrJF¯∂η∂x+G¯∂η∂r=h¯J.While in blade zone, fluid flow direction is tangential to S2 stream surface, flow surface is twisted, and flow tangential condition must be satisfied(7)vx∂ζ∂x+vr∂ζ∂r+wθ∂ζr∂θ=0.Thus, Formula (1) becomes
(8)∂∂trU¯J+∂∂ξrJF¯∂ξ∂x+G¯∂ξ∂r+H¯∂ξr∂θ+∂∂ηrJF¯∂η∂x+G¯∂η∂r+H¯∂ηr∂θ=h¯J+h¯1,where
(9)h¯1=0−r∂ζJ∂x∂p∂ζ−p∂∂ζr∂ζJ∂x−r∂ζJ∂r∂p∂ζ−p∂∂ζr∂ζJ∂r−∂ζJ∂θ∂p∂ζ−p∂∂ζ∂ζJ∂θ0.−r∂ζ/J∂x∂p/∂ζ, −r∂ζ/J∂r∂p/∂ζ, and −∂ζ/J∂θ∂p/∂ζ can be equivalently seen as components of inviscid blade force vector f⟶b in the x, r, and θ directions. Formulas (6) and (7) can be written in a universal form
(10)∂U∂t+∂F∂ξ+∂G∂η=H,where
(11)U=rU¯J,F=rJF¯∂ξ∂z+G¯∂ξ∂r+H¯∂ξr∂θ,G=rJF¯∂η∂z+G¯∂η∂r+H¯∂ηr∂θ,H=h¯J+h¯1.To avoid erroneous pressure fluctuations generated by the conservative wave model [26], Equation (8) is next transformed to nonconservative form [27]
(12)∂U~∂t+∂F~∂ξ+∂G~∂η=H,~where
(13)U~=rJρvxvrwθp,∂F~∂ξ=∂U∂U~−1∂F∂U~∂U~∂ξ,∂G~∂η=∂U∂U~−1∂G∂U~∂U~∂η,H~=∂U∂U~−1H.Since different flow conditions are assumed in blade and nonblade zones, special flow data transfer is needed between blade and nonblade zones. Li et al.’s [27] method is used to conduct the data transfer at blade leading and trailing edges, which mainly utilizes the conservation of mass, energy, and entropy, along with momentum equation. For the conservation of energy, it means the total enthalpy is invariant.


## 2.2. Numerical Schemes

Finite volume method is applied to solve the differential equations. Semi-implicit second order Crank-Nicolson’s scheme is used in the time orientation in Equation (12), and the temporal discrete difference equation is
(14)U~n+1−U~n∆t=12∂U~∂tn+1+12∂U~∂tn.Substituting Equation (12) into the above formula, the result is
(15)I+Δt∂An2∂ξ+Δt∂Bn2∂ηΔU~n=−Δt∂F~n∂ξ+∂Gn∂η−H~n,where
(16)ΔU~n=U~n+1−U~n,A=∂F~∂U~,B=∂G~∂U~.Then, approximate factorization (AF) method is implemented to reduce the computation requirement, and Equation (15) becomes
(17)I+∆t∂An2∂ξI+Δt∂Bn2∂ηΔU~n=−Δt∂F~n∂ξ+∂Gn∂η−H~n.Since eigenvalues of Jacobian matricesA and B can reflect wave propagation direction, A and B can be split with positive and negative eigenvalues
(18)A=A+−A−,(19)B=B+−B−,where
(20)A+=SξΛξ+Sξ−1,A−=SξΛξ−Sξ−1,B+=SηΛη+Sη−1,B−=SηΛη−Sη−1,Λξ+=12Λξ+Λξ,Λξ−=12Λξ−Λξ,Λη+=12Λη+Λη,Λη−=12Λη−Λη,A=SξΛξSξ−1,B=SηΛηSη−1.Λξ and Λη are diagonal matrices whose diagonal elements are the eigenvalues of A and B, while Λξ and Λη stand for diagonal matrices whose diagonal elements are the absolute values of those of Λξ and Λη. Here, forward and backward differences are defined as
(21)∆f=fi+1−fi,∇f=fi−fi−1.Substitute (18) and (19) into Equation (17), and utilize upwind scheme to discretize A+, A−, B+, and B− with local matrix freezing treatment at the same time. For the right-hand side terms, flux difference splitting (FDS) method is used. Finally, Equation (17) turns to
(22)SξI+Δt2Δξ∇Λξ+n−∇Λξ−nSξ−1SηI+Δt2Δη∇Λη+n−∇Λη−nSη−1ΔU~n=−ΔtF~i+1/2n−F~i−1/2nΔξ+G~i+1/2n−G~i−1/2nΔη−H~n.Although the approximate local linear treatment of matricesA+, A−, B+, and B− reduces the time accuracy, it is not a problem for steady computation in the present research. F~i+1/2, F~i−1/2, G~i+1/2, and G~i−1/2 are calculated by the primitive parameter fluxes U~i+1/2 and U~i−1/2. The Riemann exact solution is used to compute U~i+1/2 and U~i−1/2 at the grid interface. To realize this, the left and right values of the interface need to be reconstructed. Here, a piecewise parabolic distribution of arbitrary primitive parameter u is commenced
(23)ui+1/2n−=uin+1+φ4Δu^jn+1−φ4∇u^jn,ui−1/2n−=uin−1−φ4Δu^jn−1+φ4∇u^jn,where Δu^jn are ∇u^jn restricted increments to meet the monotonicity condition or total variation diminishing (TVD) property
(24)Δu^jn=Δujnχ∇ujn∆ujn,∇u^jn=∇ujnχ∆ujn∇ujn.The limiterχ is written as
(25)χr=0,r≤0,r,0<r<1,1,r≥1.Whenϕ is 1/3, the scheme has third order accuracy in space orientation which is a higher order Godunov’s scheme. Finally, since Equation (22) is the 5×5 block tridiagonal matrix, it can be solved by the chasing method in ξ and η orientations alternately.


## 2.3. Aerodynamic Design Method

To carry out Q3D turbine aerodynamic design, several necessary boundary conditions and design parameters are needed, including turbine inlet total pressurepin∗ and total temperature Tin∗, outlet static pressure pe, and spool rotational speed N. Geometry information of meridian flow path shape and number of blades Nb in each row are also needed. In addition, radial distribution of vθr at blade row exit should be prescribed as well. For the time-marching method, the numerical procedure of inverse design is essentially the same as direct calculation, even though the input conditions are different. Based on the CFD methodology, the computational domain and grid are in need at first. Therefore, the formulation of the physical S2 stream surface in blade zone is a key point in Q3D time-marching design work. A prevalent design process is to establish an original stream surface according to some initial conditions, and then automatically adjust surface shape during the iteration loop until vθ or vθr goal inside blade row is fulfilled. At the same time, inviscid blade is also recalculated repeatedly in the design loop, which forces flow to make a turn and influence vθ in return.A novel S2 stream surface approximation method is proposed here to decrease the iteration requirement. Since rotor specific work is only related to thevθr difference between inlet and exit, the radial distribution of vθr at blade exit is regarded as the design goal instead of the whole internal blade zone. The S2 stream surface in blade zone is formulated through distribution of θ, and θ is estimated via function of x coordinate. Noera and Satta [28] used the quadratic polynomial function of x to calculate θ of blade camber line in his SLC analysis work. Deviation angle is predicted by the Ainley and Mathieson model. Considering that blade camber line is close to streamline, Noera and Satta’s research proved the rationality of quadratic distribution of θ of streamline to some degree. To investigate the geometric characteristics of turbine blade camber line a step further, a series of camber lines of reference blades with normalized coordinates are calculated herein, as shown in Figure 2. The blade profiles come from NASA uncooled core turbine mean blade section [29], NASA high temperature engine turbine mean blade section [30], and E/TU-1 and T106 cascade section [31]. In Figure 2, camber line axial and circumferential coordinates are normalized by axial chord length Cx. The leading edges of lines are all set at the origin point in the coordinate system.Figure 2
Turbine blade camber lines.As depicted in Figure2, it can be speculated that the turbine blade camber line might be approximated by the quadratic function of x. Therefore, the distribution function of θ in blade zone is defined here similar to Noera and Satta’s, which is in the form of
(26)θx=a2x2+a1x+a0,xLE≤x≤xTE,where a0, a1, and a2 are constant coefficients. xLE and xTE denote the axial coordinates of the blade leading and trailing edges, respectively. To calculate these three constants, three additional conditions are added. First, the blade leading and trailing edge flow angles, βLE and βTE, are given, and the function should satisfy
(27)dθxLEdx=2a2xLE+a1,dθxTEdx=2a2xTE+a1.In addition, thex coordinate x0, where θ is 0 is also given. Then,
(28)θx0=a2x02+a1x0+a0=0.a0, a1, and a2 can then be obtained. Another parameter in need for Q3D design is the blade blockage factor, which reflects the influence of metal blade blocking on flow. Blade geometry information is typically used to define blockage factor, and both circumferential and normal blockage factors are often used [4, 10, 14]. Here, to accommodate the flow governing equations in the design method, ∂ζ/∂θ−1 is taken as blockage [10]. The physical meaning of ∂ζ/∂θ−1 is illustrated as follows. In orthogonal curvilinear system, the projected area Aηζ of flow unit volume on the ηζ plane is [12].
(29)Aηζ=∆VΔξ=r/J∂x/∂ξ2+∂r/∂ξ2+r∂θ/∂ξ2=1∂ζ/∂θr∂x/∂ξ∂x/∂η∂r/∂ξ∂r/∂η∂x/∂ξ2+∂r/∂ξ2+r∂θ/∂ξ2,where ΔV is the unit volume and Δξ is the edge length along the ξ axis of unit volume. Aηζ reflects the through flow area in turbine blade row passage and mass flow capability. All partial derivatives in the above formula can be calculated by numerical difference once the computational grid is given except for ∂ζ/∂θ−1. Therefore ∂ζ/∂θ−1 is regarded as the blade blockage factor, which needs to be specified manually. In nonblade zone, ∂ζ/∂θ−1 is always equal to 1. In blade zone, ∂ζ/∂θ−1 is less than 1. Blockage is distributed through an arbitrary smooth function in blade zone. Normally, the minimum blockage is located at the half Cx position. For inviscid blade force fb, many other design studies add a step of formulation and update of fb. For the present research, since the flow tangential condition is incorporated into the flow equation within blade zone, Formula (8) inherently contains the effect of fb, which no longer requires an additional special formulation anymore.After the flow field is computed, blade geometry parameters can be calculated. Incidence anglei is assigned 0, thus inlet blade angle βk1 equals inlet flow angle β1. Deviation angle δ is estimated by the empirical model in reference [32]. This model comes from the generalization of linear turbine cascade experiment results. δ result depends on blade row exit Mach number Mar2 and flow angle β2. The model can be expressed in the form of a chart as shown in Figure 3. The definitions of i and δ are
(30)i=β1−βk1,δ=β2−βk2.Figure 3
Deviation angle model.Other geometrical information of blade leading and trailing edge radius, stagger angle, thickness distribution, rotor blade tip clearance, etc., is also needed to complete blade profiling. The values of these parameters have a broader tolerance to maintain turbine cascade flow performance thanβk1 and βk2. The latter two essentially dominate mass flow rate and expansion level of cascade flow. Consequently, extra geometrical parameters are decided manually as long as the cascade passage has favorable convergence.The major Q3D design process is focused on the iteration of blade row exitvθr. When applying the Q3D design method, the radial distribution of the target swirl vθr∗ should be specified as well as the initial βLE and βTE for each blade row. The design process is illustrated in Figure 4. There are mainly two levels of loops. ε1 and ε2 denote convergence tolerance of the two loops in Figure 4. The inner loop is iteration of βLE, aiming to eliminate incidence at blade leading edge. In each step, βLE is renewed by β1:
(31)βLE′=β1,where βLE′ is the updated blade leading edge flow angle. In this way, βLE finally converges to β1. The outer loop is the iteration of blade row exit vθr. βTE reflects the degree of flow turning in blade row and has direct impact on vθr. The larger βTE is, the higher vθr is. As a consequent, βTE is corrected proportionally to the difference between vθr∗ and vθr in every step:
(32)βTE′=βTE+Kvθr∗−vθr,where βTE′ is the updated blade trailing edge flow angle, and K is a constant that is assigned 0.1 °·s/m2. To avoid an endless loop, the maximum number of iterations Imax is set. Since βLE and βTE are related to the stream surface shape, S2 surface continues to transform until convergence.Figure 4
Q3D design process.Off-design performance prediction is also important in the turbine design phase. Therefore, the Q3D method is also applied to off-design analysis with some modifications. The analysis process based on the geometry information produced by the design process is shown in Figure5. In the analysis process, the outer loop is replaced by an iteration of βTE. After Q3D design is finished, βk1 and βk2 are obtained. They are used as the initial βLE and βTE. In Figure 5, βTE′ stands for updated βTE, which is calculated by deviation angle model and βk2(33)βTE′=βk2+δ.Figure 5
Q3D analysis process.In the context of maintaining radial stage work distribution, the present design process only needs blade row exitvθr to satisfy the design target, simpler than those demanding the iteration of parameters throughout the whole blade row passage. Simultaneously, modeling and iteration of fb is no more needed, which may cause convergence problem as Jin et al. [17] pointed out. In this way, iteration requirement can be much reduced which is beneficial to the computation efficiency and robustness.


## 2.4. Loss Models

Viscous dissipative effect is modeled via empirical loss models when solving the Euler equation [33]. Generally, shock loss is thought to be accurate in the time-marching method. Extra losses, including blade profile loss, secondary flow loss, and blade tip clearance leakage loss, are supplemented. Energy loss Z models are used here. First, profile loss Zp is calculated in the following way [32]:
(34)Zp=0.0030.09kp/cosβ2+0.46cosβ1−cosβ2+0.085+0.017,where
(35)kp=−cosβ1sinβ2−β1,β1−β2≥90°,−cosβ1sinβ2−β1,β1−β2<90°,Secondary lossZs is calculated through [32]
(36)Zs=−0.0474cosβ2cosβ1tanβ1−tanβ2H/Cks+0.0118where
(37)ks=1,HC≥2,1−0.252−HC,HC<2.H stands for blade height, and C is the blade chord length. Blade tip clearance leakage loss is estimated through the Balje and Binsley model [34]
(38)Zk=0.0696tanh13dkCHCcosβ∝tanβ2−tanβ1,where dk is the rotor blade tip clearance, and β∝ is the vector mean flow angle in cascade. Next, the overall loss Zall is
(39)Zall=Zp+Zs+Zk.Empirical loss model is vital for performance prediction [35]. For viscous loss prediction in Q3D computation, it is not enough to just use the aforementioned models because these models are summarized from linear cascade flow performance. Therefore, loss parameter should be redistributed radially to reflect 3D effect. Secondary loss is mainly caused by cross flow and boundary layer separation near flow path endwall. Therefore, secondary loss at hub and tip should be higher. Tip clearance leakage loss is generated from rotor blade tip clearance leakage flow and vortex, thus loss near tip must be higher. Sullerey and Kumar [36] proposed that secondary loss could be redistributed radially in a quadratic function style. Under this assumption, secondary loss at midspan was 0, and loss increased from midspan to the endwall. Yang et al. [10] and De Jesus et al. [37] also adopted this synthesis in their Q3D computation. Based on the above work, the same quadratic distribution of Zs is used here. Furthermore, rotor blade tip clearance loss should also be redistributed. According to reference [38], redistribution of Zk is adopted from 0.7 to 1.0 span. The cubic function is used to redistribute Zk. Leakage loss at 0.7 span is 0, and increases to maximum value at 1.0 span. Finally, the radial distribution formulae of Zs and Zk are shown as follows:
(40)ZsH¯=12Z¯sH¯−0.52,0≤H¯≤1,ZkH¯=40000Z¯k81H¯−0.73,0.7≤H¯≤1,where Z¯s and Z¯k are the average secondary and tip clearance leakage losses calculated by the models. H¯ denotes relative spanwise location. Z¯s and Z¯k are separately equal to the radial definite integral of the corresponding distribution function:
(41)Z¯s=∫01ZsH¯dH¯.Z¯k=∫0.71ZkH¯dH¯.Inside blade row, loss coefficient also needs to be distributed axially or streamwise. Usually, loss coefficient predicted by the loss model is converted to an entropy increment, and then entropy gradient forms the function of viscous blade forcefv [8, 10, 13]. The entropy field in blade zone, however, needs an assumption. Persico and Rebay [13] suggested that smooth function distributions could be applicable, including linear, sinusoidal, and exponential functions. Taddei and Larocca [39], on the other hand, pointed out that the dominant contribution to entropy production came from total pressure loss when streamline curvature was neglected, and he formulated fv as a function of total pressure gradient. However, the distribution of total pressure in blade zone still requires assumption. Therefore, it is apparent that some assumptions are inevitable when modeling viscous loss. A simpler treatment is used here, which assumes that the loss coefficient is linearly distributed axially from 0 at blade leading edge to Zall at trailing edge. Further investigation to support this method is not available at present. However, loss result of the computation example could later prove some rationality. After energy loss Z is distributed at each grid cell, velocity loss coefficient φ is computed [38]:
(42)ϕ=1−Z.φ is equal to the real velocity versus isentropic velocity, and φ is used to correct the velocity field instead of adding a source term fv in flow equation.


## 2.5. Performance Validation Method

After fully 3D aerodynamic configuration of turbine is designed by Q3D method, its performance is validated through 3D viscous CFD tool. Here, NUMECA software is used. The 3D computational grid is constructed by AutoGrid software. The grid is in a single channel, and its topological structure is O-H type, as shown in Figure6. Grid y+ near wall is approximately 1, the maximum cell aspect ratio is approximately 3000, and the maximum expansion ratio is around 2.5. Orthogonal angle of the grid cell is greater than 30°. The flow solver is FINE/Turbo software, which solves N–S equation at steady condition. The temporal scheme is fourth-order Runge–Kutta’s approach, and turbulence model is the Spalart-Allmaras model. The flow data transfer at stator/rotor interface is solved by the mixed-plane approach. The same boundary condition used in Q3D computation is applied in 3D simulation.Figure 6
3D computational grid in blade row.


## 3. Research Turbines

A two-stage turbine and a single-stage high-pressure turbine were studied. Both turbine meridian flow paths were already known, as shown in Figure7. To apply the Q3D design method, a personal computer code was developed based on the above theories. The code was compiled by FORTRAN language. In Figure 7, Q3D grids are also demonstrated, and it is clear that the number of grid points is much lower than 3D grid. This dramatically improved Q3D calculation efficiency and saved considerable time. A personal computer equipped with Intel i7-6700 3.4 GHz CPU and 32 GB RAM was used to commence both Q3D and 3D calculations in the research. In each case, Q3D calculation cost less than 1 min, but 3D calculation consumed approximately 30 min. Necessary design parameters were prescribed, which are exhibited in Table 1. In Table 1, the design parameters of the two-stage turbine come from equivalent cold air condition, while the single-stage turbine parameters are from normal hot condition. The turbines were designed by the proposed Q3D method first. Fully 3D aerodynamic configurations were next constructed according to design results. Among turbine blades, the first stage stator blade stacking point was leading edge circle center, while the rest blades were stacking along the centroids. Final meridian blade profiles were a little different from initial input because of axial transportation of blade sections. 3D steady simulation was used to check the design performance to verify the effectiveness of Q3D method. Finally, characteristic curves in design rotational speed were computed by both Q3D and 3D methods.Figure 7
Turbine flow path and Q3D grid: (a) two-stage turbine; (b) single-stage turbine.
(a)(b)Table 1
Turbine design parameters.
ParameterTwo-stage turbineSingle-stage turbinepin∗/Pa1.0759×1051.3245×106Tin∗/K3781561pe/Pa3.5×1043×105N/(r/min)440713284The turbines were both designed in free-vortex style, and constant radial distribution ofvθr was given at each blade row exit. During actual design procedure, swirl goals vθr∗ at only three span locations of hub, midspan, and tip were given, in order to reduce number of iteration times and avert unreasonable twisted blade shape. Input vθr∗ of each turbine case is shown in Tables 2 and 3.Table 2
Input blade row parameters of the two-stage turbine.
Blade rowvθr∗/(m2/s)NbStator 1-73.2050Rotor 112.4561Stator 2-75.6039Rotor 29.0052Table 3
Input blade row parameters of the single-stage turbine.
Blade rowvθr∗/(m2/s)NbStator-250.3024Rotor89.2054


## 4. Results and Analysis

### 4.1. Two-Stage Turbine

The overall performance results of the two-stage turbine are demonstrated in Table4 including both Q3D and 3D results. In the table, performance parameters include mass flow rate ṁ, expansion ratio π, and isentropic efficiency ηi. It can be found that ṁ of Q3D is approximately 0.67% lower than 3D result. Q3D π is lower than 3D result by 0.022, and Q3D ηi has only 0.37 more percentage points. Among them, the accuracy of efficiency is closely related to loss prediction, which will be analyzed in detail later. Overall, the accuracy of Q3D performance results is high, and the turbine has considerably high ηi.Table 4
Overall performance results of the two-stage turbine.
ParameterQ3D3Dṁ/(kg/s)16.3616.47π2.5122.534ηi/%93.0292.65Absolute Mach numberMa and relative Mach number Mar contours in meridian plane are shown in Figures 8 and 9. Total pressure p∗ and total temperature T∗ contours are shown in Figures 10 and 11. In general, bellowing figures show that Q3D meridian results are close to 3D results with favorable expansion and compression wave capturing ability. Although a difference of shapes between approximate and real mean S2 stream surfaces does exist, the meridian contour results could verify that this difference is not that crucial, and the approximation method of stream surface is reasonable. Figure 8 shows that the maximum Ma is produced at the second stage stator exit near hub area. Turbine exit Ma is approximately 0.5, which is a low speed result that would be beneficial to downstream component performance in aeroengine application practice. In 3D result, a small low speed zone emerges at the second stage stator tip due to flow path expansion. However, there is not a noticeable low speed zone at the same location in Q3D result. The lack of viscous boundary simulation in Q3D is presumably the main cause of this error. At the same time, a higher radial Ma gradient can be seen in the second stage stator leading half zone of 3D result than Q3D result. Meanwhile, higher Mach at the first stage stator exit of 3D result can be seen compared to that of Q3D result. In Figure 9, it is found that Mar of both two stage rotors is generally under 0.8, and minimum Mar is generated at rotor tip leading edge area. However, the first stage rotor Mar of 3D result is higher than that of Q3D. Therefore, it can be inferred that the first stage flow expansion ratio of Q3D result is lower than that of 3D result, which mainly causes the lower turbine expansion ratio of Q3D result. In Figure 10, the major difference between Q3D and 3D p∗ results is that stronger radial pressure gradient in blade zone is found in Q3D contour, where p∗ at the mean radius is higher when given the same x coordinate. In Figure 11, the radial distribution of Q3D T∗ result in blade zone is also fluctuating more than 3D result, especially near the second stage rotor tip. This kind of error possibly comes from the stream surface approximation which may result in some unrealistic surface twist.Figure 8
Absolute Mach number of the two-stage turbine: (a) Q3D; (b) 3D.Figure 9
Relative Mach number of the two-stage turbine: (a) Q3D; (b) 3D.Figure 10
Total pressure of the two-stage turbine: (a) Q3D; (b) 3D.Figure 11
Total temperature of the two-stage turbine: (a) Q3D; (b) 3D.The radial distribution ofvθr at each blade row exit is shown in Figure 12. Design target vθr∗ is also demonstrated in Figure 12. Because the plus or minus sign of swirl just stands for the circumferential direction of flow, the analysis later only takes into consideration of the absolute value of swirl. Figure 12(a) shows the first stage stator vθr∗ is basically reached in Q3D result. However, the first stage stator vθr of 3D result is higher than vθr∗ by approximately 4 m2/s. However, a constant swirl distribution is still kept. The second stage stator vθr of Q3D result is higher than vθr∗ below 0.84 span, and the difference is within 2.2 m2/s. Q3D vθr result continually decreases from 0.65 to 1.0 span. The reason lies in the velocity loss correction in the upstream rotor by tip clearance loss model. The 3D result, on the other hand, is close to the target at midspan area. However, vθr near tip is higher and vθr near hub is lower. The reason could lie in the influence of casing endwall boundary layer and corner vortex near hub. In Figure 12(b), it can be seen that vθr of both stage rotors in Q3D result is close to vθr∗ from 0 to 0.5 span. However, there is a relatively large difference from target at approximately 0.85 span in both rotors, although all of vθr∗ targets are basically achieved at the prescribed span locations, i.e., 0, 0.5, and 1.0 span locations. The maximum difference between vθr and vθr∗ in the first stage rotor is 6.4 m2/s, and that in the second stage rotor is nearly 9 m2/s. This may indicate that more radial locations are needed for iteration of vθr to make sure the required vθr∗ distribution can be completely fulfilled. But as mentioned before, the more prescribed locations there are, the more twisted blade could be. Therefore, there should be a compromise in engineering practice. In 3D result, the first stage rotor vθr is even greater than Q3D result by approximately 4 m2/s from 0.5 to 0.9 span. Nevertheless, its variation trend is basically the same as Q3D result in the main flow area. The second stage rotor vθr is unexpectedly close to target from 0.1 to 0.9 span.Figure 12
Swirl of the two-stage turbine: (a) stator; (b) rotor.
(a)(b)Stator exit absolute flow angleα1 and rotor exit relative flow angle β2 are shown in Figure 13. The analysis of flow angle also only takes into consideration of the absolute value. In Figure 13, flow angle of Q3D result is generally very close to 3D result. In Figure 13(a), the first stage stator α1 of Q3D result is generally lower than 3D result. The largest error is approximately -2.5°. The remaining blade rows all reflect high accuracy of Q3D flow angle prediction results. Considering the previous Mach number comparison, it can be deduced that the difference in vθr between Q3D and 3D results mainly comes from the speed magnitude error, expect for the first stage stator. Furthermore, it can be concluded that radial variation of blade row exit flow angle does not necessarily linearly depend on swirl, such as flow angle distributions in the second stage stator and rotor. Because axial speed vx influences flow angle simultaneously. Therefore, even though free vortex design methodology is applied, the final blade shape can be twisted.Figure 13
Flow angle of the two-stage turbine: (a) stator flow angle; (b) rotor relative flow angle.
(a)(b)Stator total pressure lossY and rotor relative total pressure loss Yr are shown in Figure 14. In Figure 14, Q3D blade row total pressure losses basically also shows good accuracy, which is important for ηi prediction. In Figure 14(a), the first stage stator Y of Q3D result is very close to 3D result except near tip area. 3D Y result near tip shows a trend of reduction from 0.8 to 0.9 span, which is different from Q3D result. The reason lies in the radial migration of secondary low energy flow. In the second stage stator, this loss migration is even fiercer at tip and hub. As a result, Y of Q3D result is lower than 3D result from 0.2 to 0.85 span, and the largest error is approximately -0.06. This error indicates that the present loss model could be improved with some introduction of radial loss migration in the future. In Figure 14(b), the error of Q3D Yr resulting from 3D result mainly exists above 0.5 span. Basically, distribution of Q3D Yr result matches 3D result from 0.1 to 0.9 span quite well.Figure 14
Total pressure loss of the two-stage turbine: (a) stator total pressure loss; (b) rotor relative total pressure loss.
(a)(b)ṁ and ηi characteristic curves at design N are exhibited in Figure 15. From Figure 15, it can be seen that Q3D performance variation with π is basically close to 3D result. Both Q3D ṁ and ηi results are lower than 3D results.Figure 15
Characteristic curves of the two-stage turbine: (a) mass flow rate; (b) efficiency.
(a)(b)


### 4.2. Single-Stage Turbine

The overall performance results of the single-stage high-pressure turbine are shown in Table5. Compared with the two-stage turbine, the single-stage has much higher stage π and N. As illustrated in Table 5, Q3D ṁ result is higher than 3D result by approximately 1.25%. Q3D π is only 0.006 lower than 3D result, and Q3D ηi has just 0.03 more percentage points. Among the overall performance parameters, Q3D ṁ result has the highest error. Basically, Q3D performance results are acceptable.Table 5
Overall performance results of the single-stage turbine.
ParameterQ3D3Dṁ/(kg/s)38.1537.68π3.6573.663ηi/%90.8190.78MeridianMa and Mar contours are shown in Figures 16 and 17. Meridian p∗ and T∗ contours are shown in Figures 18 and 19. In all, Q3D meridian aerodynamic parameter contours are close to 3D results. In Figure 16, compared with 3D result, Q3D stator exit Ma is lower, and lower rotor blade zone Ma radial gradient can be seen in Q3D contour. Lower stator exit Ma of Q3D result implies lower degree of flow expansion in stator, which means smaller throat area in stator. As a consequence, Q3D ṁ result is higher than 3D. Figure 17 shows that Q3D Mar near rotor trailing edge is lower than 3D result, and rotor Mar axial gradient is higher in 3D result. Figure 18 shows that Q3D p∗ radial gradient in stator is higher than 3D result. Finally in Figure 19, the major error of Q3D T∗ lies in rotor tip area where Q3D T∗ is lower than 3D result.Figure 16
Absolute Mach number of the single-stage turbine: (a) Q3D; (b) 3D.Figure 17
Relative Mach number of the single-stage turbine: (a) Q3D; (b) 3D.Figure 18
Total pressure of the single-stage turbine: (a) Q3D; (b) 3D.Figure 19
Total temperature of the single-stage turbine: (a) Q3D; (b) 3D.The radial distribution ofvθr at blade row exit is shown in Figure 20. Figure 20(a) shows that Q3D stator exit vθr features a higher value at midspan. The absolute value of difference between Q3D vθr and vθr∗ is within 5 m2/s. For 3D stator result, vθr is close to the target from 0.1 to 0.7 span. At 0.9 span, difference from vθr∗ is nearly 6 m2/s. In Figure 20(b), both Q3D and 3D vθr results are close to the target value from 0.1 to 0.8 span. However, their vθr values at rotor tip zone are much lower than vθr∗. This may indicate that tip clearance leakage loss makes it impossible to realize the swirl target at tip.Figure 20
Swirl of the single-stage turbine: (a) stator; (b) rotor.
(a)(b)Statorα1 and rotor β2 are shown in Figure 21. In Figure 21, both stator α1 and rotor β2 of Q3D results exhibit high accuracy. From 0.1 to 0.9 span, the absolute value of Q3D stator α1 error is within 0.3°, and the absolute value of rotor β2 error is no more than 2°. Considering the Mach number difference between Q3D and 3D contours, it can be speculated that the error of Q3D result is mainly caused by the error of speed magnitude.Figure 21
Flow angle of the single-stage turbine: (a) stator flow angle; (b) rotor relative flow angle.
(a)(b)StatorY and rotor Yr are shown in Figure 22. In Figure 22(a), 3D stator Y result remains nearly constant along 0.1 to 0.9 span. However, Q3D stator Y essentially has a quadric distribution radially. This is caused by the secondary loss model in Q3D method. For the single-stage turbine, its stator tip endwall has an S shape profile, which is beneficial for reducing secondary flow strength. Therefore, radial migration of loss in 3D result is not remarkable. At the same time, stator hub boundary flow separation loss in 3D result is also depressed with high loss range is restricted below 0.1 span. In Figure 22(b), radial loss migration can be seen at rotor hub zone in 3D result, but not very strong. Both stator and rotor losses of Q3D results are lower than that of 3D results at midspan, but higher around endwalls. In the future, a more sophisticated loss migration model should be considered to make loss redistribution more accurate.Figure 22
Total pressure loss of the single-stage turbine: (a) stator total pressure loss; (b) rotor relative total pressure loss.
(a)(b)ṁ and ηi characteristic curves at design N are exhibited in Figure 23. The ṁ performance of Q3D result is generally close to 3D. For ηi characteristic curve, the change in Q3D ηi result is stronger when π is above 3. The cause of this error should be the overestimated loss by the loss model in the Q3D method.Figure 23
Characteristic curves of the single-stage turbine: (a) mass flow rate; (b) efficiency.
(a)(b)


## 4.1. Two-Stage Turbine

The overall performance results of the two-stage turbine are demonstrated in Table4 including both Q3D and 3D results. In the table, performance parameters include mass flow rate ṁ, expansion ratio π, and isentropic efficiency ηi. It can be found that ṁ of Q3D is approximately 0.67% lower than 3D result. Q3D π is lower than 3D result by 0.022, and Q3D ηi has only 0.37 more percentage points. Among them, the accuracy of efficiency is closely related to loss prediction, which will be analyzed in detail later. Overall, the accuracy of Q3D performance results is high, and the turbine has considerably high ηi.Table 4
Overall performance results of the two-stage turbine.
ParameterQ3D3Dṁ/(kg/s)16.3616.47π2.5122.534ηi/%93.0292.65Absolute Mach numberMa and relative Mach number Mar contours in meridian plane are shown in Figures 8 and 9. Total pressure p∗ and total temperature T∗ contours are shown in Figures 10 and 11. In general, bellowing figures show that Q3D meridian results are close to 3D results with favorable expansion and compression wave capturing ability. Although a difference of shapes between approximate and real mean S2 stream surfaces does exist, the meridian contour results could verify that this difference is not that crucial, and the approximation method of stream surface is reasonable. Figure 8 shows that the maximum Ma is produced at the second stage stator exit near hub area. Turbine exit Ma is approximately 0.5, which is a low speed result that would be beneficial to downstream component performance in aeroengine application practice. In 3D result, a small low speed zone emerges at the second stage stator tip due to flow path expansion. However, there is not a noticeable low speed zone at the same location in Q3D result. The lack of viscous boundary simulation in Q3D is presumably the main cause of this error. At the same time, a higher radial Ma gradient can be seen in the second stage stator leading half zone of 3D result than Q3D result. Meanwhile, higher Mach at the first stage stator exit of 3D result can be seen compared to that of Q3D result. In Figure 9, it is found that Mar of both two stage rotors is generally under 0.8, and minimum Mar is generated at rotor tip leading edge area. However, the first stage rotor Mar of 3D result is higher than that of Q3D. Therefore, it can be inferred that the first stage flow expansion ratio of Q3D result is lower than that of 3D result, which mainly causes the lower turbine expansion ratio of Q3D result. In Figure 10, the major difference between Q3D and 3D p∗ results is that stronger radial pressure gradient in blade zone is found in Q3D contour, where p∗ at the mean radius is higher when given the same x coordinate. In Figure 11, the radial distribution of Q3D T∗ result in blade zone is also fluctuating more than 3D result, especially near the second stage rotor tip. This kind of error possibly comes from the stream surface approximation which may result in some unrealistic surface twist.Figure 8
Absolute Mach number of the two-stage turbine: (a) Q3D; (b) 3D.Figure 9
Relative Mach number of the two-stage turbine: (a) Q3D; (b) 3D.Figure 10
Total pressure of the two-stage turbine: (a) Q3D; (b) 3D.Figure 11
Total temperature of the two-stage turbine: (a) Q3D; (b) 3D.The radial distribution ofvθr at each blade row exit is shown in Figure 12. Design target vθr∗ is also demonstrated in Figure 12. Because the plus or minus sign of swirl just stands for the circumferential direction of flow, the analysis later only takes into consideration of the absolute value of swirl. Figure 12(a) shows the first stage stator vθr∗ is basically reached in Q3D result. However, the first stage stator vθr of 3D result is higher than vθr∗ by approximately 4 m2/s. However, a constant swirl distribution is still kept. The second stage stator vθr of Q3D result is higher than vθr∗ below 0.84 span, and the difference is within 2.2 m2/s. Q3D vθr result continually decreases from 0.65 to 1.0 span. The reason lies in the velocity loss correction in the upstream rotor by tip clearance loss model. The 3D result, on the other hand, is close to the target at midspan area. However, vθr near tip is higher and vθr near hub is lower. The reason could lie in the influence of casing endwall boundary layer and corner vortex near hub. In Figure 12(b), it can be seen that vθr of both stage rotors in Q3D result is close to vθr∗ from 0 to 0.5 span. However, there is a relatively large difference from target at approximately 0.85 span in both rotors, although all of vθr∗ targets are basically achieved at the prescribed span locations, i.e., 0, 0.5, and 1.0 span locations. The maximum difference between vθr and vθr∗ in the first stage rotor is 6.4 m2/s, and that in the second stage rotor is nearly 9 m2/s. This may indicate that more radial locations are needed for iteration of vθr to make sure the required vθr∗ distribution can be completely fulfilled. But as mentioned before, the more prescribed locations there are, the more twisted blade could be. Therefore, there should be a compromise in engineering practice. In 3D result, the first stage rotor vθr is even greater than Q3D result by approximately 4 m2/s from 0.5 to 0.9 span. Nevertheless, its variation trend is basically the same as Q3D result in the main flow area. The second stage rotor vθr is unexpectedly close to target from 0.1 to 0.9 span.Figure 12
Swirl of the two-stage turbine: (a) stator; (b) rotor.
(a)(b)Stator exit absolute flow angleα1 and rotor exit relative flow angle β2 are shown in Figure 13. The analysis of flow angle also only takes into consideration of the absolute value. In Figure 13, flow angle of Q3D result is generally very close to 3D result. In Figure 13(a), the first stage stator α1 of Q3D result is generally lower than 3D result. The largest error is approximately -2.5°. The remaining blade rows all reflect high accuracy of Q3D flow angle prediction results. Considering the previous Mach number comparison, it can be deduced that the difference in vθr between Q3D and 3D results mainly comes from the speed magnitude error, expect for the first stage stator. Furthermore, it can be concluded that radial variation of blade row exit flow angle does not necessarily linearly depend on swirl, such as flow angle distributions in the second stage stator and rotor. Because axial speed vx influences flow angle simultaneously. Therefore, even though free vortex design methodology is applied, the final blade shape can be twisted.Figure 13
Flow angle of the two-stage turbine: (a) stator flow angle; (b) rotor relative flow angle.
(a)(b)Stator total pressure lossY and rotor relative total pressure loss Yr are shown in Figure 14. In Figure 14, Q3D blade row total pressure losses basically also shows good accuracy, which is important for ηi prediction. In Figure 14(a), the first stage stator Y of Q3D result is very close to 3D result except near tip area. 3D Y result near tip shows a trend of reduction from 0.8 to 0.9 span, which is different from Q3D result. The reason lies in the radial migration of secondary low energy flow. In the second stage stator, this loss migration is even fiercer at tip and hub. As a result, Y of Q3D result is lower than 3D result from 0.2 to 0.85 span, and the largest error is approximately -0.06. This error indicates that the present loss model could be improved with some introduction of radial loss migration in the future. In Figure 14(b), the error of Q3D Yr resulting from 3D result mainly exists above 0.5 span. Basically, distribution of Q3D Yr result matches 3D result from 0.1 to 0.9 span quite well.Figure 14
Total pressure loss of the two-stage turbine: (a) stator total pressure loss; (b) rotor relative total pressure loss.
(a)(b)ṁ and ηi characteristic curves at design N are exhibited in Figure 15. From Figure 15, it can be seen that Q3D performance variation with π is basically close to 3D result. Both Q3D ṁ and ηi results are lower than 3D results.Figure 15
Characteristic curves of the two-stage turbine: (a) mass flow rate; (b) efficiency.
(a)(b)


## 4.2. Single-Stage Turbine

The overall performance results of the single-stage high-pressure turbine are shown in Table5. Compared with the two-stage turbine, the single-stage has much higher stage π and N. As illustrated in Table 5, Q3D ṁ result is higher than 3D result by approximately 1.25%. Q3D π is only 0.006 lower than 3D result, and Q3D ηi has just 0.03 more percentage points. Among the overall performance parameters, Q3D ṁ result has the highest error. Basically, Q3D performance results are acceptable.Table 5
Overall performance results of the single-stage turbine.
ParameterQ3D3Dṁ/(kg/s)38.1537.68π3.6573.663ηi/%90.8190.78MeridianMa and Mar contours are shown in Figures 16 and 17. Meridian p∗ and T∗ contours are shown in Figures 18 and 19. In all, Q3D meridian aerodynamic parameter contours are close to 3D results. In Figure 16, compared with 3D result, Q3D stator exit Ma is lower, and lower rotor blade zone Ma radial gradient can be seen in Q3D contour. Lower stator exit Ma of Q3D result implies lower degree of flow expansion in stator, which means smaller throat area in stator. As a consequence, Q3D ṁ result is higher than 3D. Figure 17 shows that Q3D Mar near rotor trailing edge is lower than 3D result, and rotor Mar axial gradient is higher in 3D result. Figure 18 shows that Q3D p∗ radial gradient in stator is higher than 3D result. Finally in Figure 19, the major error of Q3D T∗ lies in rotor tip area where Q3D T∗ is lower than 3D result.Figure 16
Absolute Mach number of the single-stage turbine: (a) Q3D; (b) 3D.Figure 17
Relative Mach number of the single-stage turbine: (a) Q3D; (b) 3D.Figure 18
Total pressure of the single-stage turbine: (a) Q3D; (b) 3D.Figure 19
Total temperature of the single-stage turbine: (a) Q3D; (b) 3D.The radial distribution ofvθr at blade row exit is shown in Figure 20. Figure 20(a) shows that Q3D stator exit vθr features a higher value at midspan. The absolute value of difference between Q3D vθr and vθr∗ is within 5 m2/s. For 3D stator result, vθr is close to the target from 0.1 to 0.7 span. At 0.9 span, difference from vθr∗ is nearly 6 m2/s. In Figure 20(b), both Q3D and 3D vθr results are close to the target value from 0.1 to 0.8 span. However, their vθr values at rotor tip zone are much lower than vθr∗. This may indicate that tip clearance leakage loss makes it impossible to realize the swirl target at tip.Figure 20
Swirl of the single-stage turbine: (a) stator; (b) rotor.
(a)(b)Statorα1 and rotor β2 are shown in Figure 21. In Figure 21, both stator α1 and rotor β2 of Q3D results exhibit high accuracy. From 0.1 to 0.9 span, the absolute value of Q3D stator α1 error is within 0.3°, and the absolute value of rotor β2 error is no more than 2°. Considering the Mach number difference between Q3D and 3D contours, it can be speculated that the error of Q3D result is mainly caused by the error of speed magnitude.Figure 21
Flow angle of the single-stage turbine: (a) stator flow angle; (b) rotor relative flow angle.
(a)(b)StatorY and rotor Yr are shown in Figure 22. In Figure 22(a), 3D stator Y result remains nearly constant along 0.1 to 0.9 span. However, Q3D stator Y essentially has a quadric distribution radially. This is caused by the secondary loss model in Q3D method. For the single-stage turbine, its stator tip endwall has an S shape profile, which is beneficial for reducing secondary flow strength. Therefore, radial migration of loss in 3D result is not remarkable. At the same time, stator hub boundary flow separation loss in 3D result is also depressed with high loss range is restricted below 0.1 span. In Figure 22(b), radial loss migration can be seen at rotor hub zone in 3D result, but not very strong. Both stator and rotor losses of Q3D results are lower than that of 3D results at midspan, but higher around endwalls. In the future, a more sophisticated loss migration model should be considered to make loss redistribution more accurate.Figure 22
Total pressure loss of the single-stage turbine: (a) stator total pressure loss; (b) rotor relative total pressure loss.
(a)(b)ṁ and ηi characteristic curves at design N are exhibited in Figure 23. The ṁ performance of Q3D result is generally close to 3D. For ηi characteristic curve, the change in Q3D ηi result is stronger when π is above 3. The cause of this error should be the overestimated loss by the loss model in the Q3D method.Figure 23
Characteristic curves of the single-stage turbine: (a) mass flow rate; (b) efficiency.
(a)(b)


## 5. Conclusions

A quasi 3D axial-flow turbine design method based on the Euler time-marching approach is introduced here, and a corresponding person computer code is developed. The design method mainly features a novel S2 stream surface approximation method. A two-stage turbine and a single-stage high-pressure turbine were designed by the Q3D design method and simulated through 3D CFD. After comparison, conclusions can be summarized as follows:(1)
The design method is efficient and can finish a design case within 1 min. Compared with 3D simulation results at design point, mass flow rate error is no more than 1.25%, expansion ratio error is no less than -0.022, and isentropic efficiency error is no more than 0.37 percentage points(2)
The circumferential coordinateθ of S2 stream surface in blade zone can be approximated by quadratic function of axial coordinate x in both design and analysis applications. The overall performances and meridian aerodynamic parameters can basically be predicted accurately in this way during inverse design phase(3)
Blade exit flow angle has a direct impact on exit swirl, and generally a higher flow angle produces a higher swirl. Radial swirl distribution can be controlled by flow angle on S2 stream surface in Q3D design work(4)
Overall blade row loss can be calculated accurately by the energy loss model. Quadratic redistribution of secondary and cubic redistribution of tip clearance leakage loss above 0.7 span can generally reflect the loss reduction trend at midspan(5)
Loss prediction near flow path endwalls is more complex for Q3D design. In the future, radial loss migration caused by secondary flow and boundary layer separation should be taken into consideration



---

*Source: 4102723-2023-01-17.xml*