Nomenclature

$a,b,{c_i},{d_i}$

model parameters

$A,B,{C_i}$

augmented model matrices

${A_e}$

exit area

${C_v}$

nozzle exit velocity coefficient

$EGT$

exhaust gas temperature

$FAR$

fuel to air ratio

$G$

transfer function

$h$

enthalpy

$H$

fuel heating value

$I$

identity matrix

ISA

International Standard Atmosphere

$J$

moment of inertia of the rotor

$K$

percentage of pressure drop

${K_i}$

feedback gain i

$\dot m$

mass flow

$M$

Mach

$N$

spool speed

$\dot N$

spool acceleration

${N_c}$

corrected spool speed

${p_i}$

reference prefilter gain i

${P_{amb}}$

ambient pressure

${P_s}$

Sstatic pressure

${P_{std}}$

standard pressure

${P_{tin}}$

total inlet pressure

${P_{tout}}$

total outlet pressure

PI

proportional integral controller

PLA

power level angle

${r_i}$

setpoint i

SLS

standard sea-level

${t_i}$

output of memory-less function i

${T_{std}}$

standard temperature

${T_{tin}}$

total inlet temperature

${T_{tout}}$

total outlet temperature

$u$

control input

${u_r}$

control input rate

$V$

Lyapunov function

${W_f}$

fuel flow

$x$

state variable

$X$

augmented state vector

${y_i}$

increment of the output i

${z_i}$

input of memory-less function i

Greek symbol

$\delta $

non-dimensional pressure

$\eta $

efficiency

$\theta $

non-dimensional temperature

$\Lambda $

vector containing static nonlinearities

$\pi $

pressure ratio

$\tau $

torque

${\Omega _i}$

diagonal matrix i for sector definition

2.0 Thermodynamic modeling of the engine

In the current research, the AMT NL Titan 40kg thrust engine is taken into account as a small-scale turbojet engine [23]. Figure 1 illustrates the primary structure of this engine, while its specifications are detailed in Table 1. The engine features a single-stage centrifugal compressor, an annular combustion chamber, and a single-stage axial turbine, and it is used as the propulsion system for some aerial vehicles such as the DEMON UAV [Reference Savvaris, Buonanno, Jamil and Tsourdos24].

Table 1.

Specifications of the studied engine [23]

The thermodynamic model of the engine developed in this study is based on a component-level framework. It is derived from the thermodynamic relations for the various engine components, including compressor, combustion chamber, turbine and nozzle. Additionally, it incorporates performance map data for the compressor and turbine [Reference Walsh and Fletcher25, Reference Lichtsinder and Levy26]. The average properties of the fluid at the input and output of each component are taken into account, and the scaled maps for the compressor and turbine are utilised. Engine modeling is conducted within the MATLAB-Simulink environment, and some of the capabilities of the toolbox for modeling and analysis of thermodynamic systems (T-MATS) offered by NASA Glenn Research Center (GRC) are utilised [Reference Zinnecker, Chapman, Lavelle and Litt27]. The relationships governing the engine components are described in the following subsections.

2.1 Compressor

Turbine engines are equipped with a compressor that increases the pressure of the incoming air before it reaches the combustor. The algebraic equations governing the compressor are described as follows [Reference Chung, Prasad, Dhingra and Meisner28]:

(1) \begin{align} \theta & = \frac{{{T_{tin}}}}{{{T_{std}}}},\,\,\,\delta = \frac{{{P_{tin}}}}{{{P_{std}}}},\,\,\,{N_c} = \frac{N}{{\sqrt \theta }}\nonumber\\[4pt]{{\dot m}_{c,map}} & = {f_1}\left( {{N_c},\beta line} \right)\nonumber\\[4pt]\pi & = {f_2}\left( {{N_c},\beta line} \right)\nonumber\\[4pt]\eta & = {f_3}\left( {{N_c},\beta line} \right)\nonumber\\[4pt]{{\dot m}_{out}} & = \frac{{{{\dot m}_{c,map}}}}{{{{\sqrt \theta } \mathord{\left/ {\vphantom {{\sqrt \theta } \delta }} \right. } \delta }}}\nonumber\\[4pt]{h_{tout}} & = \frac{{{h_{tout,ideal}} - {h_{tin}}}}{\eta } + {h_{tin}}\nonumber\\[2pt]{P_{tout}} & = \pi \left( {{P_{tin}}} \right),\,\,\,{T_{tout}} = f\left( {{P_{tout}},{h_{tout}}} \right)\nonumber\\[2pt]\tau & = \frac{{{{\dot m}_{out}}\left( {{h_{tin}} - {h_{tout}}} \right)}}{{N\left( {\frac{{2\pi }}{{60}}} \right)}}\end{align}

