J. Aerosp. Technol. Manag., São José dos Campos, Vol.9, No 3, pp.346-356, Jul.-Sep., 2017ABSTRACT: Industry and universities around the world [627402]

J. Aerosp. Technol. Manag., São José dos Campos, Vol.9, No 3, pp.346-356, Jul.-Sep., 2017ABSTRACT: 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.Real-Time Gas Turbine Model for
Performance Simulations
Henrique Gazzetta Junior1, Cleverson Bringhenti2, João Roberto Barbosa2, Jesuíno Takachi Tomita2
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 Y ang (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 doi: 10.5028/jatm.v9i3.693
1.Empresa Brasileira de Aeronáutica – São José dos Campos/SP – Brazil. 2.Departamento de Ciência e Tecnologia Aeroespacial – Instituto Tecnológico de
Aeronáutica – Divisão de Engenharia Aeronáutica e Mecânica – São José dos Campos/SP – Brazil.
Author for correspondence: Henrique Gazzetta Junior | Avenida Cassiano Ricardo, 1.411 – Apto. 84B | CEP: 12.240-540 – São José dos Campos/SP – Brazil |
Email: henrique.gazzetta@gmail.comReceived: 19 May 2016 | Accepted: 08 Nov 2016

J. Aerosp. Technol. Manag., São José dos Campos, Vol.9, No 3, pp.346-356, Jul.-Sep., 2017347
Real-Time Gas Turbine Model for Performance Simulations
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:
• 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
Figure 1. Three shaft direct drive engine model diagram.
AMB
Air Inlet Fan air valve Compressor
Compressor
IPC

HPC

AMB AMB
IP Sha/f_t
LP Sha/f_t HP Sha/f_t HP Customer bleedFan bleed

Power
Extraction
A/f_terburner
Exhaust
Nozzle Jet pipe Turbine HPT
Turbine IPT
Turbine LPT Burner IP Customer bleed Booster 4 3 5 0 1
2 13
25 43 45 15 16
9 6 61 7 28 WB28
HPTPwrX WB15
Mixer
Spliter
Fan 12

23 18 Power
Extraction
IPTPwrX

LPTPwrX Bypass duct 17
8 Exhaust
Nozzle
WB3
HPT NGV Cooling HPT Blade Cooling IPT NGV Cooling IPT Blade Cooling LPT NGV Cooling LPT Blade Cooling HPT NGV Cooling HPT Blade Cooling IPT NGV Cooling IPT Blade Cooling LPT NGV Cooling LPT Blade Cooling WC 28F
WC 28E
WC 28D
WC 28C
WC 28B
WC 28A
WC3F
WC3 E
WC3 D
WC3 C
WC3 AWC3 B
Power Extraction

J. Aerosp. Technol. Manag., São José dos Campos, Vol.9, No 3, pp.346-356, Jul.-Sep., 2017348
Gazzetta Junior H, Bringhenti C, Barbosa JR, Tomita JT
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:
• 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 Y ang (1974), Macmillan (1974), Sellers (1975), Wittenberg (1976), Flack (1990), Stamatis et al.(1990), Ismail
and Bhinder (1991), Korakianitis and Wilson (1994), Baig and
Design point (DP)
DP input card
DP simulationOff design point (ODP)
ODP input read
ODP simulation
Simulation convergedFirst guess for iterative process
Components map reading
Components map scaling
Output of simulation results and scaled maps
Output of simulation resultsCalculate new inputs using linear systems solvers
Ye sNo
Figure 2. Engine simulation process diagram.

J. Aerosp. Technol. Manag., São José dos Campos, Vol.9, No 3, pp.346-356, Jul.-Sep., 2017349Real-Time Gas Turbine Model for Performance Simulations
Saravanamuttoo (1997), Bringhenti (1999), Saravanamuttoo
et al. (2001), Walsh and Fletcher (2004).
stAndARd AtMospheRe
Th e Standard Atmosphere defi nition implemented in the
model described in this paper is based on the U.S. Standard
Atmosphere 1976. It splits the atmosphere in 5 diff erent 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.
where: g
’ 0 = 9.80665 m/s2 is the geopotential gravity;
M0 = 28.96443 kg/mol is the air molecular weight; and
R* = 8,314.62 J/mol∙K is the universal gas constant.
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. Th e Antoine equation and constants for water are:where: P
SAT is the saturation pressure in mmHg; Twater is the
water static temperature in °C; A, B, and C are constants that
are specifi c for each substance. Th e constants for water are shown in Table 2.
index
(b)layergeopotential 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
Stratosphere20 +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
Mesosphere51 –2.8 270.65 66.9383
6 71 –2.0 214.65 3.9563Table 1. Standard atmosphere properties calculation summary table (U.S. Standard Atmosphere 1976).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.485Table 2. Constants for water saturation vapor pressure in
Antoine equation.
intAKe
Th e engine air inlet simulation was implemented following
the MIL-E-5007D, which describes the pressure recovery
factors for subsonic, supersonic, and hypersonic fl ows. Once
the engine air inlet does no thermodynamic work and the fl ow
is considered adiabatic, the stagnation temperature through the
duct remains constant. Air mass fl ow and chemical composition
also remain the same. Th e stagnation pressure downstream the air inlet is calculated as follows:
Subsonic fl ight (Mach < 1)
Supersonic fl ight (1 ≤ Mach < 5)(1)
(3)(4)
(5)
(6)(2) b b b H ALT LTT  (1)








bLRMg
b b bb
bH ALT LTTPP*0'
0
(2)







bb
TRH ALT Mg
bePP*0'
0
(3)

waterTCBA
SATP
10 (4)

IN OUT T T P RAMREC P   (5)

 35.1)1 (75.01   MN P RAMREC P
IN OUT T T (6)




 935800
4MNP RAMREC P
IN OUT T T (7)

INOUT
PtPtCPR (8)









1111

INOUT
c INOUT
PtPt
TtTt (9)

IN OUT Comp h h w  (10)
b b b H ALT LTT  (1)








bLRMg
b b bb
bH ALT LTTPP*0'
0
(2)







bb
TRH ALT Mg
bePP*0'
0
(3)

waterTCBA
SATP
10 (4)

IN OUT T T P RAMREC P   (5)

 35.1)1 (75.01   MN P RAMREC P
IN OUT T T (6)




 935800
4MNP RAMREC P
IN OUT T T (7)

INOUT
PtPtCPR (8)









1111

INOUT
c INOUT
PtPt
TtTt (9)

IN OUT Comp h h w  (10)
b b b H ALT LTT  (1)








bLRMg
b b bb
bH ALT LTTPP*0'
0
(2)







bb
TRH ALT Mg
bePP*0'
0
(3)

waterTCBA
SATP
10 (4)

IN OUT T T P RAMREC P   (5)

 35.1)1 (75.01   MN P RAMREC P
IN OUT T T (6)




 935800
4MNP RAMREC P
IN OUT T T (7)

INOUT
PtPtCPR (8)









1111

INOUT
c INOUT
PtPt
TtTt (9)

IN OUT Comp h h w  (10)
b b b H ALT LTT  (1)








bLRMg
b b bb
bH ALT LTTPP*0'
0
(2)







bb
TRH ALT Mg
bePP*0'
0
(3)

waterTCBA
SATP
10 (4)

IN OUT T T P RAMREC P   (5)

 35.1)1 (75.01   MN P RAMREC P
IN OUT T T (6)




 935800
4MNP RAMREC P
IN OUT T T (7)

INOUT
PtPtCPR (8)









1111

INOUT
c INOUT
PtPt
TtTt (9)

IN OUT Comp h h w  (10)
b b b H ALT LTT  (1)








bLRMg
b b bb
bH ALT LTTPP*0'
0
(2)







bb
TRH ALT Mg
bePP*0'
0
(3)

waterTCBA
SATP
10 (4)

IN OUT T T P RAMREC P   (5)

 35.1)1 (75.01   MN P RAMREC P
IN OUT T T (6)




 935800
4MNP RAMREC P
IN OUT T T (7)

INOUT
PtPtCPR (8)









1111

INOUT
c INOUT
PtPt
TtTt (9)

IN OUT Comp h h w  (10)
b b b H ALT LTT  (1)








bLRMg
b b bb
bH ALT LTTPP*0'
0
(2)







bb
TRH ALT Mg
bePP*0'
0
(3)

waterTCBA
SATP
10 (4)

IN OUT T T P RAMREC P   (5)

 35.1)1 (75.01   MN P RAMREC P
IN OUT T T (6)




 935800
4MNP RAMREC P
IN OUT T T (7)

INOUT
PtPtCPR (8)









1111

INOUT
c INOUT
PtPt
TtTt (9)

IN OUT Comp h h w  (10)
b b b H ALT LTT  (1)








bLRMg
b b bb
bH ALT LTTPP*0'
0
(2)







bb
TRH ALT Mg
bePP*0'
0
(3)

waterTCBA
SATP
10 (4)

IN OUT T T P RAMREC P   (5)

 35.1)1 (75.01   MN P RAMREC P
IN OUT T T (6)




 935800
4MNP RAMREC P
IN OUT T T (7)

INOUT
PtPtCPR (8)









1111

INOUT
c INOUT
PtPt
TtTt (9)

IN OUT Comp h h w  (10)

Hypersonic fl ight (Mach ≥ 5)
(7)

J. Aerosp. Technol. Manag., São José dos Campos, Vol.9, No 3, pp.346-356, Jul.-Sep., 2017350
Gazzetta Junior H, Bringhenti C, Barbosa JR, Tomita JT
where PtIN is the inlet stagnation pressure in Pa; PtOUT is the
outlet stagnation pressure in Pa; MN is the Mach Number;
RAMREC the engine air inlet pressure recovery (PtOUT/PtIN).
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:And the following equation is proposed for ER>1 considering
the air limiting the combustion:
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: b b b H ALT LTT  (1)








bLRMg
b b bb
bH ALT LTTPP*0'
0
(2)







bb
TRH ALT Mg
bePP*0'
0
(3)

waterTCBA
SATP
10 (4)

IN OUT T T P RAMREC P   (5)

 35.1)1 (75.01   MN P RAMREC P
IN OUT T T (6)




 935800
4MNP RAMREC P
IN OUT T T (7)

INOUT
PtPtCPR (8)









1111

INOUT
c INOUT
PtPt
TtTt (9)

IN OUT Comp h h w  (10)
b b b H ALT LTT  (1)








bLRMg
b b bb
bH ALT LTTPP*0'
0
(2)







bb
TRH ALT Mg
bePP*0'
0
(3)

waterTCBA
SATP
10 (4)

IN OUT T T P RAMREC P   (5)

 35.1)1 (75.01   MN P RAMREC P
IN OUT T T (6)




 935800
4MNP RAMREC P
IN OUT T T (7)

INOUT
PtPtCPR (8)









1111

INOUT
c INOUT
PtPt
TtTt (9)

IN OUT Comp h h w  (10)
b b b H ALT LTT  (1)








bLRMg
b b bb
bH ALT LTTPP*0'
0
(2)







bb
TRH ALT Mg
bePP*0'
0
(3)

waterTCBA
SATP
10 (4)

IN OUT T T P RAMREC P   (5)

 35.1)1 (75.01   MN P RAMREC P
IN OUT T T (6)




 935800
4MNP RAMREC P
IN OUT T T (7)

INOUT
PtPtCPR (8)









1111

INOUT
c INOUT
PtPt
TtTt (9)

IN OUT Comp h h w  (10)
(8)
(9)
(13)
(14)(10)The increase in the stagnation temperature due to work
added to the airflow is calculated by:
and the thermodynamic specific work is calculated by:
where: CPR is the compressor pressure ratio; TtIN is the inlet
stagnation temperature in K; TtOUT 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; wComp is the compressor
specific work in W/kg; hIN is the inlet stagnation specific enthalpy
in J/kg; hOUT 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):(12)
(11)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:
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)
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  
 
