Run time Assessment for Gas Turbine Performance Simulation

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. Correspondence author: Henrique Gazzetta Junior | Avenida Cassiano Ricardo, 1.411 – Apto. 84B | CEP: 12.240-540 – São José dos Campos/SP – Brazil | E-mail: henrique.gazzetta@gmail.com Received: May 17, 2016 | Accepted: Jul. 10, 2017 Section Editor: Dimitrios Pavlou ABSTRACT: This article describes the run time characteristics of a gas turbine performance simulation using different solvers and components off-design performance database formats. Two different nonlinear systems of equation solvers, NewtonRaphson’s and Broyden’s, and two different formats of compressor and turbine off-design performance database (maps), tabulated values and fi tted surface equations, were compared. Based on the results it is then possible to trade off and select the most appropriate combination of solver and component map type for the gas turbine performance simulation for real-time application.


INtRodUCtIoN
Th e simulation activities have been more used in the industry in earlier phases of the aircraft design.It is due to the recognition of value in identifying potential issues or even exploring the systems characteristics in early phases of the design.Th e necessary corrections or modifi cations obtained by the simulation results are far less costly to be implemented in earlier design phases, when there are no manufactured parts or signed contracts with suppliers, in comparison with the changes required in later phases, when changes in real parts or even in already closed interfaces may be required.Th e knowledge of the aircraft engine performance and transient characteristics is important to support optimization studies when the aircraft and the engine are still open for modifi cations.A gas turbine performance prediction tool is key in this process.A highfi delity real-time engine transient performance calculation tool could support steady state and transient engineering analysis, system integration studies or even enable tests with pilot in the loop in earlier phases of the aircraft development.Th erefore, a real-time gas turbine performance code was specially developed to attend those necessities.Th is article aims to compare the simulation time using two diff erent nonlinear systems of equation solver methods: Newton-Raphson's and Broyden's.Th ese two methods were implemented into the developed computer program in order to verify which one gives the best performance for real-time application.
xx/xx 02/13 methodology MODEL DESCRIPTION A brand new engine simulation model was developed to simulate gas turbine performance in real time.The model has the ability to use different solvers and different compressor and turbine maps types.In addition to the engine thermodynamic cycle parameters, the model consider also the run time.In this work, a two-shaft direct drive gas turbine was chosen to be studied.The engine architecture is shown in Fig. 1 following the proposed nomenclature from SAE AS755 (2004).The computational code was developed in a modular structure, each module representing an engine component with its own characteristics, and then linked to each other to get the whole engine performance.This type of structure gives to the program flexibility to simulate several engine configurations without changing the program structure.The components modeling and other model features are described in more details in Gazzetta (2017) and Gazzetta et al. (2017).The model simulation main process diagram is shown in Fig. 2.

AMB Air Inlet Fan
Fan air valve

DESIGN POINT
At design point, a full characterization of the ambient condition, air intake losses, turbomachine, fuel and propelling nozzles are required for the thermodynamic cycle calculation.Typical data required for that are shown in All those information are used to set the design condition and scale the component maps for off-design simulations.At off-design conditions, only the ambient condition and one single power setting parameter are required as input data for the thermodynamic cycle calculation.All the air intake, turbomachine and propelling nozzles characterization are defined by off-design characteristics mapped in tables and maps.As a power setting it may be used burner exit temperature, shaft speed, fuel flow or thrust.An iterative procedure is used to set the turbomachine operation to match the input power set.The iterative procedure is described in the off-design section.

