Real-Time Gas Turbine Model for Performance Simulations

ABSTRACT

Industry and universities around the world invest time and money to develop digital computer programs to predict gas turbine performance. This study aims to demonstrate a brand new digital model developed with the ability to simulate gas turbine real time high fidelity performance. The model herein described run faster than 30ms per point, which is compatible with a high-definition video refresh rate: 30 frames per second. This user-friendly model, built in Visual Basic in modular structure, can be easily configured to simulate almost all the existing gas turbine architectures (single, 2 or 3 shaft engines mixed or unmixed flows). In addition, its real time capability enables simulations with the pilot in the loop at earlier design phases when their feedback may lead to design changes for improvements or corrections. In this paper, besides the model description, it is presented the model run time capability as well as a comparison of the simulated performance with a commercial gas turbine tool for single, 2 and 3 shaft engine architecture.

KEYWORDS:
Propulsion; Gas turbines; Aircraft engines; Performance; Computer simulation

INTRODUCTION

SAE AIR4548 defines a real-time digital engine model as a mathematical performance computer model whose outputs are generated at a rate compatible with the response of the physical system that it represents and with the time requirements of the simulation loop where it is inserted. The early developed models were relatively simple using analog devices and they were firstly used in hardware and software development for aircraft and engine control systems. As the model complexity increased to meet more demanding requirements, analog models became too costly and difficult to use. The early mathematical models, developed to make simulations less expensive, were simply a digital implementation of the analog models and, as digital computers capabilities increased and costs reduced, the engine digital models became very popular. As listed in Bringhenti (1999), some efforts in engine analog, digital or mixed simulation development can be acknowledged through the years notably by Mckinney (1967), Koenig and Fishback (1972), Fishback and Koenig (1972), Szuk (1974), Palmer and Yang (1974), Macmillan (1974), Sellers (1975), Wittenberg (1976), Flack (1990), Stamatis et al. (1990), Ismail and Bhinder (1991), Korakianitis and Wilson (1994), Baig and Saravanamuttoo (1997), Bringhenti (1999), Grönstedt (2000), Saravanamuttoo et al. (2001), Walsh and Fletcher (2004) and ASME 95-GT-147.

Nowadays, it is wide spread in the aeronautic industry the usage of simulation models for engine or aircraft development and its systems. Most of those models can run steady state simulations only and represents specific engine architecture; others are capable of simulating also the transient states with variable geometry and are flexible to represent almost all types of gas turbines as per Bringhenti (1999, 2003), Grönstedt (2000), and Silva (2011). However, due to the number of iterations and map data readings required in the gas turbine engine simulation process, a long time is required to output the simulation results, what is not compatible with a real-time application.

The real-time engine simulation tool can enhance the simulation activities at early phases of a product design, identifying potential improvements or issues at early phases of an aircraft design when there is room for changes or even step back at virtually no cost.

In addition, a digital real-time engine model could be used for development and testing of control systems, flight simulators, and engine integration with airframe in several aspects.

METHODOLOGY

A brand new engine model was generated to provide high fidelity and real-time gas turbine performance simulation. The model is representative of a three shaft engine, which is the most complex engine architecture. Other existing jet engines architectures (single and 2 shaft engines, mixed and unmixed flows) can be simulated by activating or deactivating components or entire shafts, by defining pressure ratios and efficiencies equal to 1. Additionally, several bleed configurations were modeled in order to give the user the ability to configure the bleed port extraction position, the amount of bleed extraction and the destination of the air being bled from the compressors: outboard bleed (for engine operability, aircraft air conditioning, pressurization, and anti-ice) as well as turbine cooling. In the case of the air being bled for turbine cooling purposes the user can select where exactly the cooling flow will be inserted in the cycle: stators or rotors of the turbine stages. At last, the model can deal with power extraction from all the shafts for aircraft systems. The schematics in Fig. 1 shows the engine model architecture with the airflow paths, power extractions, and power links (components mechanically linked through the shaft) following the proposed nomenclature from SAE AS755. This diagram represents the most complex engine architure to be simulated.The model was built based on blocks with will calculate each engine module separately. The blocks developed for this model are:

Figure 1
Three shaft direct drive engine model diagram.

• AMB (Standard Atmosphere): this block reads the Altitude, Flight velocity or Mach Number, Ambient temperature or deviation from standard day and air humidity and calculates the engine air inlet properties based on the U.S. Standard Atmosphere 1976, Antoine (1888), and Gordon (1982).

• Air Inlet: this block reads the ambient properties calculated by the Standard Atmosphere block, the input pressure recovery factor and calculates the air intake performance based on the MIL-E-5007D.

• Splitter: the mass flow splitter block is used in several different places in the model, such as bypass and bleed extractions, and its basic function is to split the inlet flow in two outlet flows with the same gas properties.