Ar OYN OHYCOOH CO Ar N O CHY
 

 

    
    
0447068.041727587.320015228.010015228.0 0447068.0 727587.3
22 2 22 2 2 2
(11)

 
Ar OYY
YN OHYYCOYOH CO Ar N O CHY
 








 











   
     
0447068.042
4111727587.340015228.0
4110015228.0 0447068.0 727587.3
22 2 22 2 2 2
(12)

 
ccIN IN OUT OUT h m h mLFHV WF  (13)










1 11


INOUT
t
INOUT
PtPt
TtTt (14)

in outh hmW  (15)


11'


k kk k
kxxxf xfxf (16)

1 1 '  k k k k k xf xf xxxf (17)

1 1   k k k k k F xf xf xxxJ (18)

polygon e within this point then the2 Ifπ θi (19)
polygon theofout is point then the0 If θi (20)  
 
Ar OYN OHYCOOH CO Ar N O CHY
 

 

    
    
0447068.041727587.320015228.010015228.0 0447068.0 727587.3
22 2 22 2 2 2
(11)

 
Ar OYY
YN OHYYCOYOH CO Ar N O CHY
 








 











   
     
0447068.042
4111727587.340015228.0
4110015228.0 0447068.0 727587.3
22 2 22 2 2 2
(12)

 
ccIN IN OUT OUT h m h mLFHV WF  (13)