OFF-DESIGN
For the two shaft engine architecture chosen to be studied in this paper the nonlinear system of equation is composed by six equations and six variables (Equations: LP shaft work balance; LP shaft mass flow balance; HP 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; HP compressor pressure ratio; HP turbine pressure ratio; LP turbine pressure ratio and fuel flow).The equations are described in Eqs.1-6.
Nozzle vs. inlet mass flow balance (Conservation of mass): where: W 2 is the core inlet mass flow; W 8 is the core nozzle mass flow; WB 23 is the bleed air from booster; WB 3 is the bleed air from high pressure compressor and WF is the engine fuel flow.Low-pressure shaft power balance (Conservation of energy): (1) xx/xx 04/13 where: Pwr FAN is the power required by the fan; Pwr Booster is the power required by the booster; Pwr LPT is the power produced by the low-pressure turbine; HPX LP is the shaft power extracted from the LP shaft and η LPShaft is the low-pressure shaft mechanical efficiency.High-pressure shaft power balance (Conservation of energy): where: Pwr HPC is the power required by the high-pressure compressor; Pwr HPT is the power produced by the high-pressure turbine; HPX HP is the shaft power extracted from the HP shaft and η HPShaft is the high-pressure shaft mechanical efficiency.
Power setting constraint: one of the following constraints will be selected based on the user input.The user can run the offdesign simulation to match fuel flow, net thrust, burner exit temperature and shaft speeds: where: WF is the engine calculated fuel flow; WF Tgt is the engine target fuel flow; FN is the engine calculated net thrust; FN Tgt is the engine target net thrust; T4 is the burner calculated exit temperature; T4 Tgt is the burner target temperature; N1 and N2 are the calculated low-and high-pressure shaft speeds respectively; N1 Tgt and N2 Tgt are the low and high-pressure target shafts speeds respectively.High-pressure turbine mass flow balance (Conservation of mass): where W HPT is the calculated high-pressure turbine mass flow and W 4 is the burner exit mass flow.Low-pressure turbine mass flow balance (Conservation of mass): where: W LPT is the calculated low-pressure turbine mass flow and W 4 is the burner exit mass flow. (2) In order to find the most appropriate nonlinear system of equations solver for the gas turbine performance simulation two different methods were tested: Newton-Raphson and Broyden.
Newton-Raphson's method is widely know and used to solve nonlinear system of equations.It is based on the idea of linear approximation to calculate the next iteration steps: where: f is the function whose zeros are being searched; x is the free variable; 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 Broyden's method (Broyden 1965) is a generalization of the secant method to nonlinear systems.The secant method replaces the Newton's method derivative by a finite difference: where: f is the function whose zeros are being searched; x is the free variable; k is the iteration number.
Broyden 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 ; x is the free variable; k is the iteration number.Thus, it is not necessary to calculate the Jacobian and all its derivatives of the Newton's method, therefore this method is timesaving at a cost of lower convergence rate.
It was also compared two different types of compressor and turbine map data: tabulated and fitted surface equation.
The tabulated map format is very well known and widely used by all gas turbine performance simulation tools.In this format, the map data is organized in a table format and the simulation shall interpolate the tables to find the off-design operating condition.The format of the tabulated map is shown in Tables 2, 3 and 4. Map data was obtained from Bringhenti (1999).
In order to completely avoid the interpolation task in the maps it is proposed to use surfaces equations instead of tables to represent the component off-design characteristics.The 3D surface equations are generated from the original tabulated map data and shall be representative of the points in the table.The expression derived to represent the maps and their general format is where: x, y and f (x,y) represent any combination of parameters that better matches the tabulated data in the component map and N is the equation order.The parameters C i,j in Eq. 12 are calculated based on the least squares methodology applied for polynomials of degree N as in Miller (2006), Dai et al. (2007) and Chernov and Ma (2011).This format of surface equation was chosen because the polynomial equations can be fitted in scattered data by using least squares methodology and because it gives flexibility in the trade between fitting accuracy and equation complexity simply by changing the order of the equation.It was selected ninth order for this assessment due to better fitting to the map data.
For the fitted surface equation, it was chosen to calculate compressor or turbine speed given pressure ratio and corrected mass flow.This methodology provides a very good matching between the tabulated data and calculated by the fitting surface equation, as shown in Figs. 3 and 4.  Figure 3 shows the delta between fan speed calculated by the fitted surface and the tabulated data.Only the tabulated map knots were considered in the comparison, with no interpolation.The maximum and standard deviation around 0.4% and 0.1% respectively, in addition to the average centered in 0%, mean that the chosen surface equation is a good representation of the tabulated data.
Figure 4 shows the same comparison of the Fig. 3 identifying where in the fan map the differences are located.The importance of this chart is to show that there are no high deviation spots that could affect the off-design calculation.
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 one engine power set, burner exit temperature in this assessment.The same set of data in the same sequence was used in the two solver methods and map types.Table 5 summarizes the chosen values used to simulate different engine operational conditions and Table 6 the combinations of solvers and map types.

