Thermodynamical Material Networks for Modeling, Planning, and Control of Circular Material Flows

Waste production, carbon dioxide atmospheric accumulation, and dependence on finite natural resources are expressions of the unsustainability of the current industrial networks that supply fuels, energy, and manufacturing products. In particular, circular manufacturing supply chains and carbon control networks are urgently needed. To model and design these and, in general, any material networks, we propose to generalize the approach used for traditional networks such as water and thermal power systems by using compartmental dynamical thermodynamics and graph theory. The key idea is that the thermodynamic compartments and their connections can be added, removed or modified as needed to achieve a circular material flow. The design methodology is explained and its application is illustrated through examples. In addition, we provide a physics-based definition of circularity and, by implementing a nonlinear compartmental control, we strengthen the connection between"Industry 4.0"and"Sustainability". The paper source code is publicly available.


I. INTRODUCTION
While the human population is predicted to reach 8 billion by 2024 and 10 billion by 2056 [1], modern society strives to provide the needed services and products on a largescale.Any services and products, from the simple piece of paper to the complex graphical processing unit, require the availability of raw materials and the production of energy.As welfare and economic growth rely on material uses, a long-term sustainable management of finite natural resources is increasing in importance [2].
A natural resource particularly at risk is climate stability, which is mostly being altered by the atmospheric carbon dioxide concentration [3].The world emitted 6 billion tonnes of CO 2 in 1950, 22 billion tonnes in 1990, 36 billion tonnes in 2019 (i.e., 6 times the emissions of 1950) and the annual emissions have yet to reach their peak [4].Therefore, an effective control of this material is necessary to respect the global warming target of 1.5 degrees Celsius compared to the pre-industrial levels established with the Paris Agreement in 2015 [5].
Along with carbon dioxide, other materials requiring a more efficient management are those accumulating on lands and seas as litter or marine debris such as plastic.For example, the mass of plastic in the Great Pacific Garbage Patch was estimated to be approximately 80,000 tonnes and the mass of plastic entering the ocean each year is 1.15 to 2.41 million metric tonnes [6].In terms of the life-cycle of a material, the status of "waste" is at the final stage, hence waste accumulations are issues related to the end of the life-cycle.Similarly, today there are also increasing concerns at the beginning of the material life-cycle, i.e., at the stage of material extraction.Indeed, there are several materials classified as "critical" by the European Union [7] and the United States [8] whose supply is particularly at risk.Those materials are irreplaceable in clean technologies such as solar panels, wind turbines, electric vehicles and are also used in modern technologies such as smartphones.
To address the issues located both at the beginning and at the end of the life-cycle of materials, the paradigm of "circular economy" has gained attention over the last few years.Currently there are multiple definitions of circular economy [9].In this paper, we focus on the flows and accumulations of materials (e.g., carbon dioxide, gold, plastic, biomethane), and hence for us the adjective "circular" means "closed flow of material".As a consequence, the expression "circular economy" is equivalent to "economy based on closed flows of materials", the expression "measuring the circularity of a supply network" is equivalent to "measuring to what extent the flow of material in a supply network is closed" and the expression "circulating a material" means "closing the flow of a material".For example, hydraulic engineers seek to circulate the water by minimizing the leakages in the water network.
To enhance the modeling, planning, and control of circular material flows, initially we looked at the advanced and mature water industry and asked the question: Given that water is just a particular type of material, can we develop a capability of managing other materials as effective as the one we have with water?Then, we observed two key aspects of the design of water networks.Namely, 1) they are designed to be closed in order to minimize leakages of material and a mathematical framework arXiv:2111.10693v3[math.DS] 14 Nov 2022 that effectively depicts the network architecture is graph theory [10]; and 2) their modeling is based on the first principle of thermodynamics and the mass conservation equation [11].
Given the generality of both thermodynamics [12] and graph theory [13], in this paper we propose to extend the modeling approach of water supply networks to model the flow of any material leveraging compartmental dynamical thermodynamics [14] and graph theory [13].
The main contributions of the paper are the following.
1) We provide physics-based foundations of material circularity to add clarity to the topic [9] (see Section III).2) We propose a systematic methodology to model, plan, and control circular flows of materials (see Section IV). 3) We illustrate the use of graphs to measure material circularity (see Section V). 4) We illustrate the use of feedback control systems into the design of material flows (see Section V).By doing this, we strengthen the link between techniques typical of industrial automation and the holistic perspective of industrial ecology required to design closed-loop flows [15].
The paper is organized as follows: Section II covers related works, Section III defines the key concepts, Section IV details the proposed methodology, Section V provides two examples, and finally Section VI concludes.
Throughout the paper, vectors and matrices are indicated with bold letters, whereas sets are indicated with calligraphic letters.

