© AFM, EDP Sciences 2019

1 Introduction

A turbofan engine includes various complicated dynamics and behaviors such as the gas flow dynamics in the compressor and turbine, combustion dynamics, the heat transfer dynamics between the gas and the body, the inertia of the shafts, losses, delays and environmental factors [1–3]. Among the dynamics existing in the engine, those resulting from gas dynamics (temperature and pressure), those of sensors and actuators and ignition dynamics are much faster than the others, which are often ignored, but are considered in very accurate and complete models [4,5]. The heat transfer dynamics between the gas and the body (especially in turbines) is also usually ignored in simple models. The engine shafts dynamics are the most important dynamics, which are commonly used, in transient models for controller design [6]. However, there are no specific state equations for the nonlinear dynamics of the engine. Because the characteristic curves of the engine components are generally required for thermodynamic modeling of the engine, this model is highly iterative [7]. Therefore, the thermodynamic model should be linearized at desired operating points to generate linear state-space models for the engine controller design [8].

Many studies have been conducted about the control of turbofan engines using various control methods [9–11]. Recently, model-based controllers are extensively considered by researchers [12–14]. Since, in these controllers, the model parameters are directly present in the designed control signal, the accuracy of the model is of great importance. Sugiyama [15] presented a method to extract linear model matrices and used a series of correction factors to develop a linear model to entirely cover the flight envelope. Also, Chung et al. [16] presented an analytical method to linearize the engine equations in each time step, with the purpose of extracting a precise linear model to be used for design of controllers, especially model-based controllers. However, the linearization of thermodynamic models has commonly been performed using the perturbation method without presenting an analysis about the accuracy of the linear model [17–20]. Careful studies show that no comprehensive study has been published in the literature about the accuracy of different numerical linearization methods.

In this research, the accuracy of different numerical linearization methods that can be used in linearizing the thermodynamic model of a turbofan engine is investigated. This analysis is performed by defining different error indices in the single-input and double-input cases.

2 The engine model

A two-spool high-bypass-ratio turbofan engine is studied in this research. Figure 1 shows a schematic of the engine components and the station numbering of the engine gas path in proper sequence. Its thermodynamic nonlinear structure is similar to TMATS offered by NASA's Glenn Research Center (GRC) [21]. The engine performance modeling is based on the matching of the components operating points. It means that each component operating point should be matched to the others, which is highly iterative. Therefore, it requires several matching guesses of the operating point on the component maps and an equal number of matching constraints. In a turbofan engine dynamic model, the values of inlet air mass flow, bypass ratio and components beta are commonly considered as the matching guesses. Turbofan engine bypass ratio is the ratio of the air mass flow passing through the bypass duct to the air mass flow entering the core. Beta lines are hypothetical lines and the use of beta values is just a mathematical scheme and therefore has no physical interpretation. In Figure 2, beta lines are illustrated on HPC and HPT maps.

The guesses are updated through an iteration process until the matching constraints are satisfied. Generally, this iteration is achieved via matrix solution. In matrix iteration the overall interaction is distinguished and the error equations are solved simultaneously. For this purpose, Jacobian calculation technique is needed to evaluate matrix of partial derivatives, which contains the partial derivatives of the errors in each matching constraint with respect to each matching guess. Also, a numerical method is required in order to update all matching guesses in each iteration. This process is commonly enhanced by Newton-Raphson iterative solving technique. Figure 3 depicts the flow diagram of this methodology. In dynamic modeling, the values of the fuel flow rate and other control variables, which are calculated via the engine control system, are applied as inputs to the model. The above procedure is repeated in each time step.

The details of the design point of the engine are presented in Table 1. The design point of this engine is obtained by assuming sea level static (SLS) and international standard atmosphere (ISA) conditions.

