An Energy Graph-Based Approach to Fault Diagnosis of a Transcritical CO 2 Heat Pump

: The objective of this paper is to describe an energy-based approach to visualize, identify, and monitor faults that may occur in a water-to-water transcritical CO 2 heat pump system. A representation using energy attributes allows the abstraction of all physical phenomena present during operation into a compact and easily interpretable form. The use of a linear graph representation, with heat pump components represented as nodes and energy interactions as links, is investigated. Node signature matrices are used to present the energy information in a compact mathematical form. The resulting node signature matrix is referred to as an attributed graph and is populated in such a way as to retain the structural information, i.e., where the attribute points to in the physical system. To generate the energy and exergy information for the compilation of the attributed graphs, a descriptive thermal–fluid model of the heat pump system is developed. The thermal–fluid model is based on the specifications of and validated to the actual behavioral characteristics of a physical transcritical CO 2 heat pump test facility. As a first step to graph-matching, cost matrices are generated to represent a characteristic residual between a normal system node signature matrix and a faulty system node signature matrix. The variation in the eigenvalues and eigenvectors of the characteristic cost matrices from normal conditions to a fault condition was used for fault characterization. Three faults, namely refrigerant leakage, compressor failure and gas cooler fouling, were considered. The paper only aims to introduce an approach, with the scope limited to illustration at one operating point and considers only three relatively large faults. The results of the proposed method show promise and warrant further work to evaluate its sensitivity and robustness for small faults.


Introduction
A heat pump system extracts heat from natural sources such as air, surface water or solar energy and use it in a thermodynamic cycle. Heat pump technology is widely used in many applications that warrant cooling, heating or both. What distinguishes the considered CO 2 heat pump in this paper from chillers and other heat pumps is the fact that it is a water-to-water system and produces both hot and cold water as useful products.
According to Zogg [1], when considering the ecological and economic impact of refrigerating machines, air-conditioners, and heat pumps, the main issue is to maximize the efficiency to minimize energy consumption. This is motivated due to the global drive towards mitigating CO 2 emissions. called graph-matching, and system knowledge in terms of exergy and energy flow. This methodology is evaluated on a validated heat pump system model. This paper is outlined as follows: Section 2 discusses the heat pump test facility used to obtain practical system data. In Section 3 a system-level model of the heat pump system is derived and discussed. The model is simulated in the EES environment and validated with practical data. Section 4 introduces the energy graph-based characterization approach using energy and exergy information and then Section 5 follows with the energy graph-based fault detection and isolation methodology. The results are discussed in Section 6 followed by conclusions in Section 7. Figure 1 illustrates the water-to-water transcritical CO 2 (R744) heat pump system with its main cyclic points and components shown. The typical thermodynamic behavior of the transcritical refrigerant loop is shown in Figure 2 on a T-s diagram. Figure 1. Illustration of the transcritical heat pump system with its cyclic points and components [22]. The transcritical R744 cycle extracts heat from the water stream that flows through the evaporator component. The R744 that flows through the evaporator is therefore heated and the corresponding water stream is cooled. The R744 that enters the evaporator from the expansion valve (EV) is in a two-phase state. The R744 is heated in the evaporator until it reaches a saturated vapor state-cycle point (1). Thereafter it is heated further until it reaches a superheated vapor state-cycle point (2). Superheating is performed to increase the efficiency of the transcritical cycle and to avoid compressing a two-phase stream. Compression of a two-phase mixture that contains liquid droplets will result in compressor damage. The compressor component receives the superheated gas that leaves the evaporator and compresses the gas to a higher pressure and temperature condition as indicated in Figure 2 between cycle point (2) and (3). The R744 that is discharged from the compressor is in the supercritical gaseous state. The gas cooler component receives the high-temperature supercritical gas and transfers heat from the supercritical gas to the water stream. The R744 that flows through the gas cooler is cooled and the corresponding water stream is heated. In a transcritical heat pump, the heat exchanger component that rejects heat from the refrigerant is called a gas cooler since the R744 that flows through it remains in a supercritical state during the heat rejection process. The process of condensation, as found in a typical Freon heat pump cycle, does not take place.

System Description
The R744 in the gas cooler component may be in either the supercritical gaseous or in the sub-cooled liquid state-cycle point (23). For proper system operation, the EV should receive R744 that is in the sub-cooled liquid state. Within the expansion valve, R744 undergoes an isenthalpic throttling process-cycle point (23) to (24). The throttling process causes the R744 to undergo a significant pressure drop. The temperature of the R744 also decreases with the reduction in its pressure. The throttling process effectively lowers the temperature so that evaporation of the R744 may again take place, via heat transfer from the water stream to the R744 in the evaporator component-cycle point (24) to (2). The actual heat pump test facility viewed from the front and the back is illustrated in Figure 3. Table 1 lists the corresponding components found in the heat pump system with their respective functions.   13 Evaporator centrifugal pump Circulates the water to be cooled through the evaporator 14 Gas cooler centrifugal pump Circulates the water to be heated through the gas cooler

