Turbine Engine Starting Simulation

The process of engine control development requires the models that describe engine operation and its response on a control action. The development flow required numerous models to be engaged, like component-level non-linear model, engine-level non-linear model, linear dynamic model, etc. Models made a great progress during the recent years and became reliable tools for control engineers. However, most models are derivatives from the component-level non-linear model, which in its turn consumes the component performances. Things turn different when one addresses the starting range of engine operation. The problem here is all about the missing performances of the engine components, as it is quite hard to harvest these performances in this region as the processes that happen in the engine are transient by nature. Different scientists offered different approaches to the problem of building the component level non-linear model of the sub-idle region, but the general idea is to somehow extrapolate the known performances to the sub-idle region. However, there are no known reports about a model that considers all aspects of this approach and simulates the engine starting. In this chapter, you can find an alternative view on a problem of simulation of a sub-idle operation. The proposed model belongs to a group of linear dynamic models including the static model as well as simplified static model to support the dynamic model. Instead of trying to extrapolate component performances and get the full-scale component-level model, you will see that the canonical component performances are replaced by the direct relations between parameters that are used in the control algorithms, like gas-path parameters against the RPM. As well in this chapter, you will find the exact instructions on how to create the model and an example of the one with the real test data.


Introduction
Models are known for their benefit to reduce the price for engine development replacing the numerous experiments with the simulations. This has the biggest value for the control engineers as they need the greatest amount of experimental data. Industry came to a gold standard on how to simulate the models of engine operation within the range limited with idle on one side and the maximum on the other side. These models are well structured and validated, having the well-known thermodynamic relations under the hood. But their usage is impossible or considerably limited in the sub-idle range, because of no experimentally registered performances of engine components. The other reason is that the required performances are challenging to be determined for the sub-idle region, especially the performance of a combustion chamber.
Here we have reached the point that suites the most to introduce the existing models of sub-idle operation. They can go into two baskets. In the first we find the methods that deal with the extrapolation of the component performances into the sub-idle range, while in the second, methods that approximate the experimental data with various polynomials.
We can dive deeper to the first basket and discriminate four groups of models • based on straight extrapolation of known performances (no physics considered) [1,2,13]; • "smart" extrapolation that consumes the theoretical knowledge about component performances in the sub-idle region [3][4][5][6]; • statistical extrapolation [7]; • based on fitting the generic performances to the experimental data [8,12,14].
The brightest method of the first group is based on the simplest similarity laws. The method was developed by Gaudet et al. [1]. The method is only applicable for incompressible fluid. Another method by Sexton [2] suggests obtaining new constant corrected speed lines of compressor spool using experimentally measured ones: where subscript i refers to experimentally measured parameter and i À 1 refers to extrapolated parameter. Gaudet and Gauthier improved Sexton's method by adopting it for compressible fluid. Powers in Eq. (2) are obtained by the analysis of two experimentally measured corrected constant speed lines with the lowest values of corrected rotational speed. Next, a new corrected constant speed line is calculated referring existing by the equations Powers p, q, and r stay constant for all sub-idle operating range. The second group is very similar to the first one and differs in the fact that extrapolation takes into account prior information about processes taking place during starting. This allows improving the tolerance and preventing errors.
All methods of the second group consider only separate facts from prior information (e.g. negative efficiency of compressor spool till the moment of combustion chamber lighting [4]), but there are still no known methods that are able to take into account all prior information. For example, there is a method [5,6] to extrapolate performances of compressor or turbine taking into account the continuity equation violation at negative efficiency. Authors suggest routine performance representation to be changed to the new one. The efficiency performance of compressor spool is expressed as follows: Another method, proposed by Agrawal and Yunis [7], belongs to methods of the third group. This method requires a big number of experimental data for identifying true performances of engine components from universal ones. Relations to carry out identification are the following: • Specific work • Efficiency where K ϕ , K ψ , and K are coefficients depending on a corrected rotational speed. This relation is formed from analysis of a big number of similar engines.
In addition to the above parameters, this method also allows simulation of torque, pressure, and temperature at compressor discharge.
The last group of methods can be represented by Kong's method [8]. The point of the method is getting universal component performances from well-known thermodynamic relations for gas turbine engines and prior information about the engine with its further identification at modes from idle to maximum power (thrust). Here, identification must be understood as scaling factors calculating to match theoretical and true performances. This method has validated its effectiveness for the set of turboshaft engines. But the shortage of this method was poor convergence at modes with low rotational speeds, including sub-idle modes.
To cope with the problem of poor convergence, Jones et al. had improved Kong's method, suggesting application of non-metering components and new criterion for identification quality estimation (surge parameter). Applied measures have improved identification, especially at sub-idle modes [9].
Final improvement of Kong's method was carried out by the group of scientists from Sharif University. They have suggested a new approach to scaling factor calculation. This approach provides good identification quality at all modes. All parameters necessary for compressor map identification can be calculated as follows [10]: W a cor ¼ W a cor ð Þ univ Á W a cor calc = W a cor calc ð Þ univ , N cor ¼ N cor ð Þ univ Á N cor calc = N cor calc ð Þ univ , π ¼ π univ Á π calc À 1 ð Þþ1= π calc univ À 1 ð Þþ1, η ¼ η univ Á η calc =η calc univ : Methods of the second category accept the hypothesis about similarity of starting processes in all engines. Here appears the conclusion that any starting process can be circumscribed by universal relations between control factors and measured parameters. These methods require storing rig and in-flight experimental information to form appropriate data for identification. The polynomial relations are conserved to have the same structure from engine to engine. Only coefficients of polynomial change.
Wrapping up the above overview of the methods to simulate sub-idle operation of the gas turbine, one can conclude the following • most described methods can hardly find the usage in an every day's engineering practice; • most methods address the compressor map extrapolation ignoring the problem of other compressor components, especially the combustion chamber operation; • to identify the performances from the experimental data, the latest must be of high quality, however, as the starting is a dynamic process, the collected data may be corrupted with time lag as well as poor sensor performance at low operational modes; • getting the accurate experimental data to build the high-quality performance map is not always a good value for money as it usually needs mounting extra sensors and carrying out additional experiments (for example, low-inertia thermocouple, strain gage, etc.); • the methods from the second basket have extremely limited application range as they require the same starter to be used as well as the minor changes to the starting control loop are allowed.