Generally, the state-space equations for the engine dynamic system are as follows:$$\begin{gathered} {\stackrel{\dot{}}{\mathit{N}}}_{\mathit{L}\mathit{P}}=\frac{1}{{\mathit{I}}_{\mathit{L}\mathit{P}}}\times {\mathit{f}}_1({\mathit{N}}_{\mathit{L}\mathit{P}},{\mathit{N}}_{\mathit{H}\mathit{P}},\mathit{u},\mathit{h}) \\ {\stackrel{\dot{}}{\mathit{N}}}_{\mathit{H}\mathit{P}}=\frac{1}{{\mathit{I}}_{\mathit{H}\mathit{P}}}\times {\mathit{f}}_2({\mathit{N}}_{\mathit{L}\mathit{P}},{\mathit{N}}_{\mathit{H}\mathit{P}},\mathit{u},\mathit{h}) \\ {\mathit{y}}_{\mathit{i}}={\mathit{g}}_{\mathit{i}}({\mathit{N}}_{\mathit{L}\mathit{P}},{\mathit{N}}_{\mathit{H}\mathit{P}},\mathit{u},\mathit{h}), \end{gathered}$$(1)where N_LP and N_HP represent the fan speed and core speed, respectively, ${\dot{\mathit{N}}}_{\mathit{L}\mathit{P}}\mathrm{a}\mathrm{n}\mathrm{d}{\mathit{N}}_{\mathit{H}\mathit{P}}$ are the angular accelerations of the low-pressure and high-pressure shafts and I_LP and I_HP are the moments of inertia of the moving parts of the low-pressure and high-pressure spools, respectively. Also, f₁ and f₂ are the net torques produced by low and high pressure turbines, respectively. Vector u contains control input variables. In a turbofan engine, vector u generally includes the fuel mass flow rate, W_f , and other actuators such as VBV, VSV, HPTACC and LPTACC. In addition, vector h contains the health parameters of the engine, which represents the engine health conditions at any time and is a function of engine degradation. In this research, fuel flow rate and the bleed air taken from between the two compressors are considered as the control variables and the engine is assumed to be in a healthy condition. In addition, g_i equations represent the engine outputs and according to the designer requirements, each operating parameter of the engine may be selected as an output.

Fig. 1 Modules and typical stations numbering of a two-spool turbofan engine.

Fig. 2 Characteristic curves of HPC and HPT.

Fig. 3 Matrix iteration approach for transient thermodynamic modeling of two-spool turbofan engine.

3 Linearization of the model

3.1 Linearization using the perturbation methods