System Model
The system model is developed in the Engineering Equation Solver (EES ) environment [23]. The simulation methodology follows a systems-level approach. Semi-theoretical models are used during the implementation of conservation laws to solve for the required parameters which quantify the behavioral properties of the components in the system. The laws of conservation of mass, momentum and energy (for steady incompressible flow) are each applied to every component in the system. The use of the incompressible form of the conservation of momentum law is valid due to the Mach number of the refrigerant never exceeding 0.1 throughout the system within the range of simulated system conditions. At flow Mach numbers below 0.3, the compressible and incompressible form of the law of conservation of momentum both yield the same thermo-physical property results. The law of conservation of mass for steady incompressible flow through an individual component's control volume (CV) is expressed as [24]: where the first term is the summation of all the mass flow streams entering the CV and the second term is the summation of all the mass flow streams flowing out of the CV. The conservation of mass for every system component may thus be simplified to: The law of conservation of momentum for steady (incompressible low Mach number) flow through an individual component's CV is expressed as [24]: where the first and second terms are the total and hydrostatic pressures of fluids entering the CV, the third and fourth terms are the total-and hydrostatic pressures of fluids leaving the CV and the fifth term is the pressure drop which occurs as the various interacting fluids flow through the CV. Equation (3) may be simplified by only considering the interacting fluids for each component. The conservation of momentum for the refrigerant side of the evaporator, gas cooler, EV and compressor components is given by: where the subscript, g, refers to the R744 refrigerant. The conservation of momentum for the water side of the evaporator and gas cooler components is given by: where the subscript, w, refers to the water stream. The pressure ratio of the compressor is used to find the discharge pressure of the refrigerant leaving the compressor. The pressure ratio for a compressor is defined as [22]: Finally, the law of conservation of energy (First Law of Thermodynamics) for a steady flow process of an individual component's CV is given by [24]: where the first term is the net heat transfer into the CV, the second term is the net work absorbed by the CV, the third term is the total flow energy entering the CV, the fourth term is the gravitational potential energy of the flow(s) entering the system, the fifth term represents the flow energy leaving the CV, and the sixth term represents the gravitational potential energy of the flow(s) departing the CV. The conservation of energy for each component was compiled under the assumption that the elevation head difference throughout the system is negligible. The elevation head difference for the experimental system is 0.5 m. The conservation of energy for the evaporator and gas cooler component is given by: The conservation of energy for the compressor component is given by: The expansion valve component's throttling process is assumed to be isenthalpic. Under the isenthalpic assumption, the conservation of energy for the expansion valve simplifies to: In addition to the traditional conservation laws, an exergy balance was also applied to solve for its corresponding internal irreversibility rate within each system component. An exergy balance for an individual component's CV has the form [25]: where the first term is the exergy transfer to the CV due to heat transfer, the second term is the net work absorbed by the CV, the third term is the exergy entering the CV, the fourth term is exergy leaving the CV and the fifth term is the rate at which exergy is irreversibly destroyed inside of the CV. The evaporator and gas cooler components do not have associated work terms in their respective exergy balance equations. The exergy balance for the evaporator is given by: The work term is omitted in (11) and (12) since no work is done on the working fluid in these components, only heat transfer takes place. The exergy balance for the gas cooler component is identical to the exergy balance for the evaporator component, as in (12), with the only difference being that the irreversibility rate's subscript changes to GC. The exergy balance for the compressor component is given by The exergy balance for the expansion valve component is given by: To validate the EES system model, experimental results for expansion valve (EV) openings of 30%, 40%, 50%, 60% and 70% were compared to the simulation results for the same operating conditions. The results for an EV orifice opening of 30% are discussed in detail whereafter the results for the other conditions are merely summarized. The system input parameters for the 30% EV orifice opening condition are summarized in Table 2.  Figure 4 shows the experimental and simulation results overlaid on a T-s diagram. The pressures are also indicated on the figure for the gas side. The pressure on the water side is in the order of 200 kPa. Similarly, Figure 5 shows the overlaid results and simulation results on a log. P-h diagram. Table 3 shows the key parameter values associated with the 30% EV orifice opening operating point illustrated in Figures 4 and 5 with the deviation between the simulated and test bench parameter values for the operating point.    Table 3 reveals satisfactory validation of the EES system model for the 30% EV orifice opening condition. When considering EV orifice opening results for 40%, 50%, 60% and 70% together with the 30% case, the patterns produced when plotting successive property diagrams for the simulated and measured cases, were similar. As the EV orifice opening is increased, the area encapsulated by the refrigerant cycle property plot decreases in overall size. As the EV orifice opening is increased the percentage deviation between the simulated and test bench results gradually increases as summarized in Table 4. At 60% EV orifice opening the deviation becomes unacceptably high. The reason for the higher deviations at larger EV orifice openings is two-fold; less data was available at these operating points together with the fact that the cycle is not fully developed or thermodynamically stable at the larger EV orifice openings. The EES system model shows satisfactory prediction (if the maximum deviation of 10.7% is viewed as acceptable) of the heat pump's energy characteristics for the 30%, 40% and 50% EV orifice opening operating points under the condition that the enthalpy associated with the expansion valve's throttling process is at a value smaller than the critical enthalpy of the system's refrigerant. The 30% working point is however used for the analysis results as this represents an operating point where the cycle is fully developed and where the simulation accurately represents the test bench. The modeling of the heat pump system should incorporate the interdependent behavior that is a characteristic of the series connected system components. The modeling and characterization of each component in the system and that of the transcritical cycle will be fundamentally based on the control volume approach of thermal-fluid analysis.
The conservation equations and exergy balance equations for system components will be derived under the assumption that the heat pump system and all its components are operating under steady-state conditions. In the range of heat pump operating conditions investigated in this study the flow regime of the refrigerant and the two interacting water streams is always turbulent. The Nusselt number correlations that are used in the simulation to predict convection coefficient values are thus for turbulent flow. The finite term and incompressible flow form of the conservation of momentum equation will be used for the modeling of the heat pump system components. This assumption is valid even for the compressible carbon dioxide refrigerant due to the low Mach number values associated with the carbon dioxide flow. The effect of elevation height difference between system components is assumed to be negligible in the compilation of the conservation of momentum equations for the system's components.
To ensure component interdependency, it will be assumed in the composition of the conservation and component-specific equations that the input values of the devised control volume regions are known. The output thermo-physical properties of the associated CV regions will then be predicted by solving the parameters associated with the conservation and component-specific equations. The component interactions and interdependency will be interlaced in the EES environment by setting the output of each component equal to the input of the component that chronologically succeeds it in the transcritical cycle. Implementation of this approach assures that the operational characteristics of each component in the system are directly dependent on the output of the component that precedes it. Detail on the model derivation can be found in Appendix A.

