# Modeling Techniques for a Computational Efficient Dynamic Turbofan Engine Model

Authors: Rory A. Roberts; Scott M. Eastbourn
Journal: International Journal of Aerospace Engineering
 (2014)

---

## Abstract

A transient two-stream engine model has been developed. Individual component models developed exclusively in MATLAB/Simulink including the fan, high pressure compressor, combustor, high pressure turbine, low pressure turbine, plenum volumes, and exit nozzle have been combined to investigate the behavior of a turbofan two-stream engine. Special attention has been paid to the development of transient capabilities throughout the model, increasing physics model, eliminating algebraic constraints, and reducing simulation time through enabling the use of advanced numerical solvers. The lessening of computation time is paramount for conducting future aircraft system-level design trade studies and optimization. The new engine model is simulated for a fuel perturbation and a specified mission while tracking critical parameters. These results, as well as the simulation times, are presented. The new approach significantly reduces the simulation time.


---

## Body

## 1. Introduction

Conceptual design groups have traditionally designed aircraft from a subsystem-level viewpoint. Consequently, subsystems such as the propulsion, electrical, and thermal management systems are often optimized without consideration of vehicle-level interactions, resulting in a final aircraft design that is not truly optimized. It is believed, however, that vehicle-level analysis of subsystem interactions could result in significant performance gains across the aircraft, potentially improving the overall effectiveness of future platforms. The development of a vehicle level tip-to-tail (T2T) modeling and simulation tool would allow these performance gains to be quantified in a cost effective manner [1, 2].Recent work completed by the authors focused on the development of a nonproprietary, thermal T2T aircraft model in Simulink [3]. At least some of the motivation for building a T2T model was to perform design trade studies. In order to run these trade studies effectively and efficiently, computation times should be no slower than real time. The T2T model is intended to stimulate the optimization of individual subsystems for improving overall vehicle-level performance and mitigating the thermal and power challenges of future aircraft platforms [4]. In addition, the nonproprietary nature of the model allows the tool to be distributed to various conceptual design groups and researchers. Specifically, it is foreseen that conceptual designers will use the model to conduct design trade studies, allowing the analysis of multiple design configurations and the resulting subsystem interactions in short time periods [5–7]. In order for effective trade studies to be conducted, the model must have relatively fast computation times. Previous work has demonstrated that while effective and accurate, the developed T2T model has extremely large simulation times of half real-time. As a result, the tool would fail to meet a major requirement for conducting valuable design trade studies and optimization in a practical time frame. Through further investigation, it was identified that the turbofan engine model was contributing to slow simulation times. The T2T aircraft model is a stiff system with time constants varying between milliseconds to decaseconds. The turbofan engine has many interfaces with other subsystems and must be compatible with the various time constants captured in the T2T model.Gas turbine engine modeling tools such as Numerical Propulsion System Simulation (NPSS) tool provide a top-down propulsion system approach to provide designers with a tool to incorporate the relevant factors which affect propulsion performance early in the design and analysis process [8]. NPSS was primarily a steady-state tool that has been widely adopted by government and industry over the past decades but has expanded its capabilities in dynamic simulation. NPSS can be expanded to incorporate systems beyond the main engine. The NPSS based models were inefficient when utilized in Matlab-Simulink based T2T model simulations [5] with simulation times 1/6 of real time. Surrogate models were used in the T2T model to reduce computational times [6]. Dynamic models have been developed for real time operation with unsteady effects due to mass accumulation considered by adding a plenum between each compressor and turbine stage and shaft dynamics [9–11]. These models capture some of the dynamics of the engine but assume mass flows from previous time-steps in different sections of the engine, that is, the compressor mass flow. Other work in dynamic modeling captured the plenum volume dynamics needed to predict the secondary flows within the engine for blade cooling [12].As previously mentioned, using NPSS models with Simulink is less efficient. The computational inefficiencies reside in the compiling of source code within the Matlab-Simulink framework. These subtle inefficiencies do not pose issues when executing a limited amount of trade studies or simulations but do pose a problem when performing thousands of simulations for an optimization routine. In this work, an engine model developed entirely in Matlab-Simulink environment is presented. Several techniques were attempted to increase the simulation speed of the turbofan engine model. The techniques adopted to achieve a computational efficient turbofan engine model are presented in this work. The turbofan engine model is developed without iteration loops (algebraic constraints) and all states are continuous. This approach is very important for complex system level simulations of stiff dynamic systems. By modeling all the significant states as continuous states and not steady-state approximations with discontinuities, advanced numerical stiff solvers for stiff systems may be used. Numerical stiff solvers rely on the Jacobian matrix and thus require accurate approximations for gradients of all continuous states. Stiff solvers dramatically reduced the computational time for simulating stiff systems. Also, by having the engine model constructed completely in Simulink, a complete T2T model can be compiled then executed or exported as an efficient executable limiting unnecessary callbacks.The following sections provide a detailed description of the modeling approach for the transient engine model. Simulation results are presented along with the comparison in simulation times. Key approaches in performance map interpretation and plenum volume dynamics enable the development of a computationally efficient engine model.


## 2. Fundamental Equation Development

The engine consists of several key component models presented in Figure1. These models include: Fan, High Pressure (HP) Compressor, Combustor, High Pressure (HP) Turbine, Low Pressure (LP) Turbine, Bypass Plenum Volume, Nozzle, High Pressure (HP) Shaft, and Low Pressure (LP) Shaft. Each of the turbo-machinery components in the engine share several common equations defining the physics of the model. Equations common to the fan, HP compressor, HP turbine, and LP turbine, are outlined in the following section, including mass flow, temperature, and work terms.Figure 1
Two-stream turbofan engine diagram of major components.


### 2.1. Mass Flow

