Advanced Nonlinear Modeling of Gas Turbine Dynamics

Abstract

The process of gas turbine development requires different mathematical models. In particular, physics-based nonlinear dynamic models are widely used in the development of control and diagnostic systems. The present chapter firstly reviews known works on nonlinear dynamic engine modeling centering on model applications and developments. As an important development, modeling of heating up engine components is considered. This phenomenon consists in a radial clearance change during transients that influences engine static and dynamic performances. This clearance change is usually computed by a finite element method that is critical to computer resources. The chapter secondly presents a new and more rapid simulation methodology to integrate two dynamic processes, a general engine transient and a clearance change. This allows creating a more accurate and relatively fast engine dynamic model that is easy to use in the design of control and diagnostic systems. Finally, the chapter introduces further methodology enhancement consisting in the consideration of the influence of varying metal temperature on the strains induced by mechanical loads. To validate methodology, it is applied to a particular turbofan engine, and the simulated and real engine dynamic performances are compared.

1. Introduction

The gas turbine (GT) is a powerful source of energy that has relatively low size and weight. It is a principal power plant for aviation and electric energy production and has many other successful applications. For instance, aero-derivative engines are widely used for electricity generation in offshore platforms and as marine power plants because these engines are more compact and have faster dynamics than industrial GTs. A significant growth of a GT industry has been observed in the last decades [1].

1.1. Gas turbine modeling

Along with the development of new GTs, the use of mathematical modeling and simulation in the design of these engines and their systems becomes more intensive. Creating gas turbine models has been an effective design and manufacture strategy. In addition to the development of the engine itself, GT modeling and simulation have many other applications, such as design of a control system, condition monitoring, fault diagnosis, and system identification. The latter, for example, enables simulating the performances of a particular engine by model fitting to experimental data collected in test beds or at field conditions. In this way, control and diagnostic algorithms can be improved due to a more accurate individual engine model used instead of a general model [2].

The models of technical systems, in particular GT models, can be divided, on the one hand, into linear and nonlinear and, on the other hand, data-driven (also well-known as black-box models) and physics-based (also called white-box models). In spite of wide application of simplified linear modeling and simulation of GTs, the behavior of these machines is usually nonlinear, and precise nonlinear models are unavoidable [3].

The data-driven models do not need detailed knowledge about the system to model. Instead, they use available empiric information and are determined by optimization methods or, in the case of artificial neural networks, through machine learning. Because of their simplicity, such models are widely used in GT design. A detailed description of different gas turbine data-driven models can be found, for example, in book [3].

Physics-based modeling relies on physical laws of the functioning of turbo-machines and therefore allows realistic simulation of their behavior. These models are more complex and less used. However, they contain the information difficult to draw from empiric data and are frequently employed as a basis to create simpler data-driven models. Thus, physics-based modeling may be considered as a main gas turbine mathematical modeling type.

1.2. Thermodynamic model

The above reasoning explains why a component-based nonlinear gas turbine model is considered as principle for the design of the engine itself and for developing its control and monitoring systems. It is a highly complex thermodynamic model based on the aerothermal calculations of a gas path and the description of engine’s components (compressor, combustion chamber, turbine, etc.) by nonlinear performance maps. Foundations of the thermodynamic models can be found in [2, 4].

The thermodynamic model comprises two interrelated parts, namely nonlinear static model and nonlinear dynamic model (NDM). The nonlinear static model allows investigating steady-state performances of the engine before its final creation. This model may include more than 100 algebraic and transcendent mathematical relations and, in general, presents a system of nonlinear equations (see [2]). These equations reflect the mass and energy balance between engine components during stationary operation, and the number of such equations typically varies from 5 to 15.

Once the static model has been created, the detailed nonlinear dynamic model can be developed with fewer efforts because it is similar to the static model. The principal difference is that the mentioned algebraic equations of mass and energy balance at steady states are now written in the form of differential equations at transients. The number of such equations corresponds to the number of mass and energy accumulators simulated. Since NDM is a complex and relatively slow procedure, many simplified models are constructed on its basis to be used for the aims of engine control and diagnosis. Nevertheless, along with increasing processor and developing the methods of execution time minimization, direct use of the NDM in real time is becoming possible [2].

The development and use of the thermodynamic models have started in the 1970s, in many respects, by the studies of Saravanamuttoo et al. (for example, [5]). Since that time, many improvements related to higher accuracy and more detailed engine’s component description were introduced in this model; some of them are mentioned below.

Stamatis et al. proposed in [6] the scheme of adaptive simulation by a nonlinear system identification technique that later was used for multipoint gas turbine diagnosis [7]. Ellipsoid functions were introduced in [8] for more accurate description of the components’ maps and better identification of a whole engine at steady states and transients. The authors of paper [9] developed a stage-based compressor model to be used in the thermodynamic model instead of a compressor performance map. This modified thermodynamic model allows the localization of the faulty stages of a multistage compressor and identification of the three compressor degradation mechanisms: fouling, tip clearance increase, and erosion of aerofoils. Thus, gas turbine diagnostics become more profound. Since the early 1990s, Joachim Kurzke has developed the universal program GasTurb for nonlinear physics-based gas turbine simulation [10, 11]. This commercial software allows simulating different types of engines and helps to solve various design and analysis problems. The program GasTurb has special tools to analyze, correct, and enhance the component maps contributing in this way to the accuracy of final engine simulation. Another way to improve the simulation accuracy is proposed by Volponi et al. [12]. As an engine measurement system has individual systematic measurement errors, the authors propose to compensate them by an additional data-driven model on the basis of artificial neural networks. The introduced hybrid model is constructed from a traditional thermodynamic model and this data-driven model. It is shown that the hybrid model can more accurately simulate the performance of a particular engine than the thermodynamic model itself.

The above improvements are related to a static part of the thermodynamic model or both static and dynamic parts. However, the description of engine transients has specific problems to solve, and their solution can additionally improve the dynamic part, namely, detailed nonlinear dynamic model. For gas turbine control and monitoring systems as a whole and, more importantly, for the systems of aircraft gas turbine engines, accurate and fast NDMs are in high demand [2, 13]. These detailed nonlinear models will be useful for the implementation of model predictive control and more effective diagnosis at transients where simplified Kalman filter-based techniques have often been used to date [14]. Since here, this chapter will deal only with such models.

