Port-Hamiltonian Modeling of Thermofluid Systems and Object-Oriented Implementation With Modelica I: Thermodynamic Part

In this paper, we present the physical foundations and the development of the thermodynamic part of a Modelica library with the fundamental components for modeling thermofluid systems. We have chosen Modelica because it is an object-oriented modeling language that allows an elegant design of the library, with a top-down conception that starts from very general components where we model the thermodynamic properties common to all simple substances and descend by inheritance to model the properties of each particular substance. To model the behavior of each component, we have used: classical thermodynamics to define the equilibrium states, the local equilibrium hypothesis of Classical Irreversible Thermodynamics to model the changes of state, and the port-Hamiltonian approach to obtain the equations of the system dynamics. With this formulation, we implement the thermodynamic behavior of ideal gases (including monatomic gases as a particular case), the 2073 substances defined for the CEA (Chemical Equilibrium with Applications) NASA Glenn computer program, the IAPWS Formulation 1995 for the Thermodynamic Properties of Water Substance for General and Scientific Use, and the Syltherm 800 HTF (Heat Transfer Fluid). We also define graphical symbols for each library component that facilitate modeling complex systems with simple drag-and-drop manipulations, component connection, and parameter selection. These symbols are a slightly modified version of those used in bond graphs to facilitate their reading and the representation of the structure of complex systems. We also show the modeling, simulation, and comparison for accuracy, performance, and scalability of some thermodynamic systems implemented with the Modelica Standard Library (MSL) and the proposed library.


I. INTRODUCTION
This paper is the first part of a two-part series where we will implement a library of components to model the storage and transport phenomena of mass, momentum, and energy in thermofluid systems from a port-Hamiltonian approach. In [28], we presented the general theory of port-Hamiltonian systems, its graphical representation with bond graphs, its objectoriented implementation with the Modelica language, and its application to simple mechanical systems and electrical networks. In the current paper, we present the port-Hamiltonian The associate editor coordinating the review of this manuscript and approving it for publication was Yizhang Jiang . modeling of thermodynamic systems, we implement in Modelica the main building blocks of these systems, we propose a graphical representation based on bond graphs adapted for thermodynamic systems, and we model the thermodynamic behavior of many substances of scientific and industrial interest.

A. MOTIVATION
Thermofluid systems play a central role in industries such as oil and gas extraction and distribution, water processing and distribution, or solar thermal energy production and storage. Therefore, it is important to have tools for easy and quick modeling and simulation of these systems. These models will also be helpful to inexpensively test control strategies for an efficient system management. A general thermofluid system may be composed of several subsystems with different substances in motion or at rest, mixing or chemically reacting and exchanging mass, momentum, and energy with each other and their surroundings.
The modeling of these systems is necessarily multidisciplinary and relies on disciplines such as fluid dynamics to model the movement of thermofluids, chemistry to model the reactions between them, thermodynamics to model energy exchanges, electromechanics to model fluid pumping devices and their control, etc. Added to this complexity is the fact that a complete description of the behavior of a continuous substance requires describing its time evolution at each point in space, i.e., each magnitude x that provides information about a thermofluid system will depend on the time variable t ∈ R and the spatial variable r ∈ ⊂ R 3 .
To tackle a problem of this complexity and to be able to draw engineering applicable conclusions, it is necessary to make some simplifications that make it intellectually and computationally tractable. The main simplification of this work is to consider thermofluid systems composed of connected subsystems that can be described by a set of easily characterized lumped parameters. With this simplification, we can describe the system behavior with a differentialalgebraic system of equations (DAE) which allows using the tools for modeling lumped parameter multiphysics systems described in [28]: the port-Hamiltonian formalism, the bond graph representation, and the object-oriented implementation with the Modelica language. We can summarize the rationale for this choice as follows: • The port-Hamiltonian formalism starts from the Hamiltonian formalism of classical mechanics and extends it to include non-conservative systems (i.e, systems with energy dissipation) [50]. It provides a framework for modeling physical systems based on the storage and exchange of energy through power ports that interconnect the parts of a system with each other and with their surroundings. The connection of systems with port-Hamiltonian structure forms a new system that preserves this structure, so this formalism is suitable for modeling complex systems by connecting simpler ones. The port-Hamiltonian approach has proved to be also suitable for controller designing [50,Chap.7].
• The bond graph methodology proposed by H. Paynter in the 1960s arises from modeling ideas used in electrical networks and it is based on principles similar to the port-Hamiltonian formalism, but is not formulated with the rigor of differential geometry as the latter [4], [40]. However, bond graphs define a symbology to represent the structure of complex systems and the routes of energy exchange between their components, so that they are commonly used to represent port-Hamiltonian systems.
• The Modelica language is also used to describe the behavior of complex systems by connecting simple components (composition mechanism) [35]. As an object-oriented language, it is designed to model the behavior of general systems and to extend those models to particular systems derived from the previous ones (inheritance mechanism). Tools such as OpenModelica [19] or Dymola R process symbolically these models to obtain the differential equations for simulating the system behavior.
In this article, we study the modeling of phenomena related to energy storage in thermofluid systems without mass transfer (fluids at rest). For this purpose, we employ Equilibrium Thermodynamics [7], [48], which allows us to identify states in a thermodynamic system, and also employ Classical Irreversible Thermodynamics [12] to study the time evolution of these states.
We also show the implementation in Modelica language and the graphical representation of components necessary to model these phenomena, thus extending the library pHlib for modeling multiphysics systems started in [28].
We exemplify the applicability of the procedure by showing how to model the thermodynamic behavior of some substances of industrial interest: we started with the model of general simple substances, continue with the model of ideal gases, add the models of 2073 substances described in [30] for the CEA NASA Glenn computer program, implement the model for water and steam described by the International Association for the Properties of Water and Steam (IAPWS-95 [22]), and finish with the model of the commercial Syltherm 800 HTF [15].
Finally, we also implement models of simple systems with MSL and with pHlib in order to compare the results of both simulations. The new library components have been implemented and simulated with the Dymola tool.