Each turbo-machinery model contains a performance map that determines a corrected mass flow for a given shaft speed and pressure ratio. The traditional approach of interpreting the performance maps is to create an additional independent variable known as theR-line [13, 14]. Figure 2 shows the compressor map with R-lines included. Typically the first R-line is given a value of 1.0 that coincides with the surge line. The remaining R-lines are roughly parallel to the surge line with increasing R-line value corresponding to a higher surge margin. The R-lines have no physical meaning, but are simply mathematical constructs orthogonal to corrected speed so that any value of R-line and corrected speed will have a unique compressor operating point. The engine model did not adopt the R-line methodology due to the need for additional iterations (algebraic constraint) within the model.Figure 2
Compressor map withR-lines and approximations near the surge line.The adopted approach adjusts the original compressor maps to eliminate the curl back of the speed lines near the surge line. A minimum slope is assumed for the pressure ratio versus corrected mass flow in order to accommodate two-dimensional interpolation in lookup tables as shown in Figure2. This approach introduces some error near the surge line, but it is assumed that the engine will not operate in this region for all valid designs. If operations near surge were part of the study, then the approximations near the surge line would introduce excessive uncertainty.The maps are represented by two-dimensional lookup tables that contain a predetermined matrix for the specific turbo-machine being used. Row and column vectors are also defined within the map, allowing interpolation within the matrix based on the input signals to the lookup table. These input signals are normalized speeds and pressure ratio, shown below by (1) and (2) respectively. Using these two normalized signals, the performance map interpolates within the predefined matrix an output of a normalized mass flow rate. This normalized mass flow rate is used to calculate an actual mass flow rate using (3):
(1)Pr,normalized=PoutPinPr,design,(2)Nnormalized=(NTin)(Tin,designNdesign),(3)m˙=m˙normalized(m˙designTin,designPin,design)(PinTin).


### 2.2. Temperature

Each turbo-machinery model contains a performance map that determines an efficiency for a given shaft speed and pressure ratio. In a manner similar to the mass flow rate performance map, the efficiency performance map contains matrix defining efficiencies for predetermined shaft speeds and pressure ratios. The normalized signals for pressure ratio and shaft speed are shown by (1) and (2), respectively. The efficiency term yielded from the performance map is then used to calculate the outlet temperature for the compressor and turbine models, shown by
(4)Tout,C=Tin(1+1ηC((PoutPin)((k-1)/k)-1)),(5)Tout,T=Tin(1+ηT((PoutPin)((k-1)/k)-1)).


### 2.3. Work

The power (work rate) absorbed or produced is based on the outlet mass flow rates as well as the inlet and outlet temperatures for each of the turbo-machine models. The fan and compressor models consume power (negative) while the HP and LP turbines produce power (positive). The inlet and outlet temperatures of each model are used to calculate an enthalpy value using (6). These inlet and outlet enthalpies are combined with the outlet mass flow rate to calculate the power for the compressor and turbine models, as shown by (7) and (8), respectively,
(6)h=∫0TCp(T)dT,(7)W˙C=m˙C(hin-hout),(8)W˙T=m˙T(hin-hout).


### 2.4. Plenum Volume Dynamics

Another important technique in increasing the simulation speed of the engine model was the implementation of an isentropic plenum volume dynamics using (9). Plenum volume dynamics provide a continuous solution for the pressures within the engine model [15, 16]. For simplicity and efficiency, typically plenum volumes are assumed to be steady-state with the flow rate entering the volume equal to the flow rate leaving the volume. A steady-state assumption for the plenum volumes is valid in the sense the time constants of the plenum volumes are negligible when compared to the shaft and thermal time constants. On the contrary, the engine model has higher simulation speeds with dynamic plenum volumes while adding more physics. The increase in simulation speed is attributed to the numerical methods employed within the advanced stiff solvers in Simulink. Stiff solvers utilize the Jacobian matrix, which consists of the partial derivatives of all the states. The Jacobian matrix is used to estimate the new states for the next time step [17]. If a very sensitive parameter such as operating pressures within the engine is assumed steady-state, then the pressures are not states and excluded from the Jacobian matrix. The numerical solver is unable to account for the changes in pressure and the gradients of all the states with respect to pressure. The numerical solver will have a reduced time step increasing computational time and numerical round off errors during changes in pressure:
(9)Pout=∫(m˙in-m˙out)RoutTVdt.


## 2.1. Mass Flow

Each turbo-machinery model contains a performance map that determines a corrected mass flow for a given shaft speed and pressure ratio. The traditional approach of interpreting the performance maps is to create an additional independent variable known as theR-line [13, 14]. Figure 2 shows the compressor map with R-lines included. Typically the first R-line is given a value of 1.0 that coincides with the surge line. The remaining R-lines are roughly parallel to the surge line with increasing R-line value corresponding to a higher surge margin. The R-lines have no physical meaning, but are simply mathematical constructs orthogonal to corrected speed so that any value of R-line and corrected speed will have a unique compressor operating point. The engine model did not adopt the R-line methodology due to the need for additional iterations (algebraic constraint) within the model.Figure 2
Compressor map withR-lines and approximations near the surge line.The adopted approach adjusts the original compressor maps to eliminate the curl back of the speed lines near the surge line. A minimum slope is assumed for the pressure ratio versus corrected mass flow in order to accommodate two-dimensional interpolation in lookup tables as shown in Figure2. This approach introduces some error near the surge line, but it is assumed that the engine will not operate in this region for all valid designs. If operations near surge were part of the study, then the approximations near the surge line would introduce excessive uncertainty.The maps are represented by two-dimensional lookup tables that contain a predetermined matrix for the specific turbo-machine being used. Row and column vectors are also defined within the map, allowing interpolation within the matrix based on the input signals to the lookup table. These input signals are normalized speeds and pressure ratio, shown below by (1) and (2) respectively. Using these two normalized signals, the performance map interpolates within the predefined matrix an output of a normalized mass flow rate. This normalized mass flow rate is used to calculate an actual mass flow rate using (3):
(1)Pr,normalized=PoutPinPr,design,(2)Nnormalized=(NTin)(Tin,designNdesign),(3)m˙=m˙normalized(m˙designTin,designPin,design)(PinTin).


## 2.2. Temperature

Each turbo-machinery model contains a performance map that determines an efficiency for a given shaft speed and pressure ratio. In a manner similar to the mass flow rate performance map, the efficiency performance map contains matrix defining efficiencies for predetermined shaft speeds and pressure ratios. The normalized signals for pressure ratio and shaft speed are shown by (1) and (2), respectively. The efficiency term yielded from the performance map is then used to calculate the outlet temperature for the compressor and turbine models, shown by
(4)Tout,C=Tin(1+1ηC((PoutPin)((k-1)/k)-1)),(5)Tout,T=Tin(1+ηT((PoutPin)((k-1)/k)-1)).


## 2.3. Work