Modern NDMs generally take into considerations three “accumulators”:

• mass and energy accumulation in pneumatic gas path volumes,

• mechanical energy accumulation in the rotors,

• heat accumulation in the stator and rotor heated parts (disks, blades, vanes, case elements).

The volume dynamics is very fast and it is important for controller design. The rotor dynamics lasts for aircraft engines about 10–15 s and has the largest influence on engine performance. The heat exchange dynamics may last many minutes but its direct influence on gas path variables is small because the heat interchange between gas flow and engine-heated parts (HPs) is by far smaller than total energy of the gas. These reasons explain why the models that simulate the rotor dynamics only are still used in diagnostics.

1.3. Tip clearance dynamic effect

There is also an important indirect dynamic effect of the warming-up of engine parts, and this effect has not been taken into consideration in the NDMs yet. The point is that during engine transients, the dynamics of the warming-up is different for a rotor and a stator. The rotor parts, especially massive disks, change slowly their temperature while relatively thin stator parts are warmed up faster. Consequently, the radial displacements of rotor blade tips delay from those of the corresponding casing surfaces, and tip clearances dynamically increase during engine acceleration and decrease during the deceleration.

As the result of increased tip clearance losses, the efficiency of compressor and turbine components lowers and overall engine performance significantly degrades. Sobey and Suggs 1963 demonstrate in book [15] that the 1% turbine tip clearance increase results in the 1% reduction of turbine efficiency and the 1.5–2% increase of engine-specific fuel consumption. The impact of a compressor clearance is even greater: the 1% clearance increase causes the 1.5–3% specific fuel consumption growth. As shown in [16] for the acceleration from the idle to the take-off regime, the increase rate of aircraft engine thrust can reduce twice due to the dynamic clearance increase. The maintenance results show that the corresponding thrust loss can reach from 3 to 15% and takes place from 20th to 60th second after the engine regime change. Thus, significant thrust reduction can happen during the aircraft take-off putting the flight at risk.

Since modern aircraft gas turbine engines need effective control and monitoring systems, accurate detailed nonlinear dynamic models are in increasing demand. In this way, the modeling of the above-described dynamic clearance effect must be implemented in NDMs. So far, such dynamic models use fixed component performance maps obtained at steady states for warmed-up components. The difficulty to introduce the effect of dynamic clearances consists in the fact that they depend on stress-strain state of the stator and rotor parts, and the stresses and strains have irregular distribution that varies in time. Thus, it becomes clear that accurate modeling of the dynamic clearance effect needs the application of finite element methods to the heated parts. The problem is that such calculations are very time-consuming and cannot be directly implemented in NDMs.

As mentioned above, the clearance effect can be modeled only by the analysis of stress-strain state of both stator and rotor HPs with a known initial temperature distribution and a heating rate. The clearance model that meets such conditions was considered in [17]. Unfortunately, this model does not take into account a real shape of HPs and therefore cannot ensure a high accuracy of dynamic engine simulation. The model presented by Archipov et al. in [18] already accounts for the shape but takes other strong limitation that the material properties are independent on material. The authors also make a disputable statement that aerodynamic and pressure gas forces have significant influence on the clearances only for high-power low-speed turbines of industrial power plants. Kurzke proposed NDM with the dynamic clearance model integrated [19]. However, this model does not take into account a radial disk extension and the aerodynamic action of gases on the blade and the casing. Paper [20] compares three variations of the dynamic clearance model. It was found that impulse response model is the most accurate but also time- and memory-consuming. Thus, this model cannot be directly integrated into NDM.

Our previous papers [21, 22] generally follow the ideas of the impulse response model. First, solid models of hot parts of a turbofan had been created. Second, multiple calculations were conducted using the finite element method to understand how the displacements of the hot part surfaces depend on the external temperatures and the loads applied. Third, on the basis of these numerical experiments, a simplified dynamic clearance model (SDCM) was formed. Forth, SDCM was integrated into a nonlinear dynamic model of the turbofan resulting in an enhanced nonlinear dynamic model (ENDM). The rest of this chapter is devoted to the results of the mentioned finite element method calculations, description of the SDCM and its integration into NDM, and the results of simulation by the enhanced model. A high-pressure turbine and its disk are mostly used to exemplify the proposed methodology.

2. Enhanced nonlinear dynamic model

The structure and operation of the above-mentioned ENDM are illustrated by Figure 1. Modules 1–4 constitute SDCM of an engine component, compressor or turbine. Module 5 presents an initial engine model, which, in conjunction with SDCM, presents the enhanced nonlinear model. Only one component is presented on this scheme for simplicity. The software of ENDM includes SDCMs for all the components where the dynamic clearance effect is significant. The enhanced model has been developed for a low-bypass two-spool turbofan engine of a maneuverable aircraft. All simplified relations for the clearance model were obtained through exhaustive calculations in ANSYS with the solid models of engine stator and rotor hot parts, namely disk, blade, and casing. The below description of the clearance model is given for a high-pressure turbine (HPT) of this engine as a component example.

Figure 1.

Structure of the enhanced nonlinear dynamic model.

The initial nonlinear dynamic model includes the dynamics of two engine rotors. The corresponding differential equations are solved through their integration by an iterative procedure. In each step, the corrections to state variables (rotation speeds) are obtained and all engine variables are renewed. The enhanced engine model conserves this iterative character. At an actual step, the variables computed by NDM are used in SDCM to calculate a new tip clearance and corrections to component performances. The modified performances are employed at the next step. In this way, as with real engine dynamics, ENDM uses component performances that are dynamically changed.