B. LITERATURE REVIEW
The development of the bond graph methodology for modeling multiphysics systems began in the 1960s with the work of H. Paynter [40] and, since then, this methodology has been used to model multidisciplinary engineering systems [4], [5], [24].
The theory of port-Hamiltonian systems began in the late 1990s [49]. It is an evolution of Hamiltonian mechanics [1,Chap.3] to which it adds the concepts of power ports (inspired by bond graph modeling) and Dirac structures [11] to formalize with the language of differential geometry a framework for modeling complex systems by connecting simple systems that exchange energy [16].
In [20], [28], we can find a discussion of the relationship between bond graphs and port-Hamiltonian systems, and [14], [42] have presented the derivation of port-Hamiltonian models from bond graphs.
We can find the Modelica language specification in [35], and [18] is a comprehensive guide showing the use of the language for modeling multidomain systems. References [10], [13], [54] show examples of the development of a Modelica library to model systems using the bond graph methodology. VOLUME 9, 2021 Two standard references in Equilibrium Thermodynamics and Classical Irreversible Thermodynamics are [7] and [12], respectively. The extension of the Hamiltonian formalism to include the phenomena of dissipation that appear in irreversible thermodynamics can be found in [25], [32]. In [29], [33] this formalism is expressed in terms of metriplectic geometry where the phase space of a system with dissipation is modeled with a metriplectic manifold. This geometric object is a manifold with a double structure: simplectic and metric. The conservative part of the behavior (Hamiltonian dynamics) is modeled with the simplectic structure, and the dissipative behavior (irreversible thermodynamics) with the metric structure. In Section III, we show the relationship between the metriplectic and port-hamiltonian formalisms. We will focus in this article on the modeling of lumped parameter systems, as discussed in the previous section. For a review of the literature on port-Hamiltonian modeling of distributed systems, we refer to [44].
The use of bond graphs for modeling thermodynamic systems can be found in [8,Chap.8], [45], [16,Chap.3]. The description of the MSL implementation for modeling thermofluid systems is in [17], and [9] describes the Modelica implementation of a bond graph library for modeling thermofluid systems.

C. CONTRIBUTIONS
In this paper, we develop an object-oriented implementation, with the Modelica language, of library (pHlib) which provides components for modeling thermofluid systems. As this objective is quite ambitious, we start modeling the thermodynamic properties (fluid at rest). In a later article, we will include the transport phenomena (fluid in motion). The implementation is based on a port-Hamiltonian approach to the concepts of Equilibrium Thermodynamics and Classical Irreversible Thermodynamics and is complemented by a new set of symbols derived from the traditional bond graph representation, modified for expressiveness and readability.
This approach brings together three traditions for modeling physical systems whose communities have had little connection with each other: port-based modeling with bond graphs, Hamiltonian modeling, and object-oriented modeling with Modelica. For this purpose, we have taken advantage of the fact that these three approaches share concepts such as power ports (Modelica generalizes this concept with the class connector) or the hierarchical organization (class inheritance and composition in Modelica). Although there are references that combine some of these modeling approaches (e.g., [13], [16], [20], and [42]), they do not address all three together. Moreover, our approach has an additional advantage: it is not necessary to perform a causality analysis when building the models (as required by traditional bond graph modeling) because neither the port-Hamiltonian modeling nor the tools that process the Modelica language need causality indications to automatically generate the computational equations of the models and simulate them.
Most developments based on the ISO representation of port-Hamiltonian systems rely on coordinatedependent matrix representations [50,Chap.6]. In Section III, we propose a coordinate-free definition based on differential geometry that allows us to compare the port-Hamiltonian formalism with the metriplectic one [29] and to demonstrate in a general way its ability to model dissipative systems.
The object-oriented approach of the library allows to: (i) Implement a generic component to model the thermodynamic behavior of any substance defined by its fundamental equation. (ii) Particularize, from the previous generic component, other components to model less general substances such as ideal gases, perfect gases, and monatomic gases (in this phase of the implementation, we have used the inheritance mechanism of the Modelica language and the possibility of adding new behavioral equations in each refinement step). (iii) Implement the substance library defined for the CEA NASA Glenn computer program. Our implementation takes into account the different coefficients defined for each temperature range (from one-range substances up to five-range substances). The MSL implementation takes into account only two intervals. (iv) Implement the version for scientific use (known as IAPWS-95) of the behavioral model and thermodynamic properties of water defined by the International Association for the Properties of Water and Steam. MSL implements the industrial formulation IAPWS-IF97. (v) Implement the thermodynamic behavior of substances whose fundamental equation is unknown, but some thermodynamic properties (e.g., specific heat, density, etc.) are available (we have chosen the commercial Syltherm 800 HTF as an example).
We have always worked with the variables that Thermodynamics, the port-Hamiltonian formalism, and the bond graph modeling postulate for the thermodynamic domain: • entropy S, volume V , and mass m for the state x, • temperature T , pressure p, and chemical potential g for the effort dU , and • entropy flow rateṠ, volume flow rateV , and mass flow rateṁ for the flowẋ.
Note that the handling of entropy introduces conceptual and computational complexity into the models. For this reason, in industrial practice, it is common to work with heat Q and heat flow rateQ instead of entropy and entropy flow rate. The terms true bond graphs and pseudo-bond graphs have been introduced to distinguish bond graphs using entropy flow rate from those using heat flow rate [4], [24].
Despite the difficulty of working with entropy, the use of this variable can help to detect inconsistencies in thermodynamic models by checking compliance with the second principle of thermodynamics.