A steady-state operating point of an aeroengine is a condition that the values of inputs and flight conditions are constant and power equilibrium between the components of a spool exists. According to the perturbation method, when the engine parameters deviated from their equilibrium conditions, the expansion of equation (1) can be written as follows [22]: $$\begin{gathered} \mathrm{\Delta }{\stackrel{\dot{}}{\mathit{N}}}_{\mathit{L}\mathit{P}=}\frac{1}{{\mathit{I}}_{\mathit{L}\mathit{P}}}\times (\frac{\partial {\mathit{f}}_1}{\partial {\mathit{N}}_{\mathit{L}\mathit{P}}}\mathrm{\Delta }{\mathit{N}}_{\mathit{L}\mathit{P}}+\frac{\partial {\mathit{f}}_1}{\partial {\mathit{N}}_{\mathit{H}\mathit{P}}}\mathrm{\Delta }{\mathit{N}}_{\mathit{H}\mathit{P}}+\frac{\partial {\mathit{f}}_1}{\partial {\mathit{u}}_1}\mathrm{\Delta }{\mathit{u}}_2+\frac{\partial {\mathit{f}}_1}{\partial {\mathit{u}}_2}\mathrm{\Delta }{\mathit{u}}_2+\underset{0}{\mathrm{...}+\underset{}{\frac{\partial {\mathit{f}}_1}{\partial {\mathit{h}}_1}\mathrm{\Delta }{\mathit{h}}_1+\frac{\partial {\mathit{f}}_1}{\partial {\mathit{h}}_2}\mathrm{\Delta }{\mathit{h}}_1+\mathrm{....}}}), \\ \mathrm{\Delta }{\stackrel{\dot{}}{\mathit{N}}}_{\mathit{H}\mathit{P}=}\frac{1}{{\mathit{I}}_{\mathit{H}\mathit{P}}}\times (\frac{\partial {\mathit{f}}_2}{\partial {\mathit{N}}_{\mathit{L}\mathit{P}}}\mathrm{\Delta }{\mathit{N}}_{\mathit{L}\mathit{P}}+\frac{\partial {\mathit{f}}_2}{\partial {\mathit{N}}_{\mathit{H}\mathit{P}}}\mathrm{\Delta }{\mathit{N}}_{\mathit{H}\mathit{P}}+\frac{\partial {\mathit{f}}_2}{\partial {\mathit{u}}_1}\mathrm{\Delta }{\mathit{u}}_1+\frac{\partial {\mathit{f}}_2}{\partial {\mathit{u}}_2}\mathrm{\Delta }{\mathit{u}}_2+\underset{0}{\mathrm{..}+\mathrm{.}\underset{}{\frac{\partial {\mathit{f}}_1}{\partial {\mathit{h}}_1}\mathrm{\Delta }{\mathit{h}}_1+\frac{\partial {\mathit{f}}_1}{\partial {\mathit{h}}_2}\mathrm{\Delta }{\mathit{h}}_1+\mathrm{....}}})\mathrm{,} \\ \mathrm{\Delta }{\mathit{y}}_{\mathit{i}}=\frac{\partial {\mathit{g}}_{\mathrm{i}}}{\partial {\mathit{N}}_{\mathit{L}\mathit{P}}}\mathrm{\Delta }{\mathit{N}}_{\mathit{L}\mathit{P}}+\frac{\partial {\mathit{g}}_{\mathrm{i}}}{\partial {\mathit{N}}_{\mathit{H}\mathit{P}}}\mathrm{\Delta }{\mathit{N}}_{\mathit{H}\mathit{P}}+\frac{\partial {\mathit{g}}_{\mathrm{i}}}{\partial {\mathit{u}}_1}\mathrm{\Delta }{\mathit{u}}_1+\frac{\partial {\mathit{g}}_{\mathrm{i}}}{\partial {\mathit{u}}_2}\mathrm{\Delta }{\mathit{u}}_2+\mathrm{...}+\underset{0}{\underset{}{\frac{\partial {\mathit{g}}_{\mathrm{i}}}{\partial {\mathit{h}}_1}\mathrm{\Delta }{\mathit{h}}_1+\frac{\partial {\mathit{g}}_{\mathrm{i}}}{\partial {\mathit{h}}_2}\mathrm{\Delta }{\mathit{h}}_2+\mathrm{...}}}\mathrm{.} \end{gathered}$$(2)