II. RELATED WORK
Compartmental and dynamical thermodynamics: Compartmental thermodynamics refers to the thermodynamic analysis at equilibrium of a set of connected machines.An example of compartmental thermodynamics is the simple Rankine cycle, in which a turbine, a boiler, a water pump, and a condenser are connected through pipes.Its invention dates back to 1859.Dynamical thermodynamics, instead, is an emerging topic as it studies the dynamical (i.e., non-equilibrium) behavior of systems from a thermodynamic perspective without a focus on the multi-machine nature typical of compartmental thermodynamics.Examples of works are [16] for electrical networks, [17] for electronic circuits, [18] for mechanical systems, and [19], [20] for chemical reaction networks.The combination of these two branches of thermodynamics is compartmental dynamical thermodynamics [14].
Graph theory and thermodynamics for sustainability: In the context of circular economy, some graph-based approaches have been proposed recently.Moktadir et al. [21] used a graph architecture to examine and prioritize the driving factors of sustainable manufacturing practices; in [22] and [23] a graph architecture is used to analyze the different barriers to the implementation of a circular economy in the mining industry and in a biomass supply chain, respectively.The work of Gribaudo et al. [24] proposed the use of graphs to model the production of chitin by bio-conversion of municipal waste.
The idea of using thermodynamics for ecological modeling dates back almost thirty years [25].In 2011, the authors of [26] further extended this vision by proposing thermodynamics as the science of sustainability.With this paper, we aim at clarifying the application of thermodynamic principles for sustainable and circular design.
Material flow analysis: Material flow analysis (MFA) is one of the key techniques developed and used in industrial ecology and circular economy to assess material flows and stocks in urban and natural environments [27].The holistic perspective at the core of MFA had a strong influence on the methodology proposed in this paper.While MFA is mainly based on mass balances [28], our methodology extends it by adding dynamical power balances [14] and control systems [29].
Control systems in life-cycles: Nowadays, control theory is a well-established discipline sitting between applied mathematics and engineering.Two selected works from the vast literature in the field are [30] for linear feedback control and [29] for advanced non-linear methods.Control systems are distributed across the entire life-cycle of products and services, therefore they can play a critical role in the transition from a linear to a circular economy.In this paper, we illustrate the use of control theory into the design of material flows, specifically for bio-methane production.