D. PAPER ORGANIZATION
This paper is organized as follows. In Section II, we summarize the main results on Equilibrium Thermodynamics and justify the election of Classical Irreversible Thermodynamics for modeling the time evolution of thermodynamic systems. From the concept of the fundamental equation that models the internal energy of a thermodynamic system, we can define port-Hamiltonian models for the storage elements of that internal energy. Since this fundamental equation is not always available, in many practical cases the system dynamics must be reconstructed from measurable thermodynamic properties. For this reason, in Section II we also review: i) other ways of expressing the fundamental equation by the Legendre transformation (Helmholtz free energy, enthalpy, and Gibss free energy), ii) the fundamental set of response functions (specific isochoric heat capacity, adiabatic compressibility, and adiabatic expansibity), and iii) the primary set of response functions (specific isobaric heat capacity, isothermal compressibility, and isobaric expansibity).
In Section III, we study the ISO representation of dissipative port-Hamiltonian systems to identify its energy conservative and entropy generating parts, and to relate this representation to the metriplectic formalism.
In Section IV, we explain the equations and the Modelica implementation of simple thermodynamic systems. The components that make up the implemented library are as follows: power ports, bonds, effort and flow sources, storage elements, general models of single substances, models of ideal substances (ideal gases, perfect gases, monatomic gases), empirical models for non-perfect substances, IAPWS-95 model for water, and model for Syltherm 800.
In Section V, we model simple systems combining electrical and thermodynamic parts, implement them with MSL and pHlib, simulate them, and analyze the results.
In Appendix, one can see two lists with acronyms and nomenclature.

A. SIMPLE THERMODYNAMIC SYSTEMS
A simple system is a homogeneous and isotropic mixture of r chemical substances not acted by electric, magnetic, or gravitational fields. The postulates of equilibrium thermodynamics [7, Chap.1] establish for simple systems that: (i) There exist equilibrium states characterized by the following extensive parameters: internal energy U , volume V , and the amount n j (in mole) of each substance. (ii) For all equilibrium states there exists a real valued function of the extensive parameters (called entropy S), that is continuous, differentiable, additive over the constituent subsystems, and strictly monotonically increasing with the energy. (iii) The values assumed by the extensive parameters are those that maximize the entropy over the configuration space (manifold of equilibrium states).
The strict monotonic property of the entropy function implies that it can be inverted with respect to the energy and, by postulate (ii), the energy is a single-valued, continuous and differentiable function of S, V , and each n 1 ,. . . , n r (simplify n 1...r ). Both S = S(U , V , n 1...r ) and U = U (S, V , n 1...r ) are homogeneous first-order functions that contain all thermodynamic information about a system and each one is known as the fundamental equation.
The partial derivatives of U are called intensive parameters and are zero-order homogeneous functions with a concrete meaning in thermodynamics: Each curve that takes values in the configuration space is a dense sequence of equilibrium states that lacks information about the dynamics of state change. However, a real thermodynamic process is a sequence of equilibrium and nonequilibrium states, but no general theory of non-equilibrium thermodynamics is currently available. Attempts to extend equilibrium thermodynamics yielded the elaboration of the Classical Irreversible Thermodynamics (CIT) in the first half of the 20th century [12], [36], [37], [43]. Subsequent efforts to overcome the limitations of the CIT have resulted in formulations such as: Rational Thermodynamics (RT) [47] which has contributed to the development of continuous media thermodynamics; the Extended Irreversible Thermodynamics (EIT) [23] which provides a bridge between CIT and RT; and the formalism known as GENERIC (General Equation for the Non-Equilibrium Reversible-Irreversible Coupling) [39] which develops a non-equilibrium thermodynamics with Hamiltonian structure. Nowadays, more efforts to model nonequilibrium thermodynamic processes join the program of geometrization in physics and employ the language of differential geometry to address an invariant (coordinate-free) formulation of phenomenological thermodynamics [51].
In this paper, we have chosen to use CIT because it fits easily into the port-Hamiltonian approach and it has proved to be suitable for modeling a large number of non-equilibrium thermodynamic processes, as confirmed by experimental results. The main additional postulate of CIT is the local-equilibrium hypothesis: ''small volume elements of a system may be considered to be in thermodynamic equilibrium locally, although the whole system may not be so, and the same equilibrium postulates (i-iii) can be applied to those local volumes'' [48,Chap.23].
One can use this approach to approximate continuous systems models using networks of lumped parameter subsystems (also called networks of control volumes [38], [41]). Thus, if we call R + = {t ∈ R : t > 0}, we consider the smooth manifold canonical structure of R 2+r + , the thermodynamic state of a control volume x = (S, V , n 1...r ) ∈ R 2+r + , and the energy state function U : R 2+r + → R, then we can write the following energy flow rate for thermodynamic systems [12, p.23] (see [28,Appendix] for a summary on vector bundles and braket notation) that corresponds to a storage element C = (R 2+r + , U ) (see [28,Definition 23]) whose state space is the smooth manifold R 2+r + , and with Hamiltonian U (S, V , n 1...r ). The power port of C (see [28,Definition 1] and the effort variable is the covector field If we express the fundamental equation as a function of the mass m j = n j · M mol,j of each chemical species, where M mol,j kg/mol is the molar mass of the substance j, then U = (S, V , m 1...r ) anḋ Figure 1 is the thermodynamic domain adaptation of [28, Figure 7]. Part a) represents the storage element of a simple thermodynamic system defined by the smooth manifold X , Hamiltonian U , and power port P = dU |ẋ . Part b) is the bond graph representation of this storage element.
Since energy and entropy are homogeneous first-order functions, they can be scaled from their definition for a unitary mass system: where u = U /m, s = S/m, and v = V /m are called specific (or massic [46, Section 8.9]) energy, specific entropy, and specific volume respectively, and w j = m j /m is called mass fraction. As intensive parameters are zero-order homogeneous functions, they have the same form when expressed with specific variables, i.e., By applying Euler's theorem 1 to the first-order homogeneous function U we obtain the Euler's relation [7, p.51] that in the entropy representation has the form Single-substance thermofluid systems have a special interest in thermal engineering to design energy transport and storage systems. In this case, the fundamental equation, the equations of the intensive parameters, and the equation of the energy flow rate take a simple form when expressed with specific quantities:

B. THERMODYNAMIC POTENTIALS AND THERMODYNAMIC PROPERTIES
Although the fundamental equation (in either of its two representations: energy or entropy) contains all the information of 1 If f (x 1...r ) is a homogeneous functions of degree k and has continuous first partial derivatives, then: a thermodynamic system, its theoretical deduction is beyond the scope of thermodynamics and corresponds to statistical mechanics [27], [7,Part II]. Experimentally, however, some approximations to the fundamental equation or its differential can be obtained. When experimenting with a thermodynamic system, it has been found that some variables can be more easily measured than others. Thus, for example, there are no instruments to measure the entropy (extensive variable) of a system. However, the temperature (intensive variable) can be easily measured with a thermometer. It would thus be convenient to be able to formulate the fundamental equation taking as independent variables those extensive and intensive variables that can be easily measured. The appropriate tool to perform this alternative formulation without loss of information is the Legendre transformation [2].
Let y(x 1... ...n ) be a function with partial derivatives p k = ∂y/∂x k . The Legendre transformation of y with respect to x 1... is a function of the independent variables p 1... , x +1...n , denoted y[p 1... ] and defined Note that the sign criterion used in thermodynamics for the Legendre transformation is the opposite of the one used in mechanics.