In order to calculate the approximate values of the partial derivatives, the states and inputs in the equation (2) are slightly deviated from their steady values. In each step, the dependent variable of one partial derivative is perturbed while the other variables are kept constant and the change in the angular accelerations and the required outputs are collected. By dividing these two values and repeating the process for every variable, the value of each partial derivative is obtained one by one. Ultimately, the linear time-invariant model can be presented in the following state-space form:$$\begin{array}{c}\stackrel{\dot{}}{\mathit{x}}=\mathit{A}\mathit{x}+\mathit{B}\mathit{u}\\ \mathit{y}=\mathit{C}\mathit{x}+\mathit{D}\mathit{u}\end{array},$$(3)where $x^T=\left[\begin{array}{cc}\hfill \Delta N_{LP}\hfill & \hfill \Delta N_{HP}\hfill \end{array}\right]$, $y^T=\left[\begin{array}{ccc}\hfill \Delta y_1\hfill & \hfill \Delta y_2\hfill & \hfill \cdots \hfill \end{array}\right]$ and matrices A, B, C and D contain the partial derivatives as below: $$\begin{array}{l}A=\left[\begin{array}{cc}\hfill \frac{\partial f_1}{\partial N_{LP}}\hfill & \hfill \frac{\partial f_1}{\partial N_{HP}}\hfill \\ \hfill \frac{\partial f_2}{\partial N_{LP}}\hfill & \hfill \frac{\partial f_2}{\partial N_{HP}}\hfill \end{array}\right],B=\left[\begin{array}{ccc}\hfill \frac{\partial f_1}{\partial u_1}\hfill & \hfill \frac{\partial f_1}{\partial u_2}\hfill & \hfill \cdots \hfill \\ \hfill \frac{\partial f_2}{\partial u_1}\hfill & \hfill \frac{\partial f_2}{\partial u_2}\hfill & \hfill \cdots \hfill \end{array}\right],\\ C=\left[\begin{array}{cc}\hfill \frac{\partial y_1}{\partial N_{LP}}\hfill & \hfill \frac{\partial y_1}{\partial N_{HP}}\hfill \\ \hfill \frac{\partial y_2}{\partial N_{LP}}\hfill & \hfill \frac{\partial y_2}{\partial N_{HP}}\hfill \\ \hfill ⋮\hfill & \hfill ⋮\hfill \end{array}\right],D=\left[\begin{array}{ccc}\hfill \frac{\partial y_1}{\partial u_1}\hfill & \hfill \frac{\partial y_1}{\partial u_2}\hfill & \hfill \cdots \hfill \\ \hfill \frac{\partial y_2}{\partial u_1}\hfill & \hfill \frac{\partial y_2}{\partial u_2}\hfill & \hfill \cdots \hfill \\ \hfill ⋮\hfill & \hfill ⋮\hfill & \hfill ⋮\hfill \end{array}\right].\end{array}$$(4)

Another way to derive the linearized model with perturbation scheme is central difference method [23,24]. In the central difference perturbation (CDP), the opposite values of deviations are applied to the nonlinear model and the average values of the partial derivatives in equation (4) are computed. In this approach, the components of equation (4) can be derived as follows:$$\begin{array}{l}\frac{\partial f_i}{\partial N_{LP}}=\frac{{\Delta f_i|}_{N_{LP}+\Delta N_{LP}}-{\Delta f_i|}_{N_{LP}-\Delta N_{LP}}}{2\Delta N_{LP}},\frac{\partial f_i}{\partial N_{HP}}=\frac{{\Delta f_i|}_{N_{HP}+\Delta N_{HP}}-{\Delta f_i|}_{N_{HP}-\Delta N_{HP}}}{2\Delta N_{HP}}\\ \frac{\partial y_i}{\partial N_{LP}}=\frac{{\Delta y_i|}_{N_{LP}+\Delta N_{LP}}-{\Delta y_i|}_{N_{LP}-\Delta N_{LP}}}{2\Delta N_{LP}},\frac{\partial y_i}{\partial N_{HP}}=\frac{{\Delta y_i|}_{N_{HP}+\Delta N_{HP}}-{\Delta y_i|}_{N_{HP}-\Delta N_{HP}}}{2\Delta N_{HP}}.\end{array}$$(5)

This method should significantly increase the percentage of compliance (PC) of the linear model.

The dynamic models of the engine, which are valid in all the operating points within the flight envelope and can be used for control engineers, are generally a set of linearized models, each of which is valid in a neighborhood of an operating point in certain flight conditions (Mach, altitude, weather conditions, etc.).

3.2 Linearization using evolutionary optimization algorithms

Another method that can be utilized for the linearization of an engine thermodynamic model is to adjust the parameters of the linear state-space model using evolutionary optimization algorithms [25–28]. Evolutionary algorithms are randomized heuristic search methods based on the principles of natural evolution [29]. Due to the affinity to biology, many biological metaphors are used. Genetic algorithms (GA) are certainly the largest and most popular representative of the class of evolutionary algorithms [29]. The details of the method are described in the next section.