• Compressor: this block reads the compressor characteristics, such as pressure ratio and isentropic efficiency and calculates the gas outlet properties based on the inlet properties and gas compression as described in Saravanamuttoo et al. (2001).

• Burner (Combustion Chamber): this block reads the fuel characteristics, such as lower fuel heating value and hydrogen/carbon ratio, and combustion chamber characteristics, such as pressure ratio and exit temperature or fuel flow and calculates the combustion gases properties based on inlet air properties as per Gordon(1982).

• Turbine: the turbine block calculates the gas expansion based on the turbine isentropic efficiency and inlet properties as described in Saravanamuttoo et al. (2001).

• Duct losses (Bypass duct and jet pipe): this block calculates the pressure loss through a duct given the pressure recovery factor.

• Mixer: this block calculates the resulting gas properties based on the 2 inlet gas flows. The calculation is based on the chemical composition, pressure, and temperature of each gas flow.

• Exhaust Nozzle: this block calculates the exhaust gas properties and velocity, based on the nozzle inlet gas properties and nozzle coefficients and geometry (convergent or convergent-divergent), as well as gross thrust.

Figure 2 shows the model simulation main process diagram. The flowchart represents all the engine blocks, libraries, input data and iterations necessary to simulate the engine performance. The main steps necessary to perform the simulation are:

Figure 2
Engine simulation process diagram.

• Design Point input read. This block reads all input data necessary to characterize the engine modules and calculate each block at design point.

• Calculate each engine module at component level in the sequence of the gas flow in order to reach the Design Point performance at engine level.

• Read the components maps for off-design performance simulations.

• Scale the components maps based on each module performance previously calculated at Design Point.

• Output the simulation results and the components scaled maps for off-design simulation.

• After finishing the Design Point calculation read the of design inputs, such as operating condition and power setting.

• Set the iterative process starting point. In this model the starting point can be the set equal to the last successfully converged point or a pre-defined starting point calculated based on the flight condition and power setting.

• Calculate the engine components performance and overall performance.

• Check if all the energy balances, mass flow balance and power settings are respected. If so output the calculated engine performance else a new iteration shall be performed with the new operation condition calculated by Broyden or Newton-Raphson method for non-linear system of equation solving.

MODEL DESCRIPTION

The mathematical model described herein are simplified for the sake of the reader clarity. More details can be obtained in the open literature as Mckinney (1967), Koenig and Fishback (1972), Fishback and Koenig (1972), Szuk (1974), Palmer and Yang (1974), Macmillan (1974), Sellers (1975), Wittenberg (1976), Flack (1990), Stamatis et al.(1990), Ismail and Bhinder (1991), Korakianitis and Wilson (1994), Baig and Saravanamuttoo (1997), Bringhenti (1999), Saravanamuttoo et al. (2001), Walsh and Fletcher (2004).

STANDARD ATMOSPHERE

The Standard Atmosphere definition implemented in the model described in this paper is based on the U.S. Standard Atmosphere 1976. It splits the atmosphere in 5 different levels up to 85 km grouping in the same level altitudes with similar characteristics of temperature and pressure variation as the altitude increases.

A summary of the atmosphere properties calculation for each atmosphere layer and the parameters to be used in static temperature and pressure calculation are described in Table 1.

(1) $\begin{array}{c}\mathit{T}={\mathit{T}}_{\mathit{b}}+{\mathit{L}}_{\mathit{b}}\times \left(\mathit{ALT}-{\mathit{H}}_{\mathit{b}}\right)\hfill \\ \mathit{If}\to {\mathit{L}}_{\mathit{b}}\cancel{=}0\hfill \end{array}$

(2) $\begin{array}{c}\mathit{P}={\mathit{P}}_{\mathit{b}}\times {\left[\frac{{\mathit{T}}_{\mathit{b}}}{{\mathit{T}}_{\mathit{b}}+{\mathit{L}}_{\mathit{b}}\times \left(\mathit{ALT}-{\mathit{H}}_{\mathit{b}}\right)}\right]}^{\left[\frac{{{\mathit{g}}^{\prime }}_0\times {\mathit{M}}_0}{{\mathit{R}}^{*}\times {\mathit{L}}_{\mathit{b}}}\right]}\hfill \\ \mathit{If}\to {\mathit{L}}_{\mathit{b}}=0\hfill \end{array}$

(3) $\mathit{P}={\mathit{P}}_{\mathit{b}}\times {\mathit{e}}^{\left[\frac{{{\mathit{g}}^{\prime }}_0\times {\mathit{M}}_0\times \left(\mathit{ALT}-{\mathit{H}}_{\mathit{b}}\right)}{\mathit{R}*\times {\mathit{T}}_{\mathit{b}}}\right]}$