Requirements to the starting simulation model
During the starting engine speeds up from either turned off state to an idle running or windmilling to the flight idle. Based on the source of power that drives the rotor during the acceleration, the whole starting may be decomposed into three phases.
The rotor of the engine is dynamically balanced. During the first phase, the rotor is driven by the starter only, as the combustion chamber is not lighted yet. During this phase, the starter drives the rotor from turned off state to the RPM where the combustion chamber is started. The torque balance of this phase is where M res is the total torque of resistance that is generated by compressor and turbine cascades of the rotor, by bearings and driven auxiliary units (fuel and oil pumps, electric generator, deaerator, etc.).
As soon as combustion chamber is brought into a game, the turbine starts giving the positive input to the torque balance equation. Both turbine and starter drive the rotor because the turbine does not have enough power yet to accelerate the rotor alone. The second phase lasts till the turbine reaches the power where it can drive the compressor and accelerate the rotor on its own. The torque equation for the second phase is The starting ends up with the phase where the starter is moved from the table, making the turbine to be the only one source of power that drives the compressor and makes all rotors accelerate to the idle or the flight idle (see Figure 1). The torque balance equation is The first and foremost problem that every simulation engineer faces when designing the model is the problem of the choice among the available model structures, as the selected structure must make a perfect fit to the problem to be solved by the model. The major goal of the proposed model is to perform the preliminary tuning of the ACS. For this, the model must: • provide correct relations among parameters that describe the operation of the engine (as the starting is a non-steady-state process, the model must properly describe the dynamic properties of GTE including all processes that determine engine dynamics); • be able to be integrated with the vast majority of the models that are normally used to support the ACS design, e.g. model that describes the engine operation in the range above the idle, models of governors, etc.; • properly describe the features typical for the starting, for example, parameters overlap during the starting; • require minimum amount of experimental data to be generated and adjusted; • give a quantitative assessment to the required parameters that cannot be measured directly.

Structure
Form the vast majority of model structures, the good choice to fulfill the requirements described above is the dynamic model consisting of linear dynamic model, supported with steady-state model in simplified representation (nextsimplified static model, SSM).
The structure of the model is presented in a Figure 2.
The model of the three-spool turbofan must consider the mutual interaction of engine spools.