In this approach, in order to encompass all the input frequencies and amplitudes in the desired neighborhood of the operating point, the input signal given to the thermodynamic model should be a persistent excitation. So, the input signal can be an amplitude-modulated pseudo random binary sequence (APRBS) signal or a variable-frequency sinusoidal known as chirp signal. These signals are commonly applied for identification of nonlinear systems. The cost function of the desired optimization algorithm is defined as the difference between responses of the linear and nonlinear systems to the input signal. The cost function can be defined based on the normalized root mean square error (NRMSE). This criterion is described as follows: $$\begin{gathered} \mathit{R}\mathit{M}\mathit{S}\mathit{E}=\sqrt{\frac{{{\sum }_{\mathit{i}=1}^{\mathit{N}}(\mathit{Y}-\stackrel{̂}{\mathit{Y}})}^2}{\mathit{N}}} \\ \mathit{N}\mathit{R}\mathit{M}\mathit{S}\mathit{E}=\frac{\mathit{R}\mathit{M}\mathit{S}\mathit{E}}{{\mathit{Y}}_{\max }-{\mathit{Y}}_{\min }}, \end{gathered}$$(6)where Y is the output signal from the nonlinear model, $\stackrel{̂}{\mathit{Y}}$ is the output signal from the linear model and N is the number of sampling times.

3.3 Linearization using system identification methods

Another approach to obtain a linear model is the use of system identification methods [30,31]. The basis of this approach is based on input-output data. Similar to the previous method, the input signal should be a persistent excitation. First, the input signal is applied to the thermodynamic model and the input and output data are collected. A linear model can then be fitted to these data. The important point that must surely be considered in the linearization of the engine model is omitting the mean value from the data, because the incremental values of the variables about the operating point are considered in the linear models. Moreover, if necessary, the data can be filtered to omit the very high frequency data. The linear model can be extracted in the state-space form or continuous/discrete transfer function.

4 Linearization results

The state variables, inputs and outputs are listed in Table 2. Several outputs can be considered to be present in the linear model, depending on the designer requirements. In this study, some outputs which are generally known as measurable outputs and thrust as one of the important immeasurable outputs [32] are considered. In the single-input engine control problems, the fuel flow rate is assumed as the control input.

4.1 How to implement the linearization

In the perturbation method, the state and input variables are separately perturbed from their steady-state values by 1 and 3%, respectively. In the CDP method the state variables are perturbed ± 1% and the averages are computed.

Moreover, a genetic algorithm is used for tuning the parameters of the linear state-space model. This method is selected because of the large number of elements that should be obtained. In addition, since system identification methods in MATLAB are based on least square optimization, the genetic algorithm is used to compare the effect of different optimization methods on the accuracy of the linearization.

Two types of excitation signals are generated; an APRBS signal during 3000 seconds in the interval of ± 6% of the mean value of the desired operating point fuel flow rate, Figure 4, and the chirp sinusoidal signal, Figure 5. As it is shown in Figure 5, the desired chirp signal is generated with a variable frequency from 0.1 to 2.5 Hz and amplitude of 5% of the fuel flow rate during 45 seconds. Thus, different amplitudes and frequencies are covered by these two signals, so they can be considered as persistent excitations. Another point is that the defined optimization problem is multi-objective because it defines one error for each output. Here, by adding different errors with identical weighted coefficients, the multi-objective optimization problem converts to a single-objective one. The parameters of genetic algorithm are selected as common and are shown in Table 3.

The other technique for linearization is system identification. This approach is performed using the system identification toolbox of MATLAB software [33]. This method is also tested by the same inputs; APRBS excitation signal, Figure 4, and the chirp sinusoidal signal, Figure 5. These signals are given to the thermodynamic model as input and the required outputs are collected, which are shown in Figures 6 and 7. In Figure 6, for more clarity, only 500 of 3000 seconds of the APRBS signal is displayed. After omitting the mean value from the data, the type of the linear model must be selected. In this study, a continuous state-space model of order 2 in a canonical form and with a feedthrough value equal to one (in order to generate matrix D) is selected and the noise matrix is ignored.