The power (work rate) absorbed or produced is based on the outlet mass flow rates as well as the inlet and outlet temperatures for each of the turbo-machine models. The fan and compressor models consume power (negative) while the HP and LP turbines produce power (positive). The inlet and outlet temperatures of each model are used to calculate an enthalpy value using (6). These inlet and outlet enthalpies are combined with the outlet mass flow rate to calculate the power for the compressor and turbine models, as shown by (7) and (8), respectively,
(6)h=∫0TCp(T)dT,(7)W˙C=m˙C(hin-hout),(8)W˙T=m˙T(hin-hout).


## 2.4. Plenum Volume Dynamics

Another important technique in increasing the simulation speed of the engine model was the implementation of an isentropic plenum volume dynamics using (9). Plenum volume dynamics provide a continuous solution for the pressures within the engine model [15, 16]. For simplicity and efficiency, typically plenum volumes are assumed to be steady-state with the flow rate entering the volume equal to the flow rate leaving the volume. A steady-state assumption for the plenum volumes is valid in the sense the time constants of the plenum volumes are negligible when compared to the shaft and thermal time constants. On the contrary, the engine model has higher simulation speeds with dynamic plenum volumes while adding more physics. The increase in simulation speed is attributed to the numerical methods employed within the advanced stiff solvers in Simulink. Stiff solvers utilize the Jacobian matrix, which consists of the partial derivatives of all the states. The Jacobian matrix is used to estimate the new states for the next time step [17]. If a very sensitive parameter such as operating pressures within the engine is assumed steady-state, then the pressures are not states and excluded from the Jacobian matrix. The numerical solver is unable to account for the changes in pressure and the gradients of all the states with respect to pressure. The numerical solver will have a reduced time step increasing computational time and numerical round off errors during changes in pressure:
(9)Pout=∫(m˙in-m˙out)RoutTVdt.


## 3. Overview of Model Components

The one-dimensional, turbofan engine model has incorporated a lot of detail. The model includes detailed performance maps, pressure drop, plenum volume dynamics, thermal transients, chemical reactions, parallel flow paths with mass flows dependent on local densities and pressure loses and shaft dynamics.Detailed descriptions of each of the component models as well as the unique equations used to model the appropriate physics are covered in the following sections.


### 3.1. Compressor

Located at the front of the engine in Figure1, the fan is responsible for drawing air into the engine. The LP compressor, also known as the fan, is driven by the LP shaft and compresses the air entering the engine. Some of this compressed air then enters the HP compressor (core stream) where it will be compressed even further, but the majority of the fan air enters the bypass plenum volume (bypass stream). Within the fan model, several key equations are modeled to describe the relevant physics. In addition to the common equations described in the previous section, the fan has a unique inlet pressure. For subsonic conditions, the inlet pressure is found by calculating the total pressure at the front of the aircraft, as shown by (10). Equation (10) assumes 100% pressure recovery for the inlet diffuser. It is also worth mentioning that the outlet pressure term for the fan is represented by the bypass plenum volume pressure, which is outlined in the bypass component section:
(10)Pin=Pambient+12ρambient(MkambientRambientTambient)2.Air from the fan that does not enter the bypass plenum volume is sent to the HP compressor in Figure1. The HP compressor increases the core air pressure to its largest value before it enters the combustor. The HP compressor is driven by the HP shaft, which is powered by the HP turbine. The outlet pressure is provided by the combustor and will be discussed in the combustor section. The inlet pressure is equivalent to the bypass plenum volume pressure and will be discussed in the bypass component section. The HP compressor has bleed air extracted at the exit. After the actual outlet mass flow rate is calculated based on Section 2.1, bleed air is removed to cool turbine blades (secondary flow) and power additional systems within the aircraft.


### 3.2. Combustor

The combustor in the center of Figure1 receives an air stream from the HP compressor as well as a fuel stream of JP-8. Energy balances are used to determine the temperature and composition of the outgoing air stream. This mixture is sent to the HP turbine. It is assumed that complete combustion of the JP-8 fuel occurs, yielding CO2, H2O, and N2 as the sole products of the reaction. The JP-8 combustion equation is expressed by (11) [18]. The heat of reaction and the enthalpy flow are calculated using (12) and (13), respectively. The temperature of the combustor outlet stream can be found using (14). Equation (15) yields the molar flow rate and the molar composition by (16) for combustion of the air and fuel streams. Equation (17) provides the reaction vector derived from (11):
(11)C10.3H20.5+15.425[O2+3.76N2]⟶10.3CO2+10.25H2O+57.998N2,(12)Q˙Rx=∑N˙R(hf,R)-∑N˙P(hf,P),(13)hout=Q˙in+Q˙RxN˙out,(14)Toutlet=houtCpout,(15)N˙out=N˙in+∑R,(16)VcvCdXoutdt=N˙in(Xin-Xout)-Xout∑R+R,(17)R=N˙JP-8[rJP-8rCOrCO2rH2rH2OrN2rO2]=N˙JP-8[-1010.3010.250-15.425].


### 3.3. Turbine

The high pressure (HP) turbine receives the combustor outlet mixture shown in Figure1. Power generated by the turbine is used to apply a torque to the HP shaft, which then drives the HP compressor. The plenum volume located between the HP compressor and the HP turbine inlet is modeled within the HP turbine to derive the HP turbine inlet pressure. The plenum volume primarily accounts for the volume of the combustor. The mass flow rate entering this plenum volume is known from the HP compressor model. The outlet mass flow rate of the HP turbine is specified by the performance map. With the incoming and outgoing mass flows of the plenum volume known, the dynamic pressure of the plenum volume can be calculated via (9).As air enters the HP turbine, a secondary air stream for blade cooling is added thus reducing the temperature of the core air. This secondary stream is fed by the bleed air removed at the HP compressor exit. Within the HP turbine model, a subsystem exists to calculate the flow rate of bleed air which cools the HP turbine inlet as well as the flow rate of air that continues on to the LP turbine. The bleed mass flow rate calculations are shown by (18) and (19), respectively. As previously mentioned, the HP turbine bleed flow is mixed with core air from the combustor outlet before it enters the HP turbine in order to provide cooling. Two calculations are required to determine the resulting mass flow rate as well as the temperature of the newly formed mixture that enters the HP turbine. The required calculations for the mass flow rate and temperature signals entering the HP turbine are shown by (20) and (21), respectively,
(18)m˙HPT,bleed=xHPTm˙bleed,(19)m˙LPT,bleed=(1-xHPT)m˙bleed,(20)m˙in,HPT=m˙HPT,bleed+m˙combustor,(21)Tin,HPT=THPC,bleedm˙HPT,bleedm˙in,HPT+Tcombustorm˙combustorm˙in,HPT.After core air exits the HP turbine, it enters the low pressure (LP) turbine as shown in Figure1. The LP turbine produces power that drives the LP shaft, which in turn drives the fan. The inlet pressure is found in a similar fashion to the inlet pressure of the HP turbine, as shown by (9). The LP turbine bleed air mass flow rate is already known from (19). The inlet mass flow rate, which includes the core air from the HP turbine outlet as well as the LP turbine bleed air, and the inlet temperature are found using a method equivalent to (20) and (21), respectively.


