The current trend of electrification in modern aircraft has driven a need to design and control onboard power systems that are capable of meeting strict performance requirements while maximizing overall system efficiency. Model-based control provides the opportunity to meet the increased demands on system performance, but the development of a suitable model can be a difficult and time-consuming task. Due to the strong coupling between systems, control-oriented models should capture the underlying physical behavior regardless of energy domain or time-scale. This paper seeks to simplify the process of identifying a suitable control-oriented model by defining a scalable and broadly applicable approach to generating graph-based models of thermal, electrical, and turbomachinery aircraft components and systems. Subsequently, the process of assembling component graphs into a dynamical system graph that integrates multiple energy domains is shown. A sample electrical and thermal management system is used to demonstrate the capability of a graph model at matching the complex dynamics exhibited by nonlinear and empirically based simulation models.
Introduction
Modern aircraft are becoming increasingly complex systems-of-systems encompassing multiple energy domains and spanning a wide range of timescales. Traditional modeling and control approaches for complex systems-of-systems is often limited to decentralized high-fidelity modeling with robust, low-performance proportional-integral, and logic-based control [1]. Holistic modeling, analysis, and control design is difficult due to the complexity and size of the systems, especially when dynamics evolve over a wide range of timescales and are characterized by significant interactions among electrical, thermal, hydraulic, pneumatic, and mechanical systems. While these systems may have significantly different dynamics governed by their individual energy domains, each system satisfies a set of conservation equations. Thus, these systems can be unified under the umbrella of “power flow systems” wherein each system satisfies conservation laws and the coupling between systems is characterized by the exchange of power.
With each generation of military and commercial aircraft, there is increased demand for improved performance and efficiency. To achieve these demands, it is imperative to optimize the generation, storage, distribution, and consumption of power through advanced control strategies. Extensive research efforts have focused on these aspects for many types of systems that fall within the category of power flow systems, including microgrids [2], water distribution networks [3], chemical process networks [4], hydraulic hybrid vehicles [5], and thermal energy systems [6].
As aircraft become more complex, the process of developing, analyzing, and validating control designs must be conducted in simulation prior to the application within a physical system. Due to the complexity of these systems, modular simulation toolbox frameworks are often useful for conducting individual component sizing and validation, as well as testing a wide range of system configurations and operating regimes. These toolboxes serve as an excellent tool for developing and analyzing control designs. Several examples of such toolsets include the Thermosys toolbox [7] for modeling air-conditioning and refrigeration systems, the ATTMO toolbox [8] for modeling aircraft vapor cycle systems, and the PowerFlow toolbox [9] for holistic aircraft power system modeling.
While these toolboxes are excellent for testing and designing control algorithms, they are not ideal for generating models that are useful for model-based control designs such as model predictive control [10]. Graph-based approaches to modeling power flow systems have been shown to be particularly convenient for facilitating model-based control design, as shown in Ref. [11] for building thermal dynamics, Ref. [12] for process systems, and Refs. [5], [13], and [14] for vehicle energy management systems. To prove the efficacy of these models for use in control techniques for real-world implementation, it is first essential to demonstrate that graph-based modeling approaches can accurately capture the dynamics of power flow systems. The objective of this paper is to demonstrate the validity of the graph-based modeling framework through empirical validation.
Section 2 provides a general explanation of graphs and how they are used for dynamic modeling. Sections 3–5 detail the development of individual component graph models for thermal, electrical, and air cycle system components, respectively. Each component and system graph is validated experimentally or compared to published data if experimental facilities are not available. Component graphs can be assembled into system graphs, regardless of energy domain, as detailed in Sec. 6 for a candidate aircraft power and thermal management system architecture.
Graph Based Modeling
A graph consists of a set of vertices V and a set of edges E. Figure 1 shows an example of such a graph. When modeling a system as a graph, the vertices represent the capacitative elements of a system where energy can be stored while an edge represents the transport of energy between two vertices, referred to as power flow. Edges are assigned an orientation, denoted by the directional arrows in Fig. 1, which indicates the direction of positive power flow. Power can also be added to the system through sources which are denoted with dashed vertices and edges labeled as for i ∈ {1, 2} in Fig. 1. Finally, power can be rejected to sinks which are denoted by dashed vertices and labeled as for i ∈ {1, 2}. As components are assembled into a system graph, many sources and sinks of a given component are vertices of adjacent components. On a system graph, common source vertices include heat loads, atmospheric conditions, or an infinite bus voltage, while common sink vertices include atmospheric conditions that may affect the power being dissipated by the system.
Generic Graph Formulation.
Let G = (V, E) be an oriented graph that captures the energy storage and power flow throughout a system S. The graph consists of a set of vertices and a set of edges . The orientation of each edge ej represents the positive direction of the associated power Pj from the tail vertex to the head vertex . For the ith vertex, the set of edges directed into the vertex is and the set of edges directed out of the vertex is . A source vertex is defined as a vertex with an indegree of zero, meaning there are no incoming edges, while a sink vertex is defined as a vertex with an outdegree of zero, meaning there are no outgoing edges. Let and denote the source and sink vertices, respectively, such that and . Finally, let denote the Nd dynamic vertices such that .
where captures the graph structure for Vd and represents how power is flowing to the sink vertices.
Thermal Graph Modeling
where Pin is power input to the thermal element and Pout is power output from the thermal element.
where C is referred to as the vertex capacitance and has units of J/K. For most thermal components, C = ρVCp = mCp. Sections 3.1–3.3 detail the process to determine C for each component.
Heat Exchanger.
where 0 ≤ α ≤ 1 is tuned to match heat exchanger behavior. This allows the exit temperature of fluid A to exceed the inlet temperature of fluid B and achieves the behavior discussed for Fig. 2.