Compiling of the SSM
Static performance of GTE in a sub-idle region is nothing more than an abstraction, as in a sub-idle region, the engine does not have any steady modes. The SMM improves the precision of the starting model as well as provides the unique architecture of both sub-idle and above idle models. Let the SSM be formalized as Z N HP 0 is the base value that corresponds to an idle mode.
The proposed scheme of SMM for one of parameters Y from vector in the subidle region is presented in Figure 3. The static performance is extrapolated to the sub-idle region with a linear region and parabolic region (linking region). The performance is described with the functions. They are limited on both sides of the region as • values of the functions must be equal on the left and on the right; • first derivative of the functions must be equal on the left and on the right.
Such a structure allows alternating the slope of the linear region within some range (see line 1 in Figure 3). This is a good opportunity to adjust the SSM after the experiments have been carried out.
Equation of linear region is described as Initial values of slope angle are evaluated from the next formula The 0.6-1.2 range determines the adjustment range of the linear region. If k = 1, then linear region matches with the line joining two points (0, Z 0 ), (N 2 , Z 2 ). Y 0 is equal to a default value of a parameter to be simulated (temperature, pressure, RPM, etc.).
The linking function is described as Unknown coefficients k 0 , k 1 , and k 2 are evaluated from the set of linear equations that is compiled from the limitations for the regions, described in this chapter: The last unknown to be determined from the equation set (18)-(21) is an abscissa N 1 of point joining linear function with the linking function.
Set of Eqs. (18)-(21) is non-linear. Let us use Eqs. (20, 21) to solve it and express the coefficients k 0 and k 1 to be dependent on k 2 : From Eq. (19), we get the relation between N 1 and k 2 : Putting the obtained equations into (18), we get the solution: Building up any performance, we will require turned-off parameters, idle parameters and derivative of the performance parameter against the mode parameter, e.g. RPM of HPR. An example of this map with experimental points is shown in Figure 4.
Because the corrected exhaust gas temperature T T cor has another physical nature, its relation is formed using another approach method.
The turbine temperature decreases to some value when the mode decreases. Further mode decrease results in a drastic temperature increase mainly because of the inefficient turbine operation and low airflow through the engine.
To consider this prior information, the turbine discharge temperature is modeled as follows: where Т 2 cor is the temperature of airflow at compressor discharge; W f cor ÁH u W air cor Áс p is the temperature rise caused by the amount of heat added by the fuel combustion in the combustion chamber.
A coefficient K T considers the temperature drop that is caused by heat retraction in a high-pressure turbine and a free turbine. The coefficient is calculated to provide the equality of the exhaust gas temperatures. The first and the second values of the exhaust-gas temperature are calculated using the sub-idle operation model and the above-idle operation model, respectively.
The airflow through the engine is a priori known to be proportional to the square of the rotational speed: W air cor ¼ C Á N 2 HP . The dependence Т 4 cor ¼ f N HP cor ð Þis shown in Figure 5. The left branch (1) corresponds to the Stage 1 of starting (see Figure 1) when only starter rotates the rotor; the right branch presents Stages 2 and 3 when the combustion chamber is switched in.
It follows from Eq. (10) that for proper modeling of gas temperature on must compile the SSMs of Т 2 cor and airflow W air cor .

Identification of a linear dynamic model
The above dynamic model (1) in the state space for the three-spool turbofan engine under consideration has the form of a system of linear algebraic and differential equations.  • as the argument of the static model and functions for determining the dependence of the coefficients of the linear dynamic model on the engine operational mode is N2, then N HP ¼ N HP st which is why ΔN HP ¼ 0.
• the mutual effect of the rotors is determined mainly by the interaction of the spools of the turbine; therefore, it spreads along the flow (the high-pressure rotor affects the intermediate and low-pressure rotor, and the intermediatepressure rotor affects the low-pressure rotor).
Concluding the just listed assumptions, the set of Eqs. (27) transforms to The task of synthesizing a linear dynamic model is reduced to the problem of determining the coefficients of a linear dynamic model.
The coefficients b 24 , b 21 , and d 11 are determined analytically using a priori and experimental information, followed by approximation by polynomial dependencies, the arguments of which are the rotational speed of the high-pressure rotor. The coefficient b 24 was determined analytically from the following equation: where J HP is an inertia of high pressure rotor.
Having transformed dω HP dt to π 30 Á dN HP dt , we get: The coefficient b 24 does not depend on the type or characteristics of the used starter and is constant for the entire starting.
The determination of the coefficient b 21 пр showing the effect of fuel consumption on the rate of the rotational speed change for the high-pressure rotor was carried out analytically according to the following relationship: Further, the obtained values were approximated, considering the following apriori information about this coefficient at the starting region. Let us consider the following equation: where τ HP is the time constant of the high-pressure rotor. Let us divide both parts of this equation by time constant and move ( 1 τ HP ΔN HP ) to the right part of the equation, and we get It is known that the time constant at the beginning of the rotor spinning (N HP = 0) is equal to zero. At low rotational speeds, the rotor is unstable; therefore, the value of the time constant is negative. Physically, this can be explained by the fact that the turbine in this mode of operation is inefficient, and the slope of the torque characteristic of the turbine (the dependence of the moment on rotation speed at constant fuel consumption) is greater than the slope of the torque characteristic of the compressor. As rotor speed increases, the time constant decreases until the moment when the slope of the torque characteristic of the compressor becomes equal to the slope of the torque characteristic of the turbine ( ∂М HPC ∂n HP ¼ ∂М HPT ∂n HP ). At this point, the function of the dependence of the time constant on the rotor speed experiences a gap (τ HP ¼ π 30 Á J HP ∂М HPC ∂N HP À ∂М HPТ ∂N HP ¼ ∞). Thus, at low rotational speeds, the rotor cannot be stable without an additional energy source or closed-loop control system. With a further increase in the rotational speed, the time constant decreases from infinity to the value at the idle mode.
The coefficient K W has physical sense, therefore, the nature of its change is known. The ratio of this coefficient to the time constant shows the effect of changes in fuel consumption on the rate of N HP change. The dependences of the LDM coefficients of the high-pressure rotor on its rotational speed are shown in Figure 6.
The analytical determination of the coefficients a 11 , b 11 , a 31 , a 33 , and b 31 is impossible; therefore, they are determined by the least squares method engaging the theoretical and experimental information and with the subsequent approximation of the dependence of the obtained values on a parameter that determines the engine operating mode.
So, for example, the evaluation of coefficients a 11 and b 11 is.
where θ ! is a vector of coefficients to be identified (a 11 , b 11 ), the change of which is approximated with polynomials b 11 cor ¼ P 2 i¼0 q i Á N i HP cor and a 11 cor ¼ P 2 i¼0 p i Á N i HP cor ,θ ! -vector estimate θ ! . Figure 6.
Coefficients of the high-pressure rotor dynamics.
Thus, the problem of the evaluation of coefficients a 11 , b 11 reduces to the problem of determining the coefficients of polynomials approximating the change of these coefficients in the sub-idle region: As an example, Figure 7 represents the results of a 11 approximation on the engine operational mode. They correspond to the obtained dependence a 11 cor ¼ À4:0276 Á 10 À9 N 2 HP cor j À 3:0863 Á 10 À5 N HP cor j þ 0:0716: It is known that the coefficient smoothly increases when the rotational speed goes down; therefore, the results of its identification are reliable.
For a more adequate model, the coefficient is refined according to the following relationship:

Model of thermocouple
Measuring the temperature at the sub-idle modes is challenging because the thermocouple has an extremely high time constant. Estimating the real temperature requires complementing the linear dynamic model with an additional equation where ε is the time constant of the thermocouple. The time constant experimental observation requires mounting two thermocouples at the same station: a regular thermocouple and a quick-response thermocouple. Then, the stored data are analyzed to find the temperature delay in time. Unfortunately, the performance of this experiment is notably rare. Hence, the time constant is generally unknown.
Fortunately, this constant can be estimated using the one-time special analysis of the experimental data. The time constant can be identified for the sub-idle modes by substituting T T in Eq. (38) using the last equation of the system (27): where dT TC dt and T TC are determined from experimental data. The resultant time constant of thermocouple as a function of corrected rotation speed is shown in Figure 8.

Starting simulation algorithm
Starting simulation algorithm implements the structure, which is presented in Figure 1, as the following sequence of stages: 1. Initial values of parameters are set as N HP cor ¼    (27)) for the current step of the solution procedure are calculated: 7. The starter torsion torque M start is determined. The model of starter (see Figure 2) contains characteristics of the starter, which determines the starter torsion torque as a function of the starter rotation speed N start , pressure and temperature of the air at the starter inlet P air , T air , and the air pressure at the starter nozzle vanes P NB . The moment of the starter switching off is determined according to the starting sequence diagram. After this moment, M start is set to be zero.
Points 2-9 of this algorithm are repeated up to the end of the starting sequence diagram.

Starting model verification
For adequacy of the model and its quality checking, the simulation results were compared with experimental data. For this purpose, the starting process was chosen, which was not used for the model identification and differs from other experiments by a law of the fuel supplying and by ambient conditions. The experimental fuel flow (see Figure 9a) was inputted to the model. The time moments of the combustion ignition and the starter switching off are the same as at experiments, which were used for the model identification. The results are shown in Figure 9b-f.
Comparison of the simulation results with experimental data shows that the model meets all demands to models, which are used for the ACS development. Visible difference between experimental and simulated data is explained by the fact that the real combustion chamber ignition was some earlier than it is declared in the technical requirements to ACS, and the model implements these requirements.

Conclusion
The section discussed a new method for the synthesis of the mathematical model of engine operation in the sub-idle region. As shown, this model belongs to the class of dynamic models. To ensure the unity of the structure with the model of operating modes, a simplified static model is introduced into the dynamic model, which takes into account the synthesized characteristics of the engine components in the format of direct dependencies between the parameter-argument and the parameters used for the synthesis of ACS.