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 35 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 $Piin$ for i ∈ {1, 2} in Fig. 1. Finally, power can be rejected to sinks which are denoted by dashed vertices and labeled as $vit$ 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.

Fig. 1
Fig. 1
Close modal

### 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 $V={vi:i∈{1,2,…,Nv}}$ and a set of edges $E={ei:i∈{1,2,…,Ne}}$. The orientation of each edge ej represents the positive direction of the associated power Pj from the tail vertex $vjtail$ to the head vertex $vjhead$. For the ith vertex, the set of edges directed into the vertex is $Eiin={ej:vjhead=vi}$ and the set of edges directed out of the vertex is $Eiout={ej:vjtail=vi}$. 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 $Vs∈IRNs$ and $Vt∈IRNt$ denote the source and sink vertices, respectively, such that $Vs⊂V$ and $Vt⊂V$. Finally, let $Vd∈IRNd$ denote the Nd dynamic vertices such that $Vd⊂V∖(Vs∪Vt)$.

The dynamic state xi for vertex viVd represents the stored energy of that vertex, and the dynamics satisfy conservation of energy
$Cix˙i=∑ej∈EiinPj−∑ej∈EioutPj$
(1)
where the right-hand side is a summation of all edges oriented into the vertex minus all edges oriented out of the vertex. The power along each edge is constrained to be a function of the vertex states and an input uj
$Pj=fj(xjtail,xjhead,uj)$
(2)
The incidence matrix [15] $M=[mij]∈IR(Nd+Nt)×(Ne−Ns)$ captures the structure of the graph and is given by
$mij={+1vi is the tail of ej−1vi is the head of ej0else}$
(3)
and can be partitioned
$M=[M¯M¯]$
(4)

where $M¯∈IRNd×(Ne−Ns)$ captures the graph structure for Vd and $M¯$ represents how power is flowing to the sink vertices.

The dynamical graph system S can be written as
$S:Cx˙=−M¯P+DPin$
(5)
where $C=diag([Ci])∈IRNd×Nd, x=[xi]∈IRNd$ are the states of all dynamic vertices, $P=[Pi]∈IRNe−Ns$ is a vector of each edge power in G that is not from a source vertex, and $Pin=[Piin]∈IRNs$ is a vector of Ns power flows from source vertices. The matrix $D=[dij]∈IRNd×Ns$ captures the input of power from the source vertices, and its structure is determined by
$dij={+1vi is the head of Pjin0else}$
(6)
The vector of edge power flow in S is represented as
$P=F(x,xt,u)$
(7)
where $F(x,xt,u)=[fj(xjtail,xjhead,uj)]$, resulting in the nonlinear dynamics of S
$S:Cx˙=−M¯F(x,xt,u)+DPin$
(8)

## Thermal Graph Modeling

Thermal energy components can be split into two primary categories: storage and transport. Fuel tanks, oil tanks, pressurized air cabins, and phase-change materials act as thermal energy storage elements, while heat exchangers, pumps, fans, pipes, ducts, splits, and junctions are used to transport thermal energy throughout an aircraft. The storage of thermal energy is based upon the conservation of energy, more specifically the first law of thermodynamics which states U = Q + W. Assuming no work W is being done by the system, then the internal energy U for a thermal element can be written as
$U=mCpT$
(9)
The time differential of the energy added to the system Q can be represented as the net power that is flowing into the element
$ddtQ=Pin−Pout$
(10)

where Pin is power input to the thermal element and Pout is power output from the thermal element.

Therefore, to satisfy conservation of energy, the temperature dynamic for each thermal graph element can be written as
$ddtU=dmdtCpT+mCpT˙=Pin−Pout$
(11)
The term $(dm/dt)CpT$ captures the change in energy as a result of changing mass, while $mCpT˙$ is the change in energy as a result of changing temperature. For fuel tanks in aircraft, mass can change significantly over the course of a flight, but for most other components
$dmdtCpT+mCpT˙≈mCpT˙=ρVCpT˙$
(12)
where V is the volume of the fluid. The case of changing mass is addressed in Sec. 3.3. For cases in which mass is fixed, each thermal graph vertex temperature dynamic is written as
$CT˙=Pin−Pout$
(13)

where C is referred to as the vertex capacitance and has units of J/K. For most thermal components, C = ρVCp = mCp. Sections 3.13.3 detail the process to determine C for each component.

Power in Eq. (13) is modeled by one of two methods. Power flow via mass transfer from one component to another
$P=m˙CpT$
(14)
and power flow via conductive or convective heat transfer
$P=hAΔT$
(15)
where h is the heat transfer coefficient and A the heat transfer area. Power flows can be modeled with a generic equation
$P=(a+m˙)(bT1+cT2)$
(16)