Energy Graph-Based Characterization
Graphs can be used to represent interactions of large-scale and complex systems in a structured way [26]. Some useful and descriptive attributes can be assigned to a graph which may then be called an attributed graph [27]. This type of graph also allows structural information and attributes to be associated with specific locations in the actual system being considered [28]. This structural property may contribute significantly in terms of fault isolation. Different approaches can be followed when using energy-based attributes. First, the choice of what constitutes a link and what constitutes a node needs to be determined. In this case, a node is defined as a system component e.g., a compressor or turbine, and a link is defined as a connection between two system components. In the following sections this approach will be discussed in more detail from a graph-theoretic point of view and how it relates to the heat pump system considered.

Graph-Theoretic Concepts
According to Rahman [29], a graph G can be defined as a tuple (V, E) consisting of a finite set V of vertices and a finite set E of edges. Each edge can be represented by an unordered pair of vertices. For example, and edge e between two vertices u and v may be represented as (u, v). Balakrishnan and Ranganathan [30] gives the following formal mathematical definition of a graph: "A graph is an and I G is an "incidence" relation that associates with each element of E(G) an unordered pair of elements (same or distinct) of V(G)." The symbol G is thus used to denote the graph as a whole, V(G) represents the vertices of the said graph, E(G) is the graph edges and I G is the connection relation between the vertices and edges of the graph. Figure 6 illustrates an example of a linear graph. In the context of graph theory terminology, the term vertices may be used interchangeably with the term nodes, and the term edges may be used interchangeably with the term links. The author will henceforth refer to nodes and links.

Heat Pump System Graph Representation
An attributed graph is a graph representation that is assigned a set of parameters with numerical values [17]. An attributed graph is useful for system representation since numerous graph-matching operations may be performed on the node signature matrices that are compiled from attributed graphs [29]. A node signature matrix is merely a representation of an attributed graph's attributed parameters in matrix format [31]. The investigated heat pump system has four components. Nodes will be used to represent the four components and links will be used to represent the energy that flows in and out of the components. The approach followed in this paper is similar to the approach outlined in van Schoor & Uren [19], the only difference being that in this work the heat exchanger is assigned a single node with four energy links.
The attributed graph representation of the heat pump system is illustrated in Figure 7. The nodes are indicated with symbol n and the links are indicated with symbol l. The node with the number zero as a subscript is defined as the reference node. The inclusion of a zero node enables the structuring of the graph with its attributes while maintaining the structural information of the graph [19]. Following the approach of van Schoor & Uren [19] the general node signature matrix for the graph in Figure 7 is obtained as: M ag is structured with the first column allocated for the node attributes and the subsequent columns for the link attributes. For the link attributes the first numeral in the subscript refers to the origin of the link and the second refers to the destination. A subscript, "01", indicates that the link is directed from node number zero to node number one. The links in the proposed graph representation thus encode information on the magnitude and the direction of the attributed parameters. In this manner, each row of the matrix represents attributes associated with a specific node. The links that do not appear in Figure 7 are assigned values of zero. Equation (15) may thus be simplified by populating the entries of the matrix with zeros for the positions that represent links that are not present in the graph representation. A link attribute will also be assigned if its subscript notation is opposite to the flow direction assigned. The rule of l mn = −l nm will therefore be used to determine whether to assign a link attribute a positive or negative sign. In this way, the structural information of the linear graph is encoded during matrix compilation. The reduced attributed graph node signature matrix that results for the heat pump system after eliminating redundant links, is given by: The graph representation of the heat pump will use the rates of exergy destruction of the associated components as the node attributes. The link attributes will be the energy flow rates in or out of the associated component to which the links are connected. Table 5 lists the attributes that will be assigned to the graph representation of the heat pump.  From inspection of Table 5, it is evident that every parameter assigned to the graph representation is a variable with the same unit, i.e., kilowatt. The assignment of parameters with the same unit to the node signature matrices ensures that the subsequent calculations are performed with variables that represent the same physical property. The heat pump graph representation with the EES simulation output parameters assigned as attributes is illustrated in Figure 8. The reduced attributed graph node signature matrix that results from Figure 8 is given by: With the reduced attributed graph node signature matrix defined, the next step is to illustrate how this system description can be used for fault detection and isolation.

Energy Graph-Based Fault Detection and Isolation Method
This section introduces the graph-based concepts of graph-matching and furthermore the idea of eigendecomposition as a means to use the attributed graph node signature matrix for fault detection and isolation.

Graph-Matching
Graph-matching refers to the methodologies and the algorithms that compare two or more graphs numerically to determine the degree of similarity between their attributes [17]. In the context of fault detection and diagnosis, a fault signature may be defined as a distinct set of attributes that are expressed on a numerical basis for a fault scenario that uniquely identifies a specific operating condition under which a system is operational [32]. Graph-matching will be used in this work to aid in the compilation of fault signatures for different fault scenario. The approach outlined by Jouili & Tabbone [33] will be implemented to perform the matching of two different graphs.
De la Fuente [31] defines a normal condition as: "A system operating under a condition where none of its defining characteristic properties or associated parameters deviates from the standard system condition." In the context of graphs, the healthy system graph will represent a normal graph with its corresponding normal node signature matrix. Furthermore, de la Fuente [31] defines a fault condition as: "A system operating under a condition where an unpermitted deviation of at least one characteristic property or parameter of the system from the standard condition is present." In the context of graphs, the faulty system graph will represent a graph characteristic of a specific fault with its corresponding fault node signature matrix.
It is desirable to keep a heat pump system running at its normal condition for the longest period possible during operation. Operating at the normal "healthy" system condition ensures efficient and reliable energy conversion and lowers the overall maintenance required to run the system. A fault scenario for a heat pump system is thus likely to impact the efficiency and economics of the system. This section investigates refrigerant leakage, compressor failure and gas cooler fouling.