### 3.4. Bypass Plenum Volume

The bypass model determines bypass flow rate and pressure. The bypass labeled in Figure1 is comprised of the void space around the HP compressor, combustor, HP turbine, LP turbine, and shafts. The air that bypasses the HP compressor, combustor, HP turbine, and LP turbine travels through a bypass duct and enters a mixer plenum volume at the nozzle inlet. The majority of the fan mass flow enters the bypass rather than the HP compressor. The pressure drop across the bypass plenum volume drives the amount of mass flow that bypasses the core of the engine. This mass flow rate is represented by (22) [19]. Equation (22) assumes low Mach number flow incompressible flow. The dynamic pressure of the plenum volume is found using (9):
(22)m˙bypass=CdAbypass2ρ(Pin-Pout).


### 3.5. Nozzle

The nozzle is the final component in a turbofan engine flow path shown on the right of Figure1. A converging-diverging nozzle creates the thrust needed to propel the aircraft forward. Air from the LP turbine outlet and the bypass plenum volume are combined in the mixer volume before entering the nozzle. The temperature of the mixed stream is shown by (23). The pressure of the mixer volume is found using (9). Within the actual nozzle, two cases can exist: Choked Flow or Non-Choked Flow [20]. To determine which case is occurring at a given time, a critical pressure ratio is found using (24) through (30). To determine the thrust exit velocity and mass flow rate must be calculated:
(23)Tmixer=∫QnetmVCp,outdt,(24)(PoutPin)critical=(2k+1)k/(k-1).


#### 3.5.1. Choked Flow

Choked flow occurs when the actual nozzle pressure ratio is less than the critical pressure ratio, while nonchoked flow occurs for pressure ratios larger than the critical value. When the nozzle model has determined that the flow is choked, the exit mass flow rate is shown by (25). The exit temperature is calculated using (26). The speed of sound and exit velocity are found using (27) and (28), respectively,
(25)m˙out=PmixerAthroatkRTmixer(2k+1)((k+1)/2(k-1)),(26)Tout=Tmixer1+((k-1)/2),(27)cout=kRTout,(28)vout=cout.


#### 3.5.2. Non-Choked Flow

When the flow is not choked, several terms must first be derived, including the exit Mach number, the exit temperature, and the speed of sound at the nozzle exit. The nozzle exit Mach number is found by (29). Using the exit Mach number, the exit temperature is calculated using (30). With the exit temperature known, the speed of sound at the nozzle exit is calculated as in the choked case with (27). The velocity of air exiting the nozzle can be found using the speed of sound and Mach number results, as shown by (31). The density of air exiting the nozzle is found using (32). Finally, exit mass flow of the nozzle is found using (33):
(29)Mout=(2k-1)[(PmixerPout)((k-1)/k)-1],(30)Tout=Tmixer1+Mout2((k-1)/2),(31)vout=Moutcout,(32)ρout=PoutRoutTout,(33)m˙out=MoutcoutρoutAout.


#### 3.5.3. Thrust

Once the nozzle mass flow rate and exit velocity have been established, for choked or nonchoked flow, the thrust is calculated. The thrust produced by the engine is based on the mass flow entering and exiting the engine, as well as the pressure difference between the nozzle and ambient air. The inlet mass flow rate, which is equivalent to the fan mass flow rate, is already known. The inlet velocity, however, must be calculated using (34). Using the inlet velocity, the total engine thrust is represented by (35) for a convergent nozzle:
(34)vinlet=MaircraftkambientRambientTambient,(35)Thrust=(m˙outvout-m˙invin)+Ain(Pout-Pambient).


### 3.6. Shaft

The HP shaft connects the HP Turbine and the HP Compressor. Power from the HP turbine is transferred by the HP shaft to drive the HP compressor. The HP turbine work signal represents a positive load and the HP compressor work signal represents a negative load. The HP shaft speed is the only calculation performed within the shaft model and is represented by (36). Additional auxiliary loads may be added to the shaft (not shown):
(36)NHP_shaft=30π∫W˙C+W˙HPTJHP,shaftωHP,shaft.The LP shaft connects the LP Turbine and the fan. Power from the LP turbine is transferred by the LP shaft to drive the fan. The LP turbine work signal represents a positive load and the fan work signal represents a negative load. The LP shaft speed is the only calculation performed within the model and is represented by (37). The LP shaft has additional auxiliary loads, pumps, and generator, which are inputs to the engine model:
(37)NLP,shaft=30π∫W˙Fan+W˙LPT+W˙Pumps+W˙GeneratorJLP,shaftωLP,shaft.


## 3.1. Compressor

Located at the front of the engine in Figure1, the fan is responsible for drawing air into the engine. The LP compressor, also known as the fan, is driven by the LP shaft and compresses the air entering the engine. Some of this compressed air then enters the HP compressor (core stream) where it will be compressed even further, but the majority of the fan air enters the bypass plenum volume (bypass stream). Within the fan model, several key equations are modeled to describe the relevant physics. In addition to the common equations described in the previous section, the fan has a unique inlet pressure. For subsonic conditions, the inlet pressure is found by calculating the total pressure at the front of the aircraft, as shown by (10). Equation (10) assumes 100% pressure recovery for the inlet diffuser. It is also worth mentioning that the outlet pressure term for the fan is represented by the bypass plenum volume pressure, which is outlined in the bypass component section:
(10)Pin=Pambient+12ρambient(MkambientRambientTambient)2.Air from the fan that does not enter the bypass plenum volume is sent to the HP compressor in Figure1. The HP compressor increases the core air pressure to its largest value before it enters the combustor. The HP compressor is driven by the HP shaft, which is powered by the HP turbine. The outlet pressure is provided by the combustor and will be discussed in the combustor section. The inlet pressure is equivalent to the bypass plenum volume pressure and will be discussed in the bypass component section. The HP compressor has bleed air extracted at the exit. After the actual outlet mass flow rate is calculated based on Section 2.1, bleed air is removed to cool turbine blades (secondary flow) and power additional systems within the aircraft.


## 3.2. Combustor

The combustor in the center of Figure1 receives an air stream from the HP compressor as well as a fuel stream of JP-8. Energy balances are used to determine the temperature and composition of the outgoing air stream. This mixture is sent to the HP turbine. It is assumed that complete combustion of the JP-8 fuel occurs, yielding CO2, H2O, and N2 as the sole products of the reaction. The JP-8 combustion equation is expressed by (11) [18]. The heat of reaction and the enthalpy flow are calculated using (12) and (13), respectively. The temperature of the combustor outlet stream can be found using (14). Equation (15) yields the molar flow rate and the molar composition by (16) for combustion of the air and fuel streams. Equation (17) provides the reaction vector derived from (11):
(11)C10.3H20.5+15.425[O2+3.76N2]⟶10.3CO2+10.25H2O+57.998N2,(12)Q˙Rx=∑N˙R(hf,R)-∑N˙P(hf,P),(13)hout=Q˙in+Q˙RxN˙out,(14)Toutlet=houtCpout,(15)N˙out=N˙in+∑R,(16)VcvCdXoutdt=N˙in(Xin-Xout)-Xout∑R+R,(17)R=N˙JP-8[rJP-8rCOrCO2rH2rH2OrN2rO2]=N˙JP-8[-1010.3010.250-15.425].


## 3.3. Turbine

The high pressure (HP) turbine receives the combustor outlet mixture shown in Figure1. Power generated by the turbine is used to apply a torque to the HP shaft, which then drives the HP compressor. The plenum volume located between the HP compressor and the HP turbine inlet is modeled within the HP turbine to derive the HP turbine inlet pressure. The plenum volume primarily accounts for the volume of the combustor. The mass flow rate entering this plenum volume is known from the HP compressor model. The outlet mass flow rate of the HP turbine is specified by the performance map. With the incoming and outgoing mass flows of the plenum volume known, the dynamic pressure of the plenum volume can be calculated via (9).As air enters the HP turbine, a secondary air stream for blade cooling is added thus reducing the temperature of the core air. This secondary stream is fed by the bleed air removed at the HP compressor exit. Within the HP turbine model, a subsystem exists to calculate the flow rate of bleed air which cools the HP turbine inlet as well as the flow rate of air that continues on to the LP turbine. The bleed mass flow rate calculations are shown by (18) and (19), respectively. As previously mentioned, the HP turbine bleed flow is mixed with core air from the combustor outlet before it enters the HP turbine in order to provide cooling. Two calculations are required to determine the resulting mass flow rate as well as the temperature of the newly formed mixture that enters the HP turbine. The required calculations for the mass flow rate and temperature signals entering the HP turbine are shown by (20) and (21), respectively,
(18)m˙HPT,bleed=xHPTm˙bleed,(19)m˙LPT,bleed=(1-xHPT)m˙bleed,(20)m˙in,HPT=m˙HPT,bleed+m˙combustor,(21)Tin,HPT=THPC,bleedm˙HPT,bleedm˙in,HPT+Tcombustorm˙combustorm˙in,HPT.After core air exits the HP turbine, it enters the low pressure (LP) turbine as shown in Figure1. The LP turbine produces power that drives the LP shaft, which in turn drives the fan. The inlet pressure is found in a similar fashion to the inlet pressure of the HP turbine, as shown by (9). The LP turbine bleed air mass flow rate is already known from (19). The inlet mass flow rate, which includes the core air from the HP turbine outlet as well as the LP turbine bleed air, and the inlet temperature are found using a method equivalent to (20) and (21), respectively.


## 3.4. Bypass Plenum Volume

The bypass model determines bypass flow rate and pressure. The bypass labeled in Figure1 is comprised of the void space around the HP compressor, combustor, HP turbine, LP turbine, and shafts. The air that bypasses the HP compressor, combustor, HP turbine, and LP turbine travels through a bypass duct and enters a mixer plenum volume at the nozzle inlet. The majority of the fan mass flow enters the bypass rather than the HP compressor. The pressure drop across the bypass plenum volume drives the amount of mass flow that bypasses the core of the engine. This mass flow rate is represented by (22) [19]. Equation (22) assumes low Mach number flow incompressible flow. The dynamic pressure of the plenum volume is found using (9):
(22)m˙bypass=CdAbypass2ρ(Pin-Pout).


## 3.5. Nozzle

The nozzle is the final component in a turbofan engine flow path shown on the right of Figure1. A converging-diverging nozzle creates the thrust needed to propel the aircraft forward. Air from the LP turbine outlet and the bypass plenum volume are combined in the mixer volume before entering the nozzle. The temperature of the mixed stream is shown by (23). The pressure of the mixer volume is found using (9). Within the actual nozzle, two cases can exist: Choked Flow or Non-Choked Flow [20]. To determine which case is occurring at a given time, a critical pressure ratio is found using (24) through (30). To determine the thrust exit velocity and mass flow rate must be calculated:
(23)Tmixer=∫QnetmVCp,outdt,(24)(PoutPin)critical=(2k+1)k/(k-1).


### 3.5.1. Choked Flow

Choked flow occurs when the actual nozzle pressure ratio is less than the critical pressure ratio, while nonchoked flow occurs for pressure ratios larger than the critical value. When the nozzle model has determined that the flow is choked, the exit mass flow rate is shown by (25). The exit temperature is calculated using (26). The speed of sound and exit velocity are found using (27) and (28), respectively,
(25)m˙out=PmixerAthroatkRTmixer(2k+1)((k+1)/2(k-1)),(26)Tout=Tmixer1+((k-1)/2),(27)cout=kRTout,(28)vout=cout.


### 3.5.2. Non-Choked Flow