1) THERMODYNAMIC POTENTIALS
The Legendre transforms of the fundamental equation U are also called thermodynamic potentials (by analogy with the mechanical potential energy) and some of the most used are [2]: • Helmholtz free energy F (or Helmholtz potential): • Gibbs free energy G (Gibbs potential, or free enthalpy): It is important to note that the total Legendre transform U [T , −p, g 1...r ] is identically zero, leading to the Gibbs-Duhem equation (24) which establishes a relation between the intensive parameters of a simple system. Indeed, and, although systems with r substances have r + 2 intensive parameters, only r + 1 of them are independent. It is said that those systems have r + 1 thermodynamics degree of freedom.

2) PRIMARY SET OF RESPONSE FUNCTIONS
Just as the first partial derivatives of the fundamental equation have an important physical significance (temperature, pressure, and chemical potential), the second partial derivatives for single substances, called response functions [48,Chap.10], are also of great physical interest because they characterize the thermodynamic properties of these substances (e.g., Syltherm 800, see Example 5).
The primary set of response functions is defined from the second-order derivatives of the Gibbs free energy G [26, Chap.1]: • Specific isobaric heat capacity or specific heat capacity at constant pressure is the heat flux per unit mass required to produce a unit increase in the temperature of a system maintained at constant pressure. 2 • Isothermal compressibility is the fractional decrease in volume per unit increase in pressure of a system maintained at constant temperature. The inverse of κ T is called isothermal bulk modulus β T = 1/κ T .
• Isobaric expansibity is the fractional increase in volume per unit increase in temperature of a system maintained at constant pressure.

3) FUNDAMENTAL SET OF RESPONSE FUNCTIONS
The fundamental set of response functions is defined from the second-order derivatives of the internal energy U [48, Chap.10]: • Specific isochoric heat capacity or specific heat capacity at constant volume is the heat flux per unit mass required to produce a unit increase in the temperature of a system maintained at constant volume.
• Adiabatic compressibility is the fractional decrease in volume per unit increase in pressure of a system maintained at constant entropy.
• Adiabatic expansibity is the fractional increase in volume per unit increase in temperature of a system maintained at constant entropy. Some useful expressions relating the fundamental and primary sets of response functions are as follows:

III. ISO REPRESENTATION OF DISSIPATIVE PORT-HAMILTONIAN SYSTEMS
In this section, we will use the input-state-output (ISO) representation of dissipative port-Hamiltonian systems to compare the port-Hamiltonian and metriplectic formalisms [29], [33] introduced to extend the Hamiltonian formalism to dissipative systems. Here H stands for Hamiltonian, not for enthalpy. The kernel representation of the port-Hamiltonian system (X , F, D, H ) is (see [28,Definition 6 and Remark 12]) that yields the following local matrix equation with rankF n = n and rank[F n |E m ] = n + m, Equation 35 can be written (35) can be explicitly solve forẋ Since F E + E F = 0, F −1 E is skew-symmetric and can be written with J and B skew-symmetric matrices. If there is not a direct f to e coupling, then B = 0 and the explicit equations that model the behavior of a port-Hamiltonian system arė If we split F⊕E into an open input-output port U⊕Y (with Y = U * ) and a resistive port F R ⊕ E R (energy-dissipating port, see Figure 2), the matrix F −1 E, and the equations of the system can be written If the resistive port is defined by the linear relation f R = −R e R (with R symmetric and positive definite), then where The previous equations let us introduce the following geometric and coordinate-free definition of an ISO port-Hamiltonian system: Definition 1: An ISO port-Hamiltonian system is a dynamical system defined by the 6-tuple (X , H , J , R, U, g) where: 1) The smooth manifold X is the state space.
2) The Hamiltonian smooth function H : X → R models the energy of the system. 3) The skew-symmetric bundle morphism J : T * X → T X represents the internal power-preserving interconnections. The isomorphism Hom(T * X , T X ) ∼ = Hom(T * X ⊕ T X , R) allows to identify J with a skew-symmetric bilinear form (or antisymmetric order 2 tensor). 4) The symmetric positive semidefinite bundle morphism R : T * X → T X represents the dissipation of the system. The isomorphism Hom(T * X , T X ) ∼ = Hom(T * X ⊕ T X , R) allows to identify R with a symmetric positive semidefinite bilinear form (or symmetric positive semidefinite order 2 tensor). 5) The input power variable u is a constant section of the trivial vector bundle X ⊕ U, i.e., u depends on the time t but not on the state x. 6) The output power variable y is a constant section of the trivial vector bundle X ⊕ Y = X ⊕ U * . 7) The bundle morphism g : X ⊕ U → T X describes the distribution of the power incoming into the system from the environment. Its dual morphism is g * : T * X → X ⊕ Y. 8) The behavior of the system is represented by the equations: Since J is a skew-symmetric morphism then dH | J | dH = 0, and the power balance of an ISO port-Hamiltonian system iṡ To relate the dissipative port-Hamiltonian systems to the metriplectic formalism we split the Hamiltonian into two summands H = H c − H d , with the additional requirements (known as degeneracy constraints [39,): The dynamics of a closed system is then ( [33, Eq.25], [29,Eq.18]) The first term corresponds to the symplectic structure and the second one to the metric structure (so the name metriplectic system). In addition, the following properties are fulfilled: The above equations mean that H c is the conserved energy (Ḣ c = 0), H d is the dissipated energy, andḢ d ≥ 0 implies that dissipation is an irreversible process (the metriplectic formalism mentions this term as an entropy-like function). For closed thermodynamic systems we can identify H = F (Helmholtz free energy, see Equation 18), H c = U (internal energy), H d = TS, and writė i.e.,Ṡ = dS |ẋ ≥ S 2 /T dT | R | dT ≥ 0, which means an increase of entropy in an irreversible process.

IV. MODELICA IMPLEMENTATION OF SIMPLE THERMODYNAMIC SYSTEMS
In this section, we take a further step in the development of the Modelica library pHlib started in [28]. This library consists of two main packages: Dirac and pHS. In the package Dirac, we modeled the following concepts: power port, bond, general Dirac structure, and junction structures (0-junctions, 1-junctions, transformers TF, and gyrators GY). With these elements, one can define the non-dissipative energy exchange structure of a complex system by connecting elementary junction structures. In the package pHS, we modeled the following concepts: flow and effort sources (Sf, Se), storage element (C), and resistors (R, RS). These elements are used to model the energy exchanged with the environment, the energy stored in the system, and the dissipation of energy.
In this article, we take full advantage of the fact that Modelica is an object-oriented language to extend the pHlib library to the thermodynamics domain. To do so, we implement a new package called Thermodynamics by reusing the classes defined in the packages Dirac and pHS.

A. POWER PORTS AND BONDS
In the remainder of this paper, ''thermodynamic system'' will refer to a simple single-substance thermodynamic system unless otherwise stated.
We remember that, for the smooth manifold X with tangent bundle T X and cotangent bundle T * X , the fiber product P = T X ⊕ T * X is a power port if the dual product ω|X -where the effort ω ∈ T * X is a covector field and the flow X ∈ T X is a vector field-, has physical dimensions of power [28,Definition 1]. This definition applies to thermodynamic systems, where: Listing 1 shows the Modelica implementation of thermodynamic ports.
As we can see, the effort and flow variables of type Real that appear in the general definition of the connector Dirac.Interfaces.Ports.Port have been redeclared with units of the International System (SI), so they have meaning in the thermodynamic domain and can be checked for dimensional consistency.
Listing 2 shows the extension of the class Bond (used to connect subsystem ports in a complex system) to the thermodynamic domain.
It should be noted that, in defining new bonds for the thermodynamic domain, we have extended the base types defined in Dirac.Interfaces.Bonds to include variables with names associated with each type of bond. These variables In Figure 3 we can see the symbols used for a graphical representation of bonds and other elements defined in the package Thermodynamics.

B. SOURCES
The sources of flow (Ṡ,V ,ṁ) and effort (T , p, g) have been derived from the general ones defined in pHS, and redeclared according to their physical meaning. Listing 3 shows the implementation of modulated and constant sources.

C. STORAGE ELEMENTS
The class pHS.Storage.C is the base class for defining storage elements. The application of [28, Definition 23] to thermodynamic systems is as follows: Definition 2: A thermodynamic storage element is a 1-port element defined by the pair (R 3 + , U ) where R 3 + has the structure of a smooth manifold that models the state variables x = (S, V , m), and U ∈ C ∞ (R 3 + ) is the Hamiltonian that models the system energy. The flow and effort variables of the port are the vector fieldẋ ∈ T R 3 + and the covector field dU ∈ T * R 3 + . Listing 4 shows the implementation of the class Inter-nalEnergy which derives from pHS.Storage.C according to Definition 2.  (temperature), p1 (pressure), and g1 (specific Gibbs free energy) are sources of effort; Sfr1 (entropy flow rate), Vfr1 (volume flow rate), and mfr1 (mass flow rate) are sources of flow; U1 is a storage element for internal energy; TB1 T |Ṡ , pB1 p|V , and mB1 g|ṁ are one-dimensional bonds; thB1 (T , p, g)|(Ṡ,V ,ṁ) is a three-dimensional bond; TS1 (temperature), pS1 (pressure), and VfrS1 (volume flow rate) are ideal sensor. These symbols, together with those proposed in [28] for the elementary Dirac structures, facilitate the composition of complex systems (see Figure 5) with the Modelica drag-and-drop graphical user interface.

D. MODELING THE BEHAVIOR OF SINGLE SUBSTANCES
In the implementation of the class InternalEnergy, the instance substance of the abstract class SingleSubstance models the thermodynamic behavior of any simple substance, as described in its fundamental equation. The class Inter-nalEnergy encodes the structure of the storage element from the point of view of the port-Hamiltonian approach, and the class SingleSubstance provides information about the behavior of the particular substance used to define the storage element.
Listing 5 shows the Modelica code used to implement the class SingleSubstance, which is defined in the package Substances.Partial included in the package Thermodynamics. The assertions in the equation section are intended to stop the simulation when the substance temperature, pressure, or density takes values outside the valid range established for the fundamental equation. The security ranges are: T min ≤ T ≤ T max , p min ≤ p ≤ p max , and ≥ min .
SingleSubstance is a partial class (or abstract class in the terminology of object-oriented modeling). It is a template that can be used to define derived classes for modeling particular substances or families of substances with common features. In the following, we illustrate the applicability of the developed tools by modeling some ideal thermodynamic systems.

E. IDEAL THERMODYNAMIC SYSTEMS
In statistical mechanics, some idealized systems have been defined for which it is possible to calculate their fundamental equation. Although these models may seem too simple, they are sometimes a good first approximation to actual physical systems.
The main characteristic of the fundamental equation of an ideal thermodynamic system -in either of its two forms: internal energy or entropy, and expressed in specific quantities-is the separability of its variables [48,Chap.11].

1) IDEAL GASES
For ideal gases, the fundamental equation in the specific entropy representation is [6, Appendix D] where R = R mol /M mol J/(kg·K) is the specific gas constant of the substance, R mol = 8.314472(15) J/(mol·K) is the molar gas constant [31], and (s 0 , u 0 , v 0 ) are values of a reference state. The function φ depends on the characteristics of each gas and satisfies the condition φ(u o ) = 0. The intensive parameters of an ideal gas are (see (13)) Listing 5. Thermodynamics.Substances.Partial.SingleSubstance.
From (40) and (41), we can deduce some values of the response functions common to all ideal gases. In particular, The Modelica implementation of the ideal gases behavior has resulted in the class IdealGas (see Listing 6), defined in the package Thermodynamics.Substances.Partial. We implement only equations (41) and (47), that are common to all ideal gases and depend on the temperature.