Blocks 1.1–1.4 for the stator and Blocks 2.1–2.4 for the rotor illustrate what happens with heated parts (disk, blade, and casing) and tip clearances when engine operating mode dynamically changes. After the mode change, the gas path variables calculated by NMM (Block 5) begin to vary producing the change of heat exchange boundary conditions outside of the heated parts (Blocks 1.1 and 2.1). Because of heat accumulation or loss, the temperature state (distribution of metal temperature) of HPs begins to change (Blocks 1.2 and 2.2). Elevated temperatures of HP cause its thermal expansion. As the HP temperature state has a delay relative to the engine mode change, the corresponding displacement (Blocks 1.3 and 2.3) varies with a delay as well. The knowledge of the temperature state also allows us to correctly consider the action of forces on the HP radial displacements. As shown in [21], for a disk and blade, significant displacements are caused by a centrifugal force, while a pressure force is the most influencing for a casing. These forces are considered in Blocks 1.4 and 2.4 accordingly. When an engine operating point is changing, the force applied to HP changes as well, and the force-induced displacement reacts immediately. However, the engine mode variation also means the change of the heat exchange boundary conditions resulting in other temperature state, other metal elasticity, and an additional change of the displacement. Thus, the force-induced displacement has a static component that immediately reacts on the engine mode and a dynamic component that reacts with a delay. In this way, the total radial displacements of the surfaces that form a tip clearance have a complex dynamic behavior.

The calculation of a dynamic clearance change (see Module 3), which is necessary to evaluate the change in a component performance, is based on simple relations. An actual dynamic tip clearance $\delta$ that is a function of transient time can be expressed through a clearance ${\delta }_0$ of a cold turbine and total displacements ${\mathrm{u}}_{\mathrm{C}},{\mathrm{u}}_{\mathrm{D}},{\delta }_{\mathrm{B}}$ of the casing, disk, and blade accordingly, resulting in

$$\delta ={\delta }_0+{\mathrm{u}}_{\mathrm{C}}-\left({\mathrm{u}}_{\mathrm{D}}+{\mathrm{u}}_{\mathrm{B}}\right)$$ E1

Since the HPT performance map used in NDM corresponds to the turbine parts completely warmed up at steady states, the necessary dynamic correction of the performance will depend on a difference $\mathrm{\Delta }\delta$ between a dynamic clearance $\delta$ and a static clearance ${\delta }^{\mathrm{st}}$ . Let us express the static clearance in the form of Eq. (1) but using static displacements of HPs. This yields:

$${\delta }^{\mathrm{st}}={\delta }_0+{\mathrm{u}}_{\mathrm{C}}^{\mathrm{st}}-\left({\mathrm{u}}_{\mathrm{D}}^{\mathrm{st}}+{\mathrm{u}}_{\mathrm{B}}^{\mathrm{st}}\right)$$ E2

In this way, the dynamic clearance change is written by:

$$\mathrm{\Delta }\delta =\left({\mathrm{u}}_{\mathrm{C}}-{\mathrm{u}}_{\mathrm{C}}^{\mathrm{st}}\right)-\left({\mathrm{u}}_{\mathrm{D}}-{\mathrm{u}}_{\mathrm{D}}^{\mathrm{st}}\right)-\left({\mathrm{u}}_{\mathrm{B}}-{\mathrm{u}}_{\mathrm{B}}^{\mathrm{st}}\right)$$ E3

Using the clearance change $\mathrm{\Delta }\delta$ as an input parameter, Module 4 corrects the component efficiency $\eta$ because it is known that just this performance is affected by an increased clearance. Paper [23] shows that a turbine efficiency loss $\mathrm{\Delta }\eta$ is linearly dependent on a relative clearance change:

$$\mathrm{\Delta }\overline{\mathrm{\delta }}=\raisebox{1ex}{$\mathrm{\Delta }\delta $}\!\left/ \!\raisebox{-1ex}{$\mathrm{L}$}\right.$$ E4

where L is a blade length. Figure 2 illustrates this relation.

Figure 2.

Turbine efficiency losses vs. an increasing tip clearance.

Therefore, Module 4 computes the corrected component efficiency at each point of a transient process according to an expression:

$$\eta ={\eta }^{\mathrm{st}}-\mathrm{k}\mathrm{\Delta }\overline{\mathrm{\delta }}$$ E5

where ${\eta }^{\mathrm{st}}$ is the efficiency at an equivalent steady state and k is the coefficient depending on the construction of a modeled engine and its component.

Since Module 5 employs the corrected efficiencies of all the engine components as input parameters, the engine variables simulated by NDM take into consideration the effect of dynamically varying radial clearances.

As can be seen from the above description, the equations of Modules 3 and 4 are simple and do not need additional explanations. A nonlinear dynamic model of Module 5 is not simple, but it does not need additional description because this type of gas turbine models is well described in literature. However, Modules 1 and 2 and their blocks were presented in this section in a schematic form necessary for a general understanding of the enhanced NDM. As described before, in these modules, the displacements of the heated parts are calculated. To know how these displacements depend on external factors, the displacements were simulated in ANSYS by creating the solid model of each HP and by applying the finite element method to determine the HP stress-stain state. By multiple numerical experiments of this type, the simplified relations between the displacement and external dynamic factors were formed and included in Modules 1 and 2. This ANSYS-based simulation of the HP displacements is described in the next section.

3. Finite element-based displacement simulation

3.1. Thermal boundary conditions and mechanical loads

As shown in Figure 1 and described in Section 2, the displacements of the heated parts depend on thermal boundary conditions and mechanical loads applied. Shown in Figure 3, the design scheme of the disk, which is the most complex HP, illustrates these conditions and loads.

The boundary conditions are set by an external air and gas temperature T and heat transfer coefficient α at the disk surface. Since these conditions considerably vary at different parts of the surface, it is broken down into 24 sections with constant temperatures T_i and coefficients α_i. At one engine steady state called a reference mode, the values ${\mathrm{T}}_{\mathrm{i}}^0$ and ${\mathrm{\alpha }}_{\mathrm{i}}^0$ of these parameters are known on the basis of the experimental information. In a peripheral disk part, in addition to hot gases, heat is transmitted from the blades. This additional heat transfer is taken into account by elevated values T₁ and α₁ in Section 1 of the disk surface.

As to the mechanical loads, the centrifugal force acting on the disk is a body force that is applied to each elemental volume of the disk. The centrifugal force from the rotating blades is given as a surface force by a uniform distribution ${\mathrm{\sigma }}_{\mathrm{B}}$ in Section 1. The design schemes of the blade and the casing are similar.

Figure 3.

Design scheme of disk thermal boundary conditions and mechanical loads.

3.2. Stress-strain state and the displacements of heated parts

In the finite element-based simulation, the heated parts are presented by their solid models illustrated by Figure 4. In its solid model, each HP is divided on elemental 3D simplex volumes. Each volume is presented in finite element calculations by four nodes.

Figure 4.

Solid models of the heated parts of HPT (a—disk, b—blade, c—casing).

After the determination of the temperature state of a heated part, the nonuniform distribution of its temperature t is known. In addition to the action of this temperature, the heated part undergoes the action of a surface force $\mathrm{p}\left(\mathrm{x},\mathrm{y},\mathrm{z}\right)$ and a body force $\mathrm{F}\left(\mathrm{x},\mathrm{y},\mathrm{z}\right)$ . The known temperature irregularity and the forces induce in each node a displacement $\mathbf{u}={\left[{\mathrm{u}}_{\mathrm{х}}{\mathrm{u}}_{\mathrm{y}}{\mathrm{u}}_{\mathrm{z}}\right]}^{\mathrm{т}}$ , a strain $\mathbf{\epsilon }={\left[{\mathrm{\epsilon }}_{\mathrm{х}}{\mathrm{\epsilon }}_{\mathrm{y}}{\mathrm{\epsilon }}_{\mathrm{z}}{\mathrm{\gamma }}_{\mathrm{хy}}{\mathrm{\gamma }}_{\mathrm{yz}}{\mathrm{\gamma }}_{\mathrm{zх}}\right]}^{\mathrm{т}}$ , and a stress $\mathbf{\sigma }={\left[{\mathrm{\sigma }}_{\mathrm{х}}{\mathrm{\sigma }}_{\mathrm{y}}{\mathrm{\sigma }}_{\mathrm{z}}{\mathrm{\sigma }}_{\mathrm{хy}}{\mathrm{\sigma }}_{\mathrm{yz}}{\mathrm{\sigma }}_{\mathrm{zх}}\right]}^{\mathrm{т}}$ that are described by the following linear equations of the elasticity theory (see [24]):

$$\mathbf{\epsilon }=\mathbf{Ru}\mathrm{;}$$ E6

$$\mathbf{\sigma }=\mathbf{D}\left(\mathbf{\epsilon }-\mathbf{\alpha }\mathrm{t}\right)\mathrm{;}$$ E7

$${\mathbf{R}}^{\mathrm{T}}\mathbf{\sigma }+\mathbf{F}=0$$ E8

and by the equation of boundary condition:

$$\mathbf{p}-\mathbf{C\sigma }=0$$ E9

In these equations, R is a differential matrix operator, D presents a stiffness matrix depending on material elasticity and the Poisson’s ratio, α denotes a linear expansion coefficient vector, and C stands for a rotation matrix. On the basis of Eqs. (6)–(9) of one elemental volume, a huge system of linear equations of a whole heated part is formed. The number of unknown variables in this system can be partly reduced because one volume node pertains to some adjacent elemental volumes. The system is solved by the least squares method. As a result, the displacements of the external surface of the heated element are determined separately for the action of thermal expansion and the force. In addition to the reference engine mode, the displacements of all the heated parts were determined at the idle regime.

As mentioned in Section 3.1, the thermal boundary conditions are known at the reference mode and the variables necessary to determine mechanical loads at this mode are simply calculated by NDM. Using these data, the displacements ${\mathrm{u}}_{\mathrm{t}}^0$ and ${\mathrm{u}}_{\mathrm{F}}^0$ induced at this mode by temperature and force were firstly computed in ANSYS for the disk and the other heated parts. To know how these displacements vary during engine operation, let us firstly analyze how the thermal boundary conditions depend on an engine operating mode. It will be shown that the boundary conditions can be determined through actual and reference gas path variables known from NDM.

4. Varying boundary conditions

As mentioned above, the values of the boundary parameters ${\mathrm{T}}_{\mathrm{i}}$ and ${\mathrm{\alpha }}_{\mathrm{i}}$ in the sections of the HP surface (see Figure 3 for the case of the disk) are known only for the reference mode. To have the possibility to make the finite element calculation in ANSYS at any mode, we need to know how these parameters vary along with an engine operating point.

4.1. Boundary temperatures

Oleynik has shown in his thesis [25] that the distribution of boundary temperatures around HP at a current operating mode is similar to the distribution at a reference mode. The calculations made with NDM of the engine under analysis also confirm that gas path temperatures proportionally change from one operating point to another [22]. In this way, we can state that a temperature similarity coefficient ${\mathrm{k}}_{\mathrm{T}}=\raisebox{1ex}{$\left({\mathrm{T}}_{\mathrm{i}}-{\mathrm{T}}_{\mathrm{HPC}}\right)$}\!\left/ \!\raisebox{-1ex}{$\left({\mathrm{T}}_{\mathrm{i}}^0-{\mathrm{T}}_{\mathrm{HPC}}^0\right)$}\right.$ is approximately constant and a current temperature at any section “i” of the HP surface can be expressed through this coefficient by:

$${\mathrm{T}}_{\mathrm{i}}={\mathrm{k}}_{\mathrm{T}}\cdot \left({\mathrm{T}}_{\mathrm{i}}^0-{\mathrm{T}}_{\mathrm{HPC}}^0\right)+{\mathrm{T}}_{\mathrm{HPC}}$$ E10

The similarity coefficient is determined using the gas path temperatures computed by NDM at the reference and actual engine modes.

4.2. Heat transfer coefficients

Paper [22] shows that the heat transfer coefficients ${\mathrm{\alpha }}_{\mathrm{i}}$ change proportionally when an operating mode varies. Using known relations between different criteria of gas flow, this chapter derives the following equation for a similarity coefficient:

$${\mathrm{k}}_{\mathrm{\alpha }}=\frac{\mathrm{\alpha }}{{\mathrm{\alpha }}^0}={\left(\frac{\mathrm{n}}{{\mathrm{n}}^0}\frac{{\mathrm{P}}_{\mathrm{HPC}}}{{\mathrm{P}}_{\mathrm{HPC}}^0}\right)}^{0.8}{\left(\frac{{\mathrm{T}}_{\mathrm{HPC}}}{{\mathrm{T}}_{\mathrm{HPC}}^0}\right)}^{-0.567}$$ E11