where: g' ₀ = 9.80665 m/s² is the geopotential gravity; M₀ = 28.96443 kg/mol is the air molecular weight; and R* = 8,314.62 J/mol∙K is the universal gas constant.

Table 1
Standard atmosphere properties calculation summary table (U.S. Standard Atmosphere 1976).

Index (b)	Layer	Geopotential altitude Hb (km)	Thermal gradient Lb (K/km)	Reference temperature Tb (K)	Reference pressure Pb (Pa)
0	Troposphere	0	–6.5	288.15	101,325.0000
1	Tropopause	11	0.0	216.65	22,631.9500
2	Stratosphere	20	+1.0	216.65	5,475.0960
3		32	+2.8	228.65	868.0107
4	Stratopause	47	0.0	270.65	110.9002
5	Mesosphere	51	–2.8	270.65	66.9383
6		71	–2.0	214.65	3.9563

HUMIDITY

At all altitudes it is possible to set the humidity contained in the air. For this calculation the Antoine equation (Antoine 1888) determines the saturation vapor pressure for a given temperature for pure components. The Antoine equation and constants for water are:

(4) ${\mathit{P}}_{\mathit{SAT}}={10}^{\mathit{A}-\frac{\mathit{B}}{\mathit{C}+{\mathit{T}}_{\mathit{water}}}}$

where: P_SAT is the saturation pressure in mmHg; T_water is the water static temperature in °C; A, B, and C are constants that are specific for each substance. The constants for water are shown in Table 2.

Table 2
Constants for water saturation vapor pressure in Antoine equation.

Temperature	A	B	C
From 1 to 100 °C	8.07131	1,730.63	233.426
From 100 to 374 °C	8.14019	1,810.94	244.485

INTAKE

The engine air inlet simulation was implemented following the MIL-E-5007D, which describes the pressure recovery factors for subsonic, supersonic, and hypersonic flows. Once the engine air inlet does no thermodynamic work and the flow is considered adiabatic, the stagnation temperature through the duct remains constant. Air mass flow and chemical composition also remain the same. The stagnation pressure downstream the air inlet is calculated as follows:

Subsonic flight (Mach < 1)

(5) ${\mathit{P}}_{{\mathit{T}}_{\mathit{OUT}}}=\mathit{RAMREC}\times {\mathit{P}}_{{\mathit{T}}_{\mathit{IN}}}$

Supersonic flight (1 ≤ Mach < 5)

(6) ${\mathit{P}}_{{\mathit{T}}_{\mathit{OUT}}}=\mathit{RAMREC}\times {\mathit{P}}_{{\mathit{T}}_{\mathit{IN}}}\times \left[1-0.75\times {\left(\mathit{MN}-1\right)}^{1.35}\right]$

Hypersonic flight (Mach ≥ 5)

(7) ${\mathit{P}}_{{\mathit{T}}_{\mathit{OUT}}}=\mathit{RAMREC}\times {\mathit{P}}_{{\mathit{T}}_{\mathit{IN}}}\times \left(\frac{800}{{\mathit{MN}}^4+935}\right)$

where Pt_IN is the inlet stagnation pressure in Pa; Pt_OUT is the outlet stagnation pressure in Pa; MN is the Mach Number; RAMREC the engine air inlet pressure recovery (Pt_OUT/Pt_IN).

COMPRESSOR

The axial flow compressor was implemented following the classic formulation described by Saravanamuttoo et al. (2001), Walsh and Fletcher (2004), and Kurzke (2007). The main equations in the compressor model are described as:

(8) $\mathit{CPR}=\frac{{\mathit{Pt}}_{\mathit{OUT}}}{{\mathit{Pt}}_{\mathit{IN}}}$

The increase in the stagnation temperature due to work added to the airflow is calculated by:

(9) $\frac{{\mathit{Tt}}_{\mathit{OUT}}}{{\mathit{Tt}}_{\mathit{IN}}}=1+\frac{1}{{\mathit{\eta }}_{\mathit{c}}}\times \left[{\left(\frac{{\mathit{Pt}}_{\mathit{OUT}}}{{\mathit{Pt}}_{\mathit{IN}}}\right)}^{\frac{\left(\mathit{\gamma }-1\right)}{\mathit{\gamma }}}-1\right]$

and the thermodynamic specific work is calculated by:

(10) ${\mathit{w}}_{\mathit{Comp}}={\mathit{h}}_{\mathit{OUT}}-{\mathit{h}}_{\mathit{IN}}$

where CPR is the compressor pressure ratio; Tt_IN is the inlet stagnation temperature in K; Tt_OUT is the outlet stagnation temperature in K; γ is the specific heat ratio (Cp/Cv, being Cp and Cv the specific heat at constant pressure and volume respectively); η_c is the compressor isentropic efficiency; w_Comp is the compressor specific work in W/kg; h_IN is the inlet stagnation specific enthalpy in J/kg; h_OUT is the outlet stagnation specific enthalpy in J/kg.