III. CIRCULARITY AND THERMODYNAMICAL MATERIAL NETWORKS
Consider a cube of material β, infinitesimal mass dm, density ρ, and volume dV = dm ρ as in Fig. 1.Let G(t) = [x G (t), y G (t), z G (t)] be the center of mass whose coordinates are written with respect to a fixed reference frame with origin O = [0, 0, 0] .In general, x G (t), y G (t), and z G (t) can vary with the time t.Let p(t) = G(t) − O be the position vector of the cube center of mass G.
Definition 1 (Mechanics-based circularity).The flow of β is mechanically circular if there exist t 0 ≥ 0 and t * ∈ (t 0 , ∞) such that p(t 0 ) = p(t * ), ṗ(t 0 ) = 0, t * > t 0 . ( Remark 1.As the material β in Definition 1 is fixed, chemical reactions that modify the material are excluded.Therefore, we refer to Definition 1 as the mechanics-based definition of circularity. In thermal engineering, it is standard practice to define a control volume that contains the system under study before the application of mass and energy balances.Such a standard practice underpins the design of thermodynamic cycles, e.g., the Rankine and the Brayton cycles, and also the design of hydraulic networks [11].Each control volume identifies a thermodynamic compartment.For example, a simple ideal Rankine cycle is made of eight thermodynamic compartments: a feedwater pump, a boiler, a turbine, a condenser, and four pipes connecting these four machines into a closed-loop.Now note that mass and energy balances, that is, thermodynamics, are general principles valid for any system [12], [26].Hence, we can generalize the definition of circularity based on mechanics (Definition 1) with the following thermodynamicsbased definition.Let c k be the k-th thermodynamic compartment identified by the k-th control volume and let β be the material of an infinitesimal cube as in Fig. 1.
Definition 2 (Thermodynamics-based circularity).The flow of β is thermodynamically circular if there exists an ordered sequence of compartments φ = (c 1 , . . ., c k , . . ., c 1 ) processing β, which begins and ends in c 1 .Moreover, if some c k ∈ φ chemically transforms a material β 1 into a material β 2 and there exists an ordered sequence φ processing B = {β 1 , β 2 }, then the flow of the material set B = {β 1 , β 2 } is thermodynamically circular.More generally, the flow of B = {β 1 , . . ., β q , . . ., β n β } is thermodynamically circular if there exists an ordered sequence φ processing B. A well-established formalism to represent a network of systems is graph theory [13]: examples of network design theories based on graphs are electrical networks [16], hydraulic networks [10], and multiagent systems [31].Since the system in Fig. 2 can be seen as a network of thermodynamic compartments connected through the material flow, we will use graph theory to formulate the system in Fig. 2 as a network and then state the definition of a thermodynamical material network.

Definition 3 ([13]
).A directed graph D or digraph is a graph identified by a set of n v vertices {v 1 , v 2 , . . ., v nv } and a set of n a arcs {a 1 , a 2 , . . ., a na } that connect the vertices.A digraph D in which each vertex or arc is associated with a weight is a weighted digraph.
Let c k i,j be the k-th thermodynamic compartment through which the material flows from compartment i to compartment j.
Specifically, N = R ∪ T , where R ⊆ N is the subset of compartments that store, transform, or use the target material, and T ⊆ N is the subset of compartments that move the target material between the compartments belonging to R ⊆ N .A net N is associated with its weighted mass-flow digraph M (N ), which is a weighted digraph whose vertices are the compartments c k i,j ∈ R and whose arcs are the compartments c k i,j ∈ T .A vertex also results from the intersection of 3 or more arcs.For vertex-compartments c k i,j ∈ R it holds that i = j, whereas for arc-compartments c k i,j ∈ T it holds that i = j.The weight assigned to a vertex-compartment c k i,j ∈ R is identified by the mass stock m k within the corresponding compartment, whereas the weight assigned to an arc-compartment c k i,j ∈ T is the mass flow rate ṁi,j from the vertex-compartment c i i,i ∈ R to the vertex-compartment c j j,j ∈ R. Hence, the orientation of an arc-compartment c k i,j ∈ T is given by the direction of the material flow.The superscripts k v and k a in (2) are the k-th vertex and the kth arc, respectively, while n c and n v are the total number of compartments and vertices, respectively.Since n a is the total number of arcs, it holds that  2) is a weighted digraph with arcs and vertices labeled with the corresponding compartmental nomenclature c k i,j .Figure 3 shows an example of Next, we introduce a few more definitions from graph theory.The reason will be clarified afterwards.

Definition 7 ([13]
).A directed walk in D is a finite nonnull sequence W = (v 0 , a 1 , v 1 , a 2 , . . ., a l , v l ) whose terms are alternatively vertices and arcs such that, for i = 1, 2, . . ., l, the arc a i has head v i and tail v i−1 .The integer l is the length of W , while the vertices v 0 and v l are the origin and the terminus of W , respectively.Definition 8 ( [13]).If the sequence of arcs a 1 , a 2 , . . ., a l of a directed walk W are distinct, then W is a directed trial.

Definition 9 ([13]
).A directed trial is closed if it has positive length and its origin and terminus are the same, i.e., v 0 = v l .