Now, to assess the accuracy of linearization, some indices are defined. For this purpose, another APRBS signal is generated and given to both linear and nonlinear models as input. Mean error percentage (mean EP), maximum error percentage (max EP) and the percentage of compliance (PC) is defined as presented in equation (7) [34].$$\begin{gathered} \mathit{M}\mathit{a}\mathit{x}\mathit{E}\mathit{r}\mathit{r}\mathit{o}\mathit{r}=\frac{\mathit{max}\left(|\mathit{Y}-\stackrel{̂}{\mathit{Y}}|\right)}{\mathit{max}\left(\left|\mathit{Y}\right|\right)}\times 100 \\ \mathit{M}\mathit{e}\mathit{a}\mathit{n}\mathit{E}\mathit{r}\mathit{r}\mathit{o}\mathit{r}=\frac{{\displaystyle \sum _{i=1}^{\mathit{N}}\left(|\mathit{Y}(\mathit{i})-\stackrel{̂}{\mathit{Y}}(\mathit{i})|\right)/\mathit{N}}}{\mathit{max}\left(\left|\mathit{Y}\right|\right)}\times 100 \\ \mathit{P}\mathit{C}=(1-\frac{\Vert \mathit{Y}-\stackrel{̂}{\mathit{Y}}\Vert }{\Vert \mathit{Y}-\overline{\mathit{Y}}\Vert })\times 100. \end{gathered}$$(7)

In these equations, Y is the output signal from the nonlinear model, Y shows the mean output signal from the nonlinear model, ${\displaystyle \hat{Y}}$ is the output signal from the linear model and N represents the number of sampling times.

Fig. 4 APRBS signal of the fuel input.

Fig. 5 Chirp signal of the fuel input.

Fig. 6 Output parameters of the thermodynamic model due to APRBS input using for linearization.

Fig. 7 Output parameters of the thermodynamic model due to chirp input using for linearization.

4.2 Case 1: the single-input case – 100% of the LP shaft speed

The design point in SLS and ISA conditions in 100% of the LP shaft speed is considered as the first operating point. In this case, the fuel flow rate is the only control input. For example, the matrices of the linear model obtained from the perturbation method are presented in Table 4.

4.2.1 Effect of the time step and deviation value in the perturbation method

The deviation value and time step of the solver can be regarded as two factors that affect the linearization via the perturbation method. In order to analyze this effect, the deviation value is changed from −0.05 to −0.01% and from 0.01 to 0.05% and the time step from 0.01 to 0.05 seconds. Since matrix A in the state-space form represents the inherent characteristics of the system, changes in the elements of this matrix show the impact of these two factors on the linear model. Figure 8 indicates that time step has no effect on the elements of matrix A, because without changing the initial conditions in each time step the amount of angular acceleration will not be affected by the time step. But, as expected, by changing the deviation value, the values of the elements of matrix A have changed. In addition, the values of the elements are different for the deviation values of the same size with opposite signs. These two facts indicate the nonlinearity of the system. Even with a small deviation value; e.g. ± 1%, this difference will exist, suggesting that a linear model with high precision cannot be expected even in a small neighborhood around this operating point.

In this step, a question may be raised that, which perturbation size will be suitable and more accurate. In order to answer this question, the percentage of compliances (PCs) for different perturbation sizes of the state variables and different outputs are presented in Table 5. According to Table 5, the PCs have not changed significantly for different perturbation sizes. Thus, the deviation value does not have an important effect on the accuracy of the linear model. The selection of the deviation value depends on the system range of operation and the range considered for the operation of the linear model around the operating point.

Fig. 8 Effect of the time step and deviation value in the perturbation method.

4.2.2 Analysis of the indices