2) PERFECT GASES
For perfect gases, c(T ) = c is constant, and this implies c v = d(Tc)/dT = c is also constant. This feature allows to obtain explicit expressions for the fundamental equation and the intensive parameters. Indeed, Remark 3: Another way to obtain (55) is to calculate g from its definition g = ∂U /∂m. The internal energy of an ideal gas is If we take into account equations (49) (55)).
Listing 7 shows the implementation of the class PerfectGas. Here we implement equations (49) and (51), i.e., we consider the fundamental equation in the entropy representation. Note that PerfectGas inherits equations (41) and (47) from its parent class IdealGas.

3) MONATOMIC PERFECT GASES
For monatomic perfect gases c v = 3 2 R, c p = 5 2 R, and we can write The implementation of the class MonatomicPerfectGas is shown in Listing 8.

4) EMPIRICAL MODELS FOR NON-PERFECT SUBSTANCES
The perfect gas condition (constant c p ) may be suitable for modeling the behavior of polyatomic gases or liquids in a first approximation. However, to obtain more accurate results, we must take into account that c p varies with temperature. VOLUME 9, 2021 Given the difficulty of theoretically deducing this dependence from first principles, we have to employ a phenomenological approach inferred from experimental measurements. It is usual to opt for a dimensionless polynomial representation c p (T )/R = i∈I ⊂Z a i T i because it simplifies calculations and facilitates symbolic manipulation. The number of coefficients a i comes from a tradeoff between accuracy and computational cost.
For example, the library of thermodynamic data used with the CEA NASA Glenn computer program contains data for over 2000 substances modeled with seven coefficients for c p in different temperature ranges from 200 to 20000 K [30]. The expressions to calculate the specific heat capacity, specific entropy, and specific enthalpy with these coefficients are as follows: where b 1 and b 2 are integrations constants. Each substance can have several sets of coefficients (a 1...7 , b 1,2 ) grouped by temperature intervals depending on the phase (condensed or non-condensed) adopted in each interval.
To incorporate the NASA Glenn database into our Modelica implementation, we have processed the text file containing the above-mentioned coefficients [30,Appendix D] and generated the Modelica package NASAGlennCoefficients. This package contains the partial class Partial.Substance (see Listing 9) derived from class Thermodynamics.Substances. Partial.IdealGas and implements the equations (61)-(63) taking into account the temperature intervals.
The rest of the package NASAGlennCoefficients is composed of 2073 substance models. Each model is an extension of the previous class Substance.
Example 4: As an example, we show the classes modeling the thermodynamic behavior of CO 2 (Listing 10) and air (Listing 11). The latter is a mixture of substances, but in the NASA Glenn library, it is also modeled as an ideal gas. As classes CO2 and Air extend the class Substance, they inherit equations (61)-(63), where the thermodynamic behavior of all substances modeled in the package NASAGlennCoefficients is defined. It is only necessary to initialize the parameters involved in these equations to complete the definition of the derived classes.

F. IAPWS FORMULATION 1995 FOR THE THERMODYNAMIC PROPERTIES OF ORDINARY WATER
Water is a substance widely studied by the scientific community, not only for its role in many natural processes, some of them so important as the origin of life or the evolution of Earth's climate, but also for the large number and variety of its industrial applications. For these reasons, the scientific community has endeavored to find the fundamental equation that models its thermodynamic behavior, and as a result of this effort, the International Association for the Properties of Water and Steam has published two empirical formulations of it, known as IAPWS-95 [22], [53] and IAPWS-IF97 [21], [52]. IAPWS-95 is the most accurate and is intended for general scientific research. IAPWS-IF97 is a simplification of the previous one and is intended for industrial applications.
Although it is less accurate, it has the advantage of consuming less computational resources when performing simulations.
The Modelica Standard Library implements the IAPWS-IF97 model in the package Modelica.Media.Water.Wa-terIF97. We have implemented the IAPWS-95 model in the class Thermodynamics.Substances.WaterIAPWS95. This model formulates the fundamental equation of water in terms of the specific Helmholtz free energy f (T , v) = f (T , 1/ ). The model has two dimensionless terms: φ o (δ, τ ) for the ideal-gas part and φ r (δ, τ ) for the residual part: The ideal-gas part of (64) is and we can see that it fulfills the condition of being an separate variables equation (see Section IV-E). The coefficients n o 1...8 and γ o 4...8 can be consulted in [22, Table 1]. The residual part of (64) is  Table 2].
Since (64) is a fundamental equation, it contains all the information about the thermodynamic behavior of water. Therefore, all other thermodynamic properties and quantities can be reduced from f . In our implementation of the IAPWS-95 model, we are also interested in calculating the specific entropy s, pressure p, the isobaric heat capacity c p , and the isochoric heat capacity c v :  Second order partial derivatives φ o τ τ , φ r τ τ , φ r δτ , and φ r δδ (adapted from [22, Table 5]). Table 1 shows the partial derivatives φ o τ , φ r τ , and φ r δ used in (67)-(68). The second partial derivatives φ o τ τ , φ r τ τ , φ r δτ , and φ r δδ used in (69)-(70) are shown in Table 2. Listing 12 is the implementation of the Modelica class WaterIAPWS95. This implementation has been validated with the computer-program verification test values given in [22, Tables 6 and 7].

G. EMPIRICAL MODELS FOR GENERAL SUBSTANCES
We have already seen (Figure 1) that the port-Hamiltonian model of a simple thermodynamic system requires to know the covector field dU (see (4)). Although U and dU are not directly within our reach, we can approximate them numerically from the primary set of response functions (Section II-B2). As these properties are measured and tabulated as a function of temperature and pressure, it is more convenient to work with Gibbs potential G (T , p, m) or specific Gibbs potential G/m = g(T , p) rather than internal energy U (S, V , m). Thus, we can write the following equations: ds = ∂s ∂T p dT + ∂s ∂p T dp = c p T dT − vα p dp, 3 dv = ∂v ∂T p dT + ∂v ∂p T dp = vα p dT − vκ T dp, dg = ∂g ∂T p dT + ∂g ∂p T dp = −sdT + vdp, and in the CIT approacḣ For practical considerations, it is useful to relate α p and κ T to the density = m/V = 1/v because it is easier to tabulate (T , p) than α p or κ T : Listing 13. Thermodynamics.Substances.Sylterm800.
The class Syltherm800 derived from SympleSubstance and models the thermodynamic behavior of Syltherm 800 (see Listing 13).