Cost Matrix Generation
A cost matrix, denoted by C, is a matrix that describes the matching costs between two node signature matrices [33]. The method proposed by Jouili & Tabbone [33] may be used to match the same node signature matrix to itself or a node signature matrix with a different set of element values. The matching of an initial node signature matrix that contains parameter values relating to a healthy operating condition, to itself is performed to generate the reference condition to which other operating conditions may be compared. The element entries of a cost matrix thus encode information on the similarity between two compared node signature matrices.
In the method by Jouili & Tabbone [33], the entries of a cost matrix are calculated using a version of the Heterogeneous Euclidean Overlap Method (HEOM). The entries of a cost matrix C are calculated using the HEOM such that [34]: where i denotes a row entry in the generated cost matrix, and j denotes a column entry in the generated cost matrix. The first entry of a cost matrix, for example, is thus calculated by substituting i and j with 1 such that: The complete cost matrix is compiled from the individual entries as follows: The variable a is the number of the column in the applicable node signature matrix. The delta function, δ(·) is defined as: where M R ia refers to the i-th row of the reference node signature matrix and M F ja refers to the j-th row of the node signature matrix that is compared to the reference node signature matrix. The range is given by: where max a is the largest numerical entry of the a-th column of the reference node signature matrix and min a is the smallest numerical entry in the same column.

Eigenvectors and Eigenvalues
Lay et al. [35] give a good definition of eigenvalues and eigenvectors. Eigenvectors define a set of properties of those numerical values subject to a linear transformation in n-dimensional space that remain unchanged. The eigenvectors of a linear transformation are those vectors that remain on their spanning line through the origin of the n-dimensional space after a linear transformation [35]. Each eigenvector has associated with it an eigenvalue that encodes information on the change in direction and scaling that the eigenvectors undergo during a linear transformation. The eigenvectors associated with a linear transformation may thus change direction on their linear span and increase or decrease in magnitude, but they always start at the origin of the n-dimensional space and point in a direction along the length of their span.
Eigenvectors and their associated eigenvalues are thus excellent properties that may be used to uniquely characterize, summarize, and represent the data encoded in a matrix [19]. This is the motivation for calculation of the eigenvectors and eigenvalues of the cost matrices.

Results and Discussion
This section investigates the use of eigenvectors and eigenvalues of the cost matrices generated by the method by Jouili and Tabbone [33] as a means to fault diagnosis on the heat pump system.

Analysis Assumptions
Since the paper only aims to introduce the energy graph-based approach to fault diagnosis, the scope of analysis is limited to one operating point and considering only 3 large fault conditions. The analysis assumptions are summarized as follows:

Operating Point
For the results reported in the paper the case of 30% EV orifice opening is taken as the normal conditions operating point with the system input parameters and key cycle parameters as summarized in Tables 2 and 3 respectively. A smaller orifice opening degree works well since then the thermal-fluid behaviour of the cycle is fully developed.

Fault Conditions
Three fault conditions will be investigated: working fluid leakage from the heat pump, compressor failure, and fouling inside of the gas cooler on the water side. These fault conditions are seen as the main fault categories in heat pumps [36,37].

Variation in Data
Due to parameter variation and measurement uncertainty, a variation in energy and exergy characteristics can be expected which should be seen as part of normal behavior. For the actual test bench, the measurement uncertainty was less than 1% when calculating the heat transfer on the CO 2 side [38]. In this study, variation in data is not included per se. However, a threshold for FDI representative of typical heat pump energy performance variation is chosen. From [39] the typical COP variation that can be expected is in the region of 4%. The threshold for FDI for this study is therefore set at 5%.

The Normal Case Analyzed
The node signature matrix for the chosen normal case is numerically given by: The M normal matrix will thus represent the healthy fault-free case when progressive fault conditions are considered. The normal conditions cost matrix obtained by matching the matrix in (23) to itself is given by: The corresponding eigenvalues and eigenvectors of C normal is given by (25) and (26) Please note that the eigenvalues and eigenvector components are real numbers. It is also interesting to note that C normal shows symmetry about its main diagonal and that all the entries on the main diagonal are zeros.

Fault Cases
This section illustrates the implementation of graph-matching to obtain characteristic eigenvalues and eigenvectors for the fault cases of: (1) working fluid leakage, (2) compressor failure and, (3) fouling in the gas cooler.

Working Fluid Leakage
Refrigerant leakage may occur in a system during operations due to material imperfections or vibrations, which will result in a reduction of performance and therefore critical to detect at an early stage. Working fluid leakage was simulated by gradually decreasing the R744 flow rate as input parameter. The flow rate decreased from the normal case of 0.1559 kg/s to 0.0959 kg/s in increments of 0.02 kg/s. The reduced amount of refrigerant flowing through the GC and evaporator causes a reduction in the GC pressure and evaporating pressure. This is expected since at any given time during system operation there is less refrigerant charge inventory inside the heat exchanger components. Since the evaporator and compressor are being starved of refrigerant, the gas cooler will experience the knock-on effect. Starving the gas cooler will reduce the heat load because of limited refrigerant to reject any heat. With less heat to accept (reject from the compressor) the gas cooler will be operating at a lower temperature and pressure. The simulation results correspond to the results found in the study by Zhang & Dong et al. [40] who also found that gas cooler and evaporator pressure decreases with a decrease in total refrigerant mass flow rate in their investigated system. For the normal case the mass flow rate was programmed as a function of the EV orifice opening (refer to (A52) in Appendix A).
The node signature matrix obtained for the fluid leakage fault condition (ṁ R744 = 0.0959 kg/s) is numerically given by: The leakage fault condition cost matrix obtained by matching M normal to M F, leak is given by: The corresponding eigenvalues and eigenvectors of C leak are given by (29) and (30)