The values of the defined indices in the single-input case are outlined in Table 6. First, some considerations should be explained. The linearization results of some of the outputs with different PCs are shown in Figure 9. These outputs are typically displayed in 500 seconds to make a comparison of different PCs. In this regard, if mean EP is below 2% or PC is approximately above 95% it can be called an excellent result, as shown in Figures 9a and 9b. If mean EP is below 7% and PC is approximately above 85% it can be called moderate results, as illustrated in Figures 9c and 9d, and otherwise the results are not good and can be considered as weak results, as presented in Figures 9e and 9f.

Table 6 indicates that the perturbation, GAT-chirp and SID-chirp methods have acted weakly or moderately for the most outputs. Also, the CDP method has significantly improved the ordinary perturbation method accuracy. The PCs in the SID method with APRBS input (SID-APRBS), GA tuning method with APRBS input (GAT-APRBS) and CDP method are very close for all output variables and they have achieved better results compared to other methods. Figure 10a illustrates that the PCs in these methods are excellent in the variables; the LP and HP shaft speeds and temperatures Tt25 and Tt3, moderate in the variable; temperature Tt45 and weak in the variables; pressures Pt25 and Ps3.

In addition, for most of the output variables, especially in the HP shaft speed, the mean EP and the max EP indices show that the linearization errors in the SID-APRBS, GAT-APRBS and CDP methods are much smaller than the ordinary perturbation method, as shown in Figures 10b and 10c. However, for some outputs such as Ps3 and Tt3, it can be seen that the max EPs are almost the same for all methods. Moreover, the mean EPs are much smaller than max EPs which indicates that the error is small in most of the points which leads to a substantial reduction of the mean EP.

Fig 9 Typical results of linearization of some of the outputs with different percentage of compliances.

Fig. 10 Error indices of various linearization methods – Case 1.

4.3 Case 2: the single-input case – 80% of the LP shaft speed

The second operating point is considered in SLS and ISA conditions in 80% of the LP shaft speed. The values of the defined indices in this case are listed in Table 7.

As it is shown in Table 7, similar to the first operating point (case 1), the PCs in the SID-APRBS and GAT-APRBS methods are very close for all output variables and they have achieved better results compared to the other methods. Figure 11a shows that the PCs in these methods are excellent for the variables; the HP shaft speed, temperatures Tt45 and Tt3 and pressure Ps3, moderate for the variable; the LP shaft speed and weak for the variables; pressure Pt25 and temperature Tt25. Also, the CDP method has much better results compared with the ordinary perturbation method. In this case, opposite results are obtained for some variables compared to those in case, that is to say, for the variables in which good PCs are not achieved in case 1, good results are obtained and vice versa. The ordinary perturbation method has also operated very weakly for most outputs except in temperature Tt45, and thus affects the CDP performance, which, despite of case 1, is not good enough. The GAT-chirp method has weak results in this case, too. Compared with the perturbation method, the mean EP and the max EP show that the linearization errors in the SID-APRBS and the GAT-APRBS methods are small for all of the output variables, as illustrated in Figures 11b and 11c.

Fig. 11 Error indices of various linearization methods – Case 2.

4.4 Case 3: the double-input case

In order to analyze the performance of different methods in the double-input case, the bleed air taken from the end of LPC is considered as the second input. To this end, the same conditions of the second operating point are considered, that is, SLS and ISA conditions and 80% of the LP shaft speed, in addition to 2.5% of bleed air. The linearization results in this operating point using different methods are presented in Table 8.

Figure 12a indicates that the SID-APRBS, GAT-APRBS and CDP methods have higher PCs than other methods and the accuracy of the ordinary perturbation method is moderate or weak as before. However, the accuracy of the SID method has decreased compared with the single-input cases, which has made the performance of the perturbation method comparable with it. Therefore, the performance of the CDP method is almost the same as the SID-APRBS method.