Comparison of counter- (a) and parallel-flow (b) heat exchanger temperature profiles for hot (subscript h) and cold (subscript c) flows
where α = 1 for a parallel-flow heat exchanger. The heat transfer coefficient h for both fluid flows is highly dependent upon flow conditions and fluid properties, and can be calculated using empirical correlations [16].
Note that the subscripts a and b denote the fluid and geometric properties on each side of the heat exchanger. Equations (20) and (22) together form the heat exchanger graph model (8), which is provided as Eq. (A1) in the Appendix.
Figure 4 shows a comparison of the graph-based modeling approach to experimental data from a plate heat exchanger in a parallel-flow configuration for a variety of operating conditions. Subplot (a) shows the wall temperature collected experimentally using a surface mounted temperature sensor which is centered on the outer heat exchanger plate. Figure 5 shows a small temperature gradient across the wall relative to the difference in the fluid temperatures on either side of the heat exchanger. Subplot (b) shows the difference between inlet and outlet temperatures for the hot side (side A) and the cold side (side B) of the heat exchanger. Since the hot side has lower flow rates than the cold side, a higher temperature difference is achieved between the inlet and outlet. For both sides, the graph model matches the data very well with maximum errors of less than 0.5 °C, or 5.4% error, which is small enough for model-based control where feedback can compensate for some model error. Therefore, the use of a single temperature state to model the wall is acceptable for most applications. Power flow through the heat exchanger from the hot side to the cold side is shown in subplot (c). The goal of the graph model is to capture the rate at which heat is moving from one fluid flow to another, and Fig. 4 demonstrates very close matching between the data and graph model. Physical and graph heat exchanger parameters are provided in Table 2 of the Appendix.

Comparison of parallel-flow heat exchanger graph model and experimental data for (a) wall temperature, (b) temperature difference between the inlet and outlet of each side, and (c) power flow through the heat exchanger
Liquid–liquid heat exchanger graph parameters
Parameter | Variable | Value | Units |
---|---|---|---|
Wall capacitance | Cw | 900 | J/K |
Fluid A | |||
Capacitance | Ca | 0.0145 | J/K |
Heat transfer area | As,a | 0.2015 | m2 |
Heat transfer coefficient | ha | 10,500 | W/(m2 K) |
Specific heat | Cp | 3500 | J/(kg K) |
Fluid B | |||
Capacitance | Cb | 0.0145 | J/K |
Heat transfer area | As,b | 0.2015 | m2 |
Heat transfer coefficient | hb | 7500 | W/(m2 K) |
Specific heat | Cp | 3500 | J/(kg K) |
Parameter | Variable | Value | Units |
---|---|---|---|
Wall capacitance | Cw | 900 | J/K |
Fluid A | |||
Capacitance | Ca | 0.0145 | J/K |
Heat transfer area | As,a | 0.2015 | m2 |
Heat transfer coefficient | ha | 10,500 | W/(m2 K) |
Specific heat | Cp | 3500 | J/(kg K) |
Fluid B | |||
Capacitance | Cb | 0.0145 | J/K |
Heat transfer area | As,b | 0.2015 | m2 |
Heat transfer coefficient | hb | 7500 | W/(m2 K) |
Specific heat | Cp | 3500 | J/(kg K) |
Cold Plate.
Cold plate heat exchangers are widely used throughout aircraft electronic cooling systems, largely due to their smaller form factor compared to traditional finned heat spreaders. Common cold plate designs feature copper or stainless steel tubing pressed into channels of an aluminum plate. A coolant is pumped through the tubing to transfer heat from the component attached to the surface of the plate.
Figure 6 shows the graph model for the cold plate which includes two dynamic vertices and two edges, denoted by solid lines. The dashed vertex labeled T1 is the inlet condition for the cold plate fluid flow and is determined by an upstream graph component. The outlet temperature T of the cold plate is used to calculate power flows for downstream components, represented by the dashed edge and vertex. The cold plate graph model is very similar to the heat exchanger, except one fluid flow is replaced by a single source vertex and edge, labeled q. This vertex, which could be the temperature of an electrical component, and edge represent the power flow from a heat load that is affecting the cold plate wall temperature Tw, which then transfers heat to the coolant flow and increases the coolant exit temperature T.
The fluid power flow into the cold plate from the inlet condition T1 to the outlet condition T is calculated using Eq. (16) where a = 0, b = Cp, and c = 0. This gives a power flow of the form . Convective heat transfer between Tw and T is calculated using Eq. (16) with a = hAs, b = 1, and c = −1, giving P = hAs(T − Tw), where h is the convective heat transfer coefficient and As is the heat transfer surface area between the cold plate wall and the liquid coolant.
Equations (23) and (24) form the cold plate graph model (8), provided as Eq. (A2) in the Appendix. Figure 7 shows a comparison of the graph-based modeling approach to experimental data of a cold plate. Since a cold plate can have a significant spatial distribution of temperatures depending upon how the load is mounted to the plate, as seen in Fig. 8, several temperature sensors are used to capture the temperature distribution on the cold plate, indicated by the shaded band in the top subplot of Fig. 7. The graph model determines an average wall temperature for the calculation of power flow to the fluid. In the bottom subplot of Fig. 7, the exit fluid temperatures of the graph and experiment are shown to match closely. Graph parameters for the model are provided in Table 3 of the Appendix.

Comparison of cold plate graph model and experimental data for (a) wall temperature and (b) fluid exit temperature
Cold plate graph parameters
Parameter | Variable | Value | Units |
---|---|---|---|
Fluid capacitance | CT | 93.6 | J/K |
Wall capacitance | Cw | 777 | J/K |
Heat transfer area | As | 6.72 × 10−3 | m2 |
Heat transfer coefficient | h | 8500 | W/(m2 K) |
Fluid specific heat | Cp | 3500 | J/(kg K) |
Parameter | Variable | Value | Units |
---|---|---|---|
Fluid capacitance | CT | 93.6 | J/K |
Wall capacitance | Cw | 777 | J/K |
Heat transfer area | As | 6.72 × 10−3 | m2 |
Heat transfer coefficient | h | 8500 | W/(m2 K) |
Fluid specific heat | Cp | 3500 | J/(kg K) |
Fluid Tank.
Aircraft contain many fluid tanks for the storage of fuel, oil, hydraulic fluid, and various coolants. Each of these tanks store mass and thermal energy that can be transported around the aircraft. In the case of fuel, mass is removed from the tanks for combustion in the engine. The effect is a continual decrease of fuel tank thermal capacitance, which must be given special consideration in a graph model.
Figure 9 shows the graph model for the fluid tank which includes a single dynamic vertex and two edges, denoted by solid lines. The dashed vertex T1 is the temperature of the fluid flowing into the tank, where is the return mass flow rate. The temperature of the tank T is used to calculate the power flow out of the fluid tank, which is a function of the outlet mass flow rate . The final edge and vertex capture the energy loss from the fluid tank as a result of heat transfer with ambient conditions. This loss may be negligible in some cases and in such instances can be omitted.
The fluid power flow into the fluid tank from the upstream component is calculated using Eq. (16) where a = 0, b = Cp, and c = 0. This gives power flow into the fluid tank temperature T as . The power flow out of the tank is similar, with a = 0, b = Cp, and c = 0, giving . The heat transfer between the fluid tank and the ambient conditions is a function of the heat transfer coefficient h and the surface area As between fluid in the tank and the ambient fluid, including thermal resistance of the wall. This power flow can be calculated using Eq. (16) where a = hAs, b = 1, and c = −1, giving P = hAs (T − Tamb). Note that this is an oriented graph which associates positive power flow with leaving the fluid tank, but if T < Tamb then the sign of the power flow will become negative, resulting in power flow into the tank. The orientation of the edge is only reflective of the direction of positive power flow, not the physical direction of power flow.
Figure 10 shows a comparison of the graph-based modeling approach to experimental data of a fluid tank over a range of operating conditions. The experimental tank is small enough such that temperature gradients in the fluid are negligible. The graph matches the experimentally collected data very well, with a max difference of 1 °C, or 4.3% error, shown in the figure insert. The tank graph parameters are provided in Table 4 in the Appendix.
Liquid tank graph parameters
Parameter | Variable | Value | Units |
---|---|---|---|
Tank capacitance | CT | 3800 | J/K |
Heat transfer area | As | 51.0 × 10−3 | m2 |
Heat transfer coefficient | h | 15 | W/(m2 K) |
Fluid specific heat | Cp | 3500 | J/(kg K) |
Parameter | Variable | Value | Units |
---|---|---|---|
Tank capacitance | CT | 3800 | J/K |
Heat transfer area | As | 51.0 × 10−3 | m2 |
Heat transfer coefficient | h | 15 | W/(m2 K) |
Fluid specific heat | Cp | 3500 | J/(kg K) |
Electrical Graph Modeling
Electrical systems consist of generation, distribution, transformation, and dissipation components. Generators attached to aircraft engines provide the electrical power that is transformed via rectifiers and transformers and distributed to electrical alternating current (AC) and direct current (DC) buses at varying voltages. Constant current, power, or impedance loads consume electrical power and produce heat as a byproduct of inefficiencies.
Since the purpose of the graph modeling framework is the development of models for model-based control, there is no distinguishing between AC and DC power because the timescales associated with the AC power waveform are too fast for the system behavior that the graph is intended to capture. Assuming a balanced three-phase system, the graph model will capture the behavior of the d-component in a synchronous reference frame aligned with the voltage [17,18]. The d-component comes from a dq0 transform and corresponds to the magnitude of the voltage. For DC power, the graph model captures the behavior of the DC voltage.
where x is the state of an adjacent vertex, u is an input, and a, b, and c are coefficients defining the edge power. This is slightly different from thermal components where x was always a temperature and u was always a mass flow rate.
Generator.
Aircraft generators are mounted to gearboxes on each engine or an auxiliary power unit which provide the input mechanical power to generate electrical power. Throughout an aircraft mission, the engine shaft speed can vary depending upon the throttle position. Without a generator control unit regulating the AC voltage, increases and decreases in shaft speed will affect the voltage and frequency of the AC signal.
where ω is the shaft speed of the generator, Vgen is the generator voltage, and linear coefficients α and β are identified based upon open circuit voltage of the generator for various shaft speeds. The generator dynamics are given by Eq. (26) where Ci is tuned to match the open loop response of the generator. Generator parameters are provided in Table 5 in the Appendix.
Electric system graph parameters
Parameter | Variable | Value | Units |
---|---|---|---|
Generator capacitance | Cg | 0.15 | s |
Bus 1 capacitance | Cb,1 | 0.01 | s |
Bus 2 capacitance | Cb,2 | 0.01 | s |
Linear coefficients | |||
Shaft speed | α1 | 2000 | V s |
Generator voltage | β1 | −8900 | None |
Generator voltage | α2 | 650 | None |
Bus 1 voltage | β2 | −504 | None |
Generator voltage | α3 | 420 | None |
Bus 2 voltage | β3 | −820 | None |
Efficiency | η | 0.95 | None |
Constant resistance | R | 5 | Ω |
Parameter | Variable | Value | Units |
---|---|---|---|
Generator capacitance | Cg | 0.15 | s |
Bus 1 capacitance | Cb,1 | 0.01 | s |
Bus 2 capacitance | Cb,2 | 0.01 | s |
Linear coefficients | |||
Shaft speed | α1 | 2000 | V s |
Generator voltage | β1 | −8900 | None |
Generator voltage | α2 | 650 | None |
Bus 1 voltage | β2 | −504 | None |
Generator voltage | α3 | 420 | None |
Bus 2 voltage | β3 | −820 | None |
Efficiency | η | 0.95 | None |
Constant resistance | R | 5 | Ω |
Bus.
Aircraft electrical buses distribute power throughout the aircraft to multiple loads at varying voltages. Generator voltage is stepped up or down through a transformer unit, and potentially converted to DC voltage. Since the graph model does not distinguish between AC and DC voltages, graph edge parameters are based upon transformer turn ratios in order to capture the desired voltage increase or decrease.
where N1/N2 is the turn ratio for the electrical transformer. The input u regulates the bus voltage when loads turn on or off by affecting N1/N2, and often is set by a feedback controller.
The bus dynamics are given by Eq. (26) where Ci is chosen to be at least one order of magnitude smaller than the generator. If it is assumed that the bus voltage is regulated perfectly, then Ci = 0. This maintains the bus voltage at the set point and all loads attached to the bus pass to the generator.
Loads.
Three main types of electrical loads are considered: constant power, constant current, and constant impedance. A graph model containing each of the loads is shown in Fig. 13.
where η is the electrical load efficiency.
where R is the impedance of the load.
A constant impedance load is converted solely into heat. Therefore, the edge in Fig. 13 representing the constant impedance load is connected directly to the sink vertex T which represents the temperature of electrical load. This vertex could be the source vertex in a thermal component graph.
Since the graph-based modeling approach does not need to distinguish between AC and DC power, the following holds true for constant resistance loads as well.
Validation of Graph Model.
The electrical graph is validated using nonlinear electrical system components from the PowerFlow toolset which has previously been validated [9]. The architecture pictured in Fig. 14 is modeled in PowerFlow using system sizing parameters from Ref. [19] to represent a Boeing 787 electrical system. The electrical system features a generator operating at a nominal 230 V, a 270 V DC bus, and a 115 V AC bus. Each bus connects to a set of constant impedance, power, and current loads. A graph model is developed to represent the electrical architecture, shown in Fig. 15. The graph model is a combination of the generator, bus, and load graphs from Secs. 4.1–4.3. The unlabeled sink vertices correspond to the constant impedance loads, while the other sink vertices correspond to the constant power and current loads. In the Appendix, the system graph model is given as (A4) with parameters provided in Table 5.
Validation is performed using open loop responses. Therefore, generator and bus voltages are not regulated to their set points in order to test the ability of the graph to capture transients and steady-state offsets from nominal voltages when generator shaft speed is increased and loads are turned on and off. In Fig. 16(a) the generator shaft speed is shown to change between 10,000 and 14,000 RPM. Figures 16(b) and 16(c) show the varying constant power and constant current loads input to the graph and PowerFlow models. Resistive loads are held constant at 5 Ω since the load magnitude will change as voltage changes. Peak DC loads are 54 kW, and peak AC loads are 22 kW. Figure 17(a) shows the transient voltages for the generator, where the shaft speed is the largest contributor to changing voltage. Detail (c) shows a close-up of the graph and PowerFlow model matching very well during transient events. Figure 17(b) shows close matching voltages of the 270 V DC bus and the 115 V AC bus. A close-up in detail (d) shows that there is a slight steady-state offset with the graph model, but when a load is turned off shortly before 1300 s, both models react accordingly with an open-loop bus voltage increase. This steady-state error is not prohibitive to the intended use of the graph-based model for control design since feedback can eliminate small model errors.