Compressor Failure
Compressor failure may occur when the seal and mechanisms inside a compressor wear with use [36]. Compressor failure was simulated by assigning a fixed value for isentropic efficiency and gradually decreasing it as input parameter from 0.6527 to 0.5027 in increments of 0.05. For the normal case the isentropic efficiency was however determined as a function of the compressor pressure ratio (refer to (A30) in Appendix A). The node signature matrix obtained for the compressor failure fault condition (isentropic efficiency = 0.5027) is numerically given by: Please note that the rate of irreversibility generated in the compressor (indicated in red) is high relative to the other irreversibility parameter values in the first column of the matrix (31). The effect of compressor failure on the heat pump's thermo-physical behavior is thus reflected directly in the increased magnitude of exergy destruction that takes place inside of the compressor. The compressor failure fault condition cost matrix obtained by matching M normal to M F, fail is given by: The corresponding eigenvalues and eigenvectors of C fail are given by (33) and (34)

Fouling in the Gas Cooler
Fouling may occur on the water side of the gas cooler due to precipitation of dissolved or suspended solids in the water to the tube wall. Fouling has the effect of introducing additional thermal resistance in the path of heat transfer inside the gas cooler. Fouling is simulated by keeping the required heat transfer rate constant at 35.882 kW while gradually increasing the thermal resistance on the water side of the gas cooler. Additional thermal resistance proportional to a fouling factor FF io (refer to (A48) in Appendix A) and inversely proportional to the water wall area is introduced. The fouling factor is increased from 0 m 2 ·K/kW to 1 m 2 ·K/kW in increments of 0.25 m 2 ·K/kW. The node signature matrix obtained for the fouling fault condition (FF io = 1) is numerically given by: Please note that the amount of power consumed by the compressor is relatively high for the fouling fault condition. The power consumption parameter value is 14.673 kW (indicated by the two entries highlighted in green). The fouling failure fault condition cost matrix obtained by matching M normal to M F, foul is given by: The corresponding eigenvalues and eigenvectors of C foul are given by (37) and (38)

Fault Signatures
This section explains the methodology followed to compose a fault signature for each of the three heat pump system faults considered.

Compilation of Fault Signatures
Fault signatures will be compiled for each of the three fault scenarios. The aim of a fault signature is to compile a simple visual representation that may be used for fault detection and diagnosis purposes. The eigenvalues and eigenvectors given by (25) and (26) will be used as the normal condition that characterizes ideal fault-free heat pump operation. The following steps will be followed in the compilation of the fault signatures: (1) The eigenvalues and their corresponding eigenvectors are determined for each fault case.
(2) The eigenvalues and eigenvectors of the normal operating condition and that of the three fault signatures are arranged into a table format. The eigenvalues will make up the first row of the table, and the subsequent rows will contain the eigenvectors as entries. (3) The table entries of each of the fault cases will be compared to the entries of the normal case.
The percent deviation of each corresponding element entries in the tables will be determined. (4) A margin value will be defined to determine if two entries are sufficiently different to be allocated a sign in the final fault signature in the corresponding table position. (5) A positive sign will be allocated to an element of a fault signature if the value of the entry of a fault condition minus the value of the entry in the normal condition is a positive number, a negative sign will be allocated if the value is determined to be negative, and a zero entry will be allocated if the compared values were determined to be below the specified margin of difference.

Percent Deviation Margin
From Section 6.1.3 a deviation margin of 5% will be used for this study. This allows for a heat pump system to have natural variation in operational characteristics that result in up to 5% variation in the eigenvalues and eigenvectors determined from the sensor readings. An entry that is compared to the reference entry in the normal condition that is determined to not deviate more than 5% from the reference will thus be assigned a zero value in the fault signature. In this manner, the emphasis is placed on the contribution of the eigenvalues and eigenvectors that have constituent entries that deviate more than 5% relative to the normal. Figure 9 outlines the process to determine the signs of the various elements of the fault signatures while considering the deviation margin.

Working Fluid Leakage Fault Signature
The compilation of the fault signature for the working fluid leakage fault will be shown step by step. The remaining two fault signatures are compiled similarly. The percent deviation of the eigenvalues and eigenvectors of the fluid leakage condition when compared to the corresponding values of the normal condition is given in Table 6.  The two values indicated in bold in Table 6 have percent deviation values below 5%. The two corresponding positions in the fault signature for fluid leakage will thus be filled with zeros.
To determine the rest of the element entries of the fault signature the numerical difference between the eigenvalues and eigenvectors of the fault condition and normal condition are calculated. The difference in the various eigenvalue and eigenvector elements for the fluid leakage fault condition is shown in Table 7. The fault signature for the working fluid leakage fault condition is obtained by retaining the sign of each difference entry in Table 7 and assigning a value of 1 to difference. Following this convention, a fault signature containing only +1, −1 and 0 entries result. The resulting fault signature for the working fluid leakage fault is given in Table 8.

Compressor Failure Fault Signature
The fault signature for the compressor failure fault condition was compiled using the same process as used to compile the fault signature for the working fluid leakage fault condition and the fault signature in Table 9 obtained. Table 9. Fault signature for compressor failure.
6.4.5. Gas Cooler Fouling Fault Signature Table 10 shows the fault signature obtained for the fouling fault condition.

Discussion
Since none of the fault signatures in Tables 8-10 contain only zeros, all three faults are detectable. The degree of detectability can be quantified by counting the number of positive and negative signs in each signature. In this way it is possible to evaluate the uniqueness of each individual fault signature. Table 11 summarizes the sign count for each fault signature. Table 11. Degree of detectability.