When the flow is not choked, several terms must first be derived, including the exit Mach number, the exit temperature, and the speed of sound at the nozzle exit. The nozzle exit Mach number is found by (29). Using the exit Mach number, the exit temperature is calculated using (30). With the exit temperature known, the speed of sound at the nozzle exit is calculated as in the choked case with (27). The velocity of air exiting the nozzle can be found using the speed of sound and Mach number results, as shown by (31). The density of air exiting the nozzle is found using (32). Finally, exit mass flow of the nozzle is found using (33):
(29)Mout=(2k-1)[(PmixerPout)((k-1)/k)-1],(30)Tout=Tmixer1+Mout2((k-1)/2),(31)vout=Moutcout,(32)ρout=PoutRoutTout,(33)m˙out=MoutcoutρoutAout.


### 3.5.3. Thrust

Once the nozzle mass flow rate and exit velocity have been established, for choked or nonchoked flow, the thrust is calculated. The thrust produced by the engine is based on the mass flow entering and exiting the engine, as well as the pressure difference between the nozzle and ambient air. The inlet mass flow rate, which is equivalent to the fan mass flow rate, is already known. The inlet velocity, however, must be calculated using (34). Using the inlet velocity, the total engine thrust is represented by (35) for a convergent nozzle:
(34)vinlet=MaircraftkambientRambientTambient,(35)Thrust=(m˙outvout-m˙invin)+Ain(Pout-Pambient).


## 3.5.1. Choked Flow

Choked flow occurs when the actual nozzle pressure ratio is less than the critical pressure ratio, while nonchoked flow occurs for pressure ratios larger than the critical value. When the nozzle model has determined that the flow is choked, the exit mass flow rate is shown by (25). The exit temperature is calculated using (26). The speed of sound and exit velocity are found using (27) and (28), respectively,
(25)m˙out=PmixerAthroatkRTmixer(2k+1)((k+1)/2(k-1)),(26)Tout=Tmixer1+((k-1)/2),(27)cout=kRTout,(28)vout=cout.


## 3.5.2. Non-Choked Flow

When the flow is not choked, several terms must first be derived, including the exit Mach number, the exit temperature, and the speed of sound at the nozzle exit. The nozzle exit Mach number is found by (29). Using the exit Mach number, the exit temperature is calculated using (30). With the exit temperature known, the speed of sound at the nozzle exit is calculated as in the choked case with (27). The velocity of air exiting the nozzle can be found using the speed of sound and Mach number results, as shown by (31). The density of air exiting the nozzle is found using (32). Finally, exit mass flow of the nozzle is found using (33):
(29)Mout=(2k-1)[(PmixerPout)((k-1)/k)-1],(30)Tout=Tmixer1+Mout2((k-1)/2),(31)vout=Moutcout,(32)ρout=PoutRoutTout,(33)m˙out=MoutcoutρoutAout.


## 3.5.3. Thrust

Once the nozzle mass flow rate and exit velocity have been established, for choked or nonchoked flow, the thrust is calculated. The thrust produced by the engine is based on the mass flow entering and exiting the engine, as well as the pressure difference between the nozzle and ambient air. The inlet mass flow rate, which is equivalent to the fan mass flow rate, is already known. The inlet velocity, however, must be calculated using (34). Using the inlet velocity, the total engine thrust is represented by (35) for a convergent nozzle:
(34)vinlet=MaircraftkambientRambientTambient,(35)Thrust=(m˙outvout-m˙invin)+Ain(Pout-Pambient).


## 3.6. Shaft

The HP shaft connects the HP Turbine and the HP Compressor. Power from the HP turbine is transferred by the HP shaft to drive the HP compressor. The HP turbine work signal represents a positive load and the HP compressor work signal represents a negative load. The HP shaft speed is the only calculation performed within the shaft model and is represented by (36). Additional auxiliary loads may be added to the shaft (not shown):
(36)NHP_shaft=30π∫W˙C+W˙HPTJHP,shaftωHP,shaft.The LP shaft connects the LP Turbine and the fan. Power from the LP turbine is transferred by the LP shaft to drive the fan. The LP turbine work signal represents a positive load and the fan work signal represents a negative load. The LP shaft speed is the only calculation performed within the model and is represented by (37). The LP shaft has additional auxiliary loads, pumps, and generator, which are inputs to the engine model:
(37)NLP,shaft=30π∫W˙Fan+W˙LPT+W˙Pumps+W˙GeneratorJLP,shaftωLP,shaft.


## 4. Simulation Results

Two sets of simulations were performed using the engine model. First, a fuel perturbation was performed to demonstrate and isolate the transient response behavior of the turbofan engine. Second, a full mission was simulated to compare the simulation times of previous approaches.


### 4.1. Engine Fuel Perturbation

The turbofan engine model was simulated at an altitude of 20,000 ft (6,096 m) and a speed of 0.6 Mach number. The fuel flow rate was stepped from 4.41 to 4.8 lbm/s (2 to 2.2 kg/s). Figure3 presents the step in fuel flow at time equal to 1 second on the top graph with the thrust response in the middle and shaft speeds at the bottom. There are essentially three time constants affecting the response of the engine: fuel step approximately 0 seconds, temperature and pressure approximately 0.05 seconds, and shaft speed approximately 10 seconds. The step increase in fuel flow increases the thrust output of the engine from 15.6 to 18.3 klbf. There is a rapid increase (first 0.05 second) in thrust initially following the step up in fuel, which is a result in the rapid temperature increase in the components downstream of the combustor as shown in the HP and LP turbine temperatures presented in Figure 4. In addition, the pressures throughout the engine also experience a rapid increase following the step up in fuel as shown for the turbines in Figure 4 and the fan and compressor in Figure 5. The increase in pressure also contributes to the increase in thrust.Figure 3
Engine fuel flow input, specific fuel consumption, thrust and shaft speeds responses to a step in fuel flow.Figure 4
Pressures, temperatures, and mass flows for the HP turbine and LP turbine including mass flow of the 2nd stream bypass duct.Figure 5
Pressures, temperatures, and mass flows for the fan and compressor.Initially there is sharp increase in SFC in Figure3 due to the sudden (step change) increase in fuel and delayed increase in thrust over a 0.05 seconds. The SFC decreases sharply during the first 0.05 seconds following the perturbation as the thrust increases. On a larger time scale of 10 seconds, there is gradual increase in thrust due to the response in shaft speeds, which in turn increases the flow rate through the engine and continues to decrease the SFC.Figure4 presents the turbine pressures, temperatures, and mass flows as well as bypass mass flow. The effects of the pressure and temperature on thrust and SFC were previously discussed. The total temperature of the turbines in Figure 4 increases rapidly due to the increase in fuel and temperature in the combustor but begin to decline as the air mass flow rate through the combustor and turbines increase as the shaft speeds increase. The flow rates in Figure 4 present interesting behavior to the fuel perturbation. The turbine mass flow rates decrease (very slightly) and the bypass mass flow rate decreases noticeably before returning to the previous initial value. The same trends may be seen in the fan and compressor mass flows as for the turbine mass flows in Figure 5. The mass flow rates initially drop (very slightly) in the fan and compressor with the step increase in fuel. The initial drop in mass flow rate is due to the pressures increasing as shown in Figure 5 for the fan and compressor. The drop in mass flow reverses with the increase in shaft speed shown in Figure 3. The fan and compressor maps in Figure 6 support the mass flow/pressure interactions. Figure 6 illustrates the initial decrease in normalized mass flow with increased pressure ratio (moving up and to the left) before the normalized mass flow increases with shaft speed (moving right).Figure 6
Fan and compressor map for a step in fuel flow to the engine.The bypass mass flow rate noticeably decreases in Figure4 as the flow rates rebalance between the HP compressor and bypass in the two-stream engine. Some of the decline in bypass flow is due to the initial drop in fan mass flow. Additional decline in mass flow is due to the quicker response of the compressor and HP turbine to the fuel perturbation. The HP turbine is effected initially by the fuel increase since its first downstream of the combustor. In turn, the compressor speed is increased since it is connected to the HP shaft with the HP turbine. The compressor speed increases the compressor mass flow more rapidly before the fan mass flow increases resulting in less flow in the bypass. Once the fan speed catches up with the compressor speed and fan mass flow increases, the bypass flow returns to its initial value ending with a slight increase in overall flow.If the plenum volumes were assumed to be steady-state, the initial mass flow responses to the pressure perturbations would not be captured. The engine model was able to capture both the plenum volume dynamics and the shaft dynamics.