V. EXAMPLES
The examples in this section have been processed and simulated with Dymola running in a Linux virtual machine, FIGURE 5. Implementation with the library pHlib of the two-tank system previously modeled with MSL (see Figure 4).
which was installed on an iMac computer configured with a 2.66 GHz Intel Core 2 Duo processor and 8 GB of memory.
A. COMPARING pHlib AND MSL FOR ACCURACY Figure 4 shows the MSL implementation of a system consisting of two open tanks at ambient conditions (20 • C, 1 atm) connected through a thermal conductor G 1 = 300 W/K. The first tank contains 1 m 3 of water modeled with the MSL IAPWS-IF97 implementation for industrial use [21] and it has thermal losses with its surroundings modeled with the thermal conductor G 2 = 50 W/K. The second tank also contains 1 m 3 , in this case of Syltherm 800, modeled with the MSL library from data supplied by the manufacturer [15] and it has thermal losses with its surroundings modeled with the thermal conductor G 3 = 150 W/K. The water tank is heated for 60 hours with the heat generated by the electrical resistor R 1 = 32 and then allowed to cool for another 60 hours. For this purpose, we connect a modulated voltage generator V 1 = 400 V to R 1 and make the assumptions: a) all the electrical power transformed into heat is transferred to the water tank, b) we do not consider the variations that R 1 may have with temperature. We use thermal conductors and ideal resistors to focus our study on the modeling of the thermodynamic behavior of water and Syltherm substances. Figure 5 shows the implementation of the same system with the library pHlib where the proposed port-Hamiltonian representation makes it easier to visualize the connection network for energy exchange between the different components of the system.
In this figure we can see the following components (see Figure 3 and [28, Figures 9 and 12] for a detailed explanation of the symbols): • The tank water modeled by: the storage element IAPWS95, the boundary conditions modeled by the constant temperature source T1, the constant pressure source p1 and the constant mass flow rate source mfr1, the three-port slicer slc1 that separates the 3-dimension thermodynamic port #4 into three 1-dimension ports (thermal port #1, pressure port #2, and mass port #3), the thermal conductance GS2 thermal modeling losses to the environment, the temperature sensor TS1, and the 0-junction J0_1.
• The Syltherm 800 tank, with a similar structure to the water tank, and modeled by the storage element Syltherm.
• The thermal conductance GS1 connecting both tanks.
• The modulated voltage generator V1 connected to the electrical resistor RS1 for heating the water tank.
The thermodynamic behavior of water has been modeled with the class WaterIAPWS95 (see Section IV-F) that implements the equations for general and scientific use [22], and the behavior of Syltherm with the class Sytherm800 (see Section IV-G). The thermal conductors GS 1,2,3 and the electrical resistor RS 1 have been modeled as entropy generating elements derived from the class pHlib.Dissipative.RS (see [28,Remark 26] and Listing 14).
The implementation of RS satisfaying that the equation e|f + T 2 |Ṡ 2 = 0, guarantees the conservation of energy during the heat conduction process. Thus, the incoming power P 1 = e|f is equal to the outgoing heat flow ratė U 2 = T 2 |Ṡ 2 . However, in the transformation process, there VOLUME 9, 2021 FIGURE 6. Evolution of water and Syltherm temperature during a 120-hour simulation for MSL (a) and pHlib (b) models, and evolution of the relative error (c).
will be an entropy generation to be modeled in each particular class derived from RS.
The heat flow rate between the ends of a thermal conductor of conductance G at temperatures T 1 and T 2 iṡ and the equations for calculating the entropy flow rateṠ cond generated in the conduction process are [4, Chap.2] These equations are implemented in the classes MGS and GS (see Listing 15).
We consider that all the electrical energy consumed by an ideal linear resistor of resistance R is transferred to the thermal port, so the entropy flow rate in this port iṡ In part a) of Figures 6-8, we can see the temperature, volume, and density evolution of water and Syltherm during the 120-hour simulation of the MSL model. In part b) we can see the evolution of the same system modeled with pHlib.
We can see that both models show qualitatively similar behavior. For a quantitative comparison, we have taken the MSL model as a reference and calculated the relative errors of the pHlib model with the expression In part c) of Figures 6-8, we can see the time evolution of the relative error. For water, the errors are below 0.03% in temperature and below 0.005% in volume and density. For Syltherm, the maximum error is only slightly above 0.01% in all three quantities.   thermal conductors. Figure 10 shows the model of each tank. The system parameters are the number N of tanks, the mass of water stored in each tank, its initial volume, the thermal conductance between tanks tankToTankG, and the thermal conductance with the environment tankToAmbG. The model used for water is IAPWS-95. Resistor R1 connected to the modulated voltage generator V1 heats the   first tank for 60 hours, and then we let the tanks cool for another 60 hours. This system has also been implemented with MSL and the IAPWS-IF97 water model. Both models generate 2 × N state variables. Table 3 shows the execution times when simulating each model for some values of N . The results show that simulations  with the IAPWS-IF97 model are faster than simulations with the IAPWS-95 model. These results are expected because the industrial formulation simplifies the scientific formulation to increase computational speed in a tradeoff for lower accuracy (for more details on this subject see [21], [52]).
These results might seem to suggest that IAPWS-95 models are not very suitable for real-time applications. However, note that the 120-hour simulation for the 50 tanks system took only 9.35 seconds. This system produces a model with 100 state variables, which is common in modeling complex industrial systems. For this reason, we think that improvements in hardware performance over the last two decades make it possible to use the IAPWS-95 formulation in real-time applications. The IAPWS-95 model also has the advantage of formulating the fundamental equation of water with a single continuous and differentiable expression, unlike IAPWS-IF97, which is forced to split the formulation into FIGURE 12. Evolution of water temperature in tanks 11, 13, 15, 33 and 55 during a 120-hour simulation (see Figure 11). five regions, each with its corresponding fundamental equation. This partitioning complicates the construction of models with region changes and it can even lead to convergence and stability problems at the region boundaries due to chattering phenomena (see [3] for details).
C. ARRAY OF 5 × 5 TANKS Figure 11 shows a system consisting of an array with 5 open water tanks at ambient conditions (20 • C, 1 atm) connected through thermal conductors. Each tank has 250 kg of water, 10 W/k of thermal conductance with the environment, and 250 W/K of thermal conductance with each of its neighbors. Resistors R1 and R2 heat the system for 120 hours. V1 is a voltage generator that generates 400 V during the first 59 hours, and V2 generates 400 V during the last 59 hours. Therefore, during the two central hours, the resistors do not supply heat to the system. Figure 12 shows the temperature evolution of some tanks that have been chosen to check the consistency of the model. Given the symmetry of the system, some temperatures have to converge to the same value as is the case for the temperatures of tanks 11 and 55. Other temperatures have to maintain a relationship of order, e.g., the temperature of tank 13 has to be higher than that of tank 33 because it is closer to the heat sources.

VI. CONCLUSION AND FUTURE WORKS
This work is a continuation of [28] where we presented a general purpose Modelica library for modeling multiphysics systems with a port-Hamiltonian approach. In the conclusions of that article, we proposed the extension of this library as a future research line for ''modeling of thermofluid systems and its assessment in the development of models for largescale solar thermal power plants' '. In the current article, we have extended the library with the proposed line of research and presented the added components for modeling thermodynamic systems. This is a first step in approaching the modeling of thermofluid systems. We will complete the modeling and the library in a future paper studying the transport phenomena in flowing substances.
The object-oriented design of the library (taking advantage of the resources offered by the Modelica language) has allowed us to build general substance models and to refine them in a simple way to model specific substances such as water or Syltherm 800.
The proposed symbols which accompany each component allow the construction of complex system models that are easy to read for any designer familiar with bond graphs.
When contrasting the simulation results of the two-tank model with the results obtained for the same model implemented with MSL, we have found that they coincide with relative errors lower than 0.005% in some variables, even though both libraries are designed with different methodologies. This result validates the choice of the port-Hamiltonian approach and Classical Irreversible Thermodynamics to model the behavior of thermodynamic systems. The suitability of the port-Hamiltonian formulation as a multiphysics modeling tool of great theoretical and practical value has been illustrated, according to the following aspects: the mathematical rigor of the port-Hamiltonian formalism built on differential geometry, the fundamental role played by energy in the definition of port-Hamiltonian systems, and the ease of modeling complex systems by connecting simpler ones that exchange energy with each other.
The performance measurements evaluated on the N tanks example show that the system modeled with IAPWS-IF97 (MSL library) is faster than the one modeled with IAPWS-95 (pHLib library). The pHlib library, due to its more academic port-Hamiltonian approach, will be mainly useful in scientific research. However, the response times measured for different values of N indicate that modern computers can cope with real-time simulations of large systems in industrial applications. This result encourages us to continue with the research initiated in [28] to complete the library pHlib with the modeling of transport phenomena in thermofluid systems and to study its application to the modeling of complex systems such as solar thermal power plants.

APPENDIX ACRONYMS AND NOMENCLATURE
A. ACRONYMS PEDRO J. ZUFIRIA received the Ingeniero de Telecomunicación degree from the Universidad Politécnica de Madrid (UPM), in 1986, the M.Sc. degree in ME, the M.Sc. degree in EE, the Ph.D. degree from the University of Southern California (USC), in 1989, the Doctor Ingeniero de Telecomunicación degree from MEC, in 1991, and the Licenciado en Ciencias Matemáticas degree from the Universidad Complutense de Madrid, in 1997.
From 1987 to 1990, he was a Teaching and Research Assistant with USC, in collaboration with different projects of the National Science Foundation and TRW/NASA. Since 1991, he has been holding a teaching position (currently as a Full Professor) at the Escuela Técnica Superior de Ingenieros de Telecomunicación (ETSIT), UPM, where he was the Chairman of the Departamento de Matemática Aplicada a las Tecnologías de la Información y las Comunicaciones, from 1999 to 2004. In the ETSIT, he was responsible for the Office of International Relations with other academic institutions, from 1993 to 1996. He was the Vice-Dean of the research and graduate studies, in 1999, and the Co-Director of the Orange Chair, from 2009 to 2019. Since 2019, he has been the Co-Director of Cabify Chair. He also has been the Founder and a Scientist responsible for the Neural Networks Group and currently the Dynamical Systems, Learning, and Control Research Group, where he has been involved in research projects of the European Union as well as several Spanish official and private institutions. He has authored over 100 international publications in his research field. His research interests include analysis, control, and fault diagnosis of dynamical systems, and the applications of statistics, machine learning paradigms, and complex network theory on system modeling and data processing. LUIS J. YEBRA was born in Almería, Spain, in 1971. He received the degree in telecommunications technical engineering from the Universidad de Alcalá, Alcalá de Henares, Spain, in June 1993, the degree in physics from the Universidad Nacional de Educación a Distancia (UNED), Madrid, Spain, in June 1997, and the Ph.D. degree in automatic control and industrial computing from the Department of Automatic Control and Computer Science, UNED, in May 2006. From 1999 to 2007, he was a Predoctoral Fellow with the National Research Center, CIEMAT, and a Researcher with the Department of Energy, CIEMAT in Plataforma Solar de Almería (PSA-CIEMAT). Since 2007, he has been a Tenured Scientist in dynamic modeling and automatic control with PSA-CIEMAT, where he has been mainly focused on concentrated solar thermal plants (CST). He has coauthored four books, more than 30 scientific articles published in international journals (JCR), and more than 70 international conference proceeding papers. He has supervised four Ph.D. degree students. He has participated in more than 30 publicly funded research projects, two technology transfer agreements, and the creation of a spin-off company.
Dr. Yebra has received three prizes to the best academic curriculum for his Physics degree. VOLUME 9, 2021