COMBUSTION CHAMBER

The combustion chamber model calculates the amount of burnt fuel considering the amount of air and the equivalence ratio. Equivalence ratio is the ratio between the actual fuel air ratio and stoichiometric fuel air ratio, so equivalence ratio equal to 1 means stoichiometric burn, while lower and higher values mean lean and rich burns respectively. The chemical composition of the burnt gases is determined by the following equation, for equivalence ratio (ER) ≤ 1, as proposed by Gordon (1982):

(11) $\begin{array}{c}{\mathit{CH}}_{\mathit{Y}}+\mathit{\alpha }({\mathit{O}}_2+3.727587·{\mathit{N}}_2+0.0447068·\mathit{Ar}+\hfill \\ +0.0015228·{\mathit{CO}}_2+\mathit{\beta }·{\mathit{H}}_2\mathit{O})\hfill \\ \to \left(1+0.0015228·\mathit{\alpha }\right)·{\mathit{CO}}_2+\left(\frac{\mathit{Y}}{2}+\mathit{\alpha }·\mathit{\beta }\right)·{\mathit{H}}_2\mathit{O}+\hfill \\ +3.727587·\mathit{\alpha }·{\mathit{N}}_2+\left(\mathit{\alpha }-1-\frac{\mathit{Y}}{4}\right)·{\mathit{O}}_2+0.0447068·\mathit{\alpha }·\mathit{Ar}\hfill \end{array}$

And the following equation is proposed for ER>1 considering the air limiting the combustion:

(12) $\begin{array}{c}{\mathit{CH}}_{\mathit{Y}}+\mathit{\alpha }({\mathit{O}}_2+3.727587·{\mathit{N}}_2+0.0447068·\mathit{Ar}+\hfill \\ +0.0015228·{\mathit{CO}}_2+\mathit{\beta }·{\mathit{H}}_2\mathit{O})\hfill \\ \to \left[\mathit{\alpha }·\left(\frac{1}{1+\raisebox{1ex}{$\mathit{Y}$}\!\left/ \!\raisebox{-1ex}{$4$}\right.}+0.0015228\right)\right]·{\mathit{CO}}_2+\left[\mathit{\alpha }·\left(\frac{\mathit{Y}}{4+\mathit{Y}}·\mathit{\beta }\right)\right]·{\mathit{H}}_2\mathit{O}+\hfill \\ +3.727587·\mathit{\alpha }·{\mathit{N}}_2+\hfill \\ +\left[\mathit{\alpha }·\left(1-\frac{1}{\left(1+\raisebox{1ex}{$\mathit{Y}$}\!\left/ \!\raisebox{-1ex}{$4$}\right.\right)}-\frac{2·\mathit{Y}}{\left(4+\mathit{Y}\right)}\right)\right]·{\mathit{O}}_2+0.0447068·\mathit{\alpha }·\mathit{Ar}\hfill \end{array}$

where Y is the fuel hydrogen-carbon ratio; β is the water-air mass flow ratio; α is (4+Y)/(4·ER).

The unburnt air is mixed to the combustion gases and the chemical composition of the gas leaving the burner is recalculated. Burner exit temperature can be either inputted or calculated based on the fuel flow. In both cases, the following equation is used to calculate the temperature from fuel flow or fuel flow from temperature:

(13) $\mathit{WF}\times \mathit{LFHV}=\frac{\left({\dot{\mathit{m}}}_{\mathit{OUT}}\times {\mathit{h}}_{\mathit{OUT}}\right)-\left({\dot{\mathit{m}}}_{\mathit{IN}}\times {\mathit{h}}_{\mathit{IN}}\right)}{{\mathit{\eta }}_{\mathit{cc}}}$

where WF is the fuel flow in kg/s; LFHV is the lower fuel heating value in J/kg; ṁ_IN is the mass flow at burner inlet in kg/s; ṁ_OUT is the mass flow at burnet outlet in kg/s; η_cc is the Ccombustion efficiency.

TURBINE

Turbine performance prediction is calculated as follows:

(14) $\frac{{\mathit{Tt}}_{\mathit{OUT}}}{{\mathit{Tt}}_{\mathit{IN}}}=1+{\mathit{\eta }}_{\mathit{t}}\times \left[{\left(\frac{{\mathit{Pt}}_{\mathit{OUT}}}{{\mathit{Pt}}_{\mathit{IN}}}\right)}^{\frac{\left(\mathit{\gamma }-1\right)}{\mathit{\gamma }}}-1\right]$

where η_t is the turbine isentropic efficiency.

The expansion through the turbine generates the necessary power to drive the compressor mechanically linked to the turbine by a shaft. The turbine power can be calculated as follows:

(15) $\mathit{W}=\dot{\mathit{m}}\times \left({\mathit{h}}_{\mathit{out}}-{\mathit{h}}_{\mathit{in}}\right)$

where: ṁ is the gas mass flow at turbine inlet in kg/s.

PROPELLING NOZZLE

In this model, 2 different nozzle geometries were implemented: convergent and convergent-divergent (con-di). For the con-di nozzle, 7 different flow configurations were implemented, as described by Devenport (2001) and Shapiro (1953). The pressure distribution in the nozzle for each configuration is shown in Fig. 3.

Figure 3
Pressure distribution through the nozzle (Devenport 2001). (a) Not choked at throat; (b) Just choked at throat; (c) Shock in nozzle; (d) Shock at exit; (e) Overexpanded; (f) Design condition; (g) Underexpanded.

GAS PROPERTIES

A good gas properties model is key for any thermodynamic cycle analysis. In order to keep the flexibility and accuracy of the engine performance simulations the gas properties model was developed with refined and detailed data from the Reference Fluid Thermodynamic and Transport Properties (REFPROP; Lemmon et al. 2013). All the main gases present in the air and combustion gases composition (N₂, O₂, CO₂, Ar and H₂O) were modeled separately. The gas property is so calculated depending on its chemical composition and the partial contributions of each specific gas enthalpy and molar mass. Enthalpy was modeled considering the effects of different temperatures and pressures.

OFF-DESIGN

The 3 major contributors who enabled the model to converge in few iterations and, therefore, short clock time were the powerful nonlinear system of equation solver, the maps interpolation method and the definition of the starting point of the iterative process.

NON-LINEAR SYSTEM OF EQUATION SOLVER

For the 3 shaft engine architectures the nonlinear system of equation is composed by 8 equations and 8 variables - equations: LP (low pressure) shaft work balance, LP shaft mass flow balance, IP (intermediate pressure) shaft work balance, IP shaft mass flow balance, HP (high pressure) shaft work balance, HP shaft mass flow balance, engine core mass flow balance, and fuel flow/Max cycle temperature constraint; variables: engine mass flow, fan pressure ratio, IP compressor pressure ratio, HP compressor pressure ratio, HP turbine pressure ratio, IP turbine pressure ratio, LP turbine pressure ratio and fuel flow).

The Broyden's Method (Broyden 1965) was selected from trade study that was conducted to define which system of equations solver would give the shortest clock time to find the solution.

The Broyden's method is a generalization of the secant method to nonlinear systems. The secant method replaces the Newton's method derivative by a finite difference:

(16) ${\mathit{f}}^{\prime }\left({\mathit{x}}_{\mathit{k}}\right)\approx \frac{\mathit{f}\left({\mathit{x}}_{\mathit{k}}\right)-\mathit{f}\left({\mathit{x}}_{\mathit{k}-1}\right)}{{\mathit{x}}_{\mathit{k}}-{\mathit{x}}_{\mathit{k}-1}}$

(17) ${\mathit{f}}^{\prime }\left({\mathit{x}}_{\mathit{k}}\right)\left({\mathit{x}}_{\mathit{k}}-{\mathit{x}}_{\mathit{k}-1}\right)\approx \mathit{f}\left({\mathit{x}}_{\mathit{k}}\right)-\mathit{f}\left({\mathit{x}}_{\mathit{k}-1}\right)$

where f is the function whose zeros are being searched; x is the free variable; k is the iteration number.

Broyden's gave a system of equation generalization:

(18) ${\mathit{J}}_{\mathit{F}}\left({\mathit{x}}_{\mathit{k}}\right)\left({\mathit{x}}_{\mathit{k}}-{\mathit{x}}_{\mathit{k}-1}\right)\approx \mathit{f}\left({\mathit{x}}_{\mathit{k}}\right)-\mathit{f}\left({\mathit{x}}_{\mathit{k}-1}\right)$

where J_F is the Jacobian calculated for the system of equations; F is a matrix with the solution of each equation calculated for x_k .

Thus it is not necessary to calculate the Jacobian and all its derivatives of the Newton's method in every iteration, therefore this method is time saving at a cost of slightly lower convergence rate.

MAPS INTERPOLATION METHOD

The developed computer program make use of maps for compressors and turbines for off-design calculation. The implemented method to find the operating condition and interpolate within the map values is based on linear interpolation. However, in order to improve the interpolation time, the search for the nearest points for interpolation was enhanced. Usually the map would be read from the first line to the last looking for an interval that comprises the search point. It works fine if the interpolation point is close to the table head, usually close to the design point. However, the farthest the point is from the table head more data is necessary to be read and checked, which make the interpolation slow. In order to improve the searching for the nearest points it was implemented a procedure based on Point in Polygon (PIP) concept. The procedure consists in divide the map in four quadrants and check if the interpolation point is within one of the quadrants. The check is done by checking the sum of the angles between the interpolation point and the quadrant vertices. If the sum is 2π it means that it is in the quadrant and if it is 0, it is not. Once the quadrant that contains the interpolation point is found the same procedure is repeated reducing the quadrant size until the quadrant is formed only by the 4 nearest points, when it is ready for the interpolation. The closest three points defines a plane that comprises the interpolation point and therefore the interpolation within the plane can be calculated. The plane interpolation was implemented in order to avoid bilinear interpolation issues where the mass flow is constant and the interpolation in mass flow axis would lead to a division by 0.