1 11


INOUT
t
INOUT
PtPt
TtTt (14)

in outh hmW  (15)


11'


k kk k
kxxxf xfxf (16)

1 1 '  k k k k k xf xf xxxf (17)

1 1   k k k k k F xf xf xxxJ (18)

polygon e within this point then the2 Ifπ θi (19)
polygon theofout is point then the0 If θi (20)  
 
Ar OYN OHYCOOH CO Ar N O CHY
 

 

    
    
0447068.041727587.320015228.010015228.0 0447068.0 727587.3
22 2 22 2 2 2
(11)

 
Ar OYY
YN OHYYCOYOH CO Ar N O CHY
 








 











   
     
0447068.042
4111727587.340015228.0
4110015228.0 0447068.0 727587.3
22 2 22 2 2 2
(12)

 
ccIN IN OUT OUT h m h mLFHV WF  (13)










1 11


INOUT
t
INOUT
PtPt
TtTt (14)

in outh hmW  (15)


11'


k kk k
kxxxf xfxf (16)

1 1 '  k k k k k xf xf xxxf (17)

1 1   k k k k k F xf xf xxxJ (18)

polygon e within this point then the2 Ifπ θi (19)
polygon theofout is point then the0 If θi (20)  
 
Ar OYN OHYCOOH CO Ar N O CHY
 

 

    
    