Definition 10 ([13]
).A closed directed trial whose origin and internal vertices are distinct is a directed cycle φ.
Summarizing, a directed cycle φ is a directed walk W (Definition 7) in which the arcs are distinct, the origin and the internal vertices are distinct, the origin and the terminus are the same, and l > 0.
The reason for introducing these definitions is that it is now apparent that the requirement for material circularity (1) translates into requiring that the mass-flow digraph M (N ) must be a directed cycle φ.
Remark 2. Consider the network (2) with a mass-flow digraph M (N ).Then, the flow of the set of materials The mass-flow matrix Γ(N ) associated with the network (2) is given by whose entries along the diagonal are the weights of the vertex-compartments c k i,j ∈ R (i.e., mass stocks) and whose off diagonal entries are the weights of the arc-compartments c k i,j ∈ T (i.e., mass flow rates).Hence, Γ(N ) is a square matrix of size n v × n v with nonnegative real entries, i.e., Γ ∈ R nv×nv + .
The mass conservation principle [11] establishes the relationship between the entries of Γ(N ), namely, which can be further written in terms of the entries of Γ(N ) as or, equivalently, in vector form as Therefore, Remark 3. In line with the standard nomenclature adopted in thermal engineering, the mass accumulation or depletion in vertex-compartments is denoted as d dt m and not as ṁ [11].Indeed, the latter indicates a mass flow and involves a mass transfer between two vertex-compartments.The first form will be referred to as mass accumulation-depletion, whereas the second form will be referred to as mass flow rate.The SI unit is kg/s for both quantities.
Leveraging the mass-flow matrix (3), we now define the circularity indicator λ.
where n φ is the number of directed cycles in M , Q = {γ i,j | ṁi,j does not belong to any directed cycle}, (9) and is the cycle mean of φ, with The next section proposes a methodology for the design of circular material flows using the definitions given above.

IV. DESIGN METHODOLOGY
The proposed methodology is outlined in Fig. 4 and involves three main steps; its output is a TMN, its goal is designing the flow of the material set B = {β 1 , . . ., β q , . . ., β n β }, and it has two optional extensions at Steps 2 and 3 indicated by the black arrows.Specifically, the first step is the choice of the material set B to be circulated.Then, the network (2) is defined and depicted as needed using a compartmental diagram, a compartmental digraph and/or a mass-flow digraph (an example is in Fig. 3).Moreover, its circularity can be measured by computing the indicator (8) as indicated with the black arrow on the left-hand side, which requires to preliminarily define the mass-flow matrix (3).The third and last step consists of applying to each compartment c k i,j ∈ N the dynamical form of the compartmental energy balance, that is, and/or the compartmental mass balance (4).Equation ( 12) is a power balance in which E k , K k , U k , and P k are the total energy, the kinetic energy, the internal energy, and the potential energy of the compartment k, respectively, and Qk and Ẇk are the heat flow and the work flow exchanged by the compartment k with the surroundings.
Once the dynamics of the compartments are defined, it is possible to implement a compartmental control system as indicated by the black arrow in Step 3 if the power balance (12) and the mass balance (4) are written in state-space form where x k ∈ R n k is the state vector of the k-th compartment, u k ∈ R z k is the control input of the k-th compartment, and

B. Subsystem of Bio-Methane Supply Chain
This example demonstrates the application of the design methodology on a subsystem of a biomethane supply chain.Specifically, the subsystem involves three stages of the biomass life-cycle (see Fig. 7): the biomass hub, the truck to transport the biomass, and the anaerobic digestion plant to covert the biomass into biogas.Moreover, the anaerobic digestion plant is divided into two sub-compartments: P r is the plant reservoir of biomass and P d is the plant digester.We assume that G(t) = [x G (t), 0, 0] (i.e., the material motion is along x only), the position of the hub is x h , O ≡ x h , x p is the position of the chemical plant, and that the sizes of the hub and the plant are negligible compared to H = x p − x h .The first step of the methodology requires to choose the set of materials of interest, i.e., B; in this case, B = {β 1 }, where β 1 is the biomass (the details of its chemical composition are not considered in this example).The second step requires to define the net (2); for this example, N = {c 1  1,1 , c 2 2,2 , c 3 1,2 }, n c = 3, n v = 2, and n a = 1.The thermodynamic modeling of each compartment as required at Step 3 is as follows.
1) c 1 1,1 (Biomass hub): The process of exiting the solid biomass from the hub can be modeled more naturally as a discrete-time system rather than a continuous-time system since the material output flow is carried on in batches instead of as a continuous flow (as it would be with fluids).The discrete-time mass balance for the hub yields where Z + is the set of nonnegative integers, m 1 is the mass stock inside the hub, m l is the truckload, n l is the loading time, and δ n l (n) is the Kronecker delta, that is, 2) c 3 1,2 (Truck): For n l ≤ n ≤ n u , the biomass with mass m l is on the truck, it is transported to the plant, and it enters the plant at the unloading time n u .In this example, we model the truck as the three-wheel vehicle shown in Fig. 8 [32].The equations of motion of the vehicle can be derived using Lagrange's equations of motion ( [33]) where G Fig. 8: A three-wheel vehicle model of the truck, that is, compartment c 3 1,2 (modified from [32]).
is the Lagrangian function of the mechanical system, q ∈ R d is the vector of generalized coordinates, q ∈ R d is the vector of generalized velocities, and ξ ∈ R d is the vector of generalized forces associated with the generalized coordinates q.Before we use the Lagrangian formulation (20) to define the dynamical equations of the three-wheel vehicle, we demonstrate that the Lagrangian formulation ( 20) can be derived from the dynamical form of the energy balance (12).This is a key result for this paper as it shows that Lagrangian mechanics respects Step 3 of the proposed methodology, which requires that we model each compartment using the dynamical form of the energy and/or the mass balances.
Proposition 1. Lagrange's equations of motion given by ( 20) can be derived from the dynamical form of the energy balance (12).
Proof.For simplicity of exposition, consider the dynamical form of the energy balance (12) without the subscript k specifying the k-th compartment, that is, Since the heat flow Q and the internal energy U are neglected in solid mechanics, (22) reduces to Note that in solid mechanics the potential energy P corresponds to the conservative work done by the gravitational force, and hence ( [33]), Taking the partial derivative on both sides of ( 23) with respect to q and transposing the resulting equation yields d dt It now follows from (24) that the term on the left-hand side of ( 25) can be written in terms of the Lagrangian function (21), and hence, (25) can be rewritten as Now, note that the term on the left-hand side of ( 26) is equal to the term on the left-hand side of the Lagrangian formulation (20), and hence, Therefore, ( 26) can be written as (20).
For the three-wheel vehicle (20) (which can be derived from the power balance (22) as shown in Proposition 1 and thus it respects Step 3 of our methodology) yields [32] where τ 1 and τ 2 are the control torques of the engines M 1 and M 2 , respectively, and where l 2 , and I z is the principal moment of inertia of the vehicle with respect to its z-axis (note the dependence of the truck dynamics (28) on the truckload m l ).Once θ1 and θ2 are determined by integration of (28), ψ and θ3 are given by where with δ ± = sin ψ 2 ± d l , α = a+b l , and ρ = r d .As mentioned above, here we assume that the material motion is along the x-axis only, and hence, ψ(t) = 0, θ 1 (t) = θ 2 (t) = θ(t), and τ 1 (t) = τ 2 (t) = τ (t).Moreover, x G (t) = x O2 (t) + a, where x O2 (t) = rθ(t), and hence, ẍG (t) = ẍO2 (t) and ẍG (t) = r θ(t).Thus, the two equations of the system (28) are identical and give the dynamics of G which is independent of I z as the motion is purely translational.Moreover, (30) gives θ3 (t) = θ(t) and ψ(t) = 0.
3) c 2 2,2 (Chemical plant): As showed in Fig. 7, the anaerobic digestion plant is divided into P r and P d .The mass balance of P r in discrete-time yields (33) where n d is the time in which the reaction in P d begins and m r is the mass inside P r .
For n ≥ n d , the sub-compartment P r becomes a continuous-time system which supplies a continuous flow to the digester P d .The digester P d is modeled as a continuous stirred tank reactor.Specifically, the anaerobic digestion occurring inside P d is a four-state dynamical system, which results from the mass balance of the species involved and it assumes a two-stage reaction.Namely, first, the organic substrate S 1 (t) is degraded into volatile fatty acids S 2 (t) by acidogenic bacteria X 1 (t), and then the methanogenic bacteria X 2 (t) consume the volatile fatty acids to produce methane CH 4 and carbon dioxide CO 2 [34].The set of four ordinary differential equations of the anaerobic digestion resulting from the mass balance are given by [34] Ẋ1 where D j (t)| j=1,2,3,4 is the dilution rate for the j-th state, and α, S 1in , S 2in , k 1 , k 2 , k 3 , µ 1max , µ 2max , K S1 , K S2 , and K I2 are specific constants detailed in [34], [35].The input material flow to P d supplied by P r is where ρ b and V d are the biomass density and the digester working volume, respectively, while the biomethane flow q M produced by the anaerobic digestion is q M (S 2 (t), X 2 (t)) = k 6 µ 2 (S 2 (t))X 2 (t) [35].Hence, for n ≥ n d , the mass balance of P r yields 4) Compartmental control: Here we introduce two compartmental controllers following the black arrow in Step 3 of our methodology (Fig. 4): one for c 3  1,2 and one for c 2 2,2 .Note that this can be done at this stage because the dynamics of the compartments are known.
The motion of the truck must satisfy x G (t u ) = H + a and ẋG (t u ) = 0, where t u is the unloading time.In this way, the truck reaches the plant at the unloading time t u with null speed.A simple open-loop control law that satisfies the requirements is where To regulate the digester to the desired working point, that is, the equilibrium "SS6" in [34], a nonlinear control system was designed and implemented.Specifically, we implemented the control law of Haddad et al. [36] because, for dynamical systems of the form where f : R n k → R n k and G : R n k → R n k ×z k are continuous functions with f (0) = 0, it guarantees the stabilization of the zero solution x(t) ≡ 0 of (44) in a finite-time, where 0 ∈ R n k is a vector of zeros.In our case, . We designed the controller by choosing the functions L 2 ( x), R −1 2 ( x), and V ( x) (in the notation of [36]) as and where x(t) is the state vector translated to have the desired working point corresponding to the zero solution.This choice requires that G( x(t)) −1 exists, but it has the following benefits: first, the condition (33) of [36] simplifies significantly so that it can be easily checked for systems with complex expressions of f ( x(t)) and G( x(t)) since now it only depends on V ( x); second, since f (0) = 0, the condition (34) of [36] is satisfied regardless of the expressions of f ( x(t)) and G( x(t)); and third, the dynamics of the closed-loop reduce to the gradient system which is easy to interpret: the state dynamics ẋ(t) is proportional to the negative gradient V ( x) = ∂V ( x) ∂ x .Several simulation results are shown in Fig. 9, which were achieved considering the values in Table I.

VI. CONCLUSION
The long-term unsustainability of the current take-makedispose economy requires to redesign the material flows across natural and built environments.This paper establishes the foundations of a thermodynamics-based material flow design towards circularity by integrating the well-established design approach of thermodynamic cycles (e.g., Rankine cycle) with graph theory, dynamical systems, and control.Our examples demonstrate both the theoretical efficacy and the applicability of the proposed methodology.
Future work will consist of designing TMNs by starting with few compartments and targeting materials whose circularity has high priority, e.g., atmospheric carbon dioxide and critical raw materials [7], [8].With respect to the existing literature, the goal is to complement MFA with dynamical power balances, graph theory, and control systems.See [34], [35] See [34], [35] x 0 3.43 g L , −0.13 g L , 7.

Figure 2
Figure 2 is an example for n β = 2.A well-established formalism to represent a network of systems is graph theory[13]: examples of network design theories based on graphs are electrical networks[16], hydraulic networks[10], and multiagent systems[31].Since the system in Fig.2can be seen as a network of thermodynamic compartments connected through the material flow, we will use graph theory to formulate the system in Fig.2as a network and then state the definition of a thermodynamical material network.
, and n a = 1 depicted using a compartmental diagram (top), a compartmental digraph (middle), and a mass-flow digraph M (N ) (bottom).

Fig. 4 :
Fig. 4: Proposed methodology to design a circular flow of B.

Fig. 7 :
Fig. 7: The subsystem under study in Example V-B: (a) physical representation and (b) compartmental diagram.

1 dayFig. 9 :
Fig. 9: Simulation results of the plant digester P d and the truck.