where ${f_1}$ , ${f_2}$ and ${f_3}$ are functions or lookup tables that represent the compressor map. Gas properties such as enthalpy and temperature are determined using gas tables.

2.2 Combustion chamber

In the combustion chamber, injected fuel is combined with compressed air and after combustion, high temperature gas is produced. The hot gas drives the turbine and create thrust as it passes through the nozzle. The pressure drop in the combustion chamber is assumed to be constant. The relationships related to this component are expressed as follows:

(2) \begin{align} {P_{tout}} & = \left( {1 - K} \right){P_{tin}}\nonumber\\[3pt]{{\dot m}_{gas,out}} & = {{\dot m}_{air,in}} + {{\dot m}_f}\nonumber\\[4pt]{h_{tout}} & = \frac{{{\eta _b}{{\dot m}_f}H + {{\dot m}_{in}}{h_{tin}}}}{{{{\dot m}_{out}}}}\nonumber\\[4pt]FAR & = \frac{{{{\dot m}_f}}}{{{{\dot m}_{air,in}}}}\end{align}

where K = 0.05 is the percentage of pressure drop in the combustor and H is the heating value of the fuel.

2.3 Turbine

Gas expansion in the turbine generates work of the rotor and drives the compressor. The turbine processes can be represented as follows:

(3) \begin{align} \theta & = \frac{{{T_{tin}}}}{{{T_{std}}}},\,\,\,\delta = \frac{{{P_{tin}}}}{{{P_{std}}}},\,\,\,{N_c} = \frac{N}{{\sqrt \theta }}\nonumber\\[4pt]{{\dot m}_{out}} & = {g_1}\left( {{N_c},\beta line} \right)\nonumber\\[3pt]\pi & = {g_2}\left( {{N_c},\beta line} \right)\nonumber\\[3pt]\eta & = {g_3}\left( {{N_c},\beta line} \right)\nonumber\\[3pt]{h_{tout\_g}} & = {h_{tin}} - ({h_{tin}} - {h_{tout\_g,ideal}})\eta \nonumber\\[3pt]{P_{tout}} & = \pi \left( {{P_{tin}}} \right),\,\,\,{T_{tout}} = f\left( {{P_{tout}},{h_{tout}}} \right)\nonumber\\[3pt]{h_{tout}} & = \frac{1}{{{{\dot m}_{out}}}}\left[ {{h_{tin}}{{\dot m}_{tin}} + \eta ({h_{tin}} - {h_{tout\_g,ideal}}){{\dot m}_{in}}} \right]\nonumber\\[5pt]\tau & = \frac{{{{\dot m}_{out}}({h_{tin}} - {h_{tout\_g}})}}{{N\left( {\frac{{2\pi }}{{60}}} \right)}} \end{align}

similar to the compressor, ${g_1}$ , ${g_2}$ and ${g_3}$ are a set of functions or lookup tables that define the turbine map.

2.4 Nozzle

The engine nozzle has a constant outlet area and is a converging type. In addition to the characteristics of the fluid entering the nozzle, the ambient static pressure at the nozzle exit affects the nozzle flow. Flow characteristics according to choking and non-choking conditions of the flow are specified by comparing the nozzle pressure ratio to the critical pressure.

(4) \begin{align} & \pi = \dfrac{{{P_{tin}}}}{{{P_{amb}}}}\nonumber\\[4pt]& {P_s} = g({h_{tin}},{P_{tout}},M = 1)\nonumber\\[3pt]& {\pi _{critical}} = \dfrac{{{P_{tin}}}}{{{P_s}}}\nonumber\\[3pt]& {\rm{If \, }}\pi \gt {\pi _{critical}}\nonumber\\[3pt]& \left[ {{P_s},u} \right] = {g_1}({h_{tin}},{P_{tout}},M = 1)\nonumber\\[3pt]& {\rm{else}}\nonumber\\[3pt]& \left[ {{P_s},u} \right] = {g_2}({h_{tin}},{P_{tout}},{P_s} = {P_{amb}})\nonumber\\[3pt]& Thrust = {{\dot m}_{out}}u{C_v} + ({P_{out}} - {P_{amb}}){A_e},\,\,\,{C_v} = 1 \end{align}