0447068.041727587.320015228.010015228.0 0447068.0 727587.3
22 2 22 2 2 2
(11)

 
Ar OYY
YN OHYYCOYOH CO Ar N O CHY
 








 











   
     
0447068.042
4111727587.340015228.0
4110015228.0 0447068.0 727587.3
22 2 22 2 2 2
(12)

 
ccIN IN OUT OUT h m h mLFHV WF  (13)










1 11


INOUT
t
INOUT
PtPt
TtTt (14)

in outh hmW  (15)


11'


k kk k
kxxxf xfxf (16)

1 1 '  k k k k k xf xf xxxf (17)

1 1   k k k k k F xf xf xxxJ (18)

polygon e within this point then the2 Ifπ θi (19)
polygon theofout is point then the0 If θi (20)  
 
Ar OYN OHYCOOH CO Ar N O CHY
 

 

    
    
0447068.041727587.320015228.010015228.0 0447068.0 727587.3
22 2 22 2 2 2
(11)

 
Ar OYY
YN OHYYCOYOH CO Ar N O CHY
 








 











   
     
0447068.042
4111727587.340015228.0
4110015228.0 0447068.0 727587.3
22 2 22 2 2 2
(12)

 
ccIN IN OUT OUT h m h mLFHV WF  (13)










1 11


INOUT
t
INOUT
PtPt
TtTt (14)

in outh hmW  (15)


11'


k kk k
kxxxf xfxf (16)

1 1 '  k k k k k xf xf xxxf (17)

1 1   k k k k k F xf xf xxxJ (18)

polygon e within this point then the2 Ifπ θi (19)
polygon theofout is point then the0 If θi (20)  
 
Ar OYN OHYCOOH CO Ar N O CHY
 

 

    
    
0447068.041727587.320015228.010015228.0 0447068.0 727587.3
22 2 22 2 2 2
(11)

 
Ar OYY
YN OHYYCOYOH CO Ar N O CHY
 








 











   
     
0447068.042
4111727587.340015228.0
4110015228.0 0447068.0 727587.3
22 2 22 2 2 2
(12)

 
ccIN IN OUT OUT h m h mLFHV WF  (13)










1 11


INOUT
t
INOUT
PtPt
TtTt (14)

in outh hmW  (15)


11'


k kk k
kxxxf xfxf (16)

1 1 '  k k k k k xf xf xxxf (17)

1 1   k k k k k F xf xf xxxJ (18)

polygon e within this point then the2 Ifπ θi (19)
polygon theofout is point then the0 If θi (20)  
 
Ar OYN OHYCOOH CO Ar N O CHY
 

 

    
    
0447068.041727587.320015228.010015228.0 0447068.0 727587.3
22 2 22 2 2 2
(11)

 
Ar OYY
YN OHYYCOYOH CO Ar N O CHY
 








 











   
     
0447068.042
4111727587.340015228.0
4110015228.0 0447068.0 727587.3
22 2 22 2 2 2
(12)

 
ccIN IN OUT OUT h m h mLFHV WF  (13)










1 11


INOUT
t
INOUT
PtPt
TtTt (14)

in outh hmW  (15)


11'


k kk k
kxxxf xfxf (16)

1 1 '  k k k k k xf xf xxxf (17)

1 1   k k k k k F xf xf xxxJ (18)

polygon e within this point then the2 Ifπ θi (19)
polygon theofout is point then the0 If θi (20)  
 
Ar OYN OHYCOOH CO Ar N O CHY
 

 

    
    
0447068.041727587.320015228.010015228.0 0447068.0 727587.3
22 2 22 2 2 2
(11)

 
Ar OYY
YN OHYYCOYOH CO Ar N O CHY
 








 











   
     
