Abstract
This work presents fully coupled computational fluid dynamics (CFD) simulations and thermodynamic cycle analyses of a small-scale turbojet engine at several conditions along the equilibrium running line. The CFD simulations use a single mesh for the entire engine, from the intake to the exhaust, allowing information to travel in all directions. The CFD simulations are performed along the equilibrium running line by using the iterative Secant method to compute the fuel flow rate required to match the compressor and turbine power. The freestream pressure and temperature and shaft angular speed are the only inputs needed for the CFD simulations. To evaluate the consistency of the CFD results with thermodynamic cycle results, outputs from the CFD simulations are prescribed as inputs to the cycle model. This approach enables on-design and off-design cycle calculations to be performed without requiring turbomachinery performance maps. In contrast, traditional off-design cycle analyses require either scaling, calculating, or measuring compressor and turbine maps with boundary condition assumptions. In addition, the CFD simulations and the cycle analyses are compared with measurements of the turbojet engine. The CFD simulations, thermodynamic cycle analyses, and measurements agree in terms of total temperature and pressure at the diffuser–combustor interface, air and fuel mass flow rate, equivalence ratio, and thrust. The developed methods to perform CFD simulations from the intake to the exhaust of the turbojet engine are expected to be useful for guiding the design and development of future small-scale gas turbine engines.
1 Introduction
The design and development of modern gas turbine engines can benefit from emerging system-level computational fluid dynamics (CFD) simulations. System-level CFD simulations provide several benefits relative to traditional component-level simulations. First, the inlet and outlet boundaries for each component are coupled within the engine. Therefore, system-level CFD simulations remove assumptions regarding the boundary conditions for each component and enable more accurate computations of the pressure, temperature, and velocity profiles at component interfaces. Second, more realistic system-level compressor, combustor, and turbine efficiencies can be computed from CFD simulations using isentropic or polytropic approximations. The efficiencies are valuable for design optimization of individual components and provide feedback to reduced-order thermodynamic cycle models. Third, system-level approaches can lead to simulations of full engine flight scenarios with higher fidelity than traditional thermodynamic cycle analyses that utilize a single altitude compressor and turbine map with corrected flow and shaft angular speed to account for altitude variations. Last, system-level CFD simulations enable the extraction of results from any three-dimensional location within the engine flow path. The three-dimensional information provides insights into nonuniformities at component interfaces, supports precise comparisons with experimental data, and complements limited measurements of gas turbine engine systems.
Full engine CFD simulations are challenging for several reasons and only a few studies have been reported [1–9]. First, gas turbine engines often do not exhibit global periodicity but rather the periodicity varies from one component to another. Consequently, simulations are required to handle varying periodicity or pitch angle between components. Second, the turbomachinery geometry is altered during operation due to centrifugal loads, pressure forces, and heat transfer, yielding to radial growth and blade untwist. These complexities are not considered in this work and represent potential future research opportunities in the context of full engine simulations. Third, there are regions within the engine flow path that experience high Reynolds and Mach numbers. Therefore, compressible flow solvers are required for the compressor and turbine even though the combustor is at low Mach number. Fourth, the species transport equations must account for compression and expansion in the compressor rotor and turbine rotor, respectively. Simulating an entire engine requires a relatively large mesh cell count; therefore, global chemistry and flamelet approaches are the most cost-effective options for modeling species transport. The former has to be tuned to attain realistic combustion efficiencies whereas the latter is typically available for constant pressure and low Mach number flows.
Simulations of an entire gas turbine engine have been reported previously. Turner and coworkers [1–4] performed several multifidelity simulations of a commercial turbofan engine (GE90-94B) at takeoff conditions. Their approach used two dedicated solvers, average-passage multistage turbomachinery flow field analysis code and national combustion code, for simulating the turbomachinery and combustor, respectively. The solvers communicated circumferentially averaged quantities at the interfaces using a mixing plane model. The CFD component results and a one-dimensional thermodynamic analysis were used to generate performance maps. The numerical propulsion system simulation used these maps to set boundary conditions such as mass flow rates and rotational speeds in the CFD components for the coupled full engine simulation. This numerical strategy requires several codes to operate in tandem. Pitsch and coworkers [5] performed simulations of a modified commercial turbofan engine with an axial compressor and turbine. Their approach scaled the blades to model the engine with fixed periodicity. The full engine simulation was initiated from individual component simulation results. The turbomachinery components were simulated using an unsteady Reynolds-averaged Navier–Stokes (RANS) solver (SUmb), and the combustor was simulated using a large eddy simulations solver. The solvers communicated through an intermediary numerical tool. The NUMECA group [6,7] simulated a microgas turbine engine (KJ66) operating at one condition using multiple sectors per component with a RANS turbulence model, flamelet generated manifold combustion model, and either steady-state or nonlinear harmonic method, respectively. Although there was no comparison with experimental data, the results are encouraging because the flow field (of 19.2 × 106 cells) was computed with a single CFD solver in a timely manner (i.e., 15,000 and 42,000 CPU hours for the steady-state and nonlinear harmonic simulation) with reasonable computational resources (i.e., 288 cores). Utilizing conventional flamelet models such as flamelet generated manifold implies that the conditions at the combustor inlet (station 3.1) are known a priori. Ideally, there would be only three inputs: the freestream temperature, the freestream pressure, and the fuel flow rate or shaft angular speed. Another caveat of their approach is that fuel flow rate and shaft angular speed had to be specified a priori without matching the compressor and turbine power. Based on this review of previous full engine CFD simulations, there are several opportunities to further develop this capability.
We have previously performed CFD simulations of a small-scale gas turbine engine [8,9]. Initially, we modeled and simulated the turbojet engine from inlet to exit using the eddy break-up combustion model [8]. The parameters of the combustion model had to be tuned in an ad hoc manner to produce realistic combustion efficiencies and combustor exit temperatures at design and off-design conditions. Subsequently, we developed a penta-dimensional compressible flamelet/progress variable (FPV) approach [9,10] to account for compressibility and combustor operating conditions which are not known a priori. We conducted fully coupled full engine CFD simulations with prescribed freestream conditions, fuel flow rate, and shaft angular speed. Our full engine CFD simulations compared well with measurements from minimum to maximum power in terms of the diffuser-combustor thermodynamic state, air mass flow rate, and thrust. A logical progression is to perform the full engine CFD simulation without prescribing either the fuel flow rate or the shaft angular speed. This is critical to advance the state of the art of full engine CFD simulations from post-testing to pretesting capability when the relationship between fuel flow rate and shaft angular speed is not known a priori.
Motivated by the previous discussion, this investigation has two purposes. First, we develop a method for computing the equilibrium running line [11] by matching the power of the compressor and turbine in full engine CFD simulations such that the only inputs are the freestream static pressure and temperature and shaft angular speed. Second, we show that the full engine CFD simulation results compare well with one-dimensional thermodynamic analyses and experimental data. The specific objectives are: (1) to demonstrate that compressor-turbine power matching of the full engine CFD simulations can be obtained via the Secant root-finding method; (2) to qualitatively illustrate that the diffuser–combustor state and combustor–turbine state exhibit heterogeneous flow which supports the need for system-level design approaches; and (3) to show that the measurements, CFD results, and thermodynamic cycle analyses results agree in terms of total temperature and pressure at the diffusor–combustor interface, air and fuel mass flow rate, equivalence ratio, and thrust as function of shaft angular speed.
2 Computational Fluid Dynamics Simulations
The CFD simulations of a small-scale turbojet engine are performed using a modified version of ansyscfx [12]. This section describes details related to the engine geometry, computational domain, governing equations, combustion model, and method for matching the compressor and turbine power for calculating the equilibrium running line.
2.1 Turbojet Engine Geometry.
Figure 1 shows the geometry of the small-scale turbojet engine. The engine contains eight components: a converging intake (from station 1 to 2), radial compressor (from station 2 to 3), diffuser (from station 3 to 3.1), combustor (from station 3.1 to 4), turbine stator (from station 4 to 4.1), turbine rotor (from station 4.1 to 5), and converging nozzle (from station 5 to 8). A large plenum is included downstream of the nozzle to capture the exhaust flow and reduce boundary-induced disturbances. The station numbering is based on the Society of Automotive Engineers (SAE) Aerospace Recommended Practice (ARP) 755 [13].
2.2 Computational Domain, Meshes, and Interface Closure Models.
The computational domain is an assembly of eight individual computational domains, each corresponding to an engine component and the exhaust. Table 1 illustrates the modeling approach used for the engine simulations. Structured meshes with conformal periodicity are attained for each of the components, except for the diffuser and combustor. The diffuser and combustor use a structured mesh with nonconformal periodicity and a hexcore mesh with conformal periodicity, respectively. Moreover, all components exhibit ten prisms (normal to wall) for boundary layers. The wall y+ varies throughout the engine and for each design point simulation. The model [12] ensures that the velocity adjacent to the wall is always within the log-wall regime, regardless of mesh size. The tip clearance for the components with blades is assumed to be zero because it simplifies mesh generation without compromising the motivation and objectives stated in the Introduction. Future work should include a tip clearance model for the centrifugal rotor and turbine rotor. The volume average spatial resolution of the mesh inside the engine is approximately 0.57 mm, and the entire engine mesh is 16.5 × 106 cells. The number of sectors for each component is chosen to reduce the pitch angle ratio between components. For example, the compressor rotor periodicity is 360 deg/7, which provides a pitch angle ratio of 1 with the intake and 0.93 with the diffuser. Table 1 lists the remaining periodicities, cell count, number of sectors, and pitch ratios. The simulations use a mixing plane closure model [12] for interfaces between a rotating and stationary component containing vanes (e.g., compressor-diffuser interface and turbine stator-turbine rotor interface). The mixing plane closure model communicates steady-state circumferential-averaged values across the nonconformal mesh interface between two components. The simulations use a frozen rotor closure model [12,14] for interfaces when at least one of the adjacent components has no vanes (e.g., intake-compressor). The frozen rotor closure model communicates interpolated values across the nonconformal mesh interface between two components by changing the frame of reference while maintaining the relative position of the components. To account for differences in pitch ratio between the components, the flow passing through the interface is scaled according to the pitch ratio. The Frozen Rotor model is robust and utilizes less computer resources than the other frame change models.
Station | Components | Mesh type | Cell count | Periodicity | Sectors | Pitch ratio | Interface |
---|---|---|---|---|---|---|---|
1–2 | Intake | Hexahedral mesh with conformal periodicity | 0.10 | 360 deg/7 | 1 | ||
1.00 | Frozen rotor | ||||||
2–3 | Compressor | Hexahedral mesh with conformal periodicity | 0.73 | 360 deg/7 | 1 | ||
0.93 | Mixing plane | ||||||
3–3.1 | Diffuser | Tetrahedral mesh with nonconformal periodicity | 1.21 | 24 deg | 2 | ||
1.25 | Frozen rotor | ||||||
3.1–4 | Combustor | Hexcore mesh with nonconformal periodicity | 8.21 | 60 deg | 1 | ||
0.80 | Frozen rotor | ||||||
4–4.1 | Turbine Stator | Hexahedral mesh with conformal periodicity | 1.40 | 24 deg | 2 | ||
0.98 | Mixing plane | ||||||
4.1–5 | Turbine Rotor | Hexahedral mesh with conformal periodicity | 3.16 | 360 deg/23 | 3 | ||
0.96 | Frozen rotor | ||||||
5–8 | Nozzle | Hexahedral mesh with conformal periodicity | 0.76 | 45 deg | 1 | ||
1 | Frozen rotor | ||||||
8–9 | Exhaust | Hexahedral mesh with conformal periodicity | 0.98 | 45 deg | 1 |
Station | Components | Mesh type | Cell count | Periodicity | Sectors | Pitch ratio | Interface |
---|---|---|---|---|---|---|---|
1–2 | Intake | Hexahedral mesh with conformal periodicity | 0.10 | 360 deg/7 | 1 | ||
1.00 | Frozen rotor | ||||||
2–3 | Compressor | Hexahedral mesh with conformal periodicity | 0.73 | 360 deg/7 | 1 | ||
0.93 | Mixing plane | ||||||
3–3.1 | Diffuser | Tetrahedral mesh with nonconformal periodicity | 1.21 | 24 deg | 2 | ||
1.25 | Frozen rotor | ||||||
3.1–4 | Combustor | Hexcore mesh with nonconformal periodicity | 8.21 | 60 deg | 1 | ||
0.80 | Frozen rotor | ||||||
4–4.1 | Turbine Stator | Hexahedral mesh with conformal periodicity | 1.40 | 24 deg | 2 | ||
0.98 | Mixing plane | ||||||
4.1–5 | Turbine Rotor | Hexahedral mesh with conformal periodicity | 3.16 | 360 deg/23 | 3 | ||
0.96 | Frozen rotor | ||||||
5–8 | Nozzle | Hexahedral mesh with conformal periodicity | 0.76 | 45 deg | 1 | ||
1 | Frozen rotor | ||||||
8–9 | Exhaust | Hexahedral mesh with conformal periodicity | 0.98 | 45 deg | 1 |
For each component, the table lists the mesh type, cell count (in millions), periodicity, number of sectors, pitch ratio, and interface closure model.
Equations | Φ | ΓΦ | SΦ |
---|---|---|---|
Continuity | |||
Momentum | |||
Rothalpy | |||
Turbulent kinetic energy | |||
Eddy dissipation rate | |||
Species mass fraction | |||
Mixture fraction | |||
Mixture fraction variance | |||
Progress variable |
Equations | Φ | ΓΦ | SΦ |
---|---|---|---|
Continuity | |||
Momentum | |||
Rothalpy | |||
Turbulent kinetic energy | |||
Eddy dissipation rate | |||
Species mass fraction | |||
Mixture fraction | |||
Mixture fraction variance | |||
Progress variable |
2.3 Governing Equations.
Depending on the value of , this equation represents the continuity, momentum, rothalpy,3 turbulence kinetic energy, eddy dissipation rate, species mass fractions, mixture fraction, mixture fraction variance, and progress variable transport equations as shown in Table 2. This table also describes the transport coefficient () and the source term () for each equation. Density is dependent on local pressure and temperature through the ideal gas equation of state. The Reynolds stresses in the momentum equations are related to the mean velocity gradients through the Boussinesq hypothesis [15]. The turbulence kinetic energy and eddy dissipation rate Prandtl numbers are σk = 1.0 and σkε = 1.3, respectively. The other turbulence model constants are and .
The CFD simulations of the engine involve a hybrid of species transport equations () and low-dimensional manifold (Z, Z”2, Y) equations at different locations throughout the engine. The species transport equations are solved in the intake, compressor, diffuser, turbine, nozzle, and exhaust. The low-dimensional manifold equations are solved in the combustor. The rothalpy§ and species transport equations are not used in the combustor domain as indicated by their linearized source terms. The linearization () of the source terms and with large negative coefficients has the effect of fixing values in the combustor domain, resulting in and [16]. Note that n + 1 and n indicate the current and previous outer iteration, respectively. The manifold equations are solved in the combustor domain for interpolating the probability density function (PDF) table and obtaining temperature and species, which are then used to compute the density, thermodynamic properties, viscosity, and thermal conductivity. The Schmidt numbers for the species transport and low-dimensional manifold equations are 0.9. Molecular and turbulent Lewis number are assumed to be unity for all species.
A conventional distilled Jet-A fuel [17–19] is injected through point sources. The chemical kinetic model has been validated from 403 K to 720 K, from 1 to 14 atm, and from ϕ = 0.65 to 1.50. Scalable wall functions [12] are used to model the boundary layers near solid surfaces. The combustor cooling is modeled using source and sink terms [9] on the inner and outer liner surface, respectively, to avoid the computational cost of meshing through the cooling holes. The diffusion terms are discretized using central differencing, and the advection terms are discretized using second-order upwind scheme. The governing equations are solved with the pseudo-transient coupled pressure-based solver. The simulations were executed using an SGI ICE X system (Silicon Graphics, Inc., Milpitas, CA) with 16 Xeon Phi nodes (total of 576 cores). Each operating condition requires approximately 18,000 CPU hours.
2.4 Multipoint Compressible Flamelet/Progress Variable Combustion Model.
Flamelet models have been used in many canonical flame studies and combustor applications [20,21]. However, traditional flamelet models are only applicable to a single operating condition with a prescribed combustor inlet temperature ( and pressure (. This is problematic when simulating the engine across a range of engine operating conditions because the combustor inlet quantities at station 3.1 vary with shaft angular speed ().4 One solution to this problem involves generating several three-dimensional flamelet tables at different and conditions that encompass the engine operation from minimum to maximum power as schematically shown in Fig. 2. The range of and can be obtained from many sources including estimates based on thermodynamic cycle analysis. Each table [(Z, Z”2, Y)] is generated using previously developed tools [10], and each table may have a different number of discretization points. A trilinear interpolation among the combustion manifold values (Z, Z”2, Y) is performed to obtain the local properties and composition. If and are within the range of the computed tables, then an additional bilinear interpolation is performed involving four tables in the and space (i.e., effective pentalinear). If either or is outside of the range of the computed tables, then only an additional linear interpolation (i.e., effective tetralinear) is performed involving two adjacent tables. If both and are outside the range of computed tables, then the nearest table approximation is employed (i.e., effective trilinear). In practice, the flamelet tables are computed over a broad range of pressures and temperatures that encompass the entire engine operating space and allow the pentalinear interpolation approach to be used.
2.5 Method for Matching Compressor and Turbine Power.
The compressor–turbine power matching () is obtained by integrating the enthalpy fluxes over the (stationary) interface areas between the appropriate turbomachinery components. The shaft mechanical efficiency () is assumed to be unity. The subscript “m” indicates the iteration number. In this study, the initial value of is set to that of the experimental value and the second guess is set to . The compressor–turbine matching procedure implicitly ensures that the simulated conditions occur on the engine equilibrium running line. The convergence criterion for the Secant method is when the compressor and turbine power match to within less than 1%.
3 Thermodynamic Cycle Analyses
This section describes the on-design and off-design thermodynamic cycle analyses and the method for comparing the CFD results with the cycle results to evaluate and verify the CFD simulations. Traditional thermodynamic cycle analyses of gas turbine engines involve performing system-level on-design calculations at the selected design condition followed by off-design calculations at other conditions of interest. The on-design cycle analysis is typically performed at either sea-level static, takeoff, or cruise conditions. The cycle inputs are described as component-level (uncoupled) inputs because individual component performance values are used in the cycle model. An on-design cycle analysis for a single-spool engine requires twelve inputs accounting for the freestream pressure and temperature; air mass flow rate; shaft angular speed or thrust; component pressure ratios or pressure losses; combustor exit temperature; and compressor, combustor, and turbine efficiencies. The on-design cycle analysis results in the pressures, temperatures, and cross-sectional areas at each component interface for the engine to satisfy specific performance criteria (e.g., thrust and specific fuel consumption) at the selected on-design condition. In this work, the thermodynamic cycle analyses are performed using pyCycle [22,23].
Traditional off-design thermodynamic cycle analyses require additional component-level (uncoupled) compressor and turbine maps to determine performance along the equilibrium running line of the engine. These maps are obtained from either experiments, CFD simulations, or a combination of both. Compressor and turbine maps assume that the exit pressure distribution and inlet temperature distribution, respectively, are similar for the components regardless of whether the measurements or simulations are performed in a coupled or uncoupled manner. The turbomachinery maps are nonorthogonal and four-dimensional containing the relationship between corrected mass flow rate, pressure ratio, corrected shaft speed, and isentropic efficiency. The corrected quantities are used to account for different altitudes in the cycle analysis procedure. Since the equilibrium running line is not known a priori, many conditions need to be represented in the maps spanning from the surge to the choke lines. Consequently, many of the regions in the maps are not accessed by the off-design cycle analyses.
The mass flow rate, pressure ratio, component efficiency, and shaft angular speed from the on-design cycle analysis are utilized to scale the compressor and turbine maps within the engine operability (e.g., operating pressure ratio and mass flow rate ranges). If compressor and turbine maps are generated based on a specific engine design, it is reasonable to expect more accurate results than scaling existing maps. Figure 3 illustrates a typical single-spool turbojet off-design cycle model. The off-design cycle analysis proceeds as follows: The shaft angular speed, air mass flow rate (), and equivalence ratio () are varied to lower the residual equations () of conservation of energy (i.e., ) and mass (i.e., ) as well as the target thrust (i.e., ). For each solver iteration, a value of and are corrected to the given altitude pressure and temperature to access the scaled compressor map pressure ratios and efficiencies. This map outputs the pressure ratio () and efficiency () required to compute the compressor power () while accounting for irreversible losses through the output enthalpies and total temperatures. For the combustor, the equivalence ratio () and combustor efficiency () are inputs and the fuel mass flow rate () is an output. The exit temperature () is calculated by using NASA chemical equilibrium and applications (CEA) [24,25]. Then, and are corrected to the altitude pressure and temperature conditions to access the scaled turbine map. This map outputs the and efficiency () required to compute the turbine power () while accounting for irreversible losses through the output enthalpies and total temperatures. The analyses converge when the residual equations () are below a prescribed convergence threshold.
To evaluate the consistency of the CFD results with thermodynamic cycle results and to demonstrate that the former can provide system-level (coupled) inputs to the latter, we prescribe outputs from the CFD simulations as inputs to the thermodynamic cycle model for a range of conditions. If the full engine CFD simulations were not at the equilibrium running line, then the comparison with steady-state thermodynamic cycle models would not be appropriate because the CFD simulations would not be matching the compressor and turbine power. Figure 4 illustrates the single-spool turbojet off-design cycle model with inputs provided from CFD simulations. This off-design cycle analysis schematic is similar to a traditional on-design cycle analysis. The main differences between the off-design nontraditional cycle schematic (cf. Fig. 4) with the respect to the traditional off-design cycle schematic (cf. Fig. 3) include the following: (1) the compressor and turbine maps are no longer needed because the full engine CFD simulation results along the equilibrium running line contain thermodynamic state information; (2) the compressor and turbine components require efficiency [] and inputs; (3) the component areas () are replaced with Mach numbers [()] inputs, except for the nozzle which requires a velocity coefficient (); (4) air mass flow rate [] is an input to the intake component; (5) the combustor total temperature ] is an input; and (6) the shaft angular speed is prescribed. The conservation of energy is enforced by matching the compressor and turbine power [. The combustor exit temperature is prescribed to match the CFD simulations [. Both and are varied to satisfy these residual equations (). In contrast with the traditional off-design cycle model, thrust () and nozzle throat area () are now outputs of the model. The velocity coefficient () and gross thrust coefficient () are set to unity leading to ideal predictions of and . Additional inputs to the individual components, such as efficiencies and pressure ratios, are taken from the CFD simulations. The calculation time of the thermodynamic cycle analysis is in the order of seconds.
4 Experimental Methods
The small-scale turbojet engine (i.e., JetCat P400) is studied experimentally over a range of operating conditions from minimum to maximum power at freestream conditions corresponding to an altitude of 230 m, ISA + 9K, and Mach 0.0. Station 3.1 is instrumented for measuring total pressure and total temperature at a two separate locations. Total pressure is measured using Kiel-type probes. Total temperature is measured with sheathed K-type thermocouples with an accuracy of approximately ±2 K. Data is recorded at 10 Hz and shaft angular speed ( is held constant during data acquisition. The thrust () measurements are performed using a thrust stand. Air mass flow rate () is measured with a Pro-M 92 intrusive flow meter placed upstream of the engine inlet. Further details of the experimental setup and test facility are found elsewhere [26,27].
Two separate tests are conducted to measure several quantities of interest. During Test A, the total pressure () and total temperature () at the diffuser–combustor interface, thrust (), and shaft angular speed ( are measured as a function of the fuel mass flow rate (). During Test B, the air mass flow rate () and shaft angular speed are measured as a function of the fuel mass flow rate (). For Test A, the engine operation is conducted without an air flow meter, whereas, for Test B, an intrusive air flow meter is attached upstream of the engine intake to measure . It is worth mentioning that during the tests for a given , (Test B) > (Test A) and (Test B) < (Test A). The flow meter induces additional pressure losses that raises and reduces for Test B relative to Test A. It is not certain whether from Test B is equal to from Test A. In addition, note that includes the fuel for combustion and lubrication, which represents approximately 5% of the total fuel flow rate at maximum power. Nonetheless, the CFD results are compared with the experimental results for , , , and in Test A and in Test B.
5 Results and Discussion
5.1 CFD Simulation Results on the Engine Equilibrium Running Line.
Previous fully coupled turbojet engine CFD simulations [8,9] required prescribing the fuel mass flow () as a function of the shaft angular speed. The functional relationship was obtained from experimental data. The limitations of this approach are that the CFD simulations are only useful for post-testing evaluation and the compressor-turbine power matching are not satisfied (). These limitations are overcome in this work by using the iterative Secant method to compute the fuel flow rate required to match the compressor and turbine power () as described in Sec. 2.5. Matching the compressor and turbine power also ensures that the simulated conditions occur on the engine equilibrium running line.
Figure 5 shows the normalized difference between the compressor and turbine power as a function of fuel flow rate. The first iteration labeled as “1” are the results obtained previously [9], and the power matching is not satisfied with the prescribed fuel flow rate (). The difference is as high as 15% for = 80,000 rpm. Both the first and second iterations represent an initial guess for the fuel flow rate. The third and fourth iterations are obtained from the Secant method. A maximum of four CFD simulations per shaft angular speed are needed to converge to the equilibrium running line to within 1%. The initial guess for the fuel mass flow rate () overestimated the fuel requirement for cycle matching at 45,000, 60,000, 80,000, and 98,000 rpm, whereas it slightly underestimated the fuel requirement at 30,000 rpm. When , has to be decreased because the power extracted from the flow by the turbine exceeds the power consumed by the compressor. On the other hand, when has to be increased because the power extracted from the flow by the turbine is less than the power consumed by the compressor. These results are consistent with established knowledge in compressor-turbine matching [11].
The qualitative assessment of the CFD simulation results is now presented. Figures 6–8 show the centerplane distributions of static pressure (), static temperature (), and Mach number (). Contrary to component-level approaches in which the inlet pressure and temperature to the combustor at station 3.1 are prescribed, both and are outputs of the simulations. The centerplane distributions indicate that the static pressure monotonically increases in the engine with shaft angular speed (). The static pressure is nearly constant in the combustor, but slightly larger distributions are observed with increasing . For > 30,000 rpm, the temperature distributions show a larger region of hot gases upstream of the primary and secondary dilution zones. In addition, the static temperature at the exit of the converging nozzle () becomes increasingly nonuniform at larger . A plausible explanation is that the combustor exit temperature is more nonuniform at larger , and the nonuniformity persists through the turbine and nozzle. Moreover, is higher for = 30,000 rpm than it is for = 98,000 rpm, which indicates a lower turbine rotor () efficiency for the former than for the latter condition. The combustor is characterized by subsonic flow for all , whereas the flow through the upstream and downstream components becomes increasingly compressible with increasing . Although sonic flow is observed in some localized regions at ≥ 80,000 rpm, none of the components are choked because sonic flow does not occur uniformly across the cross section of the engine. The Mach number exiting the nozzle increases monotonically with , which suggests that thrust () also increases monotonically.
Traditional thermodynamic cycle analyses with component-level inputs make assumptions about the inlet and outlet boundary conditions such as uniform flow distributions. For this reason, the axial and circumferential Mach number at the compressor diffuser–combustor interface (station 3.1) is examined in Fig. 9, and the static temperature and swirl angle distributions at the combustor-turbine stator interface (station 4) are shown in Fig. 10. The mass-averaged Mach number at station 3.1 increases rapidly from 0.17 to 0.36 as the shaft angular speed increases from 30,000 to 60,000 rpm. The mass-averaged Mach number at station 3.1 remains nearly constant at 0.4 for > 60,000 rpm. The axial Mach number () distributions indicate that the flow at the compressor diffuser-combustor interface is nonuniform with minimal flow immediately downstream of the diffuser vane trailing edges and high flow rates in between the vanes. The axial Mach number reaches approximately 0.5 on the suction side and remains near zero on the pressure side of the diffuser vane when > 60,000 rpm. The circumferential Mach number () distributions in Fig. 9 suggest that the diffuser becomes increasingly deficient at straightening the flow as indicated by the larger circumferential Mach numbers with increasing shaft speed. Even though < , reaches approximately 0.25 when > 60,000 rpm. This corresponds to exit flow angles of approximately 26 deg at the diffusor–combustor interface. In component-level combustor research and development, the combustor is typically isolated in a plenum and the effects of residual swirl from the diffuser are neglected.
The temperature distributions at the combustor exit (station 4) indicate that the temperature varies in both the radial and circumferential directions. The temperature variation increases at higher equivalence ratios () corresponding to the minimum and maximum power. The combustor exit temperature differs from a combustor simulation performed outside of the engine. The combustor has no features that promote a net exit swirl angle. Therefore, if this combustor were analyzed independent of the engine, the swirl angle would be zero. That is, the exit velocity vectors would be perpendicular to the exit plane and the CFD exit flow would be symmetric with respect to the combustor sector centerplane. However, the results in Fig. 10 indicate that there is a maximum exit flow angle of approximately 20 deg for all of the shaft angular speeds. The spatially averaged swirl angle increases with from approximately 8 deg at minimum power to 14 deg at maximum power.
5.2 Comparison Among CFD Simulations, Cycle Analyses, and Measurements.
A quantitative comparison between the CFD simulations and thermodynamic cycle analyses with system-level inputs is presented Table 3, Figs. 11, and 12. Note that only three inputs are necessary for the CFD simulations, whereas twelve inputs are required for the thermodynamic cycle model. Both models need the freestream temperature and pressure as well as the shaft angular speed as inputs. The cycle model additionally requires the following inputs: air mass flow rate (); compressor, combustor, and turbine efficiencies (i.e., and ); combustor exit temperature (); and compressor, diffuser, combustor, and turbine stator pressure ratios (i.e., and ). Instead of providing estimates for these values, the fully coupled engine CFD simulation results are used to provide inputs to the thermodynamic cycle model. This approach enables us to evaluate the consistency between the CFD results and thermodynamic cycle results.
Table 3 indicates that the CFD and cycle results qualitatively and quantitatively agree at most stations throughout the engine from minimum to maximum power. The efficiencies are computed using total enthalpies for the turbomachinery components and sensible enthalpies for the combustor. With increasing , both and increase, while decreases. The positive correlation of with is consistent with the monotonic decrease of with increasing (cf. Fig. 7). Moreover, the negative correlation of with is associated with the monotonic increase in , , , and with . Although and increase monotonically with , their different rates lead to a nonmonotonic equivalence ratio () and hence . The combustion efficiency () decreases nonmonotonically from minimum to maximum power. This is associated with residence time decreasing nonmonotonically from 3.0 ms at minimum power to 2.5 ms at maximum power. The ability to accurately compute the combustion efficiency is limited by the chemistry and turbulent combustion closure models. Furthermore, both CFD and cycle analysis show that , ,, , and monotonically increase while , , and decrease with increasing . In other words, the diffuser, combustor, and turbine stator become increasingly more irreversible at greater since pressure losses are equal to . There is quantitative agreement between the CFD and cycle model in terms of most output values.
Figure 11 illustrates the temperature-entropy (T-s) and pressure-specific volume (P-v) diagrams for the CFD and cycle results. There is agreement between the CFD and cycle results with larger differences observed downstream of the combustor exit (station 4). A plausible explanation is that the turbine stator and nozzle add irreversible losses to the CFD while these components are isentropic in the cycle analysis. In other words, the total enthalpies are different at these components between CFD and cycle analysis. Both models show that with increasing the area within the T-s and P-v diagram increases, indicating greater thrust (). The slopes (in the T-s diagram) between stations 2 and 3 and those between station 4.1 and 5 indicate the irreversible losses in the components. Note that the equal area under process 2-3 and process 4.1-5 indicate that the compressor and turbine power match. The heat addition indicated by the line segment of 3.1-4 suggests that there is slightly more heat added to the cycle model than to the CFD simulations for all conditions. There are minimal differences between the P-v diagrams of the CFD simulations and cycle model because the latter uses pressure ratio inputs from the former (cf. Fig. 4). Both models reveal the total pressure losses associated with diffuser component at larger as indicated by the line segment of process 3-3.1 exhibiting < and < . This causes a triangular shape in the P-v diagrams not seen in the ideal Brayton cycle.
Figure 12 plots the compressor () or turbine () power as function of for the CFD simulations and cycle analyses. The power is similar to within 1% for the CFD and cycle results because the compressor and turbine power are matched using the previously described Secant method. Power increases monotonically and nonlinearly with in both the CFD and cycle results.
In Fig. 13, the CFD and cycle results are compared with experimental data in terms of total pressure (), total temperature (), air mass flow (), fuel flow rate (), global equivalence ratio (), and thrust () as a function of shaft angular speed (). The experimental and are measured between the diffuser and combustor (station 3.1), but the precise location of the measurements is not quantified. Therefore, the flow nonuniformity is quantified over the entire computational interface at station 3.1 in terms of mass-flow-averaged, minimum, and maximum quantities. The minimum and maximum values are indicated by the red shaded areas surrounding the CFD results.
There is qualitative and quantitative agreement between the CFD, cycle, and experimental results. The CFD results for versus and versus agree with the experiments across most conditions. Larger discrepancies for are observed at minimum power (=30,000 rpm), but the result remains within the nonuniformity range indicated by the shaded color in Fig. 13. Larger discrepancies for are observed at the maximum power condition (=98,000 rpm), but the results are within the nonuniformity range as shown by the shaded color in Fig. 13. The computed air and fuel mass flow rate agree with the experiments to within 13% and 7%, respectively. Some of these differences are attributed to the measurement systematic uncertainties as well as the fact that the measured fuel flow rate includes bearing lubrication, which could be as high as 5% at maximum power, as discussed in Sec. 4. The measured and computed equivalence ratios agree at midpower, but larger discrepancies are found at the minimum and maximum power conditions (∼19%). Finally, the measured and computed thrusts agree for most conditions expect for the 80,000 rpm condition (∼25%). Note that the differences between measurements and computations are caused by a combination of model limitations and systematic experimental errors discussed in Sec. 4. Overall, there is good agreement between the measurements, CFD simulations, and cycle analyses. The results demonstrate the potential capability of fully coupled engine CFD simulations along the equilibrium running line to guide the design and to enhance the performance of small-scale turbojet engines.
6 Conclusions
This work develops, demonstrates, and evaluates methods for performing fully coupled CFD simulations along the equilibrium running line of a small-scale turbojet engine from intake to exhaust. The work is motivated by the need for emerging system-level approaches to guide the design, development, and optimization of gas turbine engines.
Two methods are developed and used to enable the CFD simulations of the entire turbojet engine. First, multiple CFD simulations are used in conjunction with the iterative Secant root-finding method to compute the fuel flow rate required to match the compressor and turbine power. Second, a multipoint compressible flamelet/progress variable approach is used to model the combustion without a priori knowledge of the combustor inlet pressure and temperature. The simulations use existing turbulence models and turbomachinery interface models. The only inputs required for the fully coupled CFD simulations are the freestream pressure and temperature and shaft angular speed. In addition, the CFD results provide thermodynamic values such as pressure ratios and component efficiencies without requiring turbomachinery performance maps.
The CFD results indicate that the interfaces between components exhibit heterogeneous flow in both the radial and circumferential directions as expected. The ability of fully coupled CFD simulations to capture the interactions between components leads to more realistic results in comparison to simulations of individual components. Outputs from the CFD simulations are prescribed as inputs to the thermodynamic cycle model to evaluate the consistency of the CFD results with thermodynamic cycle results. The CFD results agree with the thermodynamic cycle results in terms total temperature-entropy and total pressure-specific volume diagrams. In addition, the CFD and cycle results are compared with experimental data. There is agreement between the CFD simulations, cycle analyses, and measurements in terms of diffuser-combustor total temperature and pressure, air and fuel flow rate, equivalence ratio, and thrust as a function of shaft angular speed from minimum to maximum power. The methods and results collectively demonstrate the capability of fully coupled CFD simulations to assist with the design, development, and optimization of future small-scale gas turbine engines on a system-level basis.
Acknowledgment
We thank Eric Hendricks from NASA Glenn Research Center for discussing and providing insights into the pyCycle thermodynamic calculation. We also thank Robert Olding from UDRI for assisting with programming the flamelet code and Nicholas Grannan from ISSI for performing the measurements. Special thanks to Scott Stouffer from UDRI as well as Timothy Erdmann, Joshua Sykes, and Timothy Gallagher from ISSI for insightful discussions. We thank the Department of Defense High Performance Computing Modernization Program for providing the computational resources. This material is based on research sponsored by the U.S. Air Force Research Laboratory under Agreement No. FA8650-15-D-2505. Any opinions, findings and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the United States Air Force. The U.S. Government is authorized to reproduce and distribute reprints for Governmental purposes notwithstanding any copyright notation thereon.
Funding Data
Air Force Research Laboratory (Grant No. FA8650-15-D-2505; Funder ID: 10.13039/100006602).
Nomenclature
- =
engine cross-sectional area (m2)
- =
specific heat capacity (J/kg-K)
- =
eddy dissipation rate model coefficient 1
- =
eddy dissipation rate model coefficient 2
- =
total enthalpy (J/m3)
- =
rotational total enthalpy or rothalpy (J/m3)
- k =
turbulence kinetic energy (m2/s2)
- =
mass flow rate (kg/s)
- =
Mach number
- =
shaft angular speed (rpm)
- =
pressure (Pa)
- =
turbulence Prandtl number
- =
cycle analysis residual equation
- =
temperature (K)
- =
specific entropy (J/kg-K)
- =
generic source term
- =
fuel vapor source term (kg/m3-s)
- =
turbulence Schmidt number
- =
jth component of velocity (m/s)
- =
specific volume (kg/m3)
- =
power (W)
- =
ith Cartesian coordinate (m)
- =
progress variable
- =
mass fraction of species “k”
- =
Favre-averaged mixture fraction
- =
Favre-averaged mixture fraction variance
- =
generic diffusion term
- ε =
turbulence eddy dissipation (m2/s3)
- =
component efficiency
- =
thermal conductivity (W/m-K)
- =
molecular viscosity (kg/m-s)
- =
turbulent viscosity (kg/m-s)
- =
density (kg/m3)
- =
turbulence kinetic energy equation Prandtl number
- =
eddy dissipation rate equation Prandtl number
- Φ =
dependent variables
- =
global equivalence ratio
- =
cell domain reference frame angular velocity vector (rpm)
- =
progress variable source term (kg/m3-s)
- =
air
- =
combustor
- =
compressor
- =
combustor liner cooling
- =
fuel
- =
Cartesian coordinate direction i
- =
Cartesian coordinate direction j
- =
Cartesian coordinate direction k
- =
Secant method iteration number
- =
turbine
- =
total (or stagnation) quantity
- =
PDF table quantity
- =
static quantity
- =
circumferential direction
Footnotes
The tilde is omitted from the equations to simplify explanation of methods and discussion of results.
Rotational total enthalpy (or rothalpy) is . Effectively, rothalpy ( only alters the transient and advection terms of Eq. (1) of the rotors and the total enthalpy ( is effectively solved for the stationary components.
For the CFD model, the shaft angular speed () is equal to the magnitude of the prescribed compressor and turbine rotor domain reference frame rotation. The cell domain reference frame rotation is . The reference frame rotation is nonzero for the compressor and turbine rotors and zero for the remaining components. Because the engine axis is aligned with the Z axis then it follows that .