2.5 Shaft

Using the torque supplied by the turbine and the torque used by the compressor, the shaft dynamics can be written as follows:

(5) \begin{align} \dot N = \dfrac{{\Delta \tau }}{{J\left( {\dfrac{{2\pi }}{{60}}} \right)}} \end{align}

where $\dot N$ is the shaft angular acceleration, $\Delta \tau $ represents the difference between the torque generated by the turbine and the torque demanded by the compressor, and J is the moment of inertia of the rotor. In order to obtain the thermodynamic model by matrix iteration technique, a number of engine parameters are guessed and the same number of errors are defined. For this objective, in this study, inlet air mass flow, beta values for the compressor and turbine are estimated and the errors related to mass conservation in the compressor, turbine and nozzle are considered. The guesses are updated at each iteration through a numerical procedure until the errors are less than a pre-defined value. The numerical process is generally performed using the iterative Newton–Raphson procedure [Reference Walsh and Fletcher25]. Figure 2 shows the iterative process of solving the equations to achieve the transient model.

Figure 2.

Iterative procedure for engine thermodynamic modeling.

In order to show the validity of the modeling, according to Figs. 3 and 4, two important outputs of the engine (i.e., the thrust value and fuel flow) in terms of relative spool speed, have been compared with the engine test data. The maximum mismatch error for thrust is 3%, and for fuel consumption is 6%. Therefore, the obtained thermodynamic model has sufficient validity for performance analysis and controller design.

Figure 3.

Engine thrust versus relative spool speed.

Figure 4.

Engine fuel flow versus relative spool speed.

Because there is no explicit nonlinear dynamic equation for gas turbine engines, the nonlinear thermodynamic model is linearised around an operating point through identification techniques for the purpose of controller design. The linear model of the system with integrated input is expressed as

(6) \begin{align} \begin{array}{l}\dot x = ax + bu\\\dot u = {u_r}\\{y_i} = {c_i}x + {d_i}u\end{array} \end{align}

where $x = \Delta N$ is the state variable (i.e., engine speed increment) in rpm, $u = \Delta {W_f}$ is the control input (i.e., fuel flow increment) in kilograms per second (kg/s), and ${y_i}$ are increments of the outputs.

3.0 Control strategy

The engine controller must provide the desired thrust and at the same time protect the limits of turbine temperature, maximum engine speed, compressor surge and combustion chamber blow-out. To achieve this objective, a selective control strategy [Reference Cai, Wang, Chen, Zheng and Zhang29, Reference Eslami and Banazadeh30] has been taken into account. The structure of the proposed controller is illustrated in Fig. 5. Considering that engine thrust cannot be measured in flight conditions [Reference Singh, Maity and Nataraj31], therefore, a control loop is considered for engine spool speed as thrust representative. The relationship between thrust and engine speed is non-linear, and based on a lookup table, the thrust command is converted to its equivalent spool speed. This control loop should be designed in such a way that the response of the controller is without overshoot so that while providing the desired speed command, the engine speed does not exceed the maximum allowed value. As the engine speed increases, it’s crucial that the temperature of the turbine remains within the acceptable range to prevent damage to the turbine blades. Measuring the turbine inlet temperature (TIT) directly is challenging due to its extreme heat. Consequently, to avoid exceeding TIT limit, it is managed indirectly through the regulation of the engine EGT. Therefore, an EGT control loop is included in the controller structure and its output is compared with the output of the previous loop and the minimum value is selected. A significant challenge in engine control is the prevention of the dangerous compressor surge phenomenon. Since the risk of compressor surge increases during the engine’s rapid acceleration, restricting the engine’s acceleration can greatly minimise the probability of this occurrence. Limiting the acceleration of the engine can be achieved by limiting the fuel flow rate, therefore, by considering the output of the speed control loop and the EGT control loop as the fuel flow rate, and by using the upper saturation limit for the fuel flow rate, it is possible to prevent rapid acceleration and compressor surge [Reference DeCastro, Litt and Frederick32, Reference Litt, Frederick and Guo33]. Another important challenge of the engine is the possibility of lean blow-out in the conditions of instantaneous throttle reduction. When the throttle is reduced abruptly, the fuel supplied by the power level angle (PLA) control loop is also immediately decreased. This leads to a significant drop in the fuel-to-air ratio and raises the risk of combustion chamber blow-out. To reduce this risk, the engine’s deceleration is controlled. In this circumstance, a lower saturation limit for the fuel flow rate is considered to prevent the instantaneous reduction of the fuel flow and the combustion chamber blow-out. The output of the saturation function after passing through an integrator determines the amount of engine fuel. State feedback regulators are used to control engine speed and EGT temperature in the proposed controller structure.