Open loop inputs for (a) generator rotational shaft speed, (b) constant power AC and DC loads, and (c) constant current AC and DC loads

Comparison of graph model and nonlinear simulation (a) generator voltage, (b) 270 V and 115 V bus voltages, (c) transient generator voltage, and (d) loads affecting bus voltage, for inputs from Fig. 17

Comparison of graph model and nonlinear simulation (a) generator voltage, (b) 270 V and 115 V bus voltages, (c) transient generator voltage, and (d) loads affecting bus voltage, for inputs from Fig. 17
Turbomachinery Graph Modeling
Gas turbine engines and air cycle machines (ACMs) are the primary forms of turbomachinery found in modern aircraft. The engine, or sometimes the auxiliary power unit, is the prime mover from which all power in the aircraft is extracted, while the air cycle machine acts as the refrigeration unit in the environmental control system. Although each of these systems exists for very different purposes, their dynamic behavior in a graph framework is very similar.
where is the angular acceleration. The right-hand side of Eq. (38) is equivalent to Newton's second law for rotational bodies, which states that where τ is the net torque acting on the body.
Air Cycle Machine.
The air cycle machine is the primary refrigeration unit in the environmental control system for modern commercial aircraft. Figure 18 shows a traditional closed-loop configuration of an ACM operating in a reverse Brayton cycle where air is compressed from (1) to (2) effectively increasing the air pressure and temperature, heat is rejected from (2) to (3), air temperature and pressure drop as it is expanded from (3) to (4), and heat is absorbed from loads between (4) and (1). In an open-loop ACM configuration, inlet air to the compressor (1) comes from another source such as the engine bleed air. Additionally, ACMs may have multiple turbines to provide extra power. In Fig. 18, a power turbine is shown to the right of the dashed line. This turbine expands air from (5) to (6) in order to provide extra power to the shaft which helps increase the refrigeration capability of the ACM. Typically, this air is bled from the compressor of the aircraft engine and is discharged to ambient.
A graph model of the ACM configuration in Fig. 18 is shown in Fig. 19. The graph consists of six temperature vertices, a shaft speed vertex, and thirteen edges, denoted by solid lines. The dashed vertex Tbleed is a source vertex for the temperature of the bleed air that is entering the power turbine. The heat that is absorbed by the ACM is represented by the dashed vertex and edge entering the vertex TK,in and the heat rejected by the ACM is represented by the dashed vertex and edge exiting the vertex TET,in.
where the right-hand side is the summation of two independent temperature terms, and the sign convention corresponds to the orientation of edges entering and leaving the ω vertex in Fig. 19. This introduces the first form of edge power equations for the ACM, which are denoted with green and blue lines in Fig. 19, and follow the form of Eq. (16) where . For edges ei, i ∈ {1, 2, 5, 6, 8, 11, 13} Eq. (16) is written with a = 0, b = Cp, and c = 0, while for edges ei, i ∈ {4, 7, 9} Eq. (16) is written with a = 0, b = 0, and c = Cp. This difference is a result of the negative sign in Eq. (42) and the corresponding orientation of edges in Fig. 19.
where α, β, and γ are identified linear coefficients for a given ACM at steady-state. Since by definition e7 = e11 and e9 = e13, in order to reach steady-state (43) will equal zero during steady-state operation.
Both mass flow rate relationships are provided as Eq. (A5) in the Appendix.
Throughout the rest of the paper, the graph shaft dynamics use Eq. (50) with ω0 included in the capacitance term.
Empirical Validation of Graph Model.
Full-scale experimental testing of ACMs is difficult due to the substantial infrastructure required to generate the ram, bleed, and boundary conditions for a flight envelope that an ACM would operate within. Furthermore, airborne testing is often infeasible due to the certification processes associated with installation of testing equipment and the time commitment required of the aircraft. For these reasons, modeling of ACMs and environmental control systems is well documented [20–24], but often not validated with flight data.
Limited work has been done with experimental ACM studies. In Ref. [25], a traditional two-wheel bootstrap ACM was tested using air compressors to condition bleed and ram air for various altitudes. Steady-state pressure, humidity-ratio, and ACM exit temperature are measured and reported. In Ref. [26], a two-wheel bootstrap ACM from an unidentified BAE Systems (Farnborough, UK) aircraft was tested for a limited set of flight conditions in order to validate a one-dimensional thermodynamic model, which was used to extrapolate ACM performance over an entire flight envelope. The data that are used to validate the steady-state behavior of the graph model come from Ref. [27] where a two-wheel bootstrap ACM was powered by a bleed-air driven power turbine, which is identical to the configuration in Sec. 5.1. The data are collected at six points along the edges of the flight envelope as shown in Fig. 20. The graph model is simulated over a mission profile that begins at point 1 and proceeds to point 6 in ascending order. The mission profile was designed such that the ACM graph remains at each flight envelope point for 500 s, allowing the system to reach steady-state. Graph source and sink vertices were defined using boundary conditions in Table 1 of Ref. [27]. Linear coefficients along edges 3, 10, and 11 in Fig. 19 were tuned such that steady-state temperatures in the graph matched the steady-state temperatures collected from Figs. 8–11 in Ref. [27].
Computational difference between graph and nonlinear model for 25 simulations
Graph | Nonlinear model | ||
---|---|---|---|
tavg (s) | St. Dev. (s) | tavg (s) | St. Dev. (s) |
25.9 | 5.26 | 277 | 19.9 |
Graph | Nonlinear model | ||
---|---|---|---|
tavg (s) | St. Dev. (s) | tavg (s) | St. Dev. (s) |
25.9 | 5.26 | 277 | 19.9 |
where Q is the heat load into the ACM, is the mass flow rate of bleed air going through the power turbine, Tbleed is the bleed air temperature, and Tamb is the ambient temperature at the exit of the power turbine. Figure 21 compares the graph model to the data presented in Ref. [27]. Inlet and outlet temperatures for each turbine and the compressor are compared in Fig. 22.
Figures 21 and 22 show close matching between steady-state conditions in the graph model and the data from Ref. [27]. The graph model calculates COP within 2.6% of each data point, while temperatures are within 11%. While there is a larger difference between the graph model and the data from Ref. [27], it is important to mention that the purpose of a control-oriented ACM model is to identify the effect of bleed air input to the amount of heat removed from the system, meaning that matching COP is the first priority. Second, without detailed knowledge of the mass flow rates or heat exchanger parameters used within [27], it is difficult to accurately estimate graph parameters. In Sec. 5.1.2, a transient validation will be presented where a high-fidelity ACM model is used for comparison. The advantage of such a model is the ability to easily identify the relationships between turbomachinery temperatures, shaft speed, and mass flow rate to help with the selection of graph parameters.
Transient Validation of Graph Model.
A high-fidelity ACM model is developed using the ATTMO toolbox [8] in the configuration of Fig. 18. The secondary side of each heat exchanger is modeled with a source fluid flow. Heat input qin is from a fuel loop that is being cooled by the ACM, while heat output qout from the ACM is rejected to engine bypass air. The inlet temperatures for each fluid flow are specified in the simulation. The graph is defined in Eqs. (A6)–(A11) with parameters specified in Table 6 of the Appendix.
Air cycle system graph parameters
Parameter | Variable | Value | Units |
---|---|---|---|
Power turbine | |||
Inlet capacitance | CPT,in | 1000 | J/K |
Outlet capacitance | CPT,out | 1000 | J/K |
Temperature coefficient | α | 254 | J/(K s) |
Mass flow rate coefficient | β | −11,250 | J/(kg s) |
Constant coefficient | γ | −43,400 | J/s |
Fluid specific heat | Cp | 1073 | J/(kg K) |
Bleed air temperature | Tbleed | 700 | K |
Compressor | |||
Inlet capacitance | CK,in | 100 | J/K |
Outlet capacitance | CK,out | 100 | J/K |
Temperature coefficient | α | 1000 | J/(K s) |
Mass flow rate coefficient | β | 3.115 × 106 | J/(kg s) |
Constant coefficient | γ | 0 | J/s |
Fluid specific heat | Cp | 1030 | J/(kg K) |
Expansion turbine | |||
Inlet capacitance | CET,in | 100 | J/K |
Outlet capacitance | CET,out | 100 | J/K |
Temperature coefficient | α | 1000 | J/(K s) |
Mass flow rate coefficient | β | −5.675 × 105 | J/(kg s) |
Constant coefficient | γ | 0 | J/s |
Fluid specific heat | Cp | 1008 | J/(kg K) |
Fuel-compressor heat exchanger | |||
Wall capacitance | CHX,f | 964 | J/K |
Outlet capacitance | Cout,f | 100 | J/K |
Heat transfer coefficient K side | hk | 17,500 | W/(m2 K) |
Heat transfer area K side | Ak | 1.86 | m2 |
Heat transfer coefficient fuel side | hf | 10,000 | W/(m2 K) |
Heat transfer area fuel side | Af | 1.86 | m2 |
Fuel specific heat | Cp | 2000 | J/(kg K) |
Air-expansion turbine heat exchanger | |||
Wall capacitance | CHX,f | 9500 | J/K |
Outlet temperature capacitance | Cout,f | 100 | J/K |
Heat transfer coefficient ET side | hET | 17,500 | W/(m2 K) |
Heat transfer area ET side | AET | 11.8 | m2 |
Heat transfer coefficient air side | ha | 17,500 | W/(m2 K) |
Heat transfer area air side | Aa | 11.8 | m2 |
Air specific heat | Cp | 1008 | J/(kg K) |
Parameter | Variable | Value | Units |
---|---|---|---|
Power turbine | |||
Inlet capacitance | CPT,in | 1000 | J/K |
Outlet capacitance | CPT,out | 1000 | J/K |
Temperature coefficient | α | 254 | J/(K s) |
Mass flow rate coefficient | β | −11,250 | J/(kg s) |
Constant coefficient | γ | −43,400 | J/s |
Fluid specific heat | Cp | 1073 | J/(kg K) |
Bleed air temperature | Tbleed | 700 | K |
Compressor | |||
Inlet capacitance | CK,in | 100 | J/K |
Outlet capacitance | CK,out | 100 | J/K |
Temperature coefficient | α | 1000 | J/(K s) |
Mass flow rate coefficient | β | 3.115 × 106 | J/(kg s) |
Constant coefficient | γ | 0 | J/s |
Fluid specific heat | Cp | 1030 | J/(kg K) |
Expansion turbine | |||
Inlet capacitance | CET,in | 100 | J/K |
Outlet capacitance | CET,out | 100 | J/K |
Temperature coefficient | α | 1000 | J/(K s) |
Mass flow rate coefficient | β | −5.675 × 105 | J/(kg s) |
Constant coefficient | γ | 0 | J/s |
Fluid specific heat | Cp | 1008 | J/(kg K) |
Fuel-compressor heat exchanger | |||
Wall capacitance | CHX,f | 964 | J/K |
Outlet capacitance | Cout,f | 100 | J/K |
Heat transfer coefficient K side | hk | 17,500 | W/(m2 K) |
Heat transfer area K side | Ak | 1.86 | m2 |
Heat transfer coefficient fuel side | hf | 10,000 | W/(m2 K) |
Heat transfer area fuel side | Af | 1.86 | m2 |
Fuel specific heat | Cp | 2000 | J/(kg K) |
Air-expansion turbine heat exchanger | |||
Wall capacitance | CHX,f | 9500 | J/K |
Outlet temperature capacitance | Cout,f | 100 | J/K |
Heat transfer coefficient ET side | hET | 17,500 | W/(m2 K) |
Heat transfer area ET side | AET | 11.8 | m2 |
Heat transfer coefficient air side | ha | 17,500 | W/(m2 K) |
Heat transfer area air side | Aa | 11.8 | m2 |
Air specific heat | Cp | 1008 | J/(kg K) |
The graph model in Fig. 19 is adapted to include the heat exchangers using the graph model described in Sec. 3.1. The new graph is shown in Fig. 23 with additional Twall and Tout vertices for each heat exchanger. Source vertices Tfuel and Tbypass are identical to the inlet temperatures specified for the high-fidelity model. Model comparison is conducted by stepping inputs and disturbances for each model. Fuel and bypass air inlet temperatures are stepped between their minimum and maximum points for ideal operation of the ACM, while the bleed air flow into the power turbine is stepped between high- and low-demand operation.
Heat absorption and rejection by the ACM are compared in Fig. 24. Subplot (a) shows the graph model matching the high-fidelity model very well in heat rejection to engine bypass, while detail (c) shows a close-up of the transient matching by the graph model. Similarly with subplot (b), the graph model matches the high-fidelity model within 4% at steady-state conditions, and detail (d) shows good transient matching even with the steady-state offset. Shaft speed of the ACM is primarily affected by changes in the mass flow rate into the power turbine. Figure 25 shows the comparison between the two models. While there is occasionally a 2% steady-state offset between the two models, the error is negligible in regard to controlling the performance of the ACM, specifically the heat rejection. More important is the role of shaft speed in calculating the ACM mass flow rate which directly affects the heat transfer capabilities of the ACM. During transients (insert of Fig. 25), it is important to capture the dynamics of the shaft so that correct mass flow rates are used for heat transfer calculations. Figure 26 shows the power produced or consumed by the turbines and compressor. For the graph model, powers are calculated using Eq. (42). A positive sign indicates power generated while a negative sign indicates power consumption. In subplot (a), the power turbine generates a significant amount of power compared to the expansion turbine due to sizing of each component. Subplot (b) shows the consumption of power by the compressor, which in steady-state is equivalent to the sum of the power and expansion turbine. Detail (c) and (d) show almost identical matching between the graph and high-fidelity model transients.