As the necessary actual and reference values of gas path variables are known from NDM, the similarity coefficient is simply calculated and the coefficients ${\mathrm{\alpha }}_{\mathrm{i}}$ at the HP surface sections are determined by:

$${\mathrm{\alpha }}_{\mathrm{i}}={\mathrm{k}}_{\mathrm{\alpha }}{\mathrm{\alpha }}_{\mathrm{i}}^0$$ E12

In this way, the distribution of the boundary variables T and ${\mathrm{\alpha }}_{\mathrm{i}}={\mathrm{k}}_{\mathrm{\alpha }}{\mathrm{\alpha }}_{\mathrm{i}}^0$ can be simply determined through the NDM gas path variables, namely, HPT rotation speed n, high pressure compressor (HPC) discharge temperature ${\mathrm{T}}_{\mathrm{HPC}}$ , and HPC discharge pressures ${\mathrm{P}}_{\mathrm{HPC}}$ .

The next challenging problem was to create the relations for calculating the HP displacements, both temperature induced and force induced, at any engine dynamic operating point. Let us begin from the displacements due to thermal expansion of the heated parts.

5. Varying thermal expansion displacements

The dynamics of the displacements caused by temperature loading is described below using a displacement transient performance. To determine this performance, the influence of a step change of boundary temperatures from 293 K (cold disk) to the distribution at the reference mode was simulated in ANSYS. It was found that, in addition to time τ, the displacement also depends on the heat transfer similarity coefficient ${\mathrm{k}}_{\mathrm{\alpha }}$ , and the displacement performance was presented as a relative function $\overline{\mathrm{u}}\left(\mathrm{\tau },{\mathrm{k}}_{\mathrm{\alpha }}\right)=\frac{\mathrm{u}\left(\mathrm{\tau },{\mathrm{k}}_{\mathrm{\alpha }}\right)-{\mathrm{u}}^0}{{\mathrm{u}}^{\mathrm{st}}-{\mathrm{u}}^0}$ illustrated by Figure 5.

Figure 5.

Transient performance of a disk displacement.

Figure 6 shows the transient performances of the blade and casing absolute displacements obtained in ANSYS by the same mode. For these heated parts, the influence of the coefficient ${\mathrm{k}}_{\mathrm{\alpha }}$ is insignificant.

Figure 6.

Transient performance of blade and casing displacements (a—blade, b—casing).

Using the disk as an example, let us now show how to consider its displacement performance in a total process of the ENDM computing. Paper [22] demonstrates that, for each value ${\mathrm{k}}_{\mathrm{\alpha }}$ , the corresponding curve in Figure 5 is accurately described by a weighted sum of two exponents and therefore can be presented by:

$$\overline{\mathrm{u}}\left(\mathrm{\tau },{\mathrm{k}}_{\mathrm{\alpha }}\right)=\sum _{\mathrm{j}=1}^2{\mathrm{k}}_{\mathrm{j}}\left({\mathrm{k}}_{\mathrm{\alpha }}\right)\left(1-{\mathrm{e}}^{\frac{-\mathrm{\tau }}{{\mathrm{T}}_{\mathrm{j}}\left({\mathrm{k}}_{\mathrm{\alpha }}\right)}}\right),\mathrm{j}=1,2$$ E13

where ${\mathrm{k}}_{\mathrm{j}}$ is a weighting coefficient and ${\mathrm{T}}_{\mathrm{j}}$ is a time constant. For each value of ${\mathrm{k}}_{\mathrm{\alpha }}$ , four parameters ${\mathrm{k}}_1$ , ${\mathrm{k}}_2$ , ${\mathrm{T}}_1$ , and ${\mathrm{T}}_2$ were determined. Figure 7 illustrates their dependency on the coefficient ${\mathrm{k}}_{\mathrm{\alpha }}$ .

Figure 7.

Weighting coefficients and time constants vs. heat transfer similarity coefficient [21] (a—weighting coefficients, b—time constants).

The two mentioned exponents present analytical solutions of linear differential equations that for absolute displacements take a form:

$${\mathrm{T}}_{\mathrm{j}}\left({\mathrm{k}}_{\mathrm{\alpha }}\right)\frac{{\mathrm{du}}_{\mathrm{j}}}{\mathrm{d\tau }}+{\mathrm{u}}_{\mathrm{j}}={\mathrm{u}}^0\left({\mathrm{k}}_{\mathrm{\alpha }}\right),\mathrm{j}=1,2$$ E14

The following two equations:

$${\mathrm{u}}_{\mathrm{j}}\left(\mathrm{\tau },{\mathrm{k}}_{\mathrm{\alpha }}\right)=\frac{{\mathrm{T}}_{\mathrm{j}}\left({\mathrm{k}}_{\mathrm{\alpha }}\right)}{{\mathrm{T}}_{\mathrm{j}}\left({\mathrm{k}}_{\mathrm{\alpha }}\right)+\mathrm{\Delta }\mathrm{\tau }}{\mathrm{u}}_{\mathrm{j}}\left(\mathrm{\tau }-\mathrm{\Delta }\mathrm{\tau }\right)+\frac{\mathrm{\Delta }\mathrm{\tau }}{{\mathrm{T}}_{\mathrm{j}}\left({\mathrm{k}}_{\mathrm{\alpha }}\right)+\mathrm{\Delta }\mathrm{\tau }}{\mathrm{u}}^0\left({\mathrm{k}}_{\mathrm{\alpha }}\right),\mathrm{j}=1,2$$ E15

are numerical solutions of these equations. Their weighted sum:

$$\mathrm{u}\left(\mathrm{\tau },{\mathrm{k}}_{\mathrm{\alpha }}\right)=\sum _{\mathrm{j}=1}^2{\mathrm{k}}_{\mathrm{j}}\left({\mathrm{k}}_{\mathrm{\alpha }}\right){\mathrm{u}}_{\mathrm{j}}\left(\mathrm{\tau },{\mathrm{k}}_{\mathrm{\alpha }}\right),\mathrm{j}=1,2$$ E16

is a final expression to numerically compute the dynamic displacement caused by thermal disk expansion. Eqs. (15) and (16) present final steps in Block 1.3 of the enhanced nonlinear dynamic model (see Figure 1). Through the coefficient ${\mathrm{k}}_{\mathrm{\alpha }}$ , the displacement calculation is adapted to an actual dynamic engine operating point. The blade and casing displacement (Blocks 1.3 and 2.3) are computed similarly.