0447068.042
4111727587.340015228.0
4110015228.0 0447068.0 727587.3
22 2 22 2 2 2
(12)

 
ccIN IN OUT OUT h m h mLFHV WF  (13)










1 11


INOUT
t
INOUT
PtPt
TtTt (14)

in outh hmW  (15)


11'


k kk k
kxxxf xfxf (16)

1 1 '  k k k k k xf xf xxxf (17)

1 1   k k k k k F xf xf xxxJ (18)

polygon e within this point then the2 Ifπ θi (19)
polygon theofout is point then the0 If θi (20)  
 
Ar OYN OHYCOOH CO Ar N O CHY
 

 

    
    
0447068.041727587.320015228.010015228.0 0447068.0 727587.3
22 2 22 2 2 2
(11)

 
Ar OYY
YN OHYYCOYOH CO Ar N O CHY
 








 











   
     
0447068.042
4111727587.340015228.0
4110015228.0 0447068.0 727587.3
22 2 22 2 2 2
(12)

 
ccIN IN OUT OUT h m h mLFHV WF  (13)










1 11


INOUT
t
INOUT
PtPt
TtTt (14)

in outh hmW  (15)


11'


k kk k
kxxxf xfxf (16)

1 1 '  k k k k k xf xf xxxf (17)

1 1   k k k k k F xf xf xxxJ (18)

polygon e within this point then the2 Ifπ θi (19)
polygon theofout is point then the0 If θi (20)  
 
Ar OYN OHYCOOH CO Ar N O CHY
 

 

    
    
0447068.041727587.320015228.010015228.0 0447068.0 727587.3
22 2 22 2 2 2
(11)

 
Ar OYY
YN OHYYCOYOH CO Ar N O CHY
 








 











   
     
0447068.042
4111727587.340015228.0
4110015228.0 0447068.0 727587.3
22 2 22 2 2 2
(12)

 
ccIN IN OUT OUT h m h mLFHV WF  (13)










1 11


INOUT
t
INOUT
PtPt
TtTt (14)

in outh hmW  (15)


11'


k kk k
kxxxf xfxf (16)

1 1 '  k k k k k xf xf xxxf (17)

1 1   k k k k k F xf xf xxxJ (18)

polygon e within this point then the2 Ifπ θi (19)
polygon theofout is point then the0 If θi (20)  
 
Ar OYN OHYCOOH CO Ar N O CHY
 

 

    
    
0447068.041727587.320015228.010015228.0 0447068.0 727587.3
22 2 22 2 2 2
(11)

 
Ar OYY
YN OHYYCOYOH CO Ar N O CHY
 








 











   
     
0447068.042
4111727587.340015228.0
4110015228.0 0447068.0 727587.3
22 2 22 2 2 2
(12)

 
ccIN IN OUT OUT h m h mLFHV WF  (13)










1 11


INOUT
t
INOUT
PtPt
TtTt (14)

in outh hmW  (15)


11'


k kk k
kxxxf xfxf (16)

1 1 '  k k k k k xf xf xxxf (17)

1 1   k k k k k F xf xf xxxJ (18)

polygon e within this point then the2 Ifπ θi (19)
polygon theofout is point then the0 If θi (20)  
 
Ar OYN OHYCOOH CO Ar N O CHY
 

 

    
    
0447068.041727587.320015228.010015228.0 0447068.0 727587.3
22 2 22 2 2 2
(11)

 
Ar OYY
YN OHYYCOYOH CO Ar N O CHY
 








 











   
     
0447068.042
4111727587.340015228.0
4110015228.0 0447068.0 727587.3
22 2 22 2 2 2
(12)

 
ccIN IN OUT OUT h m h mLFHV WF  (13)










1 11


INOUT
t
INOUT
PtPt
TtTt (14)

in outh hmW  (15)


11'


k kk k
kxxxf xfxf (16)

1 1 '  k k k k k xf xf xxxf (17)

1 1   k k k k k F xf xf xxxJ (18)

polygon e within this point then the2 Ifπ θi (19)
polygon theofout is point then the0 If θi (20)  
 
Ar OYN OHYCOOH CO Ar N O CHY
 

 

    
    
0447068.041727587.320015228.010015228.0 0447068.0 727587.3
22 2 22 2 2 2
(11)

 
Ar OYY
YN OHYYCOYOH CO Ar N O CHY
 








 











   
     
0447068.042
4111727587.340015228.0
4110015228.0 0447068.0 727587.3
22 2 22 2 2 2
(12)

 
ccIN IN OUT OUT h m h mLFHV WF  (13)










1 11


INOUT
t
INOUT
PtPt
TtTt (14)

in outh hmW  (15)


11'


k kk k
kxxxf xfxf (16)

1 1 '  k k k k k xf xf xxxf (17)

1 1   k k k k k F xf xf xxxJ (18)

polygon e within this point then the2 Ifπ θi (19)
polygon theofout is point then the0 If θi (20)  
 
Ar OYN OHYCOOH CO Ar N O CHY
 

 

    
    
0447068.041727587.320015228.010015228.0 0447068.0 727587.3
22 2 22 2 2 2
(11)

 
Ar OYY
YN OHYYCOYOH CO Ar N O CHY
 








 











   
     
0447068.042
4111727587.340015228.0
4110015228.0 0447068.0 727587.3
22 2 22 2 2 2
(12)

 
ccIN IN OUT OUT h m h mLFHV WF  (13)










1 11


INOUT
t
INOUT
PtPt
TtTt (14)

in outh hmW  (15)


11'


k kk k
kxxxf xfxf (16)

1 1 '  k k k k k xf xf xxxf (17)

1 1   k k k k k F xf xf xxxJ (18)

polygon e within this point then the2 Ifπ θi (19)
polygon theofout is point then the0 If θi (20)  
 
Ar OYN OHYCOOH CO Ar N O CHY
 

 

    
    
0447068.041727587.320015228.010015228.0 0447068.0 727587.3
22 2 22 2 2 2
(11)

 
Ar OYY
YN OHYYCOYOH CO Ar N O CHY
 








 











   
     
0447068.042
4111727587.340015228.0
4110015228.0 0447068.0 727587.3
22 2 22 2 2 2
(12)

 
ccIN IN OUT OUT h m h mLFHV WF  (13)










1 11


INOUT
t
INOUT
PtPt
TtTt (14)

in outh hmW  (15)


11'


k kk k
kxxxf xfxf (16)

1 1 '  k k k k k xf xf xxxf (17)

1 1   k k k k k F xf xf xxxJ (18)

polygon e within this point then the2 Ifπ θi (19)
polygon theofout is point then the0 If θi (20)

J. Aerosp. Technol. Manag., São José dos Campos, Vol.9, No 3, pp.346-356, Jul.-Sep., 2017351
Real-Time Gas Turbine Model for Performance Simulations
pressure distribution in the nozzle for each configuration is
shown in Fig. 3.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:
gAs pRopeR ties
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
2, O2, CO2,
Ar and H2O) 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 equA tion 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 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.
/T_hroat Exit
Distance down the nozzleP / PIN
1
0(a)
(b)
(c)
(d)
(e)
(f)
(g)
(16)
(17)
(18)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:
where JF is the Jacobian calculated for the system of equations;
F is a matrix with the solution of each equation calculated
for xk .
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 inteRpolA tion 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  
 
Ar OYN OHYCOOH CO Ar N O CHY
 

 

    
    
0447068.041727587.320015228.010015228.0 0447068.0 727587.3
22 2 22 2 2 2
(11)

 
Ar OYY
YN OHYYCOYOH CO Ar N O CHY
 








 











   
     
0447068.042
4111727587.340015228.0
4110015228.0 0447068.0 727587.3
22 2 22 2 2 2
(12)

 
ccIN IN OUT OUT h m h mLFHV WF  (13)










1 11


INOUT
t
INOUT
PtPt
TtTt (14)

in outh hmW  (15)


11'


k kk k
kxxxf xfxf (16)

1 1 '  k k k k k xf xf xxxf (17)

1 1   k k k k k F xf xf xxxJ (18)

polygon e within this point then the2 Ifπ θi (19)
polygon theofout is point then the0 If θi (20)  
 
Ar OYN OHYCOOH CO Ar N O CHY
 

 

    
    
0447068.041727587.320015228.010015228.0 0447068.0 727587.3
22 2 22 2 2 2
(11)

 
Ar OYY
YN OHYYCOYOH CO Ar N O CHY
 








 











   
     
0447068.042
4111727587.340015228.0
4110015228.0 0447068.0 727587.3
22 2 22 2 2 2
(12)

 
ccIN IN OUT OUT h m h mLFHV WF  (13)










1 11


INOUT
t
INOUT
PtPt
TtTt (14)

in outh hmW  (15)


11'


k kk k
kxxxf xfxf (16)

1 1 '  k k k k k xf xf xxxf (17)

1 1   k k k k k F xf xf xxxJ (18)

polygon e within this point then the2 Ifπ θi (19)
polygon theofout is point then the0 If θi (20)  
 