Figure 5.

The structure of the proposed controller.

In summary, it can be stated that the engine controller is responsible for meeting the required spool speed demand, which represents the thrust value, while also ensuring that the following variables remain within their limits during both steady-state and transient conditions:

(7) \begin{align} \left\{ \begin{array}{l}\Delta EGT\left( {{y_2} \le {r_2}} \right)\\[9pt]\dfrac{{d{W_f}}}{{dt}}\left( {{u_{r\min }} \le {u_r} \le {u_{r\max }}} \right)\end{array} \right. \end{align}

To formulate state feedback regulators, the linear system (1) can be expressed in the following augmented representation:

(8) \begin{align} \begin{array}{*{20}{c}}{\begin{array}{*{20}{c}}{\dot X = AX + B{u_r}} \quad {}{{y_i} = {C_i}X}\end{array}}\\[9pt]{\begin{array}{*{20}{c}}{A = \left[ {\begin{array}{*{20}{c}}a \quad {}b\\0 \quad {}0\end{array}} \right]} {}{\begin{array}{*{20}{c}}{B = \left[ {\begin{array}{*{20}{c}}0\\1\end{array}} \right]} {}{{C_i} = \left[ {\begin{array}{*{20}{c}}{{c_i}} \quad {}{{d_i}}\end{array}} \right]}\end{array}}\end{array}}\end{array} \end{align}

where $X = {\left[ {\begin{array}{*{20}{c}}x \quad {}u\end{array}} \right]^T}$ is the augmented state vector. State feedback control law for each individual control loop is regarded as follows:

(9) \begin{align} {u_{ri}} = - {K_i}X + {p_i}{r_i} \end{align}

where feedback gain ${K_i}$ is established employing state feedback methods to ensure that $A - B{K_i}$ remains stable, while each control loop exhibits satisfactory performance without any overshoot. The constant parameters of ${r_i}$ , for $i = 1,2$ consist of the setpoint for the main output and the limit value for the restricted output, both of which are determined based on the reference command and the characteristics of the system. The ${p_i}$ gains are constant reference prefilter gains aimed at ensuring that ${y_i} = {C_i}X = {r_i}$ in steady-state condition. They are established based on the following relationship.

(10) \begin{align} {p_i} = {[ - C{(A - B{K_i})^{ - 1}}B]^{ - 1}} \end{align}

The following section presents a methodology for analysing the stability of the proposed control system.

4.0 Stability investigation

To investigate the stability of the designed controller depicted in Fig. 5, the controller’s architecture is converted into the standard configuration of Lure’s system [Reference Khalil34]. To achieve this, the initial step involves reconfiguring the plant, as illustrated in Fig. 6, resulting in its transformation into an augmented form. In the new configuration, ${r'_2} = {r_2}{p_2}$ . After that, according to Fig. 7, the Min selector and saturation function are replaced with their nonlinear equivalents [Reference Imani and Montazeri-Gh35]. In order to use the techniques of Lure’s problem [Reference Khalil34], the dynamics between the outputs ${t_i}$ and inputs ${z_i}$ of the nonlinear memory-less functions are extracted. Hence, using Fig. 7, we have

(11) \begin{align} \left\{ \begin{array}{l}{z_1} = - {K_1}X + {K_2}X = ({K_2} - {K_1})X\\[4pt]{z_2} = - {K_1}X - {t_1}\end{array} \right. \end{align}

Figure 6.

Reconfiguration of the controller.

Figure 7.

New depiction of the controller.

The relationship between $X = [x\ \ u]^T$ and ${u_r}$ can be determined from the transfer functions of the augmented plant, which is expressed as follows:

(12) \begin{align} G(s) = {[\begin{array}{*{20}{c}}{{G_1}(s)} \quad {{G_2}(s)}\end{array}]^T} = {\left[\begin{array}{*{20}{c}}{\dfrac{{x(s)}}{{{u_r}(s)}}} \quad {}{\dfrac{{u(s)}}{{{u_r}(s)}}}\end{array}\right]^T} \end{align}

control input rate is established from Fig. 7 based on the following equation:

(13) \begin{align} {u_r} = - {K_1}X - {t_1} - {t_2} \end{align}

by applying Equations (12) and (13), $X = [x\ \ u]^T$ and ${u_r}$ are derived according to the values of ${t_i}$ . Substituting $X$ into Equation (11) reveals the relationship between ${t_i}$ and ${z_i}$ . Consequently, the control system illustrated in Fig. 7 can be represented as shown in Fig. 8, where matrix Q(s) expresses the dynamic relationship,

(14) \begin{align} Q(s) = \frac{{\left[ {\begin{array}{c@{\quad}c}{({K_2} - {K_1})G} & {}{({K_2} - {K_1})G}\\[4pt]1 & {}{ - {K_1}G}\end{array}} \right]}}{{{K_1}G + 1}} \end{align}

Figure 8.

Closed-loop system arranged in a standard configuration.

${K_i}$ are state feedback gains and G is the vector of transfer functions (12).

The updated configuration, referred to as Lure’s system, incorporates a dynamic multivariate block along with static nonlinearities in the feedback loop. Based on the multivariate circle criterion [Reference Khalil34], if Q(s) represents a square and proper transfer function matrix with the following minimal realisation,

(15) \begin{align} \left\{ {\begin{array}{*{20}{c}}{\dot X = \hat AX + \hat BT}\\[5pt]{Z = \hat CX + \hat DT}\end{array}} \right. \end{align}

where $T = - \Lambda (Z)$ , and $\Lambda ( \! \cdot \! )$ is a vector containing static nonlinearities which are locally Lipschitz in Z and this vector satisfies the sector condition defined by real diagonal matrices ${\Omega _1}$ and ${\Omega _2}$ (i.e., ${\left[ {\Lambda (Z) - {\Omega _1}Z} \right]^T}\left[ {\Lambda (Z) - {\Omega _2}Z} \right] \le 0$ ), then the closed-loop system is globally asymptotically stable (or absolutely stable) around the equilibrium point $X = 0$ , if $[I + {\Omega _2}Q]{[I + {\Omega _1}Q]^{ - 1}}$ is strictly positive real. According to the static non-linear functions depicted in Fig. 8, matrix ${\Omega _1}$ is defined as a zero matrix. Consequently, to ensure asymptotic stability, it suffices for matrix $[I + {\Omega _2}Q]$ to be strictly positive real. To check the strictly positive real condition, the following lemma is taken into account [Reference Haddad and Bernstien36]:

Lemma [Reference Haddad and Bernstien36]: Let N(s) be a transfer function matrix with minimal realisation matrices $(E,F,H,J)$ , then the following statements are equivalent:

1. (1) E is a Hurwitz matrix and N(s) is strongly positive real

2. (2) $J + {J^T} \gt 0$ and there exist positive definite matrices P and R such that

   \begin{align*} {E^T}P + PE + {(H - {F^T}P)^T}{(J + {J^T})^{ - 1}}(H - {F^T}P) + R = 0 \\[-25pt]\end{align*}

   It’s important to emphasise that strongly positive real implies strictly positive real [Reference Haddad and Bernstien36].

Utilising the minimum realisation of Q(s), we find that the minimal realisation matrices of $[I + {\Omega _2}Q]$ are $(\hat A,\hat B,{\Omega _2}\hat C,I + {\Omega _2}\hat D)$ . Now, according to the lemma statements, if matrix $\hat A$ is Hurwitz and $2I + {\Omega _2}\hat D + {\hat D^T}{\Omega _2}$ is positive definite, and if there exist positive definite matrices P and R such that fulfill the Riccati Equation (16), then $[I + {\Omega _2}Q]$ is strictly positive real and the multivariate circle condition is satisfied:

(16) \begin{align} {\hat A^T}P + P\hat A + {[{\Omega _2}\hat C - {\hat B^T}P]^T}{[2I + {\Omega _2}\hat D + {\hat D^T}{\Omega _2}]^{ - 1}}[{\Omega _2}\hat C - {\hat B^T}P] + R = 0 \end{align}

Consequently, the control system’s asymptotic stability around the origin $X = 0$ is assured, and $V = {X^T}PX$ is derived as a Lyapunov function for the closed-loop system [Reference Haddad and Bernstien36]. The asymptotic stability of the designed controller is investigated in the following section.