Fault Condition Plus Signs Minus Signs
Working fluid leakage 15 13 Compressor failure 10 19 Gas cooler fouling 14 14 Since the number of positive and negative signs for all three faults are high all three faults are highly detectable. When considering the isolability performance of the fault detection scheme, the uniqueness of the signatures is of interest. The isolability performance of the fault detection scheme can be evaluated by comparing the fault signature of each fault to all the others. A measure of the difference between two fault signatures can be obtained by subtracting the one signature from the other [21]. From [21] the distance between two 6 × 5 fault signatures x and y (as typical of the fault signatures in Tables 8-10) can be defined as F x (k, j) and F y (k, j) represent the entries in the respective fault signatures. By determining the distance between each of the three fault signatures and all others an indication of the degree of isolability was obtained. Table 12 summarizes the isolabilty performance.

Conclusions
The main conclusion from the results of the paper is that the approach of using energy attributed graphs for fault detection and isolation seems viable. The three fault cases considered are highly detectable and the isolability performance is good. Compilation of the fault signatures from the variation in the eigenvalues and eigenvectors from the normal case reveals that using only eigenvalues is not feasible. Using only eigenvalues would give the same fault signature for the working fluid leakage and compressor fault. Including the eigenvector variation in the fault signatures results in good isolability performance.
The results of the paper are however limited to illustration at on operating point, large fault sizes and an assumption of 4% energy performance variation. Further work is therefore first warranted in terms of sensitivity and robustness performance considering small faults and realistic variation in data respectively.
Further work is also warranted to investigate the feasibility of the proposed technique in a practical system. The challenge of measuring energy and exergy as parameters and determining the minimum set of attributes necessary for fault isolation needs to be investigated.
The derivation of a fault isolability matrix proves to be useful, further work is however warranted to relate the distances between fault signatures and the performance of the heat pump system. Proper validation against known heat pump reliability data is also an important next step. Funding: This research was funded by SASOL Ltd.

Conflicts of Interest:
The authors declare no conflicts of interest.

Appendix A. Component Model Description
The way in which each component in the heat pump system is represented using control volume (CV) regions is discussed in this section. The specific methodology used to solve for the key driving parameters present in the derived component-specific conservation and characteristic equations is also discussed. The modeling discussion proceeds in the chronological order of evaporator first, compressor second, gas cooler third and concludes with the modeling of the electronic expansion valve component.
The evaporator component is divided into two CVs on the refrigerant side and two corresponding CVs on the cooled water side, as illustrated in Figure A1. It is assumed that the mass flow rate, temperature, and pressure are known at the inlet of the R744 side of the evaporator. The amount of superheating is also known at the outlet of the evaporator component and is controlled by the EV orifice opening percentage. It is also assumed that the water side of the evaporator, the volumetric flow rate delivered by the pump, the inlet temperature, and the pressure of the water are known. The conservation of mass balance for the R744 side of the evaporator is given by: The conservation of mass for the water side of the evaporator is described by: The conservation of momentum for the R744 side of the evaporator is as follows: The total pressure drop, ∆P 0 L|E , that the R744 undergoes along the length of the evaporator is the sum of the total pressure drop that occurs in the evaporator's two-phase region and the total pressure drop that occurs in the preceding superheating region. The total pressure drop that the R744 undergoes in the two-phase region, ∆P 0 L|TP , is modeled using the Friedel method [41].
In the Friedel method the total lengthwise two-phase pressure drop, ∆P 0 L|TP , is a function of the pressure drop that would occur in the two-phase region if the whole pipe was filled with only a liquid, ∆P , multiplied by the squared two-phase Friedel multiplier, Φ tp f r . The superscript tp refers to a two-phase working fluid. The Friedel method is expressed mathematically as: The R744 liquid-based pressure drop is expressed as: where f l is the liquid-based friction factor, L TP is the length of the two-phase region in the evaporator, D H,R744 is hydraulic diameter of the homogeneous R744 flow, m tp vel is the mass velocity of the two-phase R744 mixture, and ρ b l is the liquid-based bulk density in the two-phase region. The liquid-based friction factor is calculated by use of the Blasius correlation [41,42]: where Re b l is the bulk dimensionless liquid-based Reynolds number in the two-phase region. The mass velocity of the R744 two-phase flow is simply the total R744 mass flow through the evaporator divided by the free flow area of the inner evaporator tube: After the liquid-based pressure drop has been determined; the two-phase Friedel multiplier is required to calculate the total two-phase pressure drop with relation to the liquid-based pressure drop in the two-phase region. The Friedel multiplier is expressed as a function of five dimensionless parameters as: (A8) The first dimensionless parameter is defined as: where x b is the bulk quality dryness fraction of the two-phase region and f v is the bulk vapor-based friction factor. The second dimensionless parameter is defined as: The third dimensionless parameter is defined as: The homogeneous Froude number is defined as: where g is the local gravitational acceleration at the physical location of the heat pump system and ρ b H is the bulk homogeneous density of the two-phase flow in the two-phase region. The liquid-based Weber number is defined as: where ρ b is the surface tension at the bulk conditions in the two-phase flow region.
The superheat region follows the two-phase region in the evaporator. The total pressure drop that occurs in the superheat region is expressed as: The Darcy friction factor for the superheat region is determined using the built-in EES MoodyChart function as follows: The total pressure drop on the R744 side of the evaporator is then the sum of the total pressure in the two-phase region and the total pressure drop in the superheat region: The length of the two-phase region added to the length of the superheating region must equate to the total length of the tube in the evaporator component: The focus of the simulation is on the modeling of the refrigerant side of the heat pump system. The total pressure drop on the water side of the evaporator is assumed to be negligible. Under this assumption the conservation of momentum for the water side of the evaporator is simply: The conservation of energy for the R744 side of the evaporator is: The total heat absorbed by the R744 during the evaporation process is equal to the total heat absorbed in the two-phase region and in the superheat region within the evaporator: The total enthalpy at the start of the superheat region can easily be calculated when the total pressure drop in the two-phase region is known. The total enthalpy at the start of the superheat region may thus be calculated in EES as follows: Once the total enthalpy at the start of the superheat region is known, the heat transfer in the two-phase region can be calculated by performing an energy balance for CV 1 as in Figure A1. The result is: Since the amount of superheat at the outlet of the evaporator component is controlled by the EEV operation, the heat transfer in the superheat region is calculated from the single-phase relation: The heat absorbed by the R744 during the evaporation process is supplied from the cooled water stream. The conservation of energy for the water side of the evaporator is thus given by: The exergy analysis of the evaporator is performed over the entire evaporator component. The exergy analysis thus considers all four CV regions that encapsulate the evaporator interactions simultaneously. The exergy balance over the entire evaporator is thus given by: The compressor component is represented using one CV. The compressor component with its CV indicated as illustrated in Figure A2. The compressor only has carbon dioxide as its interacting fluid with the conservation of mass for the compressor component given by: The discharge pressure of the R744 after compression is defined by the compressor's ratio of pressure: The conservation of energy for the compressor is given by: whereẆ C is the power absorbed by the compressor during the compression of the carbon dioxide gas. The compressor power drawn is the key parameter required to solve for the output flow energy delivered by the compressor. The compressor power drawn is calculated from the compressor's isentropic efficiency:Ẇ whereẆ isen is the power the compressor would draw in an ideal imaginary isentropic compression process between the same inlet and outlet thermo-physical conditions and η isen is the isentropic efficiency of the compressor with a value in the range of zero to one. The ratio of pressure across the compressor will be expressed as a function of the isentropic efficiency of the compressor at various operating conditions. Figure A3 illustrates the pressure-volume behavior of the refrigerant through the compression process inside the compression chamber of a reciprocating compressor. Four different stages are defined in the cyclic up-and-down movement performed by the compressor's pistons. At point, a in Figure A3, the suction valve opens during the movement of the piston and low-pressure refrigerant flows into the compression chamber. The piston continues to move and draw in refrigerant until it reaches its bottom dead center (BDC) as indicated by point b in Figure A3. At point b the suction valve closes and the piston starts moving in the opposite direction, compressing the refrigerant. The compression action is illustrated in Figure A3   The expression to determine the specific work required per piston to compress the refrigerant from the evaporating pressure to the gas cooler pressure is given by: where k avg is the average heat capacity ratio between the compressor's discharge and suction side, v b is the specific volume of the refrigerant at point b in Figure A3, and r P is the pressure ratio over the compressor during steady-state operation. The total specific work required by the compressor during steady operation is then the specific work required per piston times the number of pistons that compress the refrigerant: where #cyl is the number of compression cylinders that compress the refrigerant at any given instance during steady-state compressor operation. The isentropic efficiency of the compressor may thus be simulated by taking the theoretical specific isentropic work required by the compressor at the same operating conditions and dividing it by the calculated actual specific work required by the compressor: The accuracy of (A32) can be improved by including two constants into the equation that were determined from test bench data. The improved equation for the isentropic efficiency of the compressor is given by: where a1 = 1.629497 and b1 = 0.869919. The isentropic efficiency of the compressor is assumed to be known in this study and its value is used as an input parameter in the simulation of the heat pump system. This approach may be used to implicitly determine the ratio of pressure over the compressor when the isentropic efficiency of the compressor is known, as is evident from inspection of (A30). The exergy balance for the compressor is given by:

. Gas Cooler
The gas cooler is divided into 20 control volumes on the refrigerant side and 20 corresponding control volumes on the heated water side. This was done to accurately calculate the heat transfer in the gas cooler component using the ε-NTU method. The segmentation is also beneficial when plotting the supercritical heat rejection curve of the refrigerant on a T-s diagram. The gas cooler component with its control volume regions is illustrated in Figure A4. It is assumed that the mass flow rate, temperature, and pressure are known at the inlet of the R744 side of the gas cooler. It is also assumed that the volumetric flow rate, inlet temperature and pressure are known on the water side. The conservation of mass balance on the R744 side of the gas cooler is given by: The conservation of mass balance on the water side of the gas cooler is: The density of the water, ρ w , [kg/m 3 ] andq GC the volumetric flow rate of the heated water through the gas cooler has a unit [m 3 /s]. The conservation of momentum for the R744 side of the gas cooler is where ∆ P0L|GC is the total pressure drop of the R744 along the length of the gas cooler. The total pressure drop is a function of flow channel friction, the total coil length, the hydraulic diameter, and the thermo-physical properties of R744 as it flows through the gas cooler. The total pressure drop on the R744 side of the gas cooler is equal to the sum of the pressure drop over each segment on the refrigerant side. The pressure drop per increment on the refrigerant side is given by: where f W is the friction factor correlation by Wang et al. [43] for the prediction of supercritical carbon dioxide pressure drop in the turbulent flow regime, L inc is the length of the increment, D H,R744 is the hydraulic diameter, ρ b g is the bulk density per increment, and V b g is the bulk velocity per increment.
The parameter e ss is the surface roughness of the inner gas cooler tube material and Re b g is the dimensionless bulk Reynolds number of the R744 per increment. The correlation by Wang et al. has a valid Reynolds number range of: 3400 < Re b g < 2 × 10 6 . The pressure drop on the water side of the gas cooler is assumed to be negligible. Under this assumption the conservation of momentum equation for the entire water side of the gas cooler is given by: The conservation of energy for the R744 side of the gas cooler is: The energy extracted from the R744 side of the gas cooler during the cooling process is absorbed by the heated water side. The conservation of energy for the water side of the gas cooler is described by: where Q GC is the total heat transfer that takes place from the R744 to the water stream. The total heat transfer is the sum of the heat transfer per increment on the R744 side of the gas cooler. The heat transfer per increment will be solved by using the ε-NTU method. The heat transfer per increment is given by: where Q max|inc is the maximum theoretical heat transfer possible per increment and ε inc is the effectiveness of the gas cooler in transferring heat per increment in the range of 0.1. The maximum heat transfer per increment is given by: where C min is the minimum heat capacity rate between the two interacting fluid streams and ∆T max is the temperature difference between the R744 and water, at their respective inlets, per increment.
The effectiveness per increment for a counter-flow tube-in-tube type of heat exchanger is defined by Bergman et al. [44] as: where NTU is the dimensionless number of transfer units per segment of the gas cooler and c r is the dimensionless capacity ratio of the interacting fluids. Equation (A45) is valid for all possible values of the capacity ratio that are smaller than one. The capacity ratio is expressed as the ratio of the minimum to maximum heat capacity rates of the interacting heat exchanger fluids per increment: where c p b is the bulk specific heat capacity at constant pressure per increment for a given fluid stream.
The number of transfer units per gas cooler increment is given by: where U is the overall heat transfer coefficient per gas cooler increment and A s is the heat transfer surface area per increment. It is assumed that there is no possibility of fouling build-up on the R744 side of the gas cooler while it is possible for the water side to foul. Furthermore, it is also assumed that the gas cooler is well insulated on the outside of the outer tube such that there is no possibility of heat loss to the environment. Under these assumptions the product of the overall heat transfer coefficient with the heat transfer surface area per gas cooler increment may be expressed in terms of thermal resistance as [45]: (A48) The first term in (A48) is the total thermal resistance per gas cooler increment, the second term is the convection thermal resistance on the R744 side and the third term is the conduction thermal resistance per increment. The fourth term is the fouling resistance on the heated water side and the fifth term is the convection resistance of the heated water side on the gas cooler.
The Dittus-Boelter correlation [46] is used to solve for the convection coefficient, h conv , on both the R744 and the water side of the gas cooler. Although correlations exist that can predict the convection coefficient of supercritical R744 with higher accuracies such as those investigated by Venter [47], the study by van Eldik et al. [38] showed that the Dittus-Boelter correlation shows an excellent fit (for larger diameter tubes) with experimental data over a wider range of Reynolds number values. The study by [38] also concluded that the Dittus-Boelter correlation remains relatively stable when used to calculate convection coefficients based on thermo-physical conditions in the pseudo-critical region of carbon dioxide. The convection coefficient per gas cooler increment while implementing the Dittus-Boelter correlation is expressed as: where CF is a correlation factor specifically fitted to the specific test bench system to further increase the accuracy of the Dittus-Boelter correlation, k b is the bulk conduction heat transfer coefficient of the fluid, Pr b is the bulk dimensionless Prandtl number of the fluid per gas cooler increment and n is a power set to 0.3 if the fluid considered is cooled or set to 0.4 if the fluid is heated. The CF factor for the experimental data used in this study is 2.294 [-]. The CF factor is included to make the prediction of heat transfer values associated with the test bench system as accurate as possible. The CF factor reported here is only valid for this study. The CF factor used here is an average of all the individual factors determined by setting the simulation results equal to the experimental results. The CF factor is thus a form of direct numerical fitting of the simulation to the test bench results. The exergy analysis of the gas cooler is done over the entire component. The exergy balance over the entire gas cooler is given:

. Expansion Valve
The EV component is represented with one CV. The EV component with its CV is illustrated in Figure A5. The EV has only carbon dioxide as its interacting fluid. The conservation of mass for the R744 flowing through the EV component is given by: The EV controls the mass flow rate of the refrigerant through the entire heat pump system. The BITZER TM Software v6.7.0 utility [48] may be used to find the recommended mass flow rate for the BITZER TM type 4JTC-15K (40P) compressor operating at specific working points. The amount of refrigerant that was loaded into the heat pump test bench was done by an operator. The mass flow rate of the refrigerant at various locations throughout the heat pump system may thus not match the mass flow rate recommended by the BITZER TM Software utility. A better approach is a linear regression fit of the refrigerant mass flow as a function of the EV orifice opening percentage for the actual test bench. The built-in EES ® [23] linear regression utility was used to fit the mass flow rate of the refrigerant to the experimental data generated from the heat pump test bench. The fitted linear relation of mass flow as a function of the EV orifice opening is illustrated in Figure A6.
From Figure A6 it should be noted that a decrease in the flow restricting orifice opening area percentage leads to an increase in the mass flow rate of the refrigerant through the entire heat pump system. The specific linear equation that describes the relationship between the refrigerant mass flow rate and EV orifice opening of the line displayed in Figure A6  The conservation of momentum for the R744 that undergoes throttling as it passes through the EV is given by: Since the physical dimensions of the internal EV mechanisms are not readily available, a new method is required to predict the total pressure drop over the EV component. The static pressure drop over the EV is modeled as an inverse linearly proportional function of the opening percentage of the expansion valve orifice area. The suggested inverse linearly proportional graph of the static pressure drop over the EV as a function of the opening percentage of the expansion valve orifice is illustrated in Figure A7. From Figure A7 it should be noted that the static pressure drop over the expansion valve orifice tends to zero as the valve orifice is increasingly opened. This trend is suited because when the EEV orifice is 100% open it will not offer any resistance to the flow of the refrigerant and the induced pressure drop on the refrigerant will thus be zero. The specific linear equation that describes the line displayed in Figure A7 is: Once the static pressure drop over the EV is known, it can easily be converted to total pressure drop by adding the kinetic energy component of total pressure to the static pressure drop value. The throttling of the R744 through the EV is assumed to be an isenthalpic process. Under this assumption the conservation of energy for the EV is simply: Another important cycle parameter controlled by the EV is the amount of superheating that takes place in the evaporator. A linear equation is also used to describe the amount of superheating as a function of the EV orifice opening. The linear relationship between the amount of evaporator superheating and the EV orifice opening is shown in Figure A8. The specific linear equation that describes the line displayed in Figure A8 is: ∆T SH = 11.8 100 (Open%) + 3.62 .
The exergy balance over the EV is performed over the entire EV CV region. The exergy balance for the EV is given by: