2.1 Transient performance simulation method

Transient performance simulation is employed to assess gas turbine engine dynamic behaviour from one steady state to another. During transient process, the continuity of mass flow between engine components and the conservation of work on the same shaft are broken. The key impact factors on transient performance are rotor inertia, inter-component volumes and fuel control schedule.

To consider the effect of rotor inertia, rotor dynamics is applied where the unbalanced shaft power SP will drive the engine to either accelerate or decelerate based on rotor dynamics presented by Equation (1)

(1) \begin{align}SP = TW - CW = \frac{{dN}}{{dt}}{\,\rm{*}\,}I{\,\rm{*}\,}N{\,\rm{*}\,}{\left( {\frac{\pi }{{30}}} \right)^2}\end{align}

where TW is the turbine power and CW is the compressor power. The shaft rotational speed at the transient process can be determined by Euler’s method shown by Equation (2)

(2) \begin{align}N({t + \Delta t}) = \frac{{dN}}{{dt}}{\,\rm{*}\,}\Delta t + N(t)\end{align}

The inter-component volume effect may be dealt with by two methods in gas turbine transient performance simulations, one called Continuous Mass Flow (CMF) method [Reference Larrowe, Spencer and Tribus14] and the other called Inter-Component Volume (ICV) method [Reference Fawke and Saravanamuttoo15]. In this study, an improved CMF method called the Imbalanced Mass Flow (IMF) [Reference Palmer and Yan16] method is used to capture the volume packing effect. The imbalanced mass flow term (IMF) is calculated by volume dynamics, and its effect is applied to the corresponding component for mass flow rebalancing during the iterative calculation process. Such calculation is given by Equations (3) and (4).

(3) \begin{align}IMF = {W_1} - {W_2}\end{align}

(4) \begin{align}IMF = \frac{{{\rm{dm}}}}{{dt}} = \frac{{dP}}{{dt}}{\,\rm{*}\,}\frac{V}{{RT}}\end{align}

where W ₁ is in-coming flow rate and W ₂ is out-going flow rate of an inter-component volume. Due to the limited space in the paper, no further details of the IMF method are presented, and the interested readers may refer to Ref. (Reference Palmer and Yan16).

In addition, fuel control is also very important for transient performance simulations. In this paper, a simple fuel schedule where fuel flow varies with the corrected shaft speed of the engine is applied.

2.2 Geometry model

A geometry model is normally built on exact geometry of the key gas path components concerned, such as compressors, combustors and turbines of gas turbine engines. However, such information is not available at engine conceptual and preliminary design stages and therefore heat soakage and tip clearance variations have to be ignored in transient performance simulations in the past. In this study, a simplified geometry model is set up based on design point engine performance parameters and basic component design criteria to establish the necessary geometrical data of major gas path components required for transient performance simulations with the aim of capturing major heat soakage and tip clearance variation effects.

2.2.1 Axial compressors and turbines

Figure 1 represents the simplified geometry models for generic axial compressors and turbines. It is assumed that both of them consist of three parts: blades, casing and discs with constant external diameter ${D_t}$ . It is also assumed that the basic aerothermal design criterion for compressors and turbines are given in Table 1 [Reference Ramsden17].

Figure 1. Generic turbomachine geometry.

Table 1. Basic turbomachinery design criteria

The outer annulus diameter ${D_t}$ is estimated by Equation (5) using the shaft rational speed RPM at the design point:

(5) \begin{align}{D_t} = \frac{{U{B_{t.max}}}}{{RPM{\,\rm{*}\,}\pi {\,\rm{*}\,}60}}\end{align}

The flow annulus cross area ${A_{ann}}$ is calculated by Equation (6)

(6) \begin{align}{A_{ann}} = \frac{{W\sqrt T }}{{PQ}}\end{align}

where mass flow function Q is given by Equation (7)

(7) \begin{align}Q = MA{\,\rm{*}\,}\sqrt {\frac{\gamma }{R}} {\,\rm{*}\,}{\left( {1 + \frac{{\gamma - 1}}{2}{\,\rm{*}\,}{M^2}} \right)^{\left( { - \frac{{\frac{{\gamma + 1}}{2}}}{{\gamma - 1}}} \right)}}\end{align}

The hub diameter ${D_{\rm{h}}}$ of the turbomachinery component can be estimated by Equation (8)

(8) \begin{align}{D_h} = \sqrt {\frac{{{A_t} - {A_{ann}}}}{\pi }} {\,\rm{*}\,}2\end{align}

The number of stages ${n_{C,{\rm{\;}}stage}}$ of a compressor is estimated by the ceiling function Equation (9) to estimate the least number of the stages required to achieve the required overall pressure ratio for compressors

(9) \begin{align}{n_{C,{\rm{\;}}stage}} = Ceiling\left[ {{{\log }_{1.35}}PR} \right]\end{align}

where PR is the compressor pressure ratio and the value of 1.35 is the average stage pressure ratio at design point representing the state-of-the-art compressor technology such as that of the GE9X with a total pressure ratio of 27 and stage number of 11. The number of stages ${n_{T,{\rm{\;}}stage}}$ of a turbine is estimated by Equation (10) to achieve the required total enthalpy drop for the turbine through stage loading and blade tip speed.

(10) \begin{align}{n_{T,{\rm{\;}}stge}} = Ceiling\left[ {\frac{{DH}}{{SL{\,\rm{*}\,}U{B_{t,max}}^2}}} \right]\end{align}

where DH is the total enthalpy drop; SL is the stage loading; $U{B_{t,max}}$ is blade tip speed. The length L of the turbomachinery component is estimated by Equation (11).

(11) \begin{align}L = {n_{{\rm{\;}}sta}}{\,\rm{*}\,}{L_{stage}}\end{align}

Due to the existence of the interface surface between disc and blade, the effective wetted area ${A_{d,w}}$ of a disc is estimated by Equation (12).

(12) \begin{align}{A_{d,w}} = \pi {\,\rm{*}\,}{D_{h,m}}{\,\rm{*}\,}L{\,\rm{*}\,}Coe\end{align}

Consequently, the total interface ${A_{bd}}$ between the blade and the disc is estimated by

(13) \begin{align}{A_{bd}} = \pi {\,\rm{*}\,}{D_{h,m}}{\,\rm{*}\,}L{\,\rm{*}\,}({1 - Coe})\end{align}

Then the number of the blades in the intermediate stage can be estimated by Equation (14).

(14) \begin{align}{n_{b,stage}} = 2{\,\rm{*}\,}\pi {\,\rm{*}\,}\frac{{{D_t} + {\rm{\;}}{D_{h,m}}}}{{SR{\,\rm{*}\,}{D_t}{\rm{\;}} - {\rm{\;}}{D_{h,m}}/AR}}\end{align}

It is assumed that the number of blades is the same in each stage. Therefore, the total effective wetted area of the blades in the primary flow path is estimated by Equation (15).

(15) \begin{align}{A_{b,w}} = {n_b}{\,\rm{*}\,}{n_{stage}}{\,\rm{*}\,}\frac{{{{\left( {{D_t} - {D_{h,m}}} \right)}^2}}}{{AR}}\end{align}

The total volume of the blades ${V_b}$ and the casing ${V_c}$ can be represented by the corresponding surface area $A$ multiplied by the thickness $H$ by Equation (16).

(16) \begin{align}V = A{\,\rm{*}\,}H\end{align}

The disc volume ${V_d}$ is calculated by the mean cross-section area multiplied by the length as shown in Equation (17) where the space inside the disc is considered by coefficient 0.5.

(17) \begin{align}{V_d} = \pi {\,\rm{*}\,}{\left( {\frac{{{D_{h,m}}}}{2}} \right)^2}{\,\rm{*}\,}L{\,\rm{*}\,}0.5\end{align}

The metal mass of each component is estimated by Equation (18).

(18) \begin{align}{\rm{M}} = V*\rho \end{align}

As the shapes of the disc and blades are non-uniform in axial direction, the mass function per unit length M(x) of them are considered. For compressor disc, it is assumed that its inlet diameter is equal to half of its hub mean diameter and the outlet diameter is equal to one and a half of the hub mean diameter. Therefore, the linear relationship of the compressor disc diameter ${D_d}(x)$ could be expressed by Equation (19).

(19) \begin{align}{D_d}(x) = \frac{{{D_{h,m}}}}{L}x + \frac{{{D_{h,m}}}}{2}\end{align}

By substituting Equation (19) into Equations (17) and (18), the mass function of the compressor disc is given by Equation (20).

(20) \begin{align}{{\rm{M}}_d}\!(x) = \pi {\,\rm{*}\,}{\left({\frac{{{D_{h,m}}}}{{2L}}x + \frac{{{D_{h,m}}}}{4}}\right)^2}{\,\rm{*}\,}0.5{\,\rm{*}\,}\rho \end{align}

For the compressor blades, the mass distribution in axial direction could be treated the same way as that of the discs. As the tip diameter is constant, the cross-section area of the blade rows could be presented by Equation (21).

(21) \begin{align}{A_b}(x) = \pi {\,\rm{*}\,}\left[ {{{\left( {\frac{{{D_t}}}{2}} \right)}^2} - {{\left( {\frac{{{D_{h,m}}}}{{2L}}x + \frac{{{D_{h,m}}}}{4}} \right)}^2}} \right]\end{align}

Thus, the mass function of the compressor blades ${M_b}(x)$ can be represented by Equation (22).

(22) \begin{align}{{\rm{M}}_b}(x) = \pi {\,\rm{*}\,}\left[ {{{\left( {\frac{{{D_t}}}{2}} \right)}^2} - {{\left( {\frac{{{D_{h,m}}}}{{2L}}x + \frac{{{D_{h,m}}}}{4}} \right)}^2}} \right]{\,\rm{*}\,}\rho \end{align}

According to Equation (13), the interface area function ${A_{bd}}(x)$ between the blades and the disc is estimated by Equation (23).

(23) \begin{align}{A_{bd}}(x) = \pi {\,\rm{*}\,}\left( {\frac{{{D_{h,m}}}}{{}}x + \frac{{_{h,m}}}{2}} \right){\,\rm{*}\,}({1 - Coe})\end{align}

Similarly, it is also assumed that the inlet diameter of the turbine disc is equal to one and a half of its hub mean diameter and the outlet diameter of the turbine disc is equal to half of the hub mean diameter. Therefore, the linear variation of the turbine disc diameter ${D_d}(x)$ along axial direction could be expressed by Equation (24)

(24) \begin{align}{D_d}(x) = - \frac{{{D_{h,m}}}}{L}x + \frac{{3{D_{h,m}}}}{2}\end{align}

Consequently, the mass function of the turbine disc and blades are represented by Equations (25) and (26), respectively.

(25) \begin{align}{{\rm{M}}_d}(x) = \pi {\,\rm{*}\,}{\left( { - \frac{{{D_{h,m}}}}{{2L}}x + \frac{{3{D_{h,m}}}}{4}} \right)^2}{\,\rm{*}\,}0.5{\,\rm{*}\,}\rho \end{align}

(26) \begin{align}{{\rm{M}}_b}(x) = \pi {\,\rm{*}\,}\left[ {{{\left( {\frac{{{D_t}}}{2}} \right)}^2} - {{\left( { - \frac{{{D_{h,m}}}}{{2L}}x + \frac{{3{D_{h,m}}}}{4}} \right)}^2}} \right]{\,\rm{*}\,}\rho \;\end{align}

As the casing of the compressor and turbine has a constant diameter and could be regarded as being uniform, the casing mass function in axial direction could be given by Equation (27).

(27) \begin{align}{{\rm{M}}_c}(x) = \frac{{{{\rm{M}}_{\rm{c}}}}}{L}\end{align}

2.2.2 Combustor

The simplified combustor is assumed to be of annular type with double casing shown in Fig. 2. The whole combustor is divided into two parts: a diffuser and a primary zone.

Figure 2. Generically simplified combustor.

The combustor reference area at combustor inlet ${A_{ref}}$ is represented by Equation (28) [Reference Priyant Mark and Selwyn18] where the pressure drop factor $\frac{{\Delta {P_{3 - 4}}}}{{{q_{ref}}}}$ for annular combustor is assumed to be 20, and the corresponding pressure loss factor $\frac{{\Delta {P_{3 - 4}}}}{{{P_3}}}$ is assumed to be 0.06. RG is the universal gas constant.

(28) \begin{align}{A_{RB}} = {\left[ {{{RG} \over 2}{{\left( {{{{{\dot m}_3}T_3^{0.5}} \over {{P_3}}}} \right)}^2}{{\Delta {P_{3 - 4}}} \over {dp}}{{\left( {{{\Delta {P_{3 - 4}}} \over {{P_3}}}} \right)}^{ - 1}}} \right]^{0.5}} \end{align}

The distance $dR$ between the inner casing and the outer casing is estimated by Equation (29)

(29) \begin{align}dR = 0.2{\,\rm{*}\,}{R_{h,out}}\end{align}

where ${R_{h,out}}$ is the upstream compressor rear stage hub diameter. The height of the combustor inlet is calculated by Equation (30).

(30) \begin{align}{H_{RB}} = \sqrt {\frac{{{A_{RB}} + \pi {\,\rm{*}\,}{R_h}^2}}{\pi }} - {R_h}\end{align}

Hence, the height of the combustor liner ${H_L}$ is estimated by Equation (31).

(31) \begin{align}{H_L} = {H_{RB}} - 2{\,\rm{*}\,}dR\end{align}

As the modern annular combustor is designed short axially, the length of the annular combustor ${L_B}$ is roughly two times the height of its liner ${H_L}$ . Hence, the total length of the combustor is given by Equation (32).

(32) \begin{align}{L_B} = {H_L}{\,\rm{*}\,}2\end{align}

It is assumed that the length of the combustor diffuser ${L_D}$ is half of that of the liner length so ${L_D}$ is calculated by Equation (33).

(33) \begin{align}{L_B} = {H_{}}*0.5\end{align}

It is assumed that it has double casing configuration, so the total wetted area of the combustor ${A_{B,w}}$ is the sum of the wetted area of the inner casing ${A_{B,inner,w}}$ and the outer casing ${A_{B,outer,w}}$ as shown by Equation (34) where the wetted areas could be estimated by the multiplication between corresponding perimeter and length.

(34) \begin{align}{A_{B,w}}{\rm{\;}} = {A_{B,inner,w}} + {A_{B,outter,w}}\end{align}

It is also assumed that the combustor casing thickness ${H_B}$ is uniform, and therefore the volume of the combustor casing metal ${V_B}$ is calculated by the sum of the volumes of the inner casing ${V_{B,inner}}$ and the outer casing ${V_{B,outer}}$ shown by Equation (35).

(35) \begin{align}{V_B}\; = {V_{B,inner}} + {V_{B,outer}} = ({A_{B,inner,w}} + {A_{B,outter,w}})*{H_B}\end{align}

Finally, the metal mass of each part of the combustor could be determined by Equation (18).

As the combustor could be regarded as a duct, its mass function could be expressed by Equation (27) as well. Similarly, the combustor area function could be presented by Equation (36)

(36) \begin{align}{A_{B,w}}(x){\rm{\;}} = \frac{{{A_{B,inner,w}} + {A_{B,ou,w}}}}{L}\end{align}

2.3 Heat transfer model

Assumptions for the heat transfer model are made as followings:

• Conductive heat transfer is only considered between air/gas and component metals in radial direction.

• No conductive heat transfer between the components.

• No heat loss to ambient environment.

• All temperature profiles within the components are the same in radial direction and vary linearly along the axial direction.

• The metal mass profile functions defined as the mass per unit axial length vary linearly along the axial direction.

• The wetted area profile functions of all components defined as the area per unit axial length are linear along the axial direction.

• The heat transfer coefficients are linear along the axial direction.

• No heat loss in air system or from cooling flows.

The heat transfer models developed in this study is to capture the thermal effects of the following:

1. 1. Convective heat transfer between the gas flow in the main gas path and the turbomachinery and the combustor metals.

2. 2. Radiative heat transfer between the gas flow and the combustor metals.

3. 3. Conductive heat transfer between the blades and the discs in the turbomachinery components.

2.3.1 Turbomachinery heat transfer model

The left side of Fig. 3 shows a schematic of the turbomachinery components, including casing, blades and discs, and the right side of the figure shows a schematic of the heat transfer between them where convective and conductive heat transfers are considered.

Figure 3. Heat transfer model for turbomachinery.

2.3.1.1 2.3.1.1 Convective heat transfer

The convective heat transfer between air/gas flow and component metals is calculated based on the lumped capacitance method by Equation (37) with time constant.

(37) \begin{align}\frac{{dq}}{{dt}} = {{\rm{M}}_s}{\,\rm{*}\,}Cp{\,\rm{*}\,}\left( {{T_f} - {T_s}} \right)\left( {{e^{ - \frac{{\Delta t}}{\tau }}} - 1} \right)\end{align}

(38) \begin{align}\tau = {{\rm{M}}_s}{\,\rm{*}\,}\frac{{Cp}}{{h{\,\rm{*}\,}A}}\end{align}

The heat transfer coefficients used in this paper refer to Ref. (Reference Naylor2). For forced convection of flow through a duct,

(39) \begin{align}h = 0.36{\rm{\;*\;R}}{{\rm{e}}^{0.8}}P{r^{\frac{1}{3}}}{\left( {\frac{D}{L}} \right)^{\frac{1}{{18}}}}{\left( {\frac{{Visc}}{{Viscw}}} \right)^{0.14}}\frac{K}{D}{\rm{\;\;}}if{\rm{\;}}({\rm{Re}} \gt 2,100)\end{align}

(40) \begin{align}h = 1.86{\rm{\;*\;R}}{{\rm{e}}^{\frac{1}{3}}}P{r^{\frac{1}{3}}}{\left( {\frac{D}{L}} \right)^{\frac{1}{3}}}{\left( {\frac{{Visc}}{{Viscw}}} \right)^{0.14}}\frac{K}{D}{\rm{\;\;}}if{\rm{\;}}\left( {{\rm{Re}} \le 2,100} \right)\end{align}

For a forced convection of flow through a rotating disc,

(41) \begin{align}h = 0.0267{\,\rm{*}\,}\frac{{{\rm{R}}{{\rm{e}}^{0.6}}P{r^{0.8}}K}}{R}{\rm{\;}}if{\rm{\;}}\left( {Re \le 24,000} \right)\end{align}

(42) \begin{align}h = 0.616{\,\rm{*}\,}\frac{{{\rm{R}}{{\rm{e}}^{0.5}}P{r^{\frac{1}{3}}}K}}{R}{\rm{\;}}if{\rm{\;}}(Re \gt 24,000)\end{align}

where the Reynold number ${\rm{Re}}$ for a rotating disc and a duct, and Prandtl number $Pr$ are estimated by Equations (43) to (45),

(43) \begin{align}{\rm{Re}} = ({\rho {\,\rm{*}\,}\omega {\,\rm{*}\,}{R^2}})/Visc\end{align}

(44) \begin{align}{\rm{Re}} = W{\,\rm{*}\,}D/({{\,\rm{*}\,}Visc})\end{align}

(45) \begin{align}Pr = Cp{\,\rm{*}\,}Visc/K\end{align}

where $D$ is the hydraulic diameter of a tube, $L$ the length of the duct, $K$ the thermal diffusivity of the gas, ${V_{isc}}$ the bulk viscosity of the gas, ${V_{iscw}}$ the viscosity of the gas at the wall temperature, and $R$ the radius of the disc.

The linear temperature profiles ${T_{ref}}(x)$ in axial direction for the gas flow, the blades, the casing and the discs, are represented by Equation (46).

(46) \begin{align}{T_{ref}}(x) = \frac{{{T_{ref,{\rm{\;}}out}} - {\rm{\;}}{T_{ref,n}}}}{L}{\,\rm{*}\,}x + {T_{ref,in}}\end{align}

where the subscrip “ref” refers to each component concerned. The axial heat transfer coefficient profiles for component metals $h(x)$ is represented by Equation (47).

(47) \begin{align}h(x) = {h_{in}} + \frac{{{h_{out}} - {h_{in}}}}{L}{\,\rm{*}\,}x\end{align}

Substituting Equation (46) into Equation (37) and applying integration yields the total convective heat transfer between the air/gas flow and each part of component metal represented by Equation (48).

(48) \begin{align}\mathop \int \nolimits_0^L {q_{conv}}dx = \mathop \int \nolimits_0^L {M_{ref}}(x)*Cp*\left( {\frac{{{T_{f,out}} - \;{T_{in}}}}{L} - \frac{{{T_{ref,out}} - \;{T_{in}}}}{L}} \right)*\left( {{e^{ - \frac{{\Delta}}{\tau }}} - 1} \right)dx\end{align}

Where ${T_{f,out}}$ represents the outlet temperature of the air/gas flow, ${T_{ref,out}}$ the outlet temperature of the blades, the casing or the discs, and ${{\rm{M}}_{ref}}$ (x) is the metal mass profile functions defined as the mass per unit axial length and the subscript “ref” refers to the blades (“b”), the casing (“c”) or the discs (“d”) that is assumed to be linear and can be determined based on the designs in Section 2.2.

2.3.1.2 2.3.1.2 Conductive heat transfer

Figure 4 represents the conductive heat transfer between the blades and the discs estimated by the lumped capacitance method shown by Equation (49)

(49) \begin{align}\frac{{dq}}{{dt}} = {U_{total}}{\,\rm{*}\,}\left( {{T_b} - {T_d}} \right)\left( {{e^{ - \frac{{\Delta t}}{{Tau}}}} - 1} \right)\end{align}

(50) \begin{align}{U_{total}} = \frac{1}{{\dfrac{1}{{\dfrac{{{k_b}{\,\rm{*}\,}{A_{bd}}}}{{{X_b}}}}} + \dfrac{1}{{\dfrac{{{k_d}{\,\rm{*}\,}{A_{bd}}}}{{{X_d}}}}}}}\end{align}

(51) \begin{align}Tau = \frac{{{{\rm{M}}_b}*C{p_b} + {{\rm{M}}_d}*C{p_d}}}{{\dfrac{{{k_b}*{A_{bd}}}}{{{X_b}}} + \dfrac{{{k_d}*{A_{bd}}}}{{{X_d}}}}}\end{align}

where ${U_{total}}$ is the total thermal conductance of the blades and the discs, ${k_b}$ is the conductivity of the blades, ${k_d}$ is the conductivity of the discs, ${X_b}$ is the distance between the centre of the blades and the interface, ${X_d}$ is the distance between the centre of the discs and the interface, ${A_{bd}}$ is the interface area between the blades and the discs, and ${T_{au}}$ is the time constant of the conductive heat transfer.

Figure 4. Conductive heat transfer between blades and discs.

Substituting the temperature profiles of the blades ${T_b}(x)$ and the discs ${T_d}(x)$ and the area profile function of the interface ${A_{bd}}(x)$ into Equation (49) derives the total conductive heat transfer calculated by Equation (52).

(52) \begin{align}\mathop \int \nolimits_0^L {q_{cnd}}dx = \mathop \int \nolimits_0^L {A_{bd}}(x)*\left( {\frac{{{T_{b,out}} - \;{T_{in}}}}{L} - \frac{{{T_{d,out}} - \;{T_{in}}}}{L}} \right)*\frac{{{k_b}{k_d}}}{{{k_b}{X_d} + {k_d}{X_b}}}*\left( {{e^{ - \frac{{\Delta t}}{{Tau}}}} - 1} \right)dx\end{align}

where ${A_{bd}}$ (x) is the area profile functions defined as the area per unit axial length that is assumed to be linear and can be determined based on the designs in Section 2.2.

For the turbomachinery system, the heat transfer happens among four subsystems, which are the air/gas flow, the blades, the casing and the discs. Based on energy conservation, the temperature variation of the air/gas flow over time is determined by the sum of the convection between the air/gas flow and the blades, the casing and the disc, represented by Equation (53).

(53) \begin{align}\sum \mathop \int \nolimits_0^L ({q_{Conv,b}} + {q_{Conv,c}} + {q_{Conv,d}})dx = \mathop \int \nolimits_0^L {W_f}*\Delta t*C{p_f}*\left[ {{C_f}({t + \Delta t}) - {C_f}(t)} \right]xdx\end{align}

The temperature variation of the blades over time is determined by the sum of the convective heat ${q_{conv,b}}$ between the air/gas flow and the blade metal and the conductive heat ${q_{cond}}$ between the blade and the disc metals represented by Equation (54).

(54) \begin{align}\sum \mathop \int \nolimits_0^L ({q_{conv,b}} + {q_{cond}})dx = \mathop \int \nolimits_0^L {{\rm{M}}_b}(x)*Cp*\left[ {{C_b}({t + \Delta t}) - {C_b}(t)} \right]xdx\end{align}

In a similar way, the temperature variation of the discs over time is determined by the sum of the convective heat ${q_{conv,d}}$ between the air/gas flow and the disc metal and the conductive heat ${q_{cond}}$ between the blade and the disc metals represented by Equation (55).

(55) \begin{align}\sum \mathop \int \nolimits_0^L ({q_{conv,d}} + {q_{cond}})dx = \mathop \int \nolimits_0^L {{\rm{M}}_d}(x)*Cp*\left[ {{C_d}({t + \Delta t}) - {C_d}(t)} \right]xdx\end{align}

Finally, the temperature variation of the casing over time is determined by the convective heat ${q_{conv,d}}$ between the air/gas flow and the casing metal shown by Equation (56).

(56) \begin{align}\sum \mathop \int \nolimits_0^L {q_{Conv,c}}dx = \mathop \int \nolimits_0^L {{\rm{M}}_c}(x){\,\rm{*}\,}Cp{\,\rm{*}\,}\left[ {{C_c}({t + \Delta t}) - {C_c}(t)} \right]xdx\end{align}

Hence, by solving the above four Equations (52) to (56), the four temperature gradients ${\rm{\;}}{C_{ref}}$ where the subscrip ref refers to the air/gas flow, the blades, the casings and the discs, respectively, can be obtained. Correspondingly, the exit temperatures of the air/gas flow, the blades, the casings and the discs can be updated by Equation (57).

(57) \begin{align}{T_{ref,\;\;out}}({t + \Delta t}) = {C_{ref}}({t + \Delta t})*L + {T_{ref,in}}({t + \Delta t})\end{align}

(58) \begin{align}{C_{ref}} = \frac{{{T_{ref,\;out}} - \;{T_{ref,in}}}}{L}\end{align}

where ${C_{ref}}$ is temperature gradient where the ${\rm{\;}}$ subscript “ref” refers to the air/gas flow (“f”), the blades (“b”), the discs (“d”) and the casing (“c”).

2.3.2 Combustor heat transfer

Due to high gas temperature inside the combustor, the radiative heat transfer cannot be ignored. Therefore, the heat transfer model for combustors mainly considers convection and radiation. Consequently, it is reasonable to split the combustor into two sections, i.e. diffusion and primary zones due to significant temperature difference between the two parts.

It is assumed that the metal components in the combustor do not have conductive heat transfer between each other. The heat transfer in both the diffusor and the primary zone is mainly convection and radiation. The convective heat transfers for the diffusor ${q_{conv,\;D}}$ and the primary zone ${q_{conv,\;P}}$ could be calculated using Equation (48) where the subscript “ref” refers to either the diffusor or the primary zone.

The radiative heat transfer between the combustor metal and the gas flow at any axial position within the combustor is estimated by Equation (59)

(59) \begin{align}{q_{rad}} = \sigma {\,\rm{*}\,}\left( {\frac{{1 + {\varepsilon _w}}}{2}} \right){\,\rm{*}\,}{A_B}{\,\rm{*}\,}\left( {{\varepsilon _g}T_f^4 - {\alpha _g}T_B^4} \right)\left( {{e^{ - \frac{{\Delta t}}{\tau }}} - 1} \right)\end{align}

where ${\varepsilon _w}$ and ${\varepsilon _g}$ are radiative emissivity of the wall and the gas, respectively. ${\alpha _g}$ is gas absorptivity and $\sigma $ the Stefan-Boltzmann constant. The total radiative heat transfer ${q_{rad}}$ for both the diffusor ${q_{rad,\;\;}}$ and the primary zone ${q_{rad,\;P}}$ can be calculated by Equation (60).

(60) \begin{align}\mathop \int \nolimits_0^{{L_{ref}}} {q_{rad,ref}}{\rm{\;}}dx & = {\mathop \int \nolimits_0^{{L_{ref}}} {ref}}(x){\,\rm{*}\,}\sigma {\,\rm{*}\,}\left( {\frac{{1 + {\varepsilon _w}}}{2}} \right)\left[ {\varepsilon _g}\left( {{T_{f,in}} + \frac{{{T_{f,out - }}{T_{f,in}}}}{{{L_{ref}}}}x} \right)_{ref}^4 \right. \nonumber \\ & \quad \left. - {\alpha _g}\left( {{T_{B,in}} + \frac{{{T_{B,out - }}{T_{B,in}}}}{{{L_{ref}}}}x} \right)_{ref}^4 \right] \left( {{e^{ - \frac{{\Delta t}}{{{\tau _B}}}}} - 1} \right)dx\end{align}

where ${T_{f,in}}$ and ${T_{f,out}}$ are the inlet and outlet gas temperature in the diffuser or the primary zone, respectively, ${T_{B,in}}$ and ${T_{B,out}}$ are the inlet and outlet metal temperature of the diffuser or the primary zone, respectively, and ${A_{ref}}(x)$ is the surface area profile function defined as the area per unit axial length. ${A_{ref}}(x)$ is assumed to be linear and can be determined based on the designs in Section 2.2.

According to the heat transfer model of the combustor system, the heat transfer process could be divided into three sub-systems, which are the heat transfer to or from the air/gas flow, the heat transfer to or from the diffuser metal and the heat transfer to or from the primary zone metal, respectively. Based on the conservation of energy, the total heat energy transferred to or from the air/gas flow is equal to the accumulative convective heat energy and the radiative heat energy between the air/gas and the combustor metals, as shown by Equation (61).

(61) \begin{align}\sum \mathop \int \nolimits_0^{{L_B}} ({q_{conv}}{\rm{\;}} + {q_{rad}})dx = \mathop \int \nolimits_0^{{L_B}} {W_f}{\,\rm{*}\,}\Delta t{\,\rm{*}\,}C{p_f}{\,\rm{*}\,}\left[ {{C_f}({t + \Delta t}) - {C_f}(t)} \right]dx\end{align}

The total energy transferred to or from the diffuser metal is equal to the sum of convective energy ${q_{conv,D}}$ and radiative energy ${q_{rad,D}}$ , as shown by Equation (62).

(62) \begin{align}\sum \mathop \int \nolimits_0^{{L_D}} ({q_{conv,D}} + {q_{rad,D}})\;dx = \mathop \int \nolimits_0^{ref} {{\rm{M}}_D}(x)*C{p_B}*\{ {C_D}({t + \Delta t}) - {C_D}(t)]dx\end{align}

Similarly, total energy transferred to or from the primary zone metal is equal to the sum of the convective energy ${q_{conv,P}}$ and radiative energy ${q_{rad,P}}$ , as shown by Equation (63).

(63) \begin{align}\sum \mathop \int \nolimits_0^{{L_P}} ({q_{conv,P}} + {q_{rad,P}}){\rm{\;}}dx = \mathop \int \nolimits_0^{{L_P}} {M_P}{\,\rm{*}\,}C{p_B}{\,\rm{*}\,}\{ {C_P}({t + \Delta t}) - {C_P}(t)]dx\end{align}

where ${L_D}$ and ${L_P}$ stand for the length of the diffuser and the primary zone, respectively, and ${M_D}{\rm{\;}}({\rm{x}}){\rm{\;and\;}}{M_P}\left( {\rm{x}} \right){\rm{\;}}$ represent the metal mass of the diffuser and the primary zone respectively. By solving Equations (61) to (63), the temperature gradients for the gas flow ${C_f}$ , the diffusor metal ${C_D}$ and the primary zone metal ${C_P}$ can be obtained.

Once the temperature gradients for the three sub-systems become available, Equation (57) is applied to update the outlet temperature of the gas flow ${T_{f,\;out}}$ , the diffuser metal ${T_{D,\;out}}$ and the primary zone metal ${T_{P,\;out}}$ by substituting the corresponding values of the combustor.

2.4 Tip clearance model

Compressor and turbine tip clearances could vary due to centrifugal force and thermal expansion or contraction of casing, blades and discs during engine transient processes. Consequently, the changing tip clearances could change the characteristics of the turbomachinery components.

The tip clearance model is divided into two parts: the tip clearance variation estimation and map correction. It is assumed that the tip clearance variation across a turbomachinery component is uniform and that only the expansion or contraction of blades, casing and discs in radial direction is considered.

2.4.1 Tip clearance variation

The thermal expansion or contraction Δ of blades, casing and discs due to changing temperature of the components can be estimated by Equation (64) [Reference Chapman13]

(64) \begin{align}\delta ({t + \Delta t}) = {\rm{\;}}\alpha {\,\rm{*}\,}H{\,\rm{*}\,}\left[ {T({t + \Delta t}) - T(t)} \right]\end{align}

where $\alpha $ is expansion rate and H the radius of the corresponding components.

The radial expansion or contraction of discs ${\Delta _d}$ and blades ${\Delta _b}$ due to centrifugal force are given by Equations (65) and (66) [Reference Kypyuros and Melcher19].

(65) \begin{align}{\Delta _d}({t + \Delta t}) = \frac{1}{{4{\,\rm{*}\,}{E_d}}}{\,\rm{*}\,}\left( {1 - {\upsilon _d}} \right){\,\rm{*}\,}{\rho _d}{\,\rm{*}\,}R_d^3{\,\rm{*}\,}{\left( {2{\,\rm{*}\,}P{\,\rm{*}\,}\frac{{N({t + \Delta t}) - N(t)}}{{60}}} \right)^2}\end{align}

(66) \begin{align}{\Delta _b}({t + \Delta t}) = {\rho _b}*H_b^2*{R_b}/{E_b}*{\left( {2*PI*\frac{{N({t + \Delta t}) - N(t)}}{{60}}} \right)^2}\end{align}

The overall growth $\vartheta $ of the blades, casing and discs could be estimated by the summation of their thermal and mechanical growth by Equation (67)

(67) \begin{align}\vartheta = {\delta _c} - ({\delta _b} + {\Delta _b}) - \left( {{\delta _d} + {\Delta _d}} \right)\end{align}

Finally, the changing tip clearance $TC$ is determined based on the tip clearance ${\rm{TC}}\left( {\rm{t}} \right)$ at time t and its variation $\vartheta $ ,

(68) \begin{align}TC({t + \Delta t}) = TC(t) + \vartheta \end{align}

2.4.2 Map correction

According to Ref. (Reference Pilidis and Maccallum7), the characteristic map of a compressor could change due to the impact of heat transfer and tip clearance variation. Such map change could be represented by a shift of the speed lines of the map estimated by Equation (69)

(69) \begin{align}\frac{{\Delta N}}{N} = C1{\,\rm{*}\,}\frac{{\left( {{T_b} - {T_f}} \right)}}{{{T_f}}} + C2{\,\rm{*}\,}\frac{{\dot Q}}{{{W_a}{\,\rm{*}\,}{C_P}{\,\rm{*}\,}{T_f}}} + C3{\,\rm{*}\,}\frac{\vartheta }{{TC}}\end{align}

where the first term on the right side of the equation represents the effect of flow boundary layer development, the second term represents the effect of heat soakage, and the last term represents the tip clearance variation effect, and C1, C2 and C3 are corresponding coefficients that are compressor specific. These coefficients obtained from tests for low- and high-pressure compressors of a particular engine [Reference Maccallun5] are shown in Table 2 and used in this study.

Table 2. Shaft speed correction coefficients [Reference Maccallun5]

Based on the estimated shift of the compressor shaft speed, the compressor map is corrected by scaling the map using the scaling factors calculated by Equations (70) to (72) referring to Ref. [Reference Li20].

(70) \begin{align}M{F_{Scl}} = \frac{{MF\left( {N + \Delta N} \right)}}{{MF\left( N \right)}}\end{align}

(71) \begin{align}P{R_{Scl}} = \frac{{PR\left( {N + \Delta N} \right) - 1}}{{PR\left( N \right) - 1}}\end{align}

(72) \begin{align}EF{F_{Sc}} = \frac{{EFF\left( {N + \Delta N{\rm{\;}}} \right)}}{{EFF\left( N \right)}}\end{align}

The explanation of the equations and their applications are as follows:

• A new shaft speed is obtained by estimating the change in shaft speed using Equation (69).

• Two sets of mass flow, pressure ratio and efficiency are obtained using the original map at shaft speed $N$ and $\left( {N + \Delta N} \right)$ , respectively.

• The scaling factors of mass flow, pressure ratio and efficiency are calculated by Equations (70) to (72) at the transient operating point represented by the shaft speed at $N$ .

• The map is corrected by applying the scaling factors.

• Finally, the engine transient performance is updated using the corrected map.

The same map correction approach is also applied to turbines and therefore no further details are repeated.

2.5 Integration of models

The transient performance simulation system considering heat soakage and the impact of tip clearance variation has been established and illustrated in Fig. 5. This architecture of the whole system clearly shows the relationship and communication between the performance model and the heat soakage and tip clearance models.

Figure 5. Transient performance simulation system with impact of heat soakage and tip clearance variation.

3.1 Engine model

The heat soakage and tip clearance models have been integrated with TurboCycle software and applied to a model twin-spool turbojet engine shown in Fig. 6. Table 3 gives the key performance specification of the model engine at sea level static ISA condition. The LP shaft corrected relative speed is set as the engine control parameter, or handle parameter to control the engine operating conditions.

Figure 6. Schematic of model twin-spool turbojet.

Table 3. Engine performance specification at sea level static ISA condition

3.2 Engine model verification

To verify the baseline engine model set up in TurboCycle the same engine model has been created in commercial software GasTurb [Reference Kurzke21] with the same model configuration, component maps, performance specification and fuel control schedule but without the inclusion of heat soakage and tip clearance models. Then the steady state and the transient performance of both engine models have been calculated and the results are compared with each other.

Regarding the verification of the steady state performance for the model, the part-load performance of the engine model indicated by different LP shaft rotational speed is simulated by both software packages and a comparison of the simulation results of selected performance parameters including net thrust (FN), fuel flow rate (WFE), turbine entry temperature (T41) and exhuast gas temperature (T5) are shown in Fig. 7. It can be seen that the TurboCycle model produces very close part-load performance to that of the GasTurb model.

Figure 7. Comparison of steady state performance between TurboCycle and GasTurb models (a)–(d).

The verification of the engine model on simulated transient performance is conducted with an acceleration process from idle (70%) to full power (100%) using the same fuel control schedule as shown in Fig. 8.

Figure 8. Fuel schedule during acceleration.

Figure 9 shows a comparison of the simulation results of selected performance parameters including the corrected relative LP shaft speed (XN1R), air mass flow rate (W2), turbine entry temperature (T41) and net thrust (FN) of both models during the transient process. It can be seen that the simulated transient performance results of the TurboCycle model are also very close to those of the GasTurb model. Therefore, it has proven that the TurboCycle model is accurate enough as a baseline model to investigate the heat soakage and tip clearance variation effects in this study.

Figure 9. Comparison of transient performance between TurboCycle and GasTurb models (a)–(d).

3.3 Transient performance affected by heat soakage and tip clearance variation