Figure 4 shows a compressor map, as an example, and it is possible to observe how the map is divided into 4 quadrants successively until the quadrant is formed only by the 4 nearest points to the operating condition. Figure 5 shows how the angles between the operating point and the quadrant vertices shall be considered for PIP evaluation, whose possible values are described in Eqs. 19 and 20.

(19) $\mathrm{If}\sum {\mathit{\theta }}_{\mathit{i}}=2\mathit{\pi }\mathrm{then}\mathrm{the}\mathrm{point}\mathrm{is}\mathrm{within}\mathrm{the}\mathrm{polygon}$

(20) $\mathrm{If}\sum {\mathit{\theta }}_{\mathit{i}}=0\mathrm{then}\mathrm{the}\mathrm{point}\mathrm{is}\mathrm{within}\mathrm{the}\mathrm{polygon}$

Figure 4
Quadrant division example in a compressor map.

Figure 5
PIP graphic representation.

ITERATION STARTING POINT

Another extremely powerful feature of the model which improves both convergence success rate and time until the solution is the selection of the starting point close to the final solution. Obviously, the final solution is not known until the simulation is completed, but an approximation of the final solution can be estimated based on some engine parameters. In the model developed for this paper the engine parameters to start the iteration are set based on the flight condition and the a power setting parameter. The design point parameters are corrected to the off-design flight condition and then corrected to the input power setting. The power setting parameter defines the engine power such as fuel flow, burner exit temperature and shaft speed. All of them can be set as input to the model.

MODEL VERIFICATION

The developed model was compared in terms of thrust and fuel flow calculation with an existing commercial gas turbine performance model. The model for reference was GasTurb11® (Kurzke 2007) which is very known, reliable and flexible to receive the same kind of inputs necessary to set the model developed for this paper. The simulations, for all the 3 architectures, were based on a burner exit temperature sweep at ISA Sea Level Static condition and compared using the same compressors and turbine maps. Figures 6 to 11 show the comparison between the GasTurb11® and the developed model. The divergences found in thrust and fuel flow are due to differences in the combustion gas model. The gas model in GasTurb11® does not consider pressure in the enthalpy calculation while the developed model does. Also, the combustion gases composition calculation may lead to differences in the cycle calculation mainly downstream to the burner. The model could not be compared in terms of run time because no models were found in the literature with the ability to run in real time.

Figure 6
Single shaft engine thrust comparison.

Figure 7
Single shaft simulated engine fuel flow comparison.

Figure 8
Two shaft simulated engine thrust comparison.

Figure 9
Two shaft simulated engine fuel flow comparison.

Figure 10
Three shaft simulated engine thrust comparison.

Figure 11
Three shaft simulated engine fuel flow comparison.

Three different engine architectures were simulated and compared in terms of thrust and fuel flow with the engines modeled in GasTurb11® with same configuration. The architectures are the most utilized in the aeronautic industry: single, 2 and 3 shafts direct drive engines with unmixed flows and convergent nozzles. The Design Point of the models is shown in Table 3.

Table 3
Simulated engines design point for output comparison with GasTurb11®.

	Single shaft	2shaft	3shaft
Altitude (m)	0	0	0
Temperature (K)	288.15	288.15	288.15
Flight Mach Number	0	0	0
Air inlet mass flow (kg/s)	50	50	100
Air intake pressure recovery	0.99	0.99	0.99
Bypass ratio	0	5	5
Inner LP compressor pressure ratio	-	2	1.8
Inner LP compressor isentropic efficiency	-	0.88	0.88
Outer LP compressor pressure ratio	-	1.8	2
Outer LP compressor isentropic efficiency	-	0.5	0.88
IP compressor pressure ratio	-	-	2
IP compressor isentropic efficiency	-	-	0.88
HP compressor pressure ratio	12	7	5
HP compressor isentropic efficiency	0.85	0.85	0.88
Fuel heating value (kJ/kg)	43,124	43,124	43,124
Burner exit temperature (K)	1,500	1,600	1,650
Burner pressure ratio	0.97	0.97	0.98
Burner isentropic efficiency	0.9999	0.9995	0.9995
HP turbine isentropic efficiency	0.89	0.9	0.9
HP shaft mechanical efficiency	0.99	0.99	0.98
IP turbine isentropic efficiency	-	-	0.9
IP shaft mechanical efficiency	-	-	0.99
LP turbine isentropic efficiency	-	0.9	0.9
LP shaft mechanical efficiency	-	0.99	0.99