6. Varying force-induced displacements

The displacements induced in HPs by mechanical loads can be considered elastic and proportional to the load. For the disk and the blade, the main load is a centrifugal force and the displacements will be proportional to the rotation speed squared n². As the casing is mainly loaded by a pressure force, the displacement will linearly depend on the HPC pressure ${\mathrm{P}}_{\mathrm{HPC}}$ . The action of these forces has no delay and the displacement will change along with the load change.

However, since the elasticity coefficient depends on the HP temperature, the HP displacement should be simulated regarding this dependency. The temperature distribution within HP is nonuniform and dynamically changes during transient engine operation. For this reason, it will be difficult to directly simulate the elasticity change. To solve this problem, paper [22] proposes the concept of an equivalent temperature.

6.1. Equivalent temperature

The equivalent temperature t_e is defined as a temperature of a uniformly heated engine part, which has load-induced displacements equal to the displacements of HP with an actual temperature state and the same mechanical loading. Using the temperature t_e, the displacement at an actual dynamic point is written for the disk and blade by:

$${\mathrm{u}}_{\mathrm{F}}={\mathrm{u}}_{\mathrm{F}}^{\circ }\left({\mathrm{t}}_{\mathrm{e}}\right)\cdot {\left(\frac{\mathrm{n}}{{\mathrm{n}}^0}\right)}^2$$ E17

and for the casing by:

$${\mathrm{u}}_{\mathrm{F}}={\mathrm{u}}_{\mathrm{F}}^{\circ }\left({\mathrm{t}}_{\mathrm{e}}\right)\cdot \left(\frac{{\mathrm{P}}_{\mathrm{HPC}}}{{\mathrm{P}}_{\mathrm{HPC}}^0}\right)$$ E18

The displacement ${\mathrm{u}}_{\mathrm{F}}^{\circ }$ corresponds to a hypothetical situation when HP is under the constant mechanical load of the reference mode, but the HP heating conditions are varying and correspond to the actual engine operating point. A function ${\mathrm{u}}_{\mathrm{F}}^{\circ }\left({\mathrm{t}}_{\mathrm{e}}\right)$ was determined by simulating such hypothetical loading in ANSYS. Figure 8 illustrates the results of the disk displacement simulations. These results are approximated by:

$${\mathrm{u}}_{\mathrm{F}}^{\circ }\left({\mathrm{t}}_{\mathrm{e}}\right)=0.38557+8.55627\times {10}^{‐5}\cdot {\mathrm{t}}_{\mathrm{e}}+1.83458\times {10}^{‐8}\cdot {\mathrm{t}}_{\mathrm{e}}^2$$ E19

Figure 8.

Disk displacement at the reference mode vs. equivalent disk temperature.

6.2. Characteristic temperature

As follows from Eqs. (10) and (12), thermal loading on each heated part (disk, blade, and casing) depends on the temperature ${\mathrm{T}}_{\mathrm{HPC}}$ (temperature of HPC air) and the similarity coefficients ${\mathrm{k}}_{\mathrm{T}}$ and ${\mathrm{k}}_{\mathrm{\alpha }}$ . As described in Section 4, the radial displacement ${\mathrm{u}}_{\mathrm{F}}$ caused by the force depends on the temperature state of HP and, therefore, is related to the thermal loading. Thus, this relation can be written by a function ${\mathrm{u}}_{\mathrm{F}}=\mathrm{f}\left({\mathrm{T}}_{\mathrm{HPC}},{\mathrm{k}}_{\mathrm{T}},{\mathrm{k}}_{\mathrm{\alpha }}\right)$ . The three interrelated arguments make this function complex for realization. Paper [22] proposes the concept of a characteristic temperature to be used as the unique function argument. The characteristic temperature $\tilde{\mathrm{T}}$ is defined as a weighted mean of a boundary temperature T.

$$\tilde{\mathrm{T}}=\frac{\underset{\mathrm{A}}{\int }\mathrm{T}\left(\mathrm{A}\right)\mathrm{\alpha }\left(\mathrm{A}\right)\mathrm{dA}}{\underset{\mathrm{A}}{\int }\mathrm{\alpha }\left(\mathrm{A}\right)\mathrm{dA}}$$ E20

where A is the surface of a heated part.

The characteristic temperature $\tilde{\mathrm{T}}$ has an important property that temperatures t of a heated part tend to a value $\tilde{\mathrm{T}}$ when heat transfer approaches zero, i.e.:

$$\underset{{\mathrm{k}}_{\mathrm{\alpha }}\to 0}{\text{lim}}\mathrm{t}=\tilde{\mathrm{T}}$$ E21

This property allows us to determine the characteristic temperature through ANSYS simulation of the heated part with an extremely low similarity coefficient ${\mathrm{k}}_{\mathrm{\alpha }}$ . It is proven that such simulation yields low errors relatively a direct calculation of $\tilde{\mathrm{T}}$ according to Eq. (20). For example, given ${\mathrm{k}}_{\mathrm{\alpha }}\approx {10}^{-3}$ , the error was 0.01 K.

The characteristic temperature was firstly computed at the reference mode and, with the known value ${{\displaystyle \tilde{\mathrm{T}}}}^0$ , a temperature coefficient:

$$\tilde{\mathrm{\Theta }}=\frac{{{\displaystyle \tilde{\mathrm{T}}}}^0-{\mathrm{T}}_{\mathrm{HPC}}^0}{{\mathrm{T}}_{\mathrm{g}}^0-{\mathrm{T}}_{\mathrm{HPC}}^0}$$ E22

was formed, where ${\mathrm{T}}_{\mathrm{g}}$ denotes a HPT input temperature. Then, it was found that this coefficient does not depend on an operating mode and can be used to determine the characteristic temperature at any mode by a simple relation:

$$\tilde{\mathrm{T}}=\tilde{\mathrm{\Theta }}\cdot \left({\mathrm{T}}_{\mathrm{g}}-{\mathrm{T}}_{\mathrm{HPC}}\right)+{\mathrm{T}}_{\mathrm{HPC}}$$ E23

6.3. Static force-induced displacement

To determine the relation between the temperatures $\tilde{\mathrm{T}}$ and t_e, series of simulations in ANSYS have been conducted. For the disk under reference mechanical load, the thermal load parameters ${\mathrm{T}}_{\mathrm{HPC}}$ , ${\mathrm{k}}_{\mathrm{t}}$ , and ${\mathrm{k}}_{\mathrm{\alpha }}$ are varied and the displacement ${\mathrm{u}}_{\mathrm{F}}^0=\mathrm{f}\left({\mathrm{T}}_{\mathrm{HPC}},{\mathrm{k}}_{\mathrm{t}},{\mathrm{k}}_{\mathrm{\alpha }}\right)$ was determined for each combination of ${\mathrm{T}}_{\mathrm{HPC}}$ , ${\mathrm{k}}_{\mathrm{t}}$ , and ${\mathrm{k}}_{\mathrm{\alpha }}$ . The equivalent temperature t_e corresponding to each displacement was found from Eq. (19). The characteristic temperature $\tilde{\mathrm{T}}$ was calculated according to Eq. (23) using a known value ${\mathrm{T}}_{\mathrm{HPC}}$ and a gas temperature ${\mathrm{T}}_{\mathrm{g}}$ computed by NDM. By doing so, multiple pairs of t_e and $\tilde{\mathrm{T}}$ values were found. With these data, the relation is between t_e and $\tilde{\mathrm{T}}$ is described by:

$${\mathrm{t}}_{\mathrm{e}}=-16.631+1.0518\cdot \tilde{\mathrm{T}}-3.9362\times {10}^{-5}\cdot {{\displaystyle \tilde{\mathrm{T}}}}^2$$ E24

Thus, through a consecutive application of Eqs. (23), (24), (19), and (17), we can calculate a force-induced radial displacements of the disk as a function of the gas path variable ${\mathrm{T}}_{\mathrm{HPC}}$ and ${\mathrm{T}}_{\mathrm{g}}$ computed by NDM. The displacements of this enhanced algorithm as well as the original algorithm that consider constant disk elasticity were estimated by the comparison with the results of ANSYS-based simulations. Figure 9 presents the errors of both algorithms for different characteristic temperatures and consequently for different engine operating points. We can see that the original algorithm has significant errors (up to 9%), whereas for the enhanced algorithm, the errors are negligible (within 0.1%). So, the accuracy of the displacement simulation was drastically enhanced despite the simplicity of the proposed algorithm.

Figure 9.

Errors of two force-induced displacement algorithms (scored line: algorithm that considers temperature-dependent elasticity; dashed line: algorithm that uses constant elasticity).

6.4. Dynamic force-induced displacement

The equivalent temperature t_e determined in Eq. (24) as a function of $\tilde{\mathrm{T}}$ corresponds to a completely warmed-up heated part and its final static displacement. Let us call this temperature a static equivalent temperature ${\mathrm{t}}_{\mathrm{e}}^{\mathrm{st}}\left(\tilde{\mathrm{T}}\right)$ . When the boundary conditions have changed, the force-induced displacements will vary dynamically and the temperature t_e will dynamically approach ${\mathrm{t}}_{\mathrm{e}}^{\mathrm{st}}$ . As the relation between the displacement and the temperature t_e is practically linear (see Figure 8), their dynamic behavior will be similar. For this reason, the dynamics of t_e are described using the same displacement transient performances presented in Figures 5 and 6. For the disk, the algorithm to compute t_e is similar to that described in Section 5 for the thermal expansion displacements. The resulting equations to compute the equivalent temperature:

$${\mathrm{t}}_{\mathrm{ej}}\left(\mathrm{\tau },{\mathrm{k}}_{\mathrm{\alpha }}\right)=\frac{{\mathrm{T}}_{\mathrm{j}}\left({\mathrm{k}}_{\mathrm{\alpha }}\right)}{{\mathrm{T}}_{\mathrm{j}}\left({\mathrm{k}}_{\mathrm{\alpha }}\right)+\mathrm{\Delta }\mathrm{\tau }}{\mathrm{t}}_{\mathrm{ej}}\left(\mathrm{\tau }-\mathrm{\Delta }\mathrm{\tau }\right)+\frac{\mathrm{\Delta }\mathrm{\tau }}{{\mathrm{T}}_{\mathrm{j}}\left({\mathrm{k}}_{\mathrm{\alpha }}\right)+\mathrm{\Delta }\mathrm{\tau }}{\mathrm{t}}_{\mathrm{e}}^{\mathrm{st}}\left({\mathrm{k}}_{\mathrm{\alpha }}\right),\mathrm{j}=1,2$$ E25

and

$${\mathrm{t}}_{\mathrm{e}}\left(\mathrm{\tau },{\mathrm{k}}_{\mathrm{\alpha }}\right)=\sum _{\mathrm{j}=1}^2{\mathrm{k}}_{\mathrm{j}}\left({\mathrm{k}}_{\mathrm{\alpha }}\right){\mathrm{t}}_{\mathrm{ej}}\left(\mathrm{\tau },{\mathrm{k}}_{\mathrm{\alpha }}\right),\mathrm{j}=1,2$$ E26

are also similar to displacement Eqs. (15) and (16) and the same parameters ${\mathrm{T}}_{\mathrm{j}}$ and ${\mathrm{k}}_{\mathrm{j}}$ are employed. Using the dynamic value ${\mathrm{t}}_{\mathrm{e}}\left(\mathrm{\tau },{\mathrm{k}}_{\mathrm{\alpha }}\right)$ from Eq. (26) as an argument, a dynamic displacement ${\mathrm{u}}_{\mathrm{F}}^{\circ }\left({\mathrm{t}}_{\mathrm{e}}\right)$ is determined from Eq. (19) and a total force-induced disk displacement ${\mathrm{u}}_{\mathrm{F}}$ from Eq. (17). The blade and casing force-induced displacements are computed by similar algorithms. All these algorithms correspond to Blocks 1.4 and 2.4 of the engine ENDM presented in Figure 1.

7. Verification of the enhanced nonlinear dynamic model

7.1. Verification of the simplified dynamic clearance model

To verify the simplified dynamic clearance model (see Section 2), the following engine dynamics test case was prepared:

• during the time interval $\mathrm{\tau }=\mathrm{0\dots 120}\mathrm{s}$ , the turbofan engine operates at idle conditions ( ${\mathrm{k}}_{\mathrm{\alpha }}$ =0.2031, ${\mathrm{k}}_{\mathrm{t}}$ =0.4125, ${\mathrm{k}}_{\mathrm{n}}=0.5929$ ) to warm up turbine parts;

• during the time interval $\mathrm{\tau }=\mathrm{120\dots 500}\mathrm{s}$ , engine operates under the reference mode conditions ( ${\mathrm{k}}_{\mathrm{\alpha }}$ = 1.0, ${\mathrm{k}}_{\mathrm{t}}$ = 1.0, ${\mathrm{k}}_{\mathrm{n}}$ = 1.0).

Total HPT disk displacements (mechanical and temperature-induced) were computed for this test case in ANSYS and by the proposed SDCM (see Section 2). As shown in Figure 10, the simulation curves practically coincide. The maximum difference observed at the mode change moment is about 0.05 mm and then it lessens. Thus, the simplified model can be considered accurate enough and can be used within ENDM of the turbofan engine under analysis.

Figure 10.

Total disk displacement simulation (dashed line—ANSYS; solid line—SDCM).

7.2. Accuracy of the simulation of engine dynamic performance

As mentioned in the beginning of Section 2, the enhanced nonlinear dynamic model (ENDM) has been developed for a turbofan engine of a maneuverable aircraft. The main objective was to help with the synthesis and adjustment of the algorithms of an engine automatic control system. The developed ENDM is based on the original nonlinear dynamic model (NDM) and the simplified dynamic clearance models (SDCMs) created for a high-pressure turbine (HPT) and a low-pressure turbine.

To verify the accuracy of the ENDM, it was compared with original NDM and with experimental data. A test-case transient was set by a low-pressure rotor speed n_LP (control variable) profile and constant ambient conditions. The profile presents a constant speed value 8100 rpm during the first 175 s, than a linear change to 12,400 rpm during 12 s, and the same constant value up to the transient end.

Figure 11 illustrates the dynamics of the HPT radial clearance simulated by ENDM in comparison with the steady-state clearance simulation (completely warmed-up turbine parts). One can state that ENDM correctly reflects the physics of real warming-up. From the beginning of the engine acceleration, the clearance descends in 15 s because the blade is rapidly warmed up. Next, the clearance grows due to the casing warming up. Finally, the clearance descends once more as the disk begins to warm up.

Figure 11.

Dynamics of the HPT radial clearance (1—steady-state operating modes; 2—ENDM).

Figures 12 and 13 present the results of the comparison of the initial and enhanced dynamic models between each other and with experimental data for the same test-case transient. The plots of a fuel consumption variable in Figure 12 clearly show that the ENDM and experimental curves practically coincide. Both show the same fuel consumption overshoot after the control parameter change, and this overshoot gradually decreases during 150 s for both curves. This elevated fuel consumption is explained by increased turbine clearances due to the delay in disk warming-up. In contrast, the NDM curve does not manifest a visible overshoot, and the transient process is by far shorter. One can make the same conclusion analyzing the plots of a high-pressure rotor speed in Figure 13: the ENDM curve better fits experimental data, in particular, better reflects the effect of increased clearances.

Figure 12.

Fuel consumption dynamics (1—experimental data; 2—ENDM; 3—NDM).

Figure 13.

High-pressure rotor speed dynamics (1—experimental data; 2—ENDM; 3—NDM).

The thrust is the principal parameter of a turbofan. However, under the control law n_LP = const used in the experiments, it is constant as well, and the increased clearances are compensated by the additional fuel consumption observed in Figure 12.

To show the impact of the clearances dynamically changed on the thrust, the simulation of the turbofan under the control law of a constant low pressure turbine temperature was performed. Figure 14 shows the thrust simulated by both models. It can be seen that, during the first 5 s of intensive engine dynamics, both models are equal. Then, the NDM thrust remains constant, whereas the ENDM thrust begins to decrease with a maximal 7% thrust dip at the 12th second. Finally, the thrust gradually increases up to a steady-state value. Such behavior of the thrust simulated by ENDM completely corresponds to the known empirical information about the clearance influence.

Figure 14.

Thrust dynamics (1—ENDM; 2—NDM).

In this way, all the comparison results show that, first, the dynamic clearance influence is significant and cannot be neglected, and, second, the enhanced nonlinear dynamic model accurately simulates this effect and in general provides by far more realistic simulation than the initial dynamic model does.

8. Conclusion

This chapter describes a novel method to enhance a detailed physics-based nonlinear gas turbine model widely used for the aims of aircraft engine control and diagnostics. The method allows us to solve the issue of the impact of varying turbine tip clearances on the dynamic engine performance. This issue is especially important for the engines of maneuverable aircrafts.

Using the proposed method, an enhanced nonlinear dynamic model of a turbofan engine for a maneuverable aircraft has been developed on the basis of an initial nonlinear dynamic model and a simplified dynamic clearance model created with the results of the finite element simulation of turbine parts. The comparison with the initial model and experimental data confirmed a drastic improvement of the accuracy of dynamic gas turbine simulation.

Acknowledgments

This work has been carried out with the support of the National Polytechnic Institute of Mexico (research project 20181152).

Nomenclature

A

surface

k

coefficient

n

rotation speed

P

pressure

t

time, s; temperature of a heated part, K

T

temperature of air or gas, K

$\tilde{\mathrm{T}}$

characteristic temperature, K

U

radial displacement of a heated part, mm

α

heat transfer coefficient

Δδ

clearance change, mm

δ

radial clearance between rotor and stator parts, mm

ε

strain; relative error, %

η

efficiency

σ

stress

$\tilde{\mathrm{\Theta }}$

temperature coefficient

Superscripts °

reference engine operating mode

st

static

Subscripts B

blade

C

casing

D

disk

e

equivalent

F

centrifugal force

g

gas

HPC

high-pressure compressor

LP

low-pressure rotor

P

pressure

R

rotor

S

stator

T

temperature

α

heat transfer coefficient