### 4.2. Computational Time

The engine model is assessed for computational efficiency in this section. First, the engine model is simulated standalone for a specified mission presented in Figure7. The mission is 7,700 seconds with varying altitudes and Mach number. During a mission, the engine model has 11 inputs. The inputs are(1)
fuel flow,(2)
fuel temperature,(3)
air inlet temperature,(4)
air inlet pressure,(5)
air inlet composition,(6)
bleed flow demand from the fan,(7)
bleed flow demand from the compressor,(8)
shaft power extraction from the LP shaft,(9)
shaft power extraction from the HP shaft,(10)
pressure drop in the bypass heat exchangers,(11)
heat flux in the bypass heat exchangers.Throughout the mission all of the inputs are continuously changing. All of the inputs are treated as a disturbance except for the fuel flow. The inputs described as disturbances are provided by lookup tables utilizing previous data from a previous full mission T2T simulation. The fuel flow is controlled to achieve a demanded thrust from the engine via a feedback controller. The thrust demand is also in the form of a lookup table utilizing previous data. This approach isolates the computational efficiency of the engine and its controller. The data for the inputs contained in the lookup tables are from previous simulations of a T2T model. The engine model from previous work is also used to simulate the same mission with the same inputs for comparison [3]. The new engine model took 15 seconds to complete the entire mission simulation. The engine model from previous work took 29,304 seconds to simulate the same mission with less detail. The new turbofan engine model with the plenum volume dynamics and approximated performance maps near the surge line reduced the simulation time by 99.95%, while providing more detail and capturing high frequency dynamics.Figure 7
Altitude and Mach # Profile for a 2.1 hour mission.The engine model was then integrated with a T2T model. The T2T model is exercised to simulate the same predetermined mission presented in Figure7. Throughout the mission demands on shaft power extractions, bleed flows, fuel flow, inlet conditions and the bypass heat exchanger are continuously changing, but are not implemented via lookup tables. The T2T model contains dynamic models for the air vehicle, thermal management system, electrical system, and fuel system. The T2T model with the old engine model without the techniques presented completed the mission simulation in 72,650 seconds. The T2T model with the new turbofan engine model completed the mission simulation in 186 seconds. The new engine model reduced the simulation time for the T2T model by 99.75%. All simulations were carried out on the same hardware. The hardware used was an Intel Core 2 Duo Processor E8600 with 4 G of RAM. The increase in computational efficiency is significantly opening the opportunity for design optimization.


## 4.1. Engine Fuel Perturbation

The turbofan engine model was simulated at an altitude of 20,000 ft (6,096 m) and a speed of 0.6 Mach number. The fuel flow rate was stepped from 4.41 to 4.8 lbm/s (2 to 2.2 kg/s). Figure3 presents the step in fuel flow at time equal to 1 second on the top graph with the thrust response in the middle and shaft speeds at the bottom. There are essentially three time constants affecting the response of the engine: fuel step approximately 0 seconds, temperature and pressure approximately 0.05 seconds, and shaft speed approximately 10 seconds. The step increase in fuel flow increases the thrust output of the engine from 15.6 to 18.3 klbf. There is a rapid increase (first 0.05 second) in thrust initially following the step up in fuel, which is a result in the rapid temperature increase in the components downstream of the combustor as shown in the HP and LP turbine temperatures presented in Figure 4. In addition, the pressures throughout the engine also experience a rapid increase following the step up in fuel as shown for the turbines in Figure 4 and the fan and compressor in Figure 5. The increase in pressure also contributes to the increase in thrust.Figure 3
Engine fuel flow input, specific fuel consumption, thrust and shaft speeds responses to a step in fuel flow.Figure 4
Pressures, temperatures, and mass flows for the HP turbine and LP turbine including mass flow of the 2nd stream bypass duct.Figure 5
Pressures, temperatures, and mass flows for the fan and compressor.Initially there is sharp increase in SFC in Figure3 due to the sudden (step change) increase in fuel and delayed increase in thrust over a 0.05 seconds. The SFC decreases sharply during the first 0.05 seconds following the perturbation as the thrust increases. On a larger time scale of 10 seconds, there is gradual increase in thrust due to the response in shaft speeds, which in turn increases the flow rate through the engine and continues to decrease the SFC.Figure4 presents the turbine pressures, temperatures, and mass flows as well as bypass mass flow. The effects of the pressure and temperature on thrust and SFC were previously discussed. The total temperature of the turbines in Figure 4 increases rapidly due to the increase in fuel and temperature in the combustor but begin to decline as the air mass flow rate through the combustor and turbines increase as the shaft speeds increase. The flow rates in Figure 4 present interesting behavior to the fuel perturbation. The turbine mass flow rates decrease (very slightly) and the bypass mass flow rate decreases noticeably before returning to the previous initial value. The same trends may be seen in the fan and compressor mass flows as for the turbine mass flows in Figure 5. The mass flow rates initially drop (very slightly) in the fan and compressor with the step increase in fuel. The initial drop in mass flow rate is due to the pressures increasing as shown in Figure 5 for the fan and compressor. The drop in mass flow reverses with the increase in shaft speed shown in Figure 3. The fan and compressor maps in Figure 6 support the mass flow/pressure interactions. Figure 6 illustrates the initial decrease in normalized mass flow with increased pressure ratio (moving up and to the left) before the normalized mass flow increases with shaft speed (moving right).Figure 6
Fan and compressor map for a step in fuel flow to the engine.The bypass mass flow rate noticeably decreases in Figure4 as the flow rates rebalance between the HP compressor and bypass in the two-stream engine. Some of the decline in bypass flow is due to the initial drop in fan mass flow. Additional decline in mass flow is due to the quicker response of the compressor and HP turbine to the fuel perturbation. The HP turbine is effected initially by the fuel increase since its first downstream of the combustor. In turn, the compressor speed is increased since it is connected to the HP shaft with the HP turbine. The compressor speed increases the compressor mass flow more rapidly before the fan mass flow increases resulting in less flow in the bypass. Once the fan speed catches up with the compressor speed and fan mass flow increases, the bypass flow returns to its initial value ending with a slight increase in overall flow.If the plenum volumes were assumed to be steady-state, the initial mass flow responses to the pressure perturbations would not be captured. The engine model was able to capture both the plenum volume dynamics and the shaft dynamics.