In addition, the mean and max EPs in all methods are close and comparable to each other, as presented in Figures 12b and 12c. This matter could be explained due to how these methods work. The identification methods are based on mathematical rules. So, by increasing the number of tuning parameters, the accuracy of these methods decreases. But, perturbation method is based on the physics of the system and it does not relate to the level of complexity of the linearization case due to the number of inputs and outputs. Therefore, the accuracy of the perturbation methods, completely, relates to the location of the operating point on the component maps. As mentioned before, because of the complexity of engine dynamic model which is due to the effect of five characteristic curves (fan, LPC, HPC, HPT and LPT) and multiple loops with iterative numerical solutions, prediction of the behavior of the linear model at each operating point is difficult. Perhaps, by adding bleed to the system, the operating point has been transferred to a point on the maps that the performance of the system around that point is more similar to a linear operation.

Although GAT-APRBS method performance is comparable with the SID-APRBS method but it is very time-consuming. These calculations have been conducted with a computer having core 2 Duo of 2.5 GHz CPU and 4 GB of RAM using MATLAB R2013a. The runtime of the GAT-APRBS method program was about 10 hours in each case while the runtime of the SID-APRBS method was about 17 minutes. Thus, depending on the conditions, the proper numerical linearization method should be selected.

Fig. 12 Error indices of various linearization methods – Case 3.

5 Conclusion

In this paper, the accuracy of different numerical linearization methods for linearizing the dynamic model of turbofan engine was investigated. Since the explicit state equations are not available for the dynamic model of turbofan engine, the nonlinear thermodynamic model of the engine should be linearized at desired operating points for controller design. For this purpose, various methods including the perturbation methods, system identification method and tuning the parameters of the linear model using evolutionary optimization algorithms were employed. According to the analysis, it was illustrated that the ordinary perturbation method had the lowest accuracy and the SID-APRBS, GAT-APRBS and CDP methods generated higher percentage of compliances. Although, in the single-input cases, the SID-APRBS method was more reliable and had priority over the other methods, in the double-input case, the performance of this method decreased and the CDP method was more reliable. On the other hand, the percentage of compliances of the SID-APRBS and GAT-APRBS methods were not very different, but the evolutionary optimization methods are very time-consuming. The analysis indicated that the deviation value in the perturbation method has not a significant effect on the accuracy of the linear model and it should be selected due to the range of operation of the system, as well as the operation range of the linear model around the operating point.

Nomenclature

APRBS: Amplitude Modulated Pseudo Random Binary Sequence

h : Vector of health parameters

HPC: High pressure compressor

HPT: High pressure turbine

HPTACC: HPT active clearance control

I_LP , I_HP : Moment of inertia of LP and HP shafts

LPC: Low pressure compressor

LPT: Low pressure turbine

LPTACC: LPT active clearance control

N_LP , N_HP : LP and HP shaft speeds

N_LP , N_HP : LP and HP shaft angular accelerations

Pt25: HPC inlet total pressure

Ps3: HPC exit static pressure

Tt25: HPC inlet total temperature

Tt3: HPC exit total temperature

Tt45: HPT exit total temperature

u : Vector of control variables

VBV: Variable bleed valve

VSV: Variable stator vane

x : Vector of state variables

y : Vector of output variables

References

All Tables

All Figures

	Fig. 1 Modules and typical stations numbering of a two-spool turbofan engine.
In the text

	Fig. 2 Characteristic curves of HPC and HPT.
In the text

	Fig. 3 Matrix iteration approach for transient thermodynamic modeling of two-spool turbofan engine.
In the text

	Fig. 4 APRBS signal of the fuel input.
In the text

	Fig. 5 Chirp signal of the fuel input.
In the text

	Fig. 6 Output parameters of the thermodynamic model due to APRBS input using for linearization.
In the text

	Fig. 7 Output parameters of the thermodynamic model due to chirp input using for linearization.
In the text

	Fig. 8 Effect of the time step and deviation value in the perturbation method.
In the text

	Fig 9 Typical results of linearization of some of the outputs with different percentage of compliances.
In the text

	Fig. 10 Error indices of various linearization methods – Case 1.
In the text

	Fig. 11 Error indices of various linearization methods – Case 2.
In the text

	Fig. 12 Error indices of various linearization methods – Case 3.
In the text