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.

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.

Power Extraction
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.
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/s 2 is the geopotential gravity; M 0 = 28.96443kg/mol is the air molecular weight; and R* = 8,314.62J/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; T water 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.  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) (4) (4) (4)  (5) (5) (5) 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: 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: 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; 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): 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 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 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 2 , O 2 , CO 2 , Ar and H 2 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.roat Exit

non-lineAR systeM of equAtion solveR
Distance down the nozzle 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 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.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 are described in Eqs.19 and 20.

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 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          vertical axis.Differences are expected due to the gas properties model differences and premises in the 2 different simulation tools.

Corrected mass ow
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.

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 −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       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.

Run time [ms] Number of points
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.

For the 3
shaft engine architectures the nonlinear system of equation is composed by 8 equations and 8 variablesequations: 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 4 .
Figure 4. Quadrant division example in a compressor map.
i = 2π then the point is within the polygon (19) If Σθ i = 0 then the point is out of the polygon (20) 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

Figure 7 .
Figure 7. Single shaft simulated engine fuel flow comparison.

Figure 11 .
Figure 11.Three shaft simulated engine fuel flow comparison.

Figure 9 .
Figure 9. Two shaft simulated engine fuel flow comparison.
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.

Figure 15 .
Figure 15.Stopping criteria relaxing benefit in run time.

Figure 14 .
Figure 14.Iteration steps until the solution.

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

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

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