Ar OYN OHYCOOH CO Ar N O CHY
 

 

    
    
0447068.041727587.320015228.010015228.0 0447068.0 727587.3
22 2 22 2 2 2
(11)

 
Ar OYY
YN OHYYCOYOH CO Ar N O CHY
 








 











   
     
0447068.042
4111727587.340015228.0
4110015228.0 0447068.0 727587.3
22 2 22 2 2 2
(12)

 
ccIN IN OUT OUT h m h mLFHV WF  (13)










1 11


INOUT
t
INOUT
PtPt
TtTt (14)

in outh hmW  (15)


11'


k kk k
kxxxf xfxf (16)

1 1 '  k k k k k xf xf xxxf (17)

1 1   k k k k k F xf xf xxxJ (18)

polygon e within this point then the2 Ifπ θi (19)
polygon theofout is point then the0 If θi (20)

J. Aerosp. Technol. Manag., São José dos Campos, Vol.9, No 3, pp.346-356, Jul.-Sep., 2017352Gazzetta Junior H, Bringhenti C, Barbosa JR, Tomita JT
divide the map in four quadrants and check if the interpolation
point is within one of the quadrants. Th e 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. Th e closest three points defi nes a plane that comprises the interpolation point and therefore the interpolation
within the plane can be calculated. Th e plane interpolation was
implemented in order to avoid bilinear interpolation issues where the mass fl ow is constant and the interpolation in mass
fl ow 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 aredescribed in Eqs. 19 and 20.iteRAtion stARting point
Another extremely powerful feature of the model which
improves both convergence success rate and time untilthe solution is the selection of the starting point close to the
fi nal solution. Obviously, the fi nal solution is not known until
the simulation is completed, but an approximation of the fi nal
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 fl ight condition and the a power setting parameter. Th e design point parameters are
corrected to the off -design fl ight condition and then corrected
to the input power setting. Th e power setting parameter defi nes the engine power such as fuel fl ow, burner exit temperature
and shaft speed. All of them can be set as input to the
model.
Model veRificAtion
Th e developed model was compared in terms of thrust
and fuel fl ow calculation with an existing commercial gas
turbine performance model. Th e model for reference was GasTurb11® (Kurzke 2007) which is very known, reliable
and fl exible to receive the same kind of inputs necessary to set the model developed for this paper. Th e 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. Th e divergences found in thrust and fuel
fl ow are due to diff erences in the combustion gas model. Th e
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
diff erences in the cycle calculation mainly downstream to the
burner. Th e 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.
Th ree diff erent engine architectures were simulated and
compared in terms of thrust and fuel fl ow with the engines
modeled in GasTurb11® with same configuration. The
architectures are the most utilized in the aeronautic industry: single, 2 and 3 shaft s direct drive engines with unmixed fl ows
and convergent nozzles. Th e Design Point of the models is
shown in Table 3.
Figures 6 and 7 show the thrust and fuel fl ow comparison
with GasTurb11® for the turbojet architecture (one shaft direct

Pressure ratio
Corrected mass /f_low
Figure 4. Quadrant division example in a compressor map.
Figure 5. PIP graphic representation.

θ
θIf Σθi = 2π then the point is within the polygon (19)
If Σθi = 0 then the point is out of the polygon (20)

J. Aerosp. Technol. Manag., São José dos Campos, Vol.9, No 3, pp.346-356, Jul.-Sep., 2017353
Real-Time Gas Turbine Model for Performance Simulations
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 T4 sweep at sea level / static / ISA
10,00015,00020,00025,00030,00035,00040,00045,00050,000
1,200 1,250 1,300 1,350 1,400 1,450 1,500
Burner exit temperature [K]/T_hrust [N]
-2.0%-1.5%-1.0%-0.5%0.0%0.5%1.0%1.5%2.0%
Delta thrust between
the model and GasTurb11®
Model Gasturb Delta thrust [%]xxxx
/T_hrust [N]
Delta thrust between
the model and GasTurb11®
Burner exit temperature [K]
Model Gasturb Delta thrust [%]T4 sweep at sea level / static / ISA
15,00020,00025,00030,00035,00040,00045,000
1,250 1,300 1,350 1,400 1,450 1,500 1,550 1,600 1,650-3.0%-2.0%-1.0%0.0%1.0%2.0%3.0%
Burner exit temperature [K]/T_hrust [N]
Delta thrust between
the model and GasTurb11®
Model Gasturb Delta thrust [%]T4 sweep at sea level / static / ISA
02,0004,0006,0008,00010,00012,00014,00016,00018,00020,000
1,200 1,250 1,300 1,350 1,400 1,450 1,500 1,550 1,600-2.5%-2.0%-1.5%-1.0%-0.5%0.0%0.5%1.0%1.5%2.0%2.5%
Model Gasturb Delta fuel /f_low [%]T4 sweep at sea level / static / ISA
0.40.50.60.70.80.91.01.11.21.3
1,200 1,250 1,300 1,350 1,400 1,450 1,500
Burner exit temperature [K]Fuel /f_low [kg/s]
-5.0%-4.0%-3.0%-2.0%-1.0%0.0%1.0%2.0%3.0%4.0%
Delta fuel /f_low between
the model and GasTurb11®

Delta fuel /f_low between
the model and GasTurb11 ®
Burner exit temperature [K]
Model Gasturb Delta fuel /f_low [%]T4 sweep at sea level / static / ISA
0.00.10.20.30.40.50.6
1,250 1,300 1,350 1,400 1,450 1,500 1,550 1,600 1,650Fuel /f_low [kg/s]
-3.0%-2.0%-1.0%0.0%1.0%2.0%3.0%Burner exit temperature [K]
Delta fuel /f_low between
the model and GasTurb11®
Model Gasturb Delta fuel /f_low[%]T4 sweep at sea level / static / ISA
0.050.070.090.110.130.150.170.190.210.230.25
1,200 1,250 1,300 1,350 1,400 1,450 1,500 1,550 1,600Fuel /f_low [kg/s]
-5.0%-4.0%-3.0%-2.0%-1.0%0.0%1.0%2.0%3.0%4.0%5.0%

Figure 6. Single shaft engine thrust comparison.
Figure 10. Three shaft simulated engine thrust comparison.
Figure 8. Two shaft simulated engine thrust comparison.Figure 7. Single shaft simulated engine fuel flow comparison.
Figure 11. Three shaft simulated engine fuel flow comparison.Figure 9. Two shaft simulated engine fuel flow comparison.
vertical axis. Differences are expected due to the gas properties
model differences and premises in the 2 different simulation
tools.
test MA tRix foR Run tiMe ev AluA tion
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

J. Aerosp. Technol. Manag., São José dos Campos, Vol.9, No 3, pp.346-356, Jul.-Sep., 2017354
Gazzetta Junior H, Bringhenti C, Barbosa JR, Tomita JT
single
shaft2shaft 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
ratio12 7 5
HP compressor isentropic
efficiency0.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
efficiency0.89 0.9 0.9
HP shaft mechanical
efficiency0.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.99Table 3. Simulated engines design point for output
comparison with GasTurb11®.assessment. Table 4 summarizes the chosen values used to
simulate different engine operational conditions.
ReSul TS
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
−8. 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.
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

Number of iterationsNumber of points18,000
16,000
12,000
10,000
8,000
6,000
4,000
2,00014,000
0
036912151821242730333639424548515457<60
Figure 12. Run time histogram. 32
02,0004,0006,0008,00010,00012,00014,00016,00018,000
10 14 18 22 26 30 34 38 42 46 50 54 58 62 66>70Number of points

Number of points
Run time [ms]
Figure 13. Number of iterations histogram.

J. Aerosp. Technol. Manag., São José dos Campos, Vol.9, No 3, pp.346-356, Jul.-Sep., 2017355
Real-Time Gas Turbine Model for Performance Simulations
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.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)Table 4. Simulation test matrix.
1 –0.020.000.020.040.060.080.100.120.140.160.180.200.22
2 3 4 5Nozzle mass/f_low balance
HPT mass/f_low balance
LPT mass/f_low balance
6 7 8 9
Interation numberParameter value
03691215182125,000
20,000
15,000
10,000
5,000
0
2427303336394245485110–510–8
545760
Number of iterationsNumber of points
10–510–8Number of points
0 2,000 4,000 6,000 8,000 10,000 12,000
10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40 42 44 46 48 50
Run time [ms] Number of points

Figure 15. Stopping criteria relaxing benefit in run time.Figure 14. Iteration steps until the solution.Figure 16. Stopping criteria relaxing benefit in the number
of iterations.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.
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.
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.

J. Aerosp. Technol. Manag., São José dos Campos, Vol.9, No 3, pp.346-356, Jul.-Sep., 2017356
Gazzetta Junior H, Bringhenti C, Barbosa JR, Tomita JT
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
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
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
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
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.2906831Kurzke 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
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.AuThoR’S ConTRIBuTIon
Gazzetta Junior H, Bringhenti C, Barbosa JR and Tomita JT
conceived the idea; Gazzetta Junior H wrote the main text
and prepared the figures. All authors discussed the results and
commented on the manuscript.