## 4.2. Computational Time

The engine model is assessed for computational efficiency in this section. First, the engine model is simulated standalone for a specified mission presented in Figure7. The mission is 7,700 seconds with varying altitudes and Mach number. During a mission, the engine model has 11 inputs. The inputs are(1)
fuel flow,(2)
fuel temperature,(3)
air inlet temperature,(4)
air inlet pressure,(5)
air inlet composition,(6)
bleed flow demand from the fan,(7)
bleed flow demand from the compressor,(8)
shaft power extraction from the LP shaft,(9)
shaft power extraction from the HP shaft,(10)
pressure drop in the bypass heat exchangers,(11)
heat flux in the bypass heat exchangers.Throughout the mission all of the inputs are continuously changing. All of the inputs are treated as a disturbance except for the fuel flow. The inputs described as disturbances are provided by lookup tables utilizing previous data from a previous full mission T2T simulation. The fuel flow is controlled to achieve a demanded thrust from the engine via a feedback controller. The thrust demand is also in the form of a lookup table utilizing previous data. This approach isolates the computational efficiency of the engine and its controller. The data for the inputs contained in the lookup tables are from previous simulations of a T2T model. The engine model from previous work is also used to simulate the same mission with the same inputs for comparison [3]. The new engine model took 15 seconds to complete the entire mission simulation. The engine model from previous work took 29,304 seconds to simulate the same mission with less detail. The new turbofan engine model with the plenum volume dynamics and approximated performance maps near the surge line reduced the simulation time by 99.95%, while providing more detail and capturing high frequency dynamics.Figure 7
Altitude and Mach # Profile for a 2.1 hour mission.The engine model was then integrated with a T2T model. The T2T model is exercised to simulate the same predetermined mission presented in Figure7. Throughout the mission demands on shaft power extractions, bleed flows, fuel flow, inlet conditions and the bypass heat exchanger are continuously changing, but are not implemented via lookup tables. The T2T model contains dynamic models for the air vehicle, thermal management system, electrical system, and fuel system. The T2T model with the old engine model without the techniques presented completed the mission simulation in 72,650 seconds. The T2T model with the new turbofan engine model completed the mission simulation in 186 seconds. The new engine model reduced the simulation time for the T2T model by 99.75%. All simulations were carried out on the same hardware. The hardware used was an Intel Core 2 Duo Processor E8600 with 4 G of RAM. The increase in computational efficiency is significantly opening the opportunity for design optimization.


## 5. Discussion of Results

The results for the turbofan engine model built in Matlab-Simulink were presented. The engine model is capable of capturing the shaft dynamics, thermal dynamics, and the unsteady mass flow/plenum volume dynamics. By capturing the unsteady mass flow/plenum volume dynamics, two things were achieved. First, unsteady mass flow and pressure are simulated providing more detail in engine behavior. These perturbations propagate throughout the engine and integrated subsystems that utilize compressor bleed air. Second, the engine model reduces the simulation time of a vehicle level T2T model tool making it feasible to simulate thousands of simulations needed for a complex optimization routine of a large system. The ability to perform simulations in a timely manner for system level optimization including transient behavior and control architectures provides advanced capability to design highly integrated dynamic systems.The modeling techniques shown here introduce limitations. For example, simulations with the compressors operating near the surge line of the map should be aborted or considered with caution due to the approximations made in the maps to remove iteration loops. The unsteady mass flow/plenum volume dynamics captured by the model are approximated by using large control volumes within the engine components. For more detailed results, the control volumes need to be discreetly divided into smaller control volumes creating a one-dimensional distribution of mass flow and pressure for each of the engine components. This would increase the complexity and simulation time depending on the resolution of the one-dimensional grid of control volumes.


## 6. Conclusion

A dynamic, turbofan engine model was developed in the modeling and simulation environment of Matlab-Simulink. The model has been built without the aid of proprietary data, allowing the tool to be made available to multiple design and research groups. Special attention was also paid to the capturing of dynamic behavior. These dynamics not only increase the physics being captured, but even more importantly, they reduce the occurrence of algebraic constraints, leading to increased simulation speed. The new techniques of approximating the performance maps near the surge line without the use ofR-lines and implementing the plenum volume dynamics have proven very efficient, as the current engine model completes the 7700 second mission in just 15 seconds. When the engine model was integrated with the full T2T model significant gains in computational efficiency were maintained. The new modeling techniques of incorporating plenum volume dynamics and approximating the compressor maps near the surge line result in a more efficient T2T model. The T2T model with the new engine model will be capable of conducting trade studies and vehicle level design optimization.



---

*Source: 283479-2014-10-12.xml*