The heat soakage and tip clearance variation effects have been simulated by TurboCycle with the model twin-spool turbojet engine for the same acceleration process, i.e. from idle (70%) to full power (100%) using the same fuel schedule. The transient performance without the consideration of heat soakage and tip clearance variation is regarded as the baseline condition. To investigate and compare the heat soakage effect contributed by different components, the heat flux of individual components during the acceleration process are presented in Fig. 10.

Figure 10. Heat flux of individual components during acceleration.

It can be seen that at the beginning of the acceleration, the heat transfer in the combustor is the greatest due to the excessive fuelling, which can produce extremely high level of radiative and convective heat flux. The heat flux of HPT becomes dominant after three seconds because of the highest increasing temperature of gas flow passing the HPT and the third highest heat transfer coefficient shown in Figs 11 and 12. The heat flux of LPT is weaker than that of the burner and HPT but stronger than that of LPC and HPC because of the higher temperature in LPT. LPC has the lowest heat flux due to its lowest temperature and the lowest heat transfer coefficient as illustrated in Figs 11 and 12.

Figure 11. Temperature difference between metal and gas flow of individual components during acceleration.

Figure 12. Heat transfer coefficient of individual components during acceleration.

The simulation results of the transient performance with the inclusion of the tip clearance variations of the turbomachinery components are presented in Fig. 13. It can be observed that at the beginning of the transient process, the tip clearance drops significantly because of the quick expansion of the blades due to centrifugal force. Then the tip clearance increases because of the gradual expansion of the casing heated by high temperature gas flow. In addition, the tip clearance variations of the high-pressure spool components are slightly more obvious than those of the low spool components because of the shaft rotational speed difference.

Figure 13. Tip clearance variations in compressors and turbines (a)–(b).

Figure 14. Comparisons between transient performance with and without heat soakage and tip clearance variation effects (a)–(j).

To assess the impact of the heat soakage and tip clearance variation effects, the transient performance of the engine model in acceleration process with and without the heat soakage and tip clearance variation has been simulated using the same fuel schedule (Fig. 8). Figure 15 presents a comparison of the variations of selected performance parameters during the acceleration, including corrected relative LP and HP shaft speeds, net thrust, air flow rate, and pressure ratios of LPC, HPC, HPT and LPT, TET and EGT. It can be seen that the heat soakage and tip clearance variations would have significant impact on the acceleration time, and therefore would delay the response of the net thrust and other performance parameters.

Figure 15. Trajectories on compressor maps affected by heat soakage and tip clearance variation (a)–(b).

In more details, compared with the baseline condition where no heat soakage and tip clearance variation are considered, the heat soakage alone causes the following effects for the acceleration process:

• Delay the corrected relative low pressure shaft speed reaching the target value from 10s to 100s, Fig. 14(a). It therefore causes a delay of the increase of the air flow rate (Fig. 14(d)) and the LPC pressure ratio (Fig. 14(e)).

• The corrected relative high pressure shaft speed becomes slightly higher due to the lower inlet temperature to the HPC caused by the heat soakage of LPC, Fig. 14(b). Consequently, the pressure ratio of the HPC also becomes slightly higher as shown in Fig. 14(f).

• The increase of the net thrust (Fig. 14(c)) is delayed even at the end of the 100-second acceleration process, which is about 1.27% lower than that without heat soakage.

• The increase of TET (Fig. 14(i)) has been significantly delayed by about 100s. Consequently, the corresponding expansion ratios of both HPT and LPT (Fig. 14(g) and (h)) become higher due to the increase of the corrected relative shaft speeds of HPT and LPT, as shown in Fig. 14(g)–(h).

• The increase of EGT has also been delayed tremendously. As a result, the net thrust is reduced due to the lower exhaust velocity caused by the delayed EGT, as shown in Fig. 14(j).

Further inclusion of tip clearance variation effect above the heat soakage effect would result in the following changes:

• Delay the LP corrected relative shaft speed (Fig. 14(a)) reaching the target value from 10s to 48s. In addition, the final LP corrected relative shaft speed is slightly higher than the baseline value.

• The corrected relative high pressure shaft speed (Fig. 14(b) becomes lower than the baseline value because the higher low pressure shaft speed would result in the higher LPC exit temperature.

• The variation of tip clearance has little impact on the net thrust as shown in Fig. 14(c).

• The increase of the air flow rate and the LPC pressure ratio is higher than that with heat soakage due to further rise of the low-pressure shaft speed, Fig. 14(d)–(e).

• The pressure ratio of the HPC (Fig. 14(f)) becomes slightly lower than that with heat soakage and approaches the baseline value at the end of the acceleration.

• TET has a slight increase and approach its baseline value at the end of the acceleration, Fig. 14(i).

• The expansion ratios of HPT and LPT (Fig. 14(g) and (h)) have a slight increase compared to those with heat soakage only due to the slight increase of TET.

• EGT has a slight decrease and approaches its baseline value at the end of this acceleration, Fig. 14(j).

The acceleration trajectories of the LPC and the HPC at the baseline condition, with the heat soakage along, and with both the heat soakage and the tip clearance variation are shown in Fig. 15. It is obvious for the LPC that the heat soakage would lower the engine acceleration working line while adding the tip clearance variation would raise the acceleration working line slightly. The working line of the HPC is lowered by the heat soakage and tip clearance variation while the impact of adding the tip clearance variation on the working line is small.

A comparison of the most important performance parameters of the engine model at the end of the 100-second acceleration where the corrected low pressure shaft rotational speed reaches 100% with the heat soakage and the tip clearance variation compared with the baseline condition is summarised in Table 4. It can be seen that due to the heat soakage effect, TET has around -0.17% delay, EGT has about −1.61% delay, and the net thrust has about −1.3% delay. By adding extra tip clearance variation effect, TET increases by around 0.13%, EGT decreases by about 0.19%, and the net thrust increases by around 0.34%.

Table 4. Summary of impact of heat soakage and tip clearance variation at the end of 100 second of acceleration