Turbine Engine Starting Simulation

Abstract

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.

1. 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, 3];

• “smart” extrapolation that consumes the theoretical knowledge about component performances in the sub-idle region [4, 5, 6, 7];

• statistical extrapolation [8];

• based on fitting the generic performances to the experimental data [9, 10, 11].

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:

$$\frac{{\left({\mathrm{W}}_{\mathrm{a}}\right)}_{\mathrm{i}-1}}{{\left({\mathrm{W}}_{\mathrm{a}}\right)}_{\mathrm{i}}}=\left[\frac{{\mathrm{N}}_{\mathrm{i}-1}}{{\mathrm{N}}_{\mathrm{i}}}\right],\frac{{\left({\mathrm{L}}_{\mathrm{sp}}\right)}_{\mathrm{i}-1}}{{\left({\mathrm{L}}_{\mathrm{sp}}\right)}_{\mathrm{i}}}={\left[\frac{{\mathrm{N}}_{\mathrm{i}-1}}{{\mathrm{N}}_{\mathrm{i}}}\right]}^2,{\left(\mathrm{Pow}\right)}_{\mathrm{i}-1}/{\left(\mathrm{Pow}\right)}_{\mathrm{i}}={\left[{\mathrm{N}}_{\mathrm{i}-1}/{\mathrm{N}}_{\mathrm{i}}\right]}^3.$$ E1

where subscript i refers to experimentally measured parameter and $\mathrm{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

$$\frac{{\left({\mathrm{W}}_{\mathrm{a}}\right)}_{\mathrm{i}-1}}{{\left({\mathrm{W}}_{\mathrm{a}}\right)}_{\mathrm{i}}}={\left[\frac{{\mathrm{N}}_{\mathrm{i}-1}}{{\mathrm{N}}_{\mathrm{i}}}\right]}^{\mathrm{p}},\frac{{\left({\mathrm{L}}_{\mathrm{sp}}\right)}_{\mathrm{i}-1}}{{\left({\mathrm{L}}_{\mathrm{sp}}\right)}_{\mathrm{i}}}={\left[\frac{{\mathrm{N}}_{\mathrm{i}-1}}{{\mathrm{N}}_{\mathrm{i}}}\right]}^{\mathrm{q}};{\left(\mathrm{Pow}\right)}_{\mathrm{i}-1}/{\left(\mathrm{Pow}\right)}_{\mathrm{i}}={\left[\frac{{\mathrm{N}}_{\mathrm{i}-1}}{{\mathrm{N}}_{\mathrm{i}}}\right]}^{\mathrm{r}}.$$ E2

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 [5]), but there are still no known methods that are able to take into account all prior information. For example, there is a method [6, 7] 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:

$${\mathrm{\eta }}_{\mathrm{C}}^{\ast }\cdot \left[{\mathrm{L}}_{\mathrm{C}}^{\ast }/{\mathrm{N}}^2\right]=\mathrm{f}\left({\mathrm{W}}_{\mathrm{a}}\cdot {\mathrm{T}}_{\mathit{\text{in}}}^{\ast }/\mathrm{N}\cdot {\mathrm{P}}_{\mathit{\text{in}}}^{\ast }\right).$$ E3

Another method, proposed by Agrawal and Yunis [8], 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:

• Air flow

$$\frac{{\left({\mathrm{W}}_{\mathrm{a}}\cdot \sqrt{{\mathrm{T}}_{\mathrm{cor}}}/{\mathrm{P}}_{\mathrm{cor}}\right)}_{\mathrm{i}-1}}{{\left({\mathrm{W}}_{\mathrm{a}}\cdot \sqrt{{\mathrm{T}}_{\mathrm{cor}}}/{\mathrm{P}}_{\mathrm{cor}}\right)}_{\mathrm{i}}}={\mathrm{K}}_{\mathrm{\varphi }}\cdot \frac{{\left(\mathrm{N}/\sqrt{{\mathrm{T}}_{\mathrm{cor}}}\right)}_{\mathrm{i}-1}}{{\left(\mathrm{N}/\sqrt{{\mathrm{T}}_{\mathrm{cor}}}\right)}_{\mathrm{i}}};$$ E4

• Specific work

$$\frac{{\left({\mathrm{L}}_{\mathrm{sp}}/{\mathrm{T}}_{\mathrm{cor}}\right)}_{\mathrm{i}-1}}{{\left({\mathrm{L}}_{\mathrm{sp}}/{\mathrm{T}}_{\mathrm{cor}}\right)}_{\mathrm{i}}}={\mathrm{K}}_{\mathrm{\psi }}\cdot {\left[\frac{{\left(\mathrm{N}/\sqrt{{\mathrm{T}}_{\mathrm{cor}}}\right)}_{\mathrm{i}-1}}{{\left(\mathrm{N}/\sqrt{{\mathrm{T}}_{\mathrm{cor}}}\right)}_{\mathrm{i}}}\right]}^2;$$ E5

• Efficiency

$$\frac{{\left({\mathrm{\eta }}_{\mathrm{C}}\right)}_{\mathrm{i}-1}}{{\left({\mathrm{\eta }}_{\mathrm{C}}\right)}_{\mathrm{i}}}=\mathrm{K}\cdot {\mathrm{\eta }}_{\mathrm{C}\mathrm{i}}\cdot \left[\frac{{\left(\mathrm{N}/\sqrt{{\mathrm{T}}_{\mathrm{cor}}}\right)}_{\mathrm{i}-1}}{{\left(\mathrm{N}/\sqrt{{\mathrm{T}}_{\mathrm{cor}}}\right)}_{\mathrm{i}}}\right],$$ E6

where ${\mathrm{K}}_{\mathrm{\varphi }}$ , ${\mathrm{K}}_{\mathrm{\psi }}$ , 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 [9]. 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 [12].

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 [13]:

$${\mathrm{W}}_{\mathrm{a}\mathrm{cor}}={\left({\mathrm{W}}_{\mathrm{a}\mathrm{cor}}\right)}_{\mathit{\text{univ}}}\cdot {\mathrm{W}}_{\mathrm{a}\mathrm{cor}\text{calc}}/{\left({\mathrm{W}}_{\mathrm{a}\mathrm{cor}\text{calc}}\right)}_{\mathit{\text{univ}}},$$ E7

$${\mathrm{N}}_{\mathrm{cor}}={\left({\mathrm{N}}_{\mathrm{cor}}\right)}_{\mathit{\text{univ}}}\cdot {\mathrm{N}}_{\mathrm{cor}\text{calc}}/{\left({\mathrm{N}}_{\mathrm{cor}\text{calc}}\right)}_{\mathit{\text{univ}}},$$ E8

$$\mathrm{\pi }={\mathrm{\pi }}_{\mathit{\text{univ}}}\cdot \left({\mathrm{\pi }}_{\text{calc}}-1\right)+1/\left({\mathrm{\pi }}_{\text{calc univ}}-1\right)+1,$$ E9

$$\mathrm{\eta }={\mathrm{\eta }}_{\mathit{\text{univ}}}\cdot {\mathrm{\eta }}_{\text{calc}}/{\mathrm{\eta }}_{\text{calc univ}}.$$ E10

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 [14].

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.

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

$${\mathrm{M}}_{\mathit{\text{starter}}}={\mathrm{M}}_{\mathrm{res}}+{\mathrm{M}}_{\mathit{\text{inert}}},$$ E11

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

$${\mathrm{M}}_{\mathit{\text{starter}}}+{\mathrm{M}}_{\mathrm{T}}={\mathrm{M}}_{\mathrm{res}}+{\mathrm{M}}_{\mathit{\text{inert}}}.$$ E12

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

Figure 1.

Engine starting stages and moment performances.

$${\mathrm{M}}_{\mathrm{T}}={\mathrm{M}}_{\mathrm{res}}+{\mathrm{M}}_{\mathit{\text{inert}}}.$$ E13

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.

3. 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,

$$\left\{\begin{array}{l}\dot{\overrightarrow{\mathrm{X}}}=\mathrm{A}\cdot \mathrm{\Delta }\overrightarrow{\mathrm{X}}+\mathrm{B}\cdot \mathrm{\Delta }\overrightarrow{\mathrm{U}};\\ \mathrm{\Delta }\overrightarrow{\mathrm{Y}}=\mathrm{C}\cdot \mathrm{\Delta }\overrightarrow{\mathrm{Y}}+\mathrm{D}\cdot \mathrm{\Delta }\overrightarrow{\mathrm{U}}.\end{array}\right.$$ E14

supported with steady-state model in simplified representation (next—simplified static model, SSM).

The structure of the model is presented in a Figure 2.

Figure 2.

Structure of the proposed starting model of a turbofan.

The model of the three-spool turbofan must consider the mutual interaction of engine spools.

4. 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 $\overrightarrow{\mathrm{Z}}=\mathrm{f}\left(\overline{\mathrm{N}}\right)$ , where $\overrightarrow{\mathrm{Z}}={\left[{\overrightarrow{\mathrm{X}}}^{\mathrm{T}}{\overrightarrow{\mathrm{Y}}}^{\mathrm{T}}\right]}^{\mathrm{T}}$ ; $\overline{\mathrm{N}}=\frac{{\mathrm{N}}_{\mathrm{HP}\mathrm{cor}}}{{\mathrm{N}}_{\mathrm{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 sub-idle region is presented in Figure 3.

Figure 3.

Static model generation: 1—linear segment; 2—parabola; 3—working line.

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

$${\mathrm{Z}}^{\left(1\right)}={\mathrm{Z}}_0+\mathrm{k}\overline{\mathrm{N}}.$$ E15

Initial values of slope angle are evaluated from the next formula

$$\mathrm{k}=\left(0.6\dots 1.2\right)\frac{{\mathrm{Z}}_2-{\mathrm{Z}}_0}{{\overline{\mathrm{N}}}_2}.$$ E16

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₀), ( ${\overline{\mathrm{N}}}_2$ , Z₂). Y₀ is equal to a default value of a parameter to be simulated (temperature, pressure, RPM, etc.).

The linking function is described as

$${\mathrm{Z}}^{\left(2\right)}={\mathrm{k}}_0+{\mathrm{k}}_1\overline{\mathrm{N}}+{\mathrm{k}}_2{\overline{\mathrm{N}}}^2.$$ E17

Unknown coefficients k₀, k₁, and k₂ are evaluated from the set of linear equations that is compiled from the limitations for the regions, described in this chapter:

$${\mathrm{Z}}_0+\mathrm{k}{\overline{\mathrm{N}}}_1={\mathrm{k}}_0+{\mathrm{k}}_1{\overline{\mathrm{N}}}_1+{\mathrm{k}}_2{\overline{\mathrm{N}}}_1^2;$$ E18

$$\mathrm{k}={\mathrm{k}}_1+2{\mathrm{k}}_2{\overline{\mathrm{N}}}_1;$$ E19

$${\mathrm{k}}_0+{\mathrm{k}}_1{\overline{\mathrm{N}}}_2+{\mathrm{k}}_2{\overline{\mathrm{N}}}_2^2={\mathrm{Z}}_2;$$ E20

$${\mathrm{k}}_1+2{\mathrm{k}}_2{\overline{\mathrm{N}}}_2={\left(\frac{\mathrm{dZ}}{\mathrm{d}\overline{\mathrm{N}}}\right)}_2.$$ E21

The last unknown to be determined from the equation set (18)–(21) is an abscissa ${\overline{\mathrm{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₀ and k₁ to be dependent on k₂:

$${\mathrm{k}}_1={\left(\frac{\mathrm{dZ}}{\mathrm{d}\overline{\mathrm{N}}}\right)}_2-2{\mathrm{k}}_2{\overline{\mathrm{N}}}_2; {\mathrm{k}}_0={\mathrm{Z}}_2-{\left(\frac{\mathrm{dZ}}{\mathrm{d}\overline{\mathrm{N}}}\right)}_2{\overline{\mathrm{N}}}_2+{\mathrm{k}}_2{\overline{\mathrm{N}}}_2^2.$$ E22

From Eq. (19), we get the relation between ${\overline{\mathrm{N}}}_1$ and ${\mathrm{k}}_2$ :

$${\overline{\mathrm{N}}}_1=\frac{\mathrm{k}-\mathrm{D}+2{\mathrm{k}}_2{\overline{\mathrm{N}}}_2}{2{\mathrm{k}}_2}.$$ E23

Putting the obtained equations into (18), we get the solution:

$${\mathrm{k}}_2=\frac{{\left(\mathrm{k}-{\left(\frac{\mathrm{dZ}}{\mathrm{d}\overline{\mathrm{N}}}\right)}_2\right)}^2}{4\left({\mathrm{Z}}_2-{\mathrm{Z}}_0-{\mathrm{k}}_2{\overline{\mathrm{N}}}_2\right)};{\mathrm{k}}_1={\left(\frac{\mathrm{dZ}}{\mathrm{d}\overline{\mathrm{N}}}\right)}_2-\frac{{\left(\mathrm{k}-{\left(\frac{\mathrm{dZ}}{\mathrm{d}\overline{\mathrm{N}}}\right)}_2\right)}^2{\overline{\mathrm{N}}}_2}{2\left({\mathrm{Z}}_2-{\mathrm{Z}}_0-{\mathrm{k}}_2{\overline{\mathrm{N}}}_2\right)};$$ E24

$${\mathrm{k}}_0={\mathrm{Z}}_2-{\left(\frac{\mathrm{dZ}}{\mathrm{d}\overline{\mathrm{N}}}\right)}_2{\overline{\mathrm{N}}}_2+\frac{{\left(\mathrm{k}-{\left(\frac{\mathrm{dZ}}{\mathrm{d}\overline{\mathrm{N}}}\right)}_2\right)}^2{\overline{\mathrm{N}}}_2^2}{4\left({\mathrm{Z}}_2-{\mathrm{Z}}_0-{\mathrm{k}}_2{\overline{\mathrm{N}}}_2\right)};{\overline{\mathrm{N}}}_1=\frac{2\left({\mathrm{Z}}_2-{\mathrm{Z}}_0-{\mathrm{k}}_2{\overline{\mathrm{N}}}_2\right)}{\left(\mathrm{k}-{\left(\frac{\mathrm{dZ}}{\mathrm{d}\overline{\mathrm{N}}}\right)}_2\right)}+{\overline{\mathrm{N}}}_2.$$ E25

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.

Figure 4.

Relation ${\mathit{W}}_{\mathit{f}\mathit{st}\mathit{cor}}=\mathit{f}\left({\mathit{N}}_{\mathit{HP}\mathit{cor}}\right)$ of the SSM.

Because the corrected exhaust gas temperature ${\mathrm{T}}_{\mathrm{T}\mathrm{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:

$${\mathrm{Т}}_{4\mathrm{cor}}={\mathrm{Т}}_{2\mathrm{cor}}+{\mathrm{K}}_{\mathrm{T}}\frac{{\mathrm{W}}_{\mathrm{f}\mathrm{cor}}\cdot {\mathrm{H}}_{\mathrm{u}}}{{\mathrm{W}}_{\mathrm{air}\mathrm{cor}}\cdot {\mathrm{с}}_{\mathrm{p}}},$$ E26

where ${\mathrm{Т}}_{2\mathrm{cor}}$ is the temperature of airflow at compressor discharge; $\frac{{\mathrm{W}}_{\mathrm{f}\mathrm{cor}}\cdot {\mathrm{H}}_{\mathrm{u}}}{{\mathrm{W}}_{\mathrm{air}\mathrm{cor}}\cdot {\mathrm{с}}_{\mathrm{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: ${\mathrm{W}}_{\mathrm{air}\mathrm{cor}}=\mathrm{C}\cdot {\mathrm{N}}_{\mathrm{HP}}^2$ .

The dependence ${\mathrm{Т}}_{4\mathrm{cor}}=\mathrm{f}\left({\mathrm{N}}_{\mathrm{HP}\mathrm{cor}}\right)$ 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.

Figure 5.

Static model of the gas temperature.

It follows from Eq. (10) that for proper modeling of gas temperature on must compile the SSMs of ${\mathrm{Т}}_{2\mathrm{cor}}$ and airflow ${\mathrm{W}}_{\mathrm{air}\mathrm{cor}}$ .

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

$$\left\{\begin{array}{l}{\dot{\mathrm{N}}}_{\mathrm{IP}}={\mathrm{a}}_{11}\cdot \mathrm{\Delta }{\mathrm{N}}_{\mathrm{IP}}+{\mathrm{a}}_{12}\cdot \mathrm{\Delta }{\mathrm{N}}_{\mathrm{HP}}+{\mathrm{a}}_{32}\cdot \mathrm{\Delta }{\mathrm{N}}_{\mathrm{LP}}+{\mathrm{b}}_{11}\cdot \mathrm{\Delta }{\mathrm{W}}_{\mathrm{f}};\\ {\dot{\mathrm{N}}}_{\mathrm{HP}}={\mathrm{a}}_{21}\cdot \mathrm{\Delta }{\mathrm{N}}_{\mathrm{IP}}+{\mathrm{a}}_{22}\cdot \mathrm{\Delta }{\mathrm{N}}_{\mathrm{HP}}+{\mathrm{a}}_{23}\cdot \mathrm{\Delta }{\mathrm{N}}_{\mathrm{LP}}+{\mathrm{b}}_{21}\cdot \mathrm{\Delta }{\mathrm{W}}_{\mathrm{f}}+{\mathrm{b}}_{24}\cdot \left({\mathrm{М}}_{\mathit{\text{start}}}-{\mathrm{М}}_{\mathrm{res}}\right);\\ {\dot{\mathrm{N}}}_{\mathrm{LP}}={\mathrm{a}}_{31}\cdot \mathrm{\Delta }{\mathrm{N}}_{\mathrm{IP}}+{\mathrm{a}}_{32}\cdot \mathrm{\Delta }{\mathrm{N}}_{\mathrm{HP}}+{\mathrm{a}}_{33}\cdot \mathrm{\Delta }{\mathrm{N}}_{\mathrm{LP}}+{\mathrm{b}}_{31}\cdot \mathrm{\Delta }{\mathrm{W}}_{\mathrm{f}};\\ \mathrm{\Delta }{\mathrm{p}}_2={\mathrm{c}}_{11}\cdot \mathrm{\Delta }{\mathrm{N}}_{\mathrm{IP}}+{\mathrm{c}}_{12}\cdot \mathrm{\Delta }{\mathrm{N}}_{\mathrm{HP}}+{\mathrm{с}}_{13}\cdot \mathrm{\Delta }{\mathrm{N}}_{\mathrm{LP}}+{\mathrm{d}}_{11}\cdot \mathrm{\Delta }{\mathrm{W}}_{\mathrm{f}};\\ \mathrm{\Delta }{\mathrm{Т}}_4={\mathrm{c}}_{21}\cdot \mathrm{\Delta }{\mathrm{N}}_{\mathrm{IP}}+{\mathrm{c}}_{22}\cdot \mathrm{\Delta }{\mathrm{N}}_{\mathrm{HP}}+{\mathrm{с}}_{23}\cdot \mathrm{\Delta }{\mathrm{N}}_{\mathrm{LP}}+{\mathrm{d}}_{21}\cdot \mathrm{\Delta }{\mathrm{W}}_{\mathrm{f}},\end{array}\right.$$ E27

where $\mathrm{\Delta }{\mathrm{N}}_{\mathrm{IP}}={\mathrm{N}}_{\mathrm{IP}}-{\mathrm{N}}_{\mathrm{IP}\mathrm{st}}$ , $\mathrm{\Delta }{\mathrm{N}}_{\mathrm{HP}}={\mathrm{N}}_{\mathrm{HP}}-{\mathrm{N}}_{\mathrm{HP}\mathrm{st}}$ , $\mathrm{\Delta }{\mathrm{N}}_{\mathrm{LP}}={\mathrm{N}}_{\mathrm{LP}}-{\mathrm{N}}_{\mathrm{LP}\mathrm{st}}$ , $\mathrm{\Delta }{\mathrm{W}}_{\mathrm{f}}={\mathrm{W}}_{\mathrm{f}}-{\mathrm{W}}_{\mathrm{f}\mathrm{st}}$ , $\mathrm{\Delta }{\mathrm{p}}_2={\mathrm{p}}_2-{\mathrm{p}}_{2\mathrm{st}}$ and $\mathrm{\Delta }{\mathrm{Т}}_4={\mathrm{Т}}_4-{\mathrm{Т}}_{4\mathrm{st}}$ .

This system can be resolved by the Euler method.

The set of Eqs. (27) can be simplified by the introduction of the listed next assumptions:

• 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 ${\mathrm{N}}_{\mathrm{HP}}={\mathrm{N}}_{\mathrm{HP}\mathrm{st}}$ which is why $\mathrm{\Delta }{\mathrm{N}}_{\mathrm{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 intermediate-pressure rotor affects the low-pressure rotor).

Concluding the just listed assumptions, the set of Eqs. (27) transforms to

$$\left\{\begin{array}{l}{\dot{\mathrm{N}}}_{\mathrm{IP}}={\mathrm{a}}_{11}\cdot \mathrm{\Delta }{\mathrm{N}}_{\mathrm{IP}}+{\mathrm{b}}_{11}\cdot \mathrm{\Delta }{\mathrm{W}}_{\mathrm{f}};\\ {\dot{\mathrm{N}}}_{\mathrm{HP}}={\mathrm{b}}_{21}\cdot \mathrm{\Delta }{\mathrm{W}}_{\mathrm{f}}+{\mathrm{b}}_{24}\cdot \left({\mathrm{М}}_{\mathit{\text{start}}}-{\mathrm{M}}_{\mathrm{res}}\right);\\ {\dot{\mathrm{N}}}_{\mathrm{LP}}={\mathrm{a}}_{31}\cdot \mathrm{\Delta }{\mathrm{N}}_{\mathrm{IP}}+{\mathrm{a}}_{33}\cdot \mathrm{\Delta }{\mathrm{N}}_{\mathrm{HP}}+{\mathrm{b}}_{31}\cdot \mathrm{\Delta }{\mathrm{W}}_{\mathrm{f}};\\ \mathrm{\Delta }{\mathrm{p}}_2={\mathrm{d}}_{11}\cdot \mathrm{\Delta }{\mathrm{W}}_{\mathrm{f}};\\ \mathrm{\Delta }{\mathrm{Т}}_4={\mathrm{d}}_{21}\cdot \mathrm{\Delta }{\mathrm{W}}_{\mathrm{f}}.\end{array}\right.$$ E28

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₂₄, b₂₁, and d₁₁ 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₂₄ was determined analytically from the following equation:

$${\mathrm{J}}_{\mathrm{HP}}\frac{\mathrm{d}{\mathrm{\omega }}_{\mathrm{HP}}}{\mathrm{dt}}={\mathrm{M}}_{\mathit{\text{start}}}-{\mathrm{М}}_{\mathrm{HPC}}\left({\mathrm{N}}_{\mathrm{HP}}^2\right),$$ E29

where J_HP is an inertia of high pressure rotor.

Having transformed $\frac{\mathrm{d}{\mathrm{\omega }}_{\mathrm{HP}}}{\mathrm{dt}}$ to $\frac{\mathrm{\pi }}{30}\cdot \frac{{\mathrm{dN}}_{\mathrm{HP}}}{\mathrm{dt}}$ , we get:

$${\mathrm{b}}_{24}=\frac{30}{\mathrm{\pi }\cdot {\mathrm{J}}_{\mathrm{HP}}}$$ E30

The coefficient b₂₄ does not depend on the type or characteristics of the used starter and is constant for the entire starting.

The determination of the coefficient ${\mathrm{b}}_{21\mathrm{пр}}$ 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:

$${\mathrm{b}}_{21\mathrm{cor}}=\frac{{\dot{\mathrm{N}}}_{\mathrm{HP}}-{\mathrm{b}}_{24}\cdot {\mathrm{M}}_{\mathit{\text{start}}}}{\mathrm{\Delta }{\mathrm{W}}_{\mathrm{f}}}.$$ E31

Further, the obtained values were approximated, considering the following a-priori information about this coefficient at the starting region. Let us consider the following equation:

$${\mathrm{\tau }}_{\mathrm{HP}}\frac{{\mathrm{dN}}_{\mathrm{HP}}}{\mathrm{dt}}+\mathrm{\Delta }{\mathrm{N}}_{\mathrm{HP}}={\mathrm{K}}_{\mathrm{W}}\cdot \mathrm{\Delta }{\mathrm{W}}_{\mathrm{f}},$$ E32

where ${\mathrm{\tau }}_{\mathrm{HP}}$ is the time constant of the high-pressure rotor.

Let us divide both parts of this equation by time constant and move ( $\frac{1}{{\mathrm{\tau }}_{\mathrm{HP}}}\mathrm{\Delta }{\mathrm{N}}_{\mathrm{HP}}$ ) to the right part of the equation, and we get

$${\dot{\mathrm{N}}}_{\mathrm{HP}}=-\frac{1}{{\mathrm{\tau }}_{\mathrm{HP}}}\mathrm{\Delta }{\mathrm{N}}_{\mathrm{HP}}+\frac{{\mathrm{K}}_{\mathrm{W}}}{{\mathrm{\tau }}_{\mathrm{HP}}}\cdot \mathrm{\Delta }{\mathrm{W}}_{\mathrm{f}}={\mathrm{a}}_{22}\mathrm{\Delta }{\mathrm{N}}_{\mathrm{HP}}+{\mathrm{b}}_{21}\mathrm{\Delta }{\mathrm{W}}_{\mathrm{f}}.$$ E33

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 ( $\frac{\mathrm{\partial }{\mathrm{М}}_{\mathrm{HPC}}}{\mathrm{\partial }{\mathrm{n}}_{\mathrm{HP}}}=\frac{\mathrm{\partial }{\mathrm{М}}_{\mathrm{HPT}}}{\mathrm{\partial }{\mathrm{n}}_{\mathrm{HP}}}$ ). At this point, the function of the dependence of the time constant on the rotor speed experiences a gap ( ${\mathrm{\tau }}_{\mathrm{HP}}=\frac{\frac{\mathrm{\pi }}{30}·{\mathrm{J}}_{\mathrm{HP}}}{\frac{\mathrm{\partial }{\mathrm{М}}_{\mathrm{HPC}}}{\mathrm{\partial }{\mathrm{N}}_{\mathrm{HP}}}-\frac{\mathrm{\partial }{\mathrm{М}}_{\mathrm{HPТ}}}{\mathrm{\partial }{\mathrm{N}}_{\mathrm{HP}}}}=\infty$ ). 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.

Figure 6.

Coefficients of the high-pressure rotor dynamics.

The analytical determination of the coefficients a₁₁, b₁₁, a₃₁, a₃₃, and b₃₁ 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₁₁ and b₁₁ is.

$$\widehat{\overrightarrow{\mathrm{\theta }}}=argmin\mathrm{J}\left(\overrightarrow{\mathrm{\theta }}\right);$$ E34

where $\overrightarrow{\mathrm{\theta }}=\left[\begin{array}{l}{\mathrm{a}}_{11}\\ {\mathrm{b}}_{11}\end{array}\right]$ is a vector of coefficients to be identified (a₁₁, b₁₁), the change of which is approximated with polynomials ${\mathrm{b}}_{11\mathrm{cor}}={\sum }_{\mathrm{i}=0}^2 {\mathrm{q}}_{\mathrm{i}}\cdot {{\mathrm{N}}^{\mathrm{i}}}_{\mathrm{HP}\mathrm{cor}}$ and ${\mathrm{a}}_{11\mathrm{cor}}={\sum }_{\mathrm{i}=0}^2 {\mathrm{p}}_{\mathrm{i}}\cdot {{\mathrm{N}}^{\mathrm{i}}}_{\mathrm{HP}\mathrm{cor}}$ , $\widehat{\overrightarrow{\mathrm{\theta }}}$ —vector estimate $\overrightarrow{\mathrm{\theta }}$ .

Thus, the problem of the evaluation of coefficients a₁₁, b₁₁ reduces to the problem of determining the coefficients of polynomials approximating the change of these coefficients in the sub-idle region:

$$\begin{array}{l}\mathrm{J}\left(\overrightarrow{\mathrm{\theta }}\right)=\mathrm{J}\left({\mathrm{q}}_2,{\mathrm{q}}_1,{\mathrm{q}}_0,{\mathrm{p}}_2,{\mathrm{p}}_1,{\mathrm{p}}_0\right)\\ =\sum _{\mathrm{j}=1}^{\mathrm{N}}\left[{\dot{\mathrm{N}}}_{\mathrm{IP}\mathrm{cor}\mathrm{j}}\right.-\left\{\left({\mathrm{q}}_2{{\mathrm{N}}^2}_{\mathrm{HP}\mathrm{cor}\mathrm{j}}+{\mathrm{q}}_1{\mathrm{N}}_{\mathrm{HP}\mathrm{cor}\mathrm{j}}+{\mathrm{q}}_0\right)\mathrm{\Delta }{\mathrm{W}}_{\mathrm{f}\mathrm{cor}\mathrm{j}}\right.\\ +{\left.\left.\left({\mathrm{p}}_2{{\mathrm{N}}^2}_{\mathrm{HP}\mathrm{cor}\mathrm{j}}+{\mathrm{p}}_1{\mathrm{N}}_{\mathrm{HP}\mathrm{cor}\mathrm{j}}+{\mathrm{p}}_0\right)\mathrm{\Delta }{\mathrm{N}}_{\mathrm{IP}\mathrm{cor}\mathrm{j}}\right\}\right]}^2.\end{array}$$ E35

As an example, Figure 7 represents the results of a₁₁ approximation on the engine operational mode. They correspond to the obtained dependence

Figure 7.

Relation of coefficient a₁₁ on operation mode.

$${\mathrm{a}}_{11\mathrm{cor}}=-4.0276\cdot {10}^{-9}{{\mathrm{N}}^2}_{\mathrm{HP}\mathrm{cor}\mathrm{j}}-3.0863\cdot {10}^{-5}{\mathrm{N}}_{\mathrm{HP}\mathrm{cor}\mathrm{j}}+0.0716.$$ E36

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:

$${\mathrm{b}}_{11}=\frac{{\dot{\mathrm{N}}}_{\mathrm{IP}exp}-{\mathrm{a}}_{11}\cdot \left({\mathrm{N}}_{\mathrm{IP}exp}-{\mathrm{N}}_{\mathrm{IP}\mathrm{st}}\right)}{{\mathrm{W}}_{\mathrm{f}exp}-{\mathrm{W}}_{\mathrm{f}\mathrm{st}}}.$$ E37

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

$$\mathrm{\epsilon }\frac{{\mathrm{dT}}_{\mathrm{TC}}}{\mathrm{dt}}+{\mathrm{T}}_{\mathrm{TC}}={\mathrm{T}}_4,$$ E38

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):

$$\mathrm{\epsilon }=\frac{{\mathrm{c}}_{21}\cdot \mathrm{\Delta }{\mathrm{N}}_{\mathrm{IP}}+{\mathrm{c}}_{22}\cdot \mathrm{\Delta }{\mathrm{N}}_{\mathrm{HP}}+{\mathrm{с}}_{23}\cdot \mathrm{\Delta }{\mathrm{N}}_{\mathrm{LP}}+{\mathrm{d}}_{21}\cdot \mathrm{\Delta }{\mathrm{W}}_{\mathrm{f}}-{\mathrm{T}}_{\mathrm{TC}}}{\frac{{\mathrm{dT}}_{\mathrm{TC}}}{\mathrm{dt}}},$$ E39

where $\frac{{\mathrm{dT}}_{\mathrm{TC}}}{\mathrm{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.

Figure 8.

Estimated time constant of thermocouple as a function of rotation speed.

7. 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 ${\mathrm{N}}_{\mathrm{HP}\mathrm{cor}}={\mathrm{N}}_{\mathrm{IP}\mathrm{cor}}={\mathrm{N}}_{\mathrm{LP}\mathrm{cor}}=0$ , ${\mathrm{W}}_{\mathrm{f}}=0$ , ${\mathrm{P}}_2={\mathrm{P}}_{\mathit{\text{in}}}$ , ${\mathrm{T}}_4={\mathrm{Т}}_2={\mathrm{T}}_{\mathrm{TC}}={\mathrm{T}}_{\mathit{\text{in}}}$ , ${\mathrm{W}}_{\mathrm{air}\mathrm{cor}}=0$ .

2. Knowing ${\mathrm{N}}_{\mathrm{HP}\mathrm{cor}}$ , ${\mathrm{P}}_{\mathit{\text{in}}}$ , ${\mathrm{T}}_{\mathit{\text{in}}}$ (in first step, the initial values are used) and applying the static model ( ${\mathrm{W}}_{\mathrm{f}\mathrm{st}\mathrm{cor}}=\mathrm{f}\left({\mathrm{N}}_{\mathrm{HP}\mathrm{cor}}\right)$ , ${\mathrm{N}}_{\mathrm{IP}\mathrm{st}\mathrm{cor}}=\mathrm{f}\left({\mathrm{N}}_{\mathrm{HP}\mathrm{cor}}\right)$ , ${\mathrm{N}}_{\mathrm{LP}\mathrm{st}\mathrm{cor}}=\mathrm{f}\left({\mathrm{N}}_{\mathrm{HP}\mathrm{cor}}\right)$ , ${\mathrm{P}}_{2\mathrm{st}\mathrm{cor}}=\mathrm{f}\left({\mathrm{N}}_{\mathrm{HP}\mathrm{cor}}\right)$ , ${\mathrm{T}}_{2\mathrm{st}\mathrm{cor}}=\mathrm{f}\left({\mathrm{N}}_{\mathrm{HP}\mathrm{cor}}\right)$ , ${\mathrm{W}}_{\mathrm{air}\mathrm{st}\mathrm{cor}}=\mathrm{f}\left({\mathrm{N}}_{\mathrm{HP}\mathrm{cor}}\right)$ , ${\mathrm{Т}}_{4\mathrm{st}\mathrm{cor}}=\mathrm{f}\left({\mathrm{N}}_{\mathrm{HP}\mathrm{cor}}\right)$ ), static values of the simulated parameters are calculated: ${\mathrm{W}}_{\mathrm{f}\mathrm{st}\mathrm{cor}}$ , ${\mathrm{N}}_{\mathrm{IP}\mathrm{st}\mathrm{cor}}$ , ${\mathrm{N}}_{\mathrm{LP}\mathrm{st}\mathrm{cor}}$ , ${\mathrm{P}}_{2\mathrm{st}\mathrm{cor}}$ , ${\mathrm{T}}_{2\mathrm{st}\mathrm{cor}}$ , ${\mathrm{W}}_{\mathrm{air}\mathrm{st}\mathrm{cor}}$ , ${\mathrm{Т}}_{4\mathrm{st}\mathrm{cor}}$ .

3. Physical values of the parameters ${\mathrm{W}}_{\mathrm{f}\mathrm{st}}$ , ${\mathrm{N}}_{\mathrm{IP}\mathrm{st}}$ , ${\mathrm{N}}_{\mathrm{LP}\mathrm{st}}$ , ${\mathrm{P}}_{2\mathrm{st}}$ , ${\mathrm{T}}_{2\mathrm{st}}$ , ${\mathrm{W}}_{\mathrm{air}\mathrm{st}}$ , ${\mathrm{Т}}_{4\mathrm{st}}$ are calculated using equations of correction.

4. Knowing ${\mathrm{N}}_{\mathrm{HP}\mathrm{cor}}$ , ${\mathrm{P}}_{\mathit{\text{in}}}$ , ${\mathrm{T}}_{\mathit{\text{in}}}$ and using the relations ${\mathrm{a}}_{11\mathrm{cor}}=\mathrm{f}\left({\mathrm{N}}_{\mathrm{HP}\mathrm{cor}}\right)$ , ${\mathrm{b}}_{11\mathrm{cor}}=\mathrm{f}\left({\mathrm{N}}_{\mathrm{HP}\mathrm{cor}}\right)$ , ${\mathrm{b}}_{21\mathrm{cor}}=\mathrm{f}\left({\mathrm{N}}_{\mathrm{HP}\mathrm{cor}}\right)$ , ${\mathrm{a}}_{31\mathrm{cor}}=\mathrm{f}\left({\mathrm{N}}_{\mathrm{HP}\mathrm{cor}}\right)$ , ${\mathrm{a}}_{33\mathrm{cor}}=\mathrm{f}\left({\mathrm{N}}_{\mathrm{HP}\mathrm{cor}}\right)$ , ${\mathrm{b}}_{31\mathrm{cor}}=\mathrm{f}\left({\mathrm{N}}_{\mathrm{HP}\mathrm{cor}}\right)$ , ${\mathrm{d}}_{11\mathrm{cor}}=\mathrm{f}\left({\mathrm{N}}_{\mathrm{HP}\mathrm{cor}}\right)$ , ${\mathrm{d}}_{21\mathrm{cor}}=\mathrm{f}\left({\mathrm{N}}_{\mathrm{HP}\mathrm{cor}}\right)$ , the coefficients of the linear dynamic model are determined for a current step of the solution procedure: ${\mathrm{a}}_{11\mathrm{cor}}$ , ${\mathrm{b}}_{11\mathrm{cor}}$ , ${\mathrm{b}}_{21\mathrm{cor}}$ , ${\mathrm{a}}_{31\mathrm{cor}}$ , ${\mathrm{a}}_{33\mathrm{cor}}$ , ${\mathrm{b}}_{31\mathrm{cor}}$ , ${\mathrm{d}}_{11\mathrm{cor}}$ , ${\mathrm{d}}_{21\mathrm{cor}}$ .

5. Using the equations of correction, physical values of the LM coefficients are determined: ${\mathrm{a}}_{11}$ , ${\mathrm{b}}_{11}$ , ${\mathrm{b}}_{21}$ , ${\mathrm{a}}_{31}$ , ${\mathrm{a}}_{33}$ , ${\mathrm{b}}_{31}$ , ${\mathrm{d}}_{11}$ , ${\mathrm{d}}_{21}$ .

6. Deviations $\mathrm{\Delta }\mathrm{U},\mathrm{\Delta }\overrightarrow{\mathrm{X}}$ (see Eq. (27)) for the current step of the solution procedure are calculated:

$\mathrm{\Delta }\mathrm{U}=\mathrm{\Delta }{\mathrm{W}}_{\mathrm{f}}={\mathrm{W}}_{\mathrm{f}}‐{\mathrm{W}}_{\mathrm{f}\mathrm{st}}$ ; $\mathrm{\Delta }{\mathrm{x}}_1=\mathrm{\Delta }{\mathrm{N}}_{\mathrm{IP}}={\mathrm{N}}_{\mathrm{IP}}‐{\mathrm{N}}_{\mathrm{IP}\mathrm{st}}$ ; $\mathrm{\Delta }{\mathrm{x}}_2=\mathrm{\Delta }{\mathrm{N}}_{\mathrm{HP}}=0$ ; $\mathrm{\Delta }{\mathrm{x}}_3=\mathrm{\Delta }{\mathrm{N}}_{\mathrm{LP}}={\mathrm{N}}_{\mathrm{LP}}‐{\mathrm{N}}_{\mathrm{LP}\mathrm{st}}$ .

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

2. 8. Using the quasi-linear dynamic model, new values of the simulated parameters N_HP, N_IP, N_LP, P₂, T₄, T_TC are calculated:

$$\left\{\begin{array}{l}{\mathrm{N}}_{\mathrm{IP}\mathrm{i}+1}={\mathrm{N}}_{\mathrm{IP}\mathrm{i}}+{\mathrm{t}}_{\mathit{\text{step}}}\cdot \left({\mathrm{a}}_{11}\cdot \mathrm{\Delta }{\mathrm{N}}_{\mathrm{IP}}+{\mathrm{b}}_{11}\cdot \mathrm{\Delta }{\mathrm{W}}_{\mathrm{f}}\right);\\ {\mathrm{N}}_{\mathrm{HP}\mathrm{i}+1}={\mathrm{N}}_{\mathrm{HP}\mathrm{i}}+{\mathrm{t}}_{\mathit{\text{step}}}\cdot \left({\mathrm{b}}_{21}\cdot \mathrm{\Delta }{\mathrm{W}}_{\mathrm{f}}+{\mathrm{b}}_{24}\cdot {\mathrm{М}}_{\mathit{\text{start}}}\right);\\ {\mathrm{N}}_{\mathrm{LP}\mathrm{i}+1}={\mathrm{N}}_{\mathrm{LP}\mathrm{i}}+{\mathrm{t}}_{\mathit{\text{step}}}\cdot \left({\mathrm{a}}_{31}\cdot \mathrm{\Delta }{\mathrm{N}}_{\mathrm{LP}}+{\mathrm{a}}_{33}\cdot \mathrm{\Delta }{\mathrm{N}}_{\mathrm{LP}}+{\mathrm{b}}_{31}\cdot \mathrm{\Delta }{\mathrm{W}}_{\mathrm{f}}\right);\\ {\mathrm{P}}_{2\mathrm{i}+1}={\mathrm{P}}_{2\mathrm{st}\mathrm{i}+1}+{\mathrm{d}}_{11}\cdot \mathrm{\Delta }{\mathrm{W}}_{\mathrm{f}};\\ {\mathrm{Т}}_{4\mathrm{i}+1}={\mathrm{Т}}_{4\mathrm{st}\mathrm{i}+1}+{\mathrm{d}}_{21}\cdot \mathrm{\Delta }{\mathrm{W}}_{\mathrm{f}};\\ {\mathrm{T}}_{\mathrm{ТC}\mathrm{i}+1}={\mathrm{T}}_{\mathrm{ТC}\mathrm{i}+1}+\frac{{\mathrm{t}}_{\mathit{\text{step}}}}{\mathrm{\epsilon }}\left({\mathrm{T}}_{4\mathrm{i}+1}‐{\mathrm{T}}_{\mathrm{ТC}\mathrm{i}}\right).\end{array}\right.$$ E40

1. 9. ${\mathrm{N}}_{\mathrm{HP}\mathrm{cor}}$ is determined on known N_HP.

Points 2–9 of this algorithm are repeated up to the end of the starting sequence diagram.

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

Figure 9.

Starting diagrams.

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.

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

A mathematical model of the thermocouple was also presented and it was proved to become an integral part of the starting model, due to the large inertia of temperature sensors in the sub-idle.

The starting model can be used in the development of control laws for the starting, dead cranking, and in-flight starting.

Acknowledgments

This work has been carried out with the support of the Ministry of Education and Science of Ukraine (research project no. D203-3/2019-П).

Nomenclature

c_p

specific heat capacity at constant pressure condition

J

moment of inertia

L

specific work

M

torque

N

rotation speed

P

pressure

Pow

power

T

temperature

W

air (gas) flow rate

Δ

absolute deviation

ε

thermocouple time constant

η

efficiency

π

pressure ratio

$\mathrm{\tau }$

engine time constant

Indexes a

air

C

compressor

cor

corrected parameter

f

fuel

HP

high-pressure cascade

in

inlet

IP

intermediate-pressure cascade

LP

low-pressure cascade

res

resistance

st

static

T

turbine

TC

thermocouple