(a) Heat rejected to the bypass air, (b) heat absorbed from the fuel by the ACM, and (c) and (d) detail showing matching transient behavior by the graph model

(a) Power produced by the power turbine (top) and expansion turbine (bot), (b) power consumed by the compressor, and (c) and (d) detail showing matching transient behavior by the graph model
System Graph Models
The major advantage of graph models is the ability to assemble individual component graphs into a system graph representation. A candidate aircraft power system schematic consisting of electrical, thermal, and air cycle systems is shown in Fig. 27. This architecture is representative of a single engine aircraft, but could easily be doubled for larger dual-engined commercial aircraft, in which case the graph model would simply double as well. In this configuration, the engine is the prime source of power, providing mechanical power to the generator and bleed air from the compressor to the power turbine. Two electrical buses provide 115 V and 230 V power to avionic loads and a large radar load. Each electrical load is assumed to generate waste heat and is therefore cooled by fuel cold plates. The radar load has a dedicated fuel tank that helps mitigate large transient loads from the radar. Fuel from fuel tank #1 is pumped around the aircraft as a coolant for avionics, generator, engine oil, and fuel tank #2. After absorbing heat from each source, the fuel is cooled by a liquid––air heat exchanger on the ACM. The ACM rejects the heat through a heat exchanger in the bypass duct of the aircraft engine. Heat can also be rejected from fuel tank #1 through a pumped liquid loop through a liquid–air heat exchanger using ram air as the coolant.
Component graphs from Secs. 3–5 are used to create the full system graph shown in Fig. 28. Layout of the graph attempts to match the system schematic. Gray vertices represent system sinks and sources such as the bypass, ram, and bleed air conditions, constant power and constant current loads, and engine shaft speed. Sink and source vertices on individual component graphs become neighboring system vertices in the system graph. For example, vi for i = [8, 19, 4] is a cold plate graph, where v8 is the inlet temperature, v4 is the wall temperature, and v19 is the outlet temperature from the cold plate. Heat into the cold plate is e14 which is algebraically related to the generator load.
System Graph Validation.
A nonlinear model of the system in Fig. 27 is developed using PowerFlow [9] for the electrical system, ATTMO [8] for the air cycle system, and Thermosys [7] for the single phase thermal-fluid system. Edge and vertex parameters for the graph in Fig. 28 were identified using physical parameters such as mass, material properties, and average heat transfer coefficients of the nonlinear model. Both systems are simulated in open loop, and independently such that no states are shared between the nonlinear model and the graph model for the entirety of the simulation. Feedback control is used for the electrical system to maintain generator and bus voltages, unlike the analysis in Sec. 4.
In Figs. 29 and 30, a comparison between the graph model and nonlinear model shows how well the graph can match high-fidelity model temperatures, voltages, and power flows. In subplot (a)–(e) of Fig. 29, graph and nonlinear model temperatures for the most important components of the system match for most of the simulation. In several instances, there are deviations of several degrees, but even by the end of the 8000 s mission, the temperatures are nearly identical. Subplot (f)–(h) show generator voltage, 270 V bus voltage, and 115 V bus voltage. Several transients that occur during large loading events causes deviations up to 4 V between the two models, which is attributed to different closed-loop control algorithms used between the graph and nonlinear model.