ReSUltS
The run time distribution and the number of iterations until the convergence are shown, for each combination of solver and map type as per Table 6, in the histogram charts below.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 simulation run times for the cases specified in Table 5 are shown in Figs. 5 and 6.The results are disposed in histogram charts where it is shown the distribution of the number of converged points, in the ordinates, by the elapsed time until convergence, in abscissas.The points and the operating conditions evaluated are described in Table 5.
Similar charts are used to show the distribution of the number of iterations until the convergence.From Fig. 7 to Fig. 10, the histograms charts show the distribution of the number of converged points, in ordinates, by the number of iterations until the convergence, in the abscissas.
The histogram charts from Fig. 7 to Fig. 10 show that Broyden's method generally takes more iterations than Newton-Raphson's to reach the solution but takes shorter clock time.It is because Newton-Raphson's method performs the derivatives calculation several times to build the Jacobian, which is very time-consuming, while Broyden's method calculates the Jacobian just once in the first iteration.Another conclusion is that the fitted surface equations, defined from tabulated maps after fitting methods, are more computationally costly than the tabulated map due to the high order of the equations to keep a good matching with the tabulated data.In this study, it was used fitting surface equation of ninth order.
By reducing the order of the map, fitting surface equation to sixth instead of ninth order the simulation time was improved for both solvers while the number of iteration to reach the solution remains the same as shown in Figs.11 to 14.It means that the solver is doing the same iterations but spending less time on each one due to lower number of surface equation terms.Reduced order of the fitting surface equation may lead to lower fidelity at off-design simulation.
Figures 11 and 12 are histograms that show distribution of the number of converged points, in ordinates, by the run time in the abscissas.The peaks offset to the left in both Newton-Raphson's and Broyden's methods show that most of the points spent less clock time to reach the solution when the fitted surface equations were simpler (6 th order instead of 9 th ).
Figures 13 and 14 are histograms that show the distribution of the number of converged points, in ordinates, by the number of iterations until the solution, in the abscissas.Those charts show that there was no effect of the fitted surface equation order reduction (from 9 th to 6 th order) on the number of iterations until convergence, once there is no offset in the charts curves, neither to the left nor to the right.It means that, the model will do the same number of iterations but faster once the maps fitted surface are less costly due to the lower equation order.

CoNClUSIoNS
A brand new engine performance prediction model was built with the ability to use two different nonlinear system of equations solvers, Newton-Raphson's and Broyden's, and two different component map format, tabulated and fitted surface equation.Using the simulation results, all possible combinations of solvers and maps format were compared in terms of run time and number of iterations until the solution, including an additional map fitted surface equation option with a reduced equation order.Several different computational tools were found to fit polynomial surface equations in a scattered data; however, none of them was capable to fit high order equations, necessary to not degrade the map calculation accuracy.Therefore, a new tool was developed based on the existing least squares methodology expanded to high order surface equation to generate the maps equations.
Among the options tested and compared in this paper, the conclusion is that the most appropriate solver for a real-time gas turbine performance simulation is the Broyden's method combined with the tabulated map.An additional comparison with a lower order map fitted surface equation shows that, if high fidelity model is not extremely required or if the run time is more important than the model accuracy (in flight simulator, for instance), the map fitted surface equation could be considered combined with Broyden's method.

Figure 3 .
Figure 3. Histogram of fan speed delta between tabulated map data and calculated by the fitted surface equation at same pressure ratio and corrected mass flow.

Figure 4 .
Figure 4. Map of fan speed delta between tabulated map data and calculated by the fitted surface equation at same pressure ratio and corrected mass flow.

Figure 14 .
Figure 14.Maps surface equation order impact in the number of iterations for EQ-Broyden.

table 1 .
Required input for design point run.

table 2 .
Relative corrected mass flow tabulated map.

table 6 .
Combinations of solvers and map types.
Maps surface equation order impact in run time for EQ-NRaphson.
Figure 9. EQ-NRaphson number of iterations histogram.Figure 10.EQ-Broyden number of iterations histogram.Figure 12. Maps surface equation order impact in run time for EQ-Broyden.Figure 13.Maps surface equation order impact in the number of iterations for EQ-NRaphson.