For edges that do not have an associated mass flow rate, $m˙=0$ and a ≠ 0. In Secs. 3.13.3, Eqs. (13) and (16) are used to generate graph models of a heat exchanger, cold plate, and fluid tank.

### Heat Exchanger.

Plate heat exchangers are the primary type of heat exchangers found in aircraft thermal systems due to the high heat flux that can be achieved in a compact form factor. Typically, plate heat exchangers are configured such that the fluids flow counter or parallel to each other. In a counter-flow configuration, the exit temperature of the hot fluid may be less than the exit temperature of the cold fluid; however, in a parallel-flow configuration, the temperature of the two flows cannot cross, as shown in Fig. 2. A single graph model can be used to represent either flow configuration. The graph model for a heat exchanger is shown in Fig. 3. This includes three dynamic vertices and four edges, denoted by solid lines. The dashed source vertices T1,a and T1,b are the inlet conditions for each fluid flow that are determined by upstream graph components. The outlet temperatures for each fluid flow Ta and Tb are used to calculate power flows for downstream components, represented by the dashed edges and dashed vertices on the right. Between the outlet temperatures, the wall temperature Tw captures the average temperature of the heat exchanger wall and facilitates the transfer of heat from fluid flow B to fluid flow A, denoted by the orientation of the edges. The thermal power flows from the inlet conditions T1,a and T1,b to the outlet conditions Ta and Tb are calculated by Eq. (16) using a = 0, b = Cp, and c = 0. This gives power flows in the form $P=m˙CpT$ where $m˙$ is set by a pump. The power flow between the fluid flows and wall depend upon flow configuration. For parallel-flow plate heat exchangers, Eq. (16) takes the form of $m˙=0, a=hbAs,b$, b = 1, and c = −1, which gives
$P=hbAs,b(Tb−Tw)$
(17)
where hb is the heat transfer coefficient for fluid flow B and As,b is the heat transfer surface area for fluid flow B. Similarly, for heat transfer from the wall into the exit temperature of fluid A, a = haAs,a, b = 1, and c = −1, which gives
$P=haAs,a(Tw−Ta)$
(18)
where ha is the heat transfer coefficient for fluid flow A and As,a is the heat transfer surface area for fluid flow A. In counter-flow heat exchangers, c = −α in Eq. (16) for fluid flow A and B, which gives power flows in the form of
$P=hbAs,b(Tb−αTw)P=haAs,a(Tw−αTa)$
(19)

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.

Fig. 2
Fig. 2
Close modal
Fig. 3
Fig. 3
Close modal
The temperature dynamics are given by Eq. (13). For the walls C = mCp, resulting in the wall dynamics given as
$mCpT˙w=hbAb(Tb−αTw)−haAa(Tw−αTa)$
(20)

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].

The exit temperature capacitances are calculated as
$C=AcLρCp$
(21)
where Ac is the cross-sectional area of the fluid passage, and L is the length of the heat exchanger fluid flow passage. Cross-sectional area of the heat exchanger is calculated using the hydraulic diameter appropriate for the heat exchanger geometry [16]. Let Ca = Ac,aLaρaCp,a and Cb = Ac,bLbρbCp,b, and then the exit temperature dynamics are given as
$CaT˙a=haAs,a(Tw−Ta)+m˙1,aCp,aT1,a−m˙2,aCp,aTaCbT˙b=m˙1,bCp,bT1,b−hbAs,b(Tb−Tw)−m˙2,bCp,bTb$
(22)

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.

Fig. 4
Fig. 4
Close modal
Fig. 5
Fig. 5
Close modal
Table 2

Liquid–liquid heat exchanger graph parameters

ParameterVariableValueUnits
Wall capacitanceCw900J/K
Fluid A
CapacitanceCa0.0145J/K
Heat transfer areaAs,a0.2015m2
Heat transfer coefficientha10,500W/(m2 K)
Specific heatCp3500J/(kg K)
Fluid B
CapacitanceCb0.0145J/K
Heat transfer areaAs,b0.2015m2
Heat transfer coefficienthb7500W/(m2 K)
Specific heatCp3500J/(kg K)
ParameterVariableValueUnits
Wall capacitanceCw900J/K
Fluid A
CapacitanceCa0.0145J/K
Heat transfer areaAs,a0.2015m2
Heat transfer coefficientha10,500W/(m2 K)
Specific heatCp3500J/(kg K)
Fluid B
CapacitanceCb0.0145J/K
Heat transfer areaAs,b0.2015m2
Heat transfer coefficienthb7500W/(m2 K)
Specific heatCp3500J/(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.

Fig. 6
Fig. 6
Close modal

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 $P=m˙CpT$. 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.

The temperature dynamics are given by Eq. (13), where C = mCp for the cold plate, resulting in the temperature dynamic for the wall
$mCpT˙w=q−hAs(Tw−T)$
(23)
The exit temperature dynamic of the fluid in the cold plate is similar to the heat exchanger exit temperature dynamics where the cold plate capacitance C = AcLρCp, which gives the exit temperature dynamic
$AcLρCpT˙=hAs(Tw−T)+m˙1CpT1−m˙2CpT$
(24)

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.

Fig. 7
Fig. 7
Close modal
Fig. 8
Fig. 8
Close modal
Table 3

Cold plate graph parameters

ParameterVariableValueUnits
Fluid capacitanceCT93.6J/K
Wall capacitanceCw777J/K
Heat transfer areaAs6.72 × 10−3m2
Heat transfer coefficienth8500W/(m2 K)
Fluid specific heatCp3500J/(kg K)
ParameterVariableValueUnits
Fluid capacitanceCT93.6J/K
Wall capacitanceCw777J/K
Heat transfer areaAs6.72 × 10−3m2
Heat transfer coefficienth8500W/(m2 K)
Fluid specific heatCp3500J/(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 $m˙1$ 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 $m˙2$. 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.

Fig. 9
Fig. 9
Close modal

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 $P=m˙1CpT1$. The power flow out of the tank is similar, with a = 0, b = Cp, and c = 0, giving $P=m˙2CpT$. 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.

The temperature dynamics are given by Eq. (13), where thermal capacitance of the tank is C = mCp, giving the tank temperature dynamic as
$mCpT˙=m˙1CpT1−m˙2CpT−hAs(T−Tamb)$
(25)

This equation forms the fluid tank graph model (8) which is provided as Eq. (A3) in the Appendix.

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.

Fig. 10
Fig. 10
Close modal
Table 4

Liquid tank graph parameters

ParameterVariableValueUnits
Tank capacitanceCT3800J/K
Heat transfer areaAs51.0 × 10−3m2
Heat transfer coefficienth15W/(m2 K)
Fluid specific heatCp3500J/(kg K)
ParameterVariableValueUnits
Tank capacitanceCT3800J/K
Heat transfer areaAs51.0 × 10−3m2
Heat transfer coefficienth15W/(m2 K)
Fluid specific heatCp3500J/(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.

Unlike thermal component graphs, electrical energy is not stored within most electrical components, except batteries, and super-capacitors, so Ci in Eq. (1) does not represent a storage capacitance for each electrical graph vertex. Instead, Ci is selected to represent the transient behavior of the modeled component. For generators, Ci is equivalent to the direct axis short circuit transient time constant, which can be found on data sheets or calculated using methods from Ref. [18]. The Ci for a voltage bus is selected to be at least one order of magnitude smaller than the generator due to the fact that bus transients are significantly faster than generator transients. The general expression for the voltage time rate of change within electrical system components is defined as
$CiV˙=∑Pin−∑Pout$
(26)
where each P denotes the amount of power flowing in or out of a component. Sections 4.14.2 detail the selection of Ci and form of Eq. (26) for individual electrical system components. Edge power flows are given by
$P=(a+u)(bx1+cx2)$
(27)

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.

The generator graph model is shown in Fig. 11. This includes a vertex Vgen representing the voltage of the generator, a single edge oriented into the vertex that represents the power input from the shaft, and several edges oriented out of Vgen toward sink vertices. The single source vertex ω is considered a disturbance to the generator graph and would be specified by another graph model for the engine. Sink vertices are individual voltage buses or loads directly attached to the generator. Edge power flow into Vgen from the source vertex is given by Eq. (27) with a = 0, b = α, and c = β, giving
$P=u(αω+βVgen)$
(28)

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.

Fig. 11
Fig. 11
Close modal
Table 5

Electric system graph parameters

ParameterVariableValueUnits
Generator capacitanceCg0.15s
Bus 1 capacitanceCb,10.01s
Bus 2 capacitanceCb,20.01s
Linear coefficients
Shaft speedα12000V s
Generator voltageβ1−8900None
Generator voltageα2650None
Bus 1 voltageβ2−504None
Generator voltageα3420None
Bus 2 voltageβ3−820None
Efficiencyη0.95None
Constant resistanceR5Ω
ParameterVariableValueUnits
Generator capacitanceCg0.15s
Bus 1 capacitanceCb,10.01s
Bus 2 capacitanceCb,20.01s
Linear coefficients
Shaft speedα12000V s
Generator voltageβ1−8900None
Generator voltageα2650None
Bus 1 voltageβ2−504None
Generator voltageα3420None
Bus 2 voltageβ3−820None
Efficiencyη0.95None
Constant resistanceR5Ω

### 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.

The electrical bus graph model is shown in Fig. 12. This has a single vertex Vbus that is connected via a single edge to the source vertex Vgen. This assumes that the bus is connected directly to the generator. While this is the most common architecture, two buses may be connected together in which case the generator source vertex could be replaced with another bus voltage vertex. Sink vertices and the edges connecting them to Vbus represent individual loads on the bus which are discussed in Sec. 4.3. The edge between the source vertex and Vbus represents a transformer unit, and its power flow is defined as
$P=u(αVgen+βVbus)$
(29)
where u is a control input and
$αβ∝N1N2$
(30)

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.

Fig. 12
Fig. 12
Close modal

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.

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.

Fig. 13
Fig. 13
Close modal
For constant power loads, the sink vertex p sets the power for the electrical load and thermal load. The electrical edge power flow is given by Eq. (27), with u = 1, a = 0, b = η, and c = 0 which gives
$P=ηp$
(31)
where η is the efficiency of the load. The heat generated as a byproduct of inefficiencies is represented by the dashed line in Fig. 13 between vertex p and T. The edge terminates in a sink vertex T, which represents the temperature of the constant power electrical load. While p was a sink state in the electrical domain, it is represented as a source vertex in the thermal domain such that the edge parameters are a = 1, b = (1 − η), and c = 0 which gives
$P=(1−η)p$
(32)
To maintain the structure of Eq. (27), the power flow for constant current loads set the input u = I, and coefficients are defined as a = 0, b = η, and c = 0, such that the constant current power flow is given by
$P=uηVbus=ηVbusI$
(33)

where η is the electrical load efficiency.

The heat generated as a byproduct of the inefficiency is represented by the dashed line in Fig. 13 between vertex I and T. The edge terminates in a sink vertex T which represents the temperature of the constant current electrical load. This vertex could be the source vertex in a thermal component graph. While I was a sink state for the electrical domain, it is represented as a source vertex for the thermal domain such that the edge parameters are u = 1, a = 0, b = (1 − η), and c = 0 which gives
$P=(1−η)I$
(34)
Similar to current loads, to maintain the structure of Eq. (27) the input is set u = Vbus, and coefficients can be defined as a = 0, b = 1/R, and c = 0, giving
$P=uVbusR=Vbus2R$
(35)

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.14.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.

Fig. 14
Fig. 14
Close modal
Fig. 15
Fig. 15
Close modal

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.

Fig. 16
Fig. 16
Close modal
Fig. 17
Fig. 17
Close modal

## 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.

Recall that each thermal graph model in Sec. 3 was based upon the conservation of energy. For rotating machinery, angular momentum is an additional conserved quantity that must be considered for graph vertices that represent rotational shaft speeds. Angular momentum L for a shaft is defined as
$L=Iω$
(36)
where I is the moment of inertia and ω is the angular velocity. Taking the time derivative of the angular momentum yields
$dLdt=d(Iω)dt=dIdtω+dωdtI$
(37)
Assuming that the moment of inertia is constant, Eq. (37) can be reduced to
$dLdt=dωdtI=ω˙I$
(38)

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 $τ=Iω˙$ where τ is the net torque acting on the body.

Given Nc components producing positive or negative work on a shaft, the net power is expressed as
$P=ω∑iNcτi=∑iNcm˙iCp,iΔTi$
(39)
where $m˙i$ is the mass flow rate through the turbomachinery component, Cp,i is the specific heat of the working fluid, and ΔTi is the temperature difference from inlet to outlet of each turbomachinery component. For a single component, the torque can be written as
$τ=m˙CpΔTω$
(40)
and the shaft dynamics can be found by combining Eqs. (38) and (40) to give the governing dynamics for a shaft
$dωdt=1Iω∑iNcτi=1Iω∑iNcm˙iCp,iΔTi$
(41)

The dynamics for temperature vertices in the turbomachinery graphs follow Eq. (13) and the edges follow the power flow equation (16).

### 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.

Fig. 18
Fig. 18
Close modal

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.

Fig. 19
Fig. 19
Close modal
Recall (39) which requires a temperature difference in order to calculate the power input or output from the shaft. Each compressor and turbine in the graph model is represented with an inlet and outlet temperature vertex with edges that are oriented such that the inlet conditions are directed into the shaft speed vertex and the outlet conditions are out of the shaft speed vertex. Therefore, the power equation for edges connected to the shaft speed can be written as
$P=m˙CpΔT=m˙CpTin−m˙CpTout$
(42)

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 $P=m˙CpT$. 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.

The red and purple power equation in Fig. 19 is used from inlet to outlet conditions for each compressor and turbine. This linear power flow equation is of the form
$P=α(Tin−Tout)+βm˙+γ$
(43)

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.

Two mass flow rates are used for the ACM configuration in Fig. 18 and the graph in Fig. 19. Since the power turbine is isolated from the air cycle loop
$m˙PT=f(u)$
(44)
where u is a valve opening on the bleed air supply line, which will affect the power flow along e1 and effectively change the shaft speed. The mass flow rate in the air cycle loop $m˙ACM$ is not actuated and is only affected by the shaft speed which often is fitted with a quadratic relationship identified from simulation data
$m˙ACM=f(1,ω,ω2)$
(45)

Both mass flow rate relationships are provided as Eq. (A5) in the Appendix.

For heat absorbed and rejected by the ACM, the source and sink vertices are often wall temperatures for heat exchangers. Heat enters the ACM through a heat exchanger between the expansion turbine and compressor, which effectively increases the temperature at the inlet to the compressor. The heat transfer between the fluid flow and a heat exchanger wall follows the form of Eq. (16) where a = hAs, b = 1, and c = −1, which gives:
$P=hAs(Tw−TK,in)$
(46)
where h is the heat transfer coefficient for the fluid and As is the heat transfer surface area between the fluid and wall. Heat is rejected from the ACM using a heat exchanger between the compressor outlet and expansion turbine inlet; therefore, in the graph power exits the ACM from the TET,in vertex. The heat transfer between the fluid flow and a heat exchanger wall follows the form of Eq. (16) where a = hAs, b = 1, and c = −1, which gives
$P=hAs(TET,in−Tw)$
(47)
The temperature dynamics for each vertex in Fig. 19 simply follow from Eq. (13), where C = MCp = ρVCp and the air volume V is determined from the geometry for the ACM. The dynamics for each temperature are expressed as
$ρVCpT˙=∑Pin−∑Pout$
(48)
Shaft dynamics use Eq. (41) as a basis, but are linearized about a nominal shaft speed ω0 to avoid nonlinear behavior in the graph model
$ω˙=[1Iω0−1Iω02(ω−ω0)](∑Pin−∑Pout)$
(49)
Due to $ω02$ being very large relative to ω0, the linear term (ω − ω0) will be dominated by the constant term and can thus be removed with little affect on the dynamics. This simplifies the expression to the format of Eq. (1)
$Iω0ω˙=(∑Pin−∑Pout)$
(50)

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 [2024], 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].

Fig. 20
Fig. 20
Close modal
Table 1

Computational difference between graph and nonlinear model for 25 simulations

GraphNonlinear model
tavg (s)St. Dev. (s)tavg (s)St. Dev. (s)
25.95.2627719.9
GraphNonlinear model
tavg (s)St. Dev. (s)tavg (s)St. Dev. (s)
25.95.2627719.9
The first metric of comparison is the coefficient of performance (COP), which is defined as
$COP=Qm˙bleedCp(Tbleed−Tamb)$
(51)

where Q is the heat load into the ACM, $m˙bleed$ 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.

Fig. 21
Fig. 21
Close modal
Fig. 22
Fig. 22
Close modal

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.

Table 6

Air cycle system graph parameters

ParameterVariableValueUnits
Power turbine
Inlet capacitanceCPT,in1000J/K
Outlet capacitanceCPT,out1000J/K
Temperature coefficientα254J/(K s)
Mass flow rate coefficientβ−11,250J/(kg s)
Constant coefficientγ−43,400J/s
Fluid specific heatCp1073J/(kg K)
Bleed air temperatureTbleed700K
Compressor
Inlet capacitanceCK,in100J/K
Outlet capacitanceCK,out100J/K
Temperature coefficientα1000J/(K s)
Mass flow rate coefficientβ3.115 × 106J/(kg s)
Constant coefficientγ0J/s
Fluid specific heatCp1030J/(kg K)
Expansion turbine
Inlet capacitanceCET,in100J/K
Outlet capacitanceCET,out100J/K
Temperature coefficientα1000J/(K s)
Mass flow rate coefficientβ−5.675 × 105J/(kg s)
Constant coefficientγ0J/s
Fluid specific heatCp1008J/(kg K)
Fuel-compressor heat exchanger
Wall capacitanceCHX,f964J/K
Outlet capacitanceCout,f100J/K
Heat transfer coefficient K sidehk17,500W/(m2 K)
Heat transfer area K sideAk1.86m2
Heat transfer coefficient fuel sidehf10,000W/(m2 K)
Heat transfer area fuel sideAf1.86m2
Fuel specific heatCp2000J/(kg K)
Air-expansion turbine heat exchanger
Wall capacitanceCHX,f9500J/K
Outlet temperature capacitanceCout,f100J/K
Heat transfer coefficient ET sidehET17,500W/(m2 K)
Heat transfer area ET sideAET11.8m2
Heat transfer coefficient air sideha17,500W/(m2 K)
Heat transfer area air sideAa11.8m2
Air specific heatCp1008J/(kg K)
ParameterVariableValueUnits
Power turbine
Inlet capacitanceCPT,in1000J/K
Outlet capacitanceCPT,out1000J/K
Temperature coefficientα254J/(K s)
Mass flow rate coefficientβ−11,250J/(kg s)
Constant coefficientγ−43,400J/s
Fluid specific heatCp1073J/(kg K)
Bleed air temperatureTbleed700K
Compressor
Inlet capacitanceCK,in100J/K
Outlet capacitanceCK,out100J/K
Temperature coefficientα1000J/(K s)
Mass flow rate coefficientβ3.115 × 106J/(kg s)
Constant coefficientγ0J/s
Fluid specific heatCp1030J/(kg K)
Expansion turbine
Inlet capacitanceCET,in100J/K
Outlet capacitanceCET,out100J/K
Temperature coefficientα1000J/(K s)
Mass flow rate coefficientβ−5.675 × 105J/(kg s)
Constant coefficientγ0J/s
Fluid specific heatCp1008J/(kg K)
Fuel-compressor heat exchanger
Wall capacitanceCHX,f964J/K
Outlet capacitanceCout,f100J/K
Heat transfer coefficient K sidehk17,500W/(m2 K)
Heat transfer area K sideAk1.86m2
Heat transfer coefficient fuel sidehf10,000W/(m2 K)
Heat transfer area fuel sideAf1.86m2
Fuel specific heatCp2000J/(kg K)
Air-expansion turbine heat exchanger
Wall capacitanceCHX,f9500J/K
Outlet temperature capacitanceCout,f100J/K
Heat transfer coefficient ET sidehET17,500W/(m2 K)
Heat transfer area ET sideAET11.8m2
Heat transfer coefficient air sideha17,500W/(m2 K)
Heat transfer area air sideAa11.8m2
Air specific heatCp1008J/(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.

Fig. 23
Fig. 23
Close modal

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.

Fig. 24
Fig. 24
Close modal
Fig. 25
Fig. 25
Close modal
Fig. 26
Fig. 26
Close modal

## 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.

Fig. 27
Fig. 27
Close modal

Component graphs from Secs. 35 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.

Fig. 28
Fig. 28
Close modal

### 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.

Fig. 29
Fig. 29
Close modal
Fig. 30
Fig. 30
Close modal

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. 35, 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

$m˙$ =

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

ω =

### Appendix

Graph definition for a liquid–liquid heat exchanger
$[Ca000Cw000Cb]︷C[T˙aT˙wT˙b]=−[−10101−1000101]︷M¯[haAs,a(Tw−Ta)hbAs,b(Tb−Tw)m˙2,aCpTam˙2,bCpTb]︷P+[100001]︷D[m˙1,aCp,aT1,am˙1,bCp,bT1,b]︷Pin$
(A1)
Graph definition for a cold plate
$[Cw00CT]︷C[T˙wT˙]=−[10−11]︷M¯[hAs(Tw−T)m˙1CpT]︷P+[1001]︷D[qm˙1CpT1]︷Pin$
(A2)
Graph definition for a liquid tank
$[CT]︷C[T˙]=−[11]︷M¯[hAs(T−Tamb)m˙2CpT]︷P+[1]︷D[m˙1CpT1]︷Pin$
(A3)
Graph definition for the electrical system in Sec. 4, with subscripts g = gen and b = bus
$[Cg000Cb,1000Cb,2]︷C[V˙gV˙b,1V˙b,2]=−[11000000−101110000−1000111]︷M¯[α2Vg+β2Vb,1α3Vg+β3Vb,2ηPVb,12/RηVb,1ηPVb,22/RηVb,2]︷P+[100]︷D[α1ω+β1Vg]︷Pin$
(A4)
Air cycle machine mass flow rate relationships
$m˙PT=0.00537u−1.302m˙ACM=−2.46×10−7ω2+2.79×10−3ω−6.57$
(A5)
Air cycle machine graph model
$C=diag([CPT,in CPT,out Cω CK,in CK,out CET,in CET,out CHX,f Cout,f CHX,e Cout,e])$
(A6)
$x˙=[T˙PT,in T˙PT,out ω˙ T˙K,in T˙K,out T˙ET,in T˙ET,out T˙HX,f T˙out,f T˙HX,e T˙out,e]T$
(A7)
$M¯=[110000000000000000−10−10000000000001000−11−11−11000000000000001000100−1−10000000000−100−1100000000000000100−1100010000000000−100−110000000000000000001−1000000000000000001000100000000000000−1100000000000000000−1001]$
(A8)
$P=[α(TPT,in−TPT,out)+βm˙PT+γm˙PTCpTPT,inm˙PTCpTPT,outm˙ACMCpTK,inm˙ACMCpTK,outm˙ACMCpTET,inm˙ACMCpTET,outα(TK,in−TK,out)+βm˙ACM+γm˙ACMCpTK,outα(TET,in−TET,out)+βm˙ACM+γm˙ACMCpTET,outhkAHX,f(Tw−TK,in)hfAHX,f(Tout−Tw)hETAHX,e(TET,in−Tw)haAHX,e(Tw−Tout)m˙fuelCpToutm˙bypassCpTout]$
(A9)
$D=[100000000000000000010000000000001]T$
(A10)
$Pin=[m˙PTCpTbleedm˙fuelCpTfuelm˙bypassCpTbypass]$
(A11)

## References

1.
Liptak
,
B. G.
,
1995
,
Process Control
, 3rd ed.,
Chilton
2.
Guerrero
,
J. M.
,
Chandorkar
,
M.
,
Lee
,
T. L.
, and
Loh
,
P. C.
,
2013
, “
Advanced Control Architectures for Intelligent Microgrids—Part I: Decentralized and Hierarchical Control
,”
IEEE Trans. Ind. Electron.
,
60
(
4
), pp.
1254
1262
.
3.
Negenborn
,
R. R.
,
Sahin
,
A.
,
Lukszo
,
Z.
,
Schutter
,
B. D.
, and
Morari
,
M.
,
2009
, “
A Non-Iterative Cascaded Predictive Control Approach for Control of Irrigation Canals
,”
IEEE International Conference on Systems, Man and Cybernetics
(
ICSMC
), San Antonio, TX, Oct. 11–14, pp.
3552
3557
.
4.
Christofides
,
P. D.
,
Scattolini
,
R.
,
Muñoz de la Peña
,
D.
, and
Liu
,
J.
,
2013
, “
Distributed Model Predictive Control: A Tutorial Review and Future Research Directions
,”
Comput. Chem. Eng.
,
51
, pp.
21
41
.
5.
Deppen
,
T. O.
,
2013
, “Optimal Energy Use in Mobile Applications With Storage,”
Ph.D. thesis
, University of Illinois, Urbana, IL.https://www.ideals.illinois.edu/handle/2142/44290
6.
Jain
,
N.
,
Koeln
,
J. P.
,
Sundaram
,
S.
, and
Alleyne
,
A. G.
,
2014
, “
Partially Decentralized Control of Large-Scale Variable-Refrigerant-Flow Systems in Buildings
,”
J. Process Control
,
24
(
6
), pp.
798
819
.
7.
Rasmussen
,
B. P.
, and
Alleyne
,
A. G.
,
2006
, “Dynamic Modeling and Advanced Control of Air Conditioning and Refrigeration Systems,”
Ph.D. dissertation
, University of Illinois at Urbana-Champaign, Champaign, IL.https://www.ideals.illinois.edu/handle/2142/12355
8.
McCarthy
,
P.
,
Niedbalski
,
N.
,
McCarthy
,
K.
,
Walters
,
E.
,
Cory
,
J.
, and
Patnaik
,
S.
,
2016
, “
A First Principles Based Approach for Dynamic Modeling of Turbomachinery
,”
SAE Int. J. Aerosp.
,
9
(
1
), pp.
45
61
.
9.
Williams, M., Sridharan. S., Banerjee, S., Mak, C., Pauga, C., Krein, P., Alleyne, A., Jacobi, A., and D'Urso, S.,
2015
, “PowerFlow: A Toolbox for Modeling and Simulation of Aircraft Systems,”
SAE
Paper No. 2015-01-2417.
10.
Garcia
,
C. E.
,
Prett
,
D. M.
, and
Morari
,
M.
,
1989
, “
Model Predictive Control: Theory and Practice—A Survey
,”
Automatica
,
25
(
3
), pp.
335
348
.
11.
Moore
,
K. L.
,
Vincent
,
T. L.
,
Lashhab
,
F.
, and
Liu
,
C.
,
2011
, “
Dynamic Consensus Networks With Application to the Analysis of Building Thermal Processes
,”
IFAC Proc. Vol.
,
44
(
1
), pp.
3078
3083
.
12.
Preisig
,
H. A.
,
2009
, “
A Graph-Theory-Based Approach to the Analysis of Large-Scale Plants
,”
Comput. Chem. Eng.
,
33
(
3
), pp.
598
604
.
13.
Koeln
,
J. P.
,
Williams
,
M. A.
, and
Alleyne
,
A. G.
,
2015
, “
Hierarchical Control of Multi-Domain Power Flow in Mobile Systems—Part I: Framework Development and Demonstration
,”
ASME
Paper No. DSCC2015-9908.
14.
Williams
,
M. A.
,
Koeln
,
J. P.
, and
Alleyne
,
A. G.
,
2015
, “
Hierarchical Control of Multi-Domain Power Flow in Mobile Systems—Part II: Aircraft Application
,”
ASME
Paper No. DSCC2015-9904.
15.
West
,
D. B.
,
2001
,
Introduction to Graph Theory
,
Prentice Hall
16.
Incropera
,
F. P.
,
DeWitt
,
D. P.
,
Bergman
,
T. L.
, and
Lavine
,
A. S.
,
2006
,
Fundamentals of Heat and Mass Transfer
, 6th ed.,
Wiley
, Hoboken, NJ.
17.
Eremia
,
M.
, and
Bulac
,
C.
,
2013
,
Handbook of Electrical Power System Dynamics: Modeling, Stability, and Control
,
IEEE Press
,
Piscataway, NJ
.
18.
Sauer
,
P.
, and
Pai
,
M. A.
,
1998
,
Power System Dynamics and Stability
,
Prentice Hall
19.
Cao, Y.,
Williams
,
M. A.
,
Kearbey
,
B. J.
,
Smith
,
A. T.
,
Krein
,
P. T.
, and
Alleyne
,
A. G.
,
2016
, “
20x-Real Time Modeling and Simulation of More Electric Aircraft Thermally Integrated Electrical Power Systems
,” International Conference on Electrical Systems for Aircraft, Railway, Ship Propulsion and Road Vehicles & International Transportation Electrification Conference (
ESARS-ITEC
), Toulouse, France, Nov. 2–4, pp. 1–6.
20.
Ashford
,
R.
,
2004
, “
Verification and Validation of the F/A-22 Raptor Environmental Control System/Thermal Management System Software
,”
SAE
Paper No. 2004-01-2573
.
21.
Santos
,
A. P. P.
,
,
C. R.
, and
Zaparoli
,
E. L.
,
2014
, “
A Thermodynamic Study of Air Cycle Machine for Aeronautical Applications
,”
Int. J. Thermodyn.
,
17
(
3
), pp.
117
126
.
22.
Vargas
,
J. V. C.
, and
Bejan
,
A.
,
2001
, “
Integrative Thermodynamic Optimization of the Environmental Control System of an Aircraft
,”
Int. J. Heat Mass Transfer
,
44
(
20
), pp.
3907
3917
.
23.
Conceição
,
S. T.
,
Zaparoli
,
E. L.
, and
Turcio
,
W. H. L.
,
2007
, “
Thermodynamic Study of Aircraft Air Conditioning Air Cycle Machine: 3-Wheel × 4-Wheel
,”
SAE
Paper No. 2007-01-2579.
24.
Tu
,
Y.
, and
Lin
,
G. P.
,
2011
, “
Dynamic Simulation of Aircraft Environmental Control System Based on Flowmaster
,”
J. Aircr.
,
48
(
6
), pp.
2031
2041
.
25.
Zhao
,
H.
,
Hou
,
Y.
,
Zhu
,
Y.
,
Chen
,
L.
, and
Chen
,
S.
,
2009
, “
Experimental Study on the Performance of an Aircraft Environmental Control System
,”
Appl. Therm. Eng.
,
29
(
16
), pp.
3284
3288
.
26.
Childs
,
T.
,
Jones
,
A. B.
, and
Chen
,
R.
,
2015
, “
Development of a Full Scale Experimental and Simulation Tool for Environmental Control System Optimisation and Fault Detection
,”
AIAA
Paper No. 2015-1196.
27.
Matullch
,
D. S.
,
1989
, “
High-Temperature Bootstrap Compared With F15 Growth Air Cycle Air Conditioning System
,”
SAE
Paper No. 891436.