Validation of graph (a) avionics wall temperature, (b) generator and engine temperatures, (c) radar temperature, (d) fuel tank #2 temperature, (e) fuel tank #1 temperature, (f) generator voltage, (g) 270 V bus voltage, and (h) 115 V bus voltage

Validation of graph power flow for (a) fuel heat rejection along e43, (b) ACM heat rejection along e58, and (c) ram air heat rejection along e19
In subplot (a) of Fig. 30, the heat absorbed by the ACM through e43 shows close matching, especially during transients. A similar trend is seen in subplot (b) which is heat rejection from the ACM to the bypass duct air through e15. Subplot (c) shows the heat rejected by the fuel–air heat exchanger that is dedicated to cooling fuel tank #1 via e37 in the graph. There is a large deviation between the two models during the second half of the mission, which is likely caused by the fixed heat transfer coefficient used by the graph, while the nonlinear model has correlations being used throughout different operational regimes. However, the error is within 16%, which is acceptable for an independent simulation. When implemented in a model-based controller, the graph model would be updated with measured states in order to correct any errors that exist between the model and system.
What makes the graph model advantageous to the nonlinear model is the speed comparison displayed in Table 1. An 8000 s simulation was executed 25 times, and the graph model completed the simulation with an average time of 25.9 s, which was a full order of magnitude faster than the 277 s nonlinear model average simulation time. Faster computational times are beneficial for reducing idle time while models are running, increasing the potential design space that can be explored through design iterations, and help model-based controllers execute faster. The simulations presented in Table 1 were run on a standard desktop computer with an Intel i7-6700, 16 GB of RAM, running 64-bit Windows 10.
Discussion of Results.
The overall goal of a graph model is to capture system dynamics, provide a control-oriented model, and be computationally fast. From each validation figure in Secs. 3–5, it is clear that a graph model can accurately capture individual component behaviors. Figures 29 and 30 show that when these component graph models are combined into a full system graph, there is little error introduced and the graph model can accurately capture complex dynamics. Furthermore, when placed in a model-based controller state, updates would occur frequently, meaning large errors would be unlikely to happen. Additionally, the derivation of each model contains a controllable input that can be used in the development of models and controllers to determine the effect of control on the component or system. Finally, it is important for models to be computationally efficient such that controllers can execute rapidly, which is achievable as shown in the computational results of Table 1.
Conclusions
This paper presents the development and validation of dynamical graph models for thermal, electrical, and air cycle systems that are common on commercial and military aircraft. Graph models are designed to be modular and reconfigurable such that multiple system configurations and control strategies can be rapidly analyzed. Additionally, graph models are computationally efficient while accurately capturing dynamical interactions across multiple energy domains and timescales, which facilitates their implementation on embedded processors for system, subsystem, or component controllers. Each component graph model is based upon conservation of energy and is derived such that vertices correspond to dynamic states and edges represent the flow of energy. Thermal fluid experiments are used to validate the thermal graph models, while published experimental data and high-fidelity models are used to validate the electrical and air cycle system models. A complete system graph is constructed using each of the validated component graphs. The system graph is compared to a high-fidelity simulation model to show that a compiled graph is capable of matching system behavior at an order of magnitude faster computational speed.
Funding Data
- •
National Science Foundation (Grant No. DGE-1144245 and Agreement No. EEC-1449548).
Nomenclature
- A =
area, m2
- Ci =
capacitance for vi
- Cp =
specific heat, J/(kg K)
- E =
set of graph edges
- ei =
ith graph edge
- h =
heat transfer coefficient, W/(m2 K)
- I =
current, A
- L =
angular momentum, J s
- m =
mass, kg
- =
mass flow rate, kg/s
- M =
graph incidence matrix
- Nd =
number of graph dynamic vertices
- Ne =
number of graph edges
- Ns =
number of graph source vertices
- Nt =
number of graph sink vertices
- Nv =
number of graph vertices
- P =
power, W
- q =
heat, W
- R =
resistance, Ω
- t =
time, s
- T =
temperature, K
- u =
input
- V =
set of graph vertices
- vi =
ith graph vertex
- Vd =
dynamic vertices set
- Vs =
source vertices set
- Vt =
sink vertices set
- β =
scalar coefficient
- γ =
scalar coefficient
- η =
efficiency
- ρ =
density, kg/m3
- τ =
torque, N·m
- ω =
rotational speed, rad/s