Figures 6 and 7 show the thrust and fuel flow comparison with GasTurb11® for the turbojet architecture (one shaft direct drive engine). Figures 8 and 9 show the same comparison for the 2 shaft direct drive turbofan. Finally, Figs. 10 and 11 are the comparison between the models for 3 shaft direct drive turbofan engine.

In Figs. 6 and 7 the blue and red curves refer to the calculated parameters, thrust or fuel flow, by this paper's model and GasTurb11®, respectively, and the values are in the left vertical axis. The difference between the values calculated by the model described in this paper and Gasturb11® are shown by the green curve whose values are in the right vertical axis. Differences are expected due to the gas properties model differences and premises in the 2 different simulation tools.

TEST MATRIX FOR RUN TIME EVALUATION

In order to test the convergence time, the model was run at different off-design conditions to explore different component map regions. The off-design conditions were set by inputting different altitudes, Mach Numbers, temperatures and 1 engine power set, burner exit temperature in this assessment. Table 4 summarizes the chosen values used to simulate different engine operational conditions.

Table 4
Simulation test matrix.

Altitude	Mach Number	Delta standard day	Burner exit temperature
Sea level - 15,000 m (steps of 500 m)	Static - 0.8 (steps of 0.05)	−30 °C +30 °C (steps of 5 °C)	1,800 K - 1,000 K (steps of 100 K)

RESULTS

The run time distribution and the number of iterations until the convergence are shown in Figs. 12 and 13, respectively. The run times were achieved in a personal computer with Intel Core i7 920 at 2.67GHz and the solver convergence criteria was set to square root of the machine precision which was in the computer where the points were run, 10⁻⁸. The results are disposed in a histogram chart where it is shown the distribution of the number of converged points, in the ordinates, by the elapsed time until convergence (Fig. 12) or number or iterations until the convergence (Fig. 13), in abscissas. The points and the operating conditions evaluated are described in Table 4.

Figure 12
Run time histogram.

Figure 13
Number of iterations histogram.

An additional run time reducing opportunity was assessed in order to improve the model run time: iteration stopping criteria relaxing. In order to provide accuracy in the calculations the stopping criteria was chosen to be the square root of the machine precision. Figure 14 shows that the model converges very quickly to the solution and spends a lot of iterations refining the solution to meet the very tight stopping criteria. The chart shows the evolution of 3 of the equations in the nonlinear system of equations for off design calculation. When the equation goes to 0 it means it converged. It can be seen that the parameters converge very quickly to the solution, approximately 4 iterations in the example, and require another 5 iteration to refine the final solution to meet the excessively sharp stopping criteria.

Figure 14
Iteration steps until the solution.

The potential run time improvement due to the stopping criteria relaxation was assessed and the results are shown in Fig. 15. The chart shows the number of converged points in the ordinates and the run time in the abscissas. It can be seen that the peak and the average of the red columns, which represents the run time of the points with relaxed stopping criteria, are at lower run times when compared with the blue columns, which represents the points with original stopping criteria. It means that, by relaxing the stopping criteria, in general, the points converged faster, as expected.

Figure 15
Stopping criteria relaxing benefit in run time.

Similar result can be verified in Fig. 16. The chart shows the benefit that the stopping criteria relaxing brought in terms of numbers of iterations. The peak and the average moved to the left, what means that more points converged at lower number of iterations.

Figure 16
Stopping criteria relaxing benefit in the number of iterations.

CONCLUSIONS

A brand new engine performance prediction model was developed with the ability to run and reach the convergence in most of the times in less than 30ms, which is compatible with a high definition video format, whose refresh rate is 30 frames per second. The features implemented in the model to improve the run time were very effective and ensure good model performance, within the target run time. Additionally, the model did not lose accuracy and flexibility with those features. In fact, by setting the starting point close to the final solution, the convergence success rate was also improved. An additional feature which was also investigated, the relaxation in the iteration stopping criteria could improve even more the run time at a cost of some accuracy loss.

The aim of this study was to demonstrate a model with the ability to simulate the performance of single, 2 and 3 shaft gas turbines with run times compatible with real-time applications with high-fidelity accuracy. The developed model was verified using commercial gas turbine performance software.

ACKNOWLEDGEMENTS

The financial support from Empresa Brasileira de Aeronáutica (Embraer), Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq), Centro de Pesquisa e Inovação Sueco-Brasileiro (CISB), and Svenska Aeroplan AB (SAAB) is acknowledged.

REFERENCES

• Antoine C (1888) Tensions des vapeurs; nouvelle relation entre les tensions et les températures. Comptes Rendus des Séances de l'Académie des Sciences (in French). Paris.
• Baig MF, Saravanamuttoo HIH (1997) Off-design performance prediction of single-spool turbojets using gasdynamics. J Propul Power 13(6):808-810. doi: 10.2514/2.5240
  » https://doi.org/10.2514/2.5240
• Bringhenti C (1999) Análise de desempenho de turbinas a gás em regime permanente (Master's thesis). São José dos Campos: Instituto Tecnológico de Aeronáutica.
• Bringhenti C (2003) Variable geometry gas turbine performance analysis (PhD thesis). São José dos Campos: Instituto Tecnológico de Aeronáutica. In Portuguese.
• Broyden CG (1965) A class of methods for solving nonlinear simultaneous equations. Math Comp 19:577-593.
• Devenport WJ (2001) Instructions; [accessed 2015 May 15]. http://www.engapplets.vt.edu/fluids/CDnozzle/cdinfo.html
  » http://www.engapplets.vt.edu/fluids/CDnozzle/cdinfo.html
• Fishback LH, Koenig RW (1972) GENENG II: A program for calculating design and off-design performance of two and three spool turbofans with as many as three nozzles. (NASA TN D-6553).Washington: NASA.
• Flack RD (1990) Analysis and matching of gas turbine components. International Journal of Turbo and Jet Engines 7(3-4):217-226. doi: 10.1515/TJJ.1990.7.3-4.217
  » https://doi.org/10.1515/TJJ.1990.7.3-4.217
• Gordon S (1982) Thermodynamic and transport combustion properties of hydrocarbons with air. Lewis Research Center Cleveland.
• Grönstedt T (2000) Development of methods for analysis and optimization of complex jet engine systems (PhD thesis). Göteborg: Chalmers University of Technology.
• Ismail IH, Bhinder FS (1991) Simulation of aircraift gas turbine engines. J Eng Gas Turbines Power 113(1):95-99. doi: 10.1115/1.2906536
  » https://doi.org/10.1115/1.2906536
• Koenig RW, Fishback LH (1972) GENENG: A Program for calculating design and off-design performance for turbojet and turbofan engines. (NASA TN D-6552). Washington: NASA.
• Korakianitis T, Wilson DG (1994) Models for predicting the performance of Brayton-cycle engines. J Eng Gas Turbines Power 116(2):381-388. doi: 10.1115/1.2906831
  » https://doi.org/10.1115/1.2906831
• Kurzke J (2007) GasTurb11: design and off-design performance of gas turbines. Aachen: GasTurb GmbH.
• Lemmon EW, Huber ML, McLinden MO (2013) NIST Standard Reference Database 23: Reference Fluid Thermodynamic and Transport Properties-REFPROP, Version 9.1. Gaithersburg: National Institute of Standards and Technology.
• Macmillan WL (1974) Development of a modular type computer computer program for the calculation of gas turbine off-design performance. Cranfield: Cranfield Institute of Technology.
• Mckinney JS (1967) Simulation of turbofan engine - SMOTE: Description of method and balancing technique. AD-825197/AFAPL-TR-67-125. Air Force Aero Propulsion Lab. pt.1-2.
• Palmer JR, Yang CZ (1974) TURBOTRANS: A programming language for the performance simulation of arbitrary gas turbine engines with arbitrary control systems. ASME Paper 82-GT-200.
• Saravanamuttoo HIH, Rogers GFC, Cohen H (2001) Gas turbine theory. 5th edition. New York: Prentice Hall.
• Sellers JD (1975) DYNGEN: a program for calculating steady-state and transient performance of turbojet and turbofan engines. NASA TN D-7901. Washington: NASA.
• Shapiro AH (1953) The dynamics and thermodynamics of compressible fluid flow. Vol. I-II. New York: The Ronald Press Company.
• Silva FJS (2011) Gas turbines performance study under the influence of variable geometry transients (PhD thesis). São José dos Campos: Instituto Tecnológico de Aeronáutica. In Portuguese.
• Stamatis A, Mathioudakis K, Papailiou KD (1990) Adaptive simulation of gas turbine performance. J Eng Gas Turbines Power 112(2):168-175. doi: 10.1115/1.2906157
  » https://doi.org/10.1115/1.2906157
• Szuk JR (1974) HYDES: a generalized hybrid computer program for studying turbojet or turbofan engine dynamics. NASA TM X-3014.Washington: NASA.
• Walsh PP, Fletcher P (2004) Gas turbine performance. 2nd edition. Oxford: Blackwell Science.
• Wittenberg H (1976) Prediction of off-design performance of turbojet and turbofan engines. CP-242-76, AGARD.

Publication Dates

• Publication in this collection
  Jul-Sep 2017

History

• Received
  19 May 2016
• Accepted
  08 Nov 2016