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

ABSTRACT

Material flow analysis (MFA) is the main methodology to assess material flow circularity. It is essentially a data-analysis-based approach whose physical foundations consist of conservation of mass. To improve both the accuracy and the repeatability of MFA models, in this paper we leverage compartmental dynamical thermodynamics merged with graph theory and control 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. Thus, our methodology consists of designing thermodynamical material networks (TMNs). We also provide a physics-based definition of circularity and implement a nonlinear compartmental control, which has been possible since TMNs are highly dynamic models based on differential calculus (i.e. ordinary differential equations) rather than on arithmetic as is typical for MFA models. As we envision scalable and repeatable designs of TMNs, we made publicly available the paper source code.¹

KEYWORDS:

• Compartmental dynamical thermodynamics
• graph theory
• control systems
• material flow design
• industrial ecology
• circular economy

I. Introduction

As the human population is predicted to reach 8 billion by 2024 and 10 billion by 2056 (Max Roser and Ortiz-Ospina 2019), modern society strives to provide the needed services and products on a large-scale. 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 (Ellen MacArthur Foundation 2022).

A natural resource particularly at risk is climate stability, which is mostly being altered by the atmospheric carbon dioxide concentration (National Aeronautics and Space Administration 2022). The world emitted 6 billion tonnes of $\mathrm{C}{\mathrm{O}}_2$ in 1950, 22 billion tonnes in 1990, and 36 billion tonnes in 2019 (i.e. 6 times the emissions of 1950), and the annual emissions have yet to reach their peak (Ritchie and Roser 2020). Therefore, an effective control of this material is necessary to respect the global warming goal limiting the temperature increase to 1.5 degrees Celsius above the pre-industrial levels as established with the Paris Agreement in 2015 (United Nations 2015).

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 (The Ocean Cleanup 2022). 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 (European Commission 2020) and the United States (U.S. Department of Interior 2018) whose supply is particularly at risk. Those materials are currently required for clean technologies such as solar panels, wind turbines, and 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 (Kirchherr, Reike, and Hekkert 2017). In this paper, we focus on the flows and accumulations of materials (e.g. carbon dioxide, gold, plastic and 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 minimising the leakages in the water network.

Non-circular material flows cause harmful accumulations and depletions of resources. For example, one of the most challenging accumulations to be solved is the carbon dioxide in the atmosphere, which accelerates climate instability (National Aeronautics and Space Administration 2022). Plastic and paper carrier bags were addressed by Thakker and Bakshi (2021) as they form a relevant fraction of waste streams. The authors concluded by suggesting several recovery routes with different trade-off optimality performances. In Altamira-Algarra et al. (2022), the runoff streams of phosphorus were recirculated to reduce its contamination of water bodies as well as the demand of mining of phosphorus rocks. Municipal solid waste was targeted by Gutiérrez, Mendoza Fandiño, and Cabello Eras (2021), where 10 alternative routes were considered for the recovery chain such as incineration for electricity production and anaerobic digestion of the organic fraction. Sivanandhini, Unnikrishnan, and Gedam (2021) examined the lubricant supply chain considering a network made of one source of crude oil, one collection point, one refinery, one blending plant, and one distribution centre, while Mahdjoub et al. (2021) investigated the post-consumer use of waste glass as a coarse aggregate within burnt clay bricks and its recovery from local dumped stocks.

To enhance the modelling, 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 effectively 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 minimise leakages of material and a mathematical framework that effectively depicts the network architecture is graph theory (Deuerlein 2008); and

2. their modelling is based on the first principle of thermodynamics and the mass conservation equation (Kaminski and Jensen 2017).

Given the generality of both thermodynamics (Haddad 2017) and graph theory (Bondy and Murty 1976), in this paper we propose to extend the modelling approach of water supply networks to model the flow of any material leveraging compartmental dynamical thermodynamics (Haddad 2019) and graph theory (Bondy and Murty 1976).

The main contributions of the paper are the following.

1. We provide a physics-based foundation to material circularity to add mathematical rigour to the topic (Kirchherr, Reike, and Hekkert 2017) (see Sections III).

2. We extend the methodology of material flow analysis (MFA) by adding dynamical power balances and control systems to improve its accuracy and repeatability (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 (Bakshi and Paulson 2022).

The paper is organised 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

II.A. 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) behaviour of systems from a thermodynamic perspective without a focus on the multi-machine nature typical of compartmental thermodynamics. Examples of works are Freitas, Delvenne, and Esposito (2020) for electrical networks, Freitas, Delvenne, and Esposito (2021) for electronic circuits, Gay-Balmaz and Yoshimura (2018) for mechanical systems, and Avanzini et al. (2021) and Penocchio, Rao, and Esposito (2021) for chemical reaction networks. The combination of these two branches of thermodynamics is compartmental dynamical thermodynamics (Haddad 2019).

II.B. Graph theory and thermodynamics for sustainability

In the context of circular economy, some graph-based approaches have been proposed recently. Moktadir et al. (2018) used a graph architecture to examine and prioritise the driving factors of sustainable manufacturing practices; in Singh et al. (2020) and How et al. (2018) a graph architecture is used to analyse 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. (2020) proposed the use of graphs to model the production of chitin by bio-conversion of municipal waste.

The idea of using thermodynamics for ecological modelling dates back almost 30 years (Schneider and Kay 1994). In 2011, Bakshi, Gutowski, and Sekulić (2011) further extended this vision by proposing thermodynamics as the science of sustainability. In 2014, Capilla and Delgado (2014) presented a cradle-to-cradle view of Earth abiotic resources based on the second law of thermodynamics. More recently, Tan et al. (2019) used thermodynamics to assess the ability of a city to use the available resources and reduce their demand, while Sugar and Kennedy (2020) considered the thermodynamics of urban growth. Varbanov et al. (2020) proposed an exergy footprint evaluation of process systems and tested it on a waste-to-energy process fed by municipal solid waste, while Artuzo et al. (2021) developed an unsustainability index for agricultural systems based on the laws of thermodynamics. All these works contribute to link thermodynamics to sustainability, but they do not provide a practical and clear procedure to systematically design material flows. In this paper, we develop a systematic and general material flow design methodology based on thermodynamics; more specifically, our methodology is based on a dynamical systems theory of thermodynamics (Haddad 2019) merged with graph theory (Bondy and Murty 1976) and control theory (Haddad and Chellaboina 2011). The first theory captures the dynamics of the system, the second theory provides an effective mathematical formalism to study the properties and topology of network systems (in our case, the network results from the sequence of material transformation and transportation stages), and the third theory regulates the behaviour of the network; the goal is the design of networks with closed-loop material flows.

II.C. 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 (Graedel 2019). 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 (Brunner and Rechberger 2016), our methodology extends it by adding dynamical power balances (Haddad 2019) and control systems (Haddad and Chellaboina 2011).

II.D. Control systems in life-cycles

Control theory is a well-established discipline sitting between applied mathematics and engineering. Two selected works from the vast literature in the field are Doyle, Francis, and Tannenbaum (2013) for linear feedback control and Haddad and Chellaboina (2011) 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 using a biomethane production subsystem as an example.

III. Circularity and thermodynamical material networks

Consider a cube of material $\beta$, infinitesimal mass d$m$, density $\rho$, and volume $\mathrm{d}V=\frac{\mathrm{d}m}{\rho }$ as in Figure 1. Let $\mathit{G}(t)=[x_G(t),y_G(t),z_G(t){]}^{\mathrm{T}}$ be the centre of mass whose coordinates are written with respect to a fixed reference frame with origin $\mathit{O}=[0,0,0{]}^{\mathrm{T}}$. In general, $x_G(t)$, $y_G(t)$, and $z_G(t)$ can vary with time $t$. Let $\mathit{p}(t)=\mathit{G}(t)-O$ be the position vector of the cube centre of mass $\mathit{G}$.

Figure 1. Infinitesimal cube of material $\beta$ and centre of mass $\mathit{G}$.

Infinitesimal cube of material β and centre of mass G.

Definition 1

(Mechanics-based circularity). The flow of $\beta$ is mechanically circular if there exist $t_0\ge 0$ and $t^{\ast }\in (t_0,\mathrm{\infty })$ such that

(1) $p\left(t_0\right)=p\left(t^{\ast }\right),\stackrel{\cdot }{p}\left(t_0\right)\ne 0,t^{\ast }>t_0.$(1)

Remark 1. As the material $\beta$ 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 (Kaminski and Jensen 2017). Each control volume identifies a thermodynamic compartment. For example, a simple ideal Rankine cycle is made of eight thermodynamic compartments, namely, 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 (Bakshi, Gutowski, and Sekulić 2011; Haddad 2017). Hence, we can generalise the definition of circularity based on mechanics (Definition 1) with the following thermodynamics-based definition.

Let $c^k$ be the $k$-th thermodynamic compartment identified by the $k$-th control volume and let $\beta$ be the material of an infinitesimal cube as in Figure 2.

Definition 2

(Thermodynamics-based circularity). The flow of $\beta$ is thermodynamically circular if there exists an ordered sequence of compartments $\varphi =(c^1,\dots ,c^k,\dots ,c^1)$ processing $\beta$, which begins and ends in $c^1$. Moreover, if some $c^k\in \varphi$ chemically transforms a material ${\beta }_1$ into a material ${\beta }_2$ and there exists an ordered sequence $\varphi$ processing $\mathcal{B}=\{{\beta }_1,{\beta }_2\}$, then the flow of the material set $\mathcal{B}=\{{\beta }_1,{\beta }_2\}$ is thermodynamically circular. More generally, the flow of $\mathcal{B}=\{{\beta }_1,\dots ,{\beta }_q,\dots ,{\beta }_{n_{\beta }}\}$ is thermodynamically circular if there exists an ordered sequence $\varphi$ processing $\mathcal{B}$.

Figure 2 is an example for $n_{\beta }=2$.

Figure 2. Example of a circular flow of $\mathcal{B}=\{{\beta }_1,{\beta }_2\}$. The chemical transformation from ${\beta }_1$ to ${\beta }_2$ is in $c^2$.

An infinitesimal cube of material that crosses three compartments. A chemical reaction occurs in the second compartment c 2. The position of the cube is with respect to a Cartesian coordinate system.

A well-established formalism to represent a network of systems is graph theory (Bondy and Murty 1976); examples of network design theories based on graphs are electrical networks (Freitas, Delvenne, and Esposito 2020), hydraulic networks (Deuerlein 2008), and multiagent systems (Mesbahi and Egerstedt 2010). Since the system in Figure 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 Figure 2 as a network and then state the definition of a thermodynamical material network.

Definition 3

(Bondy and Murty 1976). A directed graph $D$ or digraph is a graph identified by a set of $n_v$ vertices $\{v_1,v_2,\dots ,v_{n_v}\}$ and a set of $n_a$ arcs $\{a_1,a_2,\dots ,a_{n_a}\}$ that connect the vertices. A digraph $D$ in which each vertex or arc is associated with a weight is a weighted digraph.

Let $c_{i,j}^k$ be the $k$-th thermodynamic compartment through which the material flows from compartment $i$ to compartment $j$.

Definition 4

(Thermodynamical material network). A thermodynamical material network (TMN) is a set $\mathcal{N}$ of connected thermodynamic compartments, that is,

(2) $\begin{array}{c}\mathcal{N}=\left\{c_{1,1}^1,\dots ,c_{k_v,k_v}^{k_v},\dots ,c_{n_v,n_v}^{n_v},\right.\\ \left.c_{i_{n_v+1},j_{n_v+1}}^{n_v+1},\dots ,c_{i_{n_v+k_a},j_{n_v+k_a}}^{n_v+k_a},\dots ,c_{i_{n_c},j_{n_c}}^{n_c}\right\},\end{array}$(2)

which transport, store, and transform a set of $n_{\beta }$ materials $\mathcal{B}=\{{\beta }_1,\dots ,{\beta }_q,\dots ,{\beta }_{n_{\beta }}\}$ and whose modelling is based on compartmental dynamical thermodynamics (Haddad 2019).

Specifically, $\mathcal{N}=\mathcal{R}\cup \mathcal{T}$, where $\mathcal{R}\subseteq \mathcal{N}$ is the subset of compartments that store, transform, or use the target material, and $\mathcal{T}\subseteq \mathcal{N}$ is the subset of compartments that move the target material between the compartments belonging to $\mathcal{R}\subseteq \mathcal{N}$. A net $\mathcal{N}$ is associated with its weighted mass-flow digraph $M(\mathcal{N})$, which is a weighted digraph whose vertices are the compartments $c_{i,j}^k\in \mathcal{R}$ and whose arcs are the compartments $c_{i,j}^k\in \mathcal{T}$. A vertex also results from the intersection of 3 or more arcs. For vertex-compartments $c_{i,j}^k\in \mathcal{R}$ it holds that $i=j$, whereas for arc-compartments $c_{i,j}^k\in \mathcal{T}$ it holds that $i\ne j$. The weight assigned to a vertex-compartment $c_{i,j}^k\in \mathcal{R}$ is identified by the mass stock $m_k$ within the corresponding compartment, whereas the weight assigned to an arc-compartment $c_{i,j}^k\in \mathcal{T}$ is the mass flow rate ${\dot{m}}_{i,j}$ from the vertex-compartment $c_{i,i}^i\in \mathcal{R}$ to the vertex-compartment $c_{j,j}^j\in \mathcal{R}$. Hence, the orientation of an arc-compartment $c_{i,j}^k\in \mathcal{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 $k$-th 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 $n_c=n_v+n_a$.

Definition 5

(Compartmental diagram). The compartmental diagram of the network (2) depicts the thermodynamic compartments $c_{i,j}^k$ and the arrows of material flow along with the material class ${\beta }_q\in \mathcal{B}$.

Definition 6

(Compartmental digraph). The compartmental digraph of the network (2) is a weighted digraph with arcs and vertices labelled with the corresponding compartmental nomenclature $c_{i,j}^k$.

Figure 3 shows an example of $\mathcal{N}=\left\{c_{1,1}^1,c_{2,2}^2,c_{1,2}^3\right\}$ with $\mathcal{B}=\{{\beta }_1,{\beta }_2\}$, $n_c=3$, $n_v=2$, and $n_a=1$ depicted using a compartmental diagram (top), a compartmental digraph (middle), and a mass-flow digraph $M(\mathcal{N})$ (bottom).

Figure 3. Graphical representations of $\mathcal{N}=\left\{c_{1,1}^1,c_{2,2}^2,c_{1,2}^3\right\}$: (a) compartmental diagram, (b) compartmental digraph, and (c) mass-flow digraph.

The three graphical representations of a TMN: compartmental diagram at the top, compartmental digraph in the middle, and mass-flow digraph at the bottom.

Next, we introduce a few more definitions from graph theory. The reason will be clarified afterwards.

Definition 7

(Bondy and Murty 1976). A directed walk in $D$ is a finite non-null sequence $W=(v_0,a_1,v_1,a_2,\dots ,a_l,v_l)$ whose terms are alternatively vertices and arcs such that, for $i=1,2,\dots ,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

(Bondy and Murty 1976). If the sequence of arcs $a_1,a_2,\dots ,a_l$ of a directed walk $W$ are distinct, then $W$ is a directed trial.

Definition 9

(Bondy and Murty 1976). 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

(Bondy and Murty 1976). A closed directed trial whose origin and internal vertices are distinct is a directed cycle $\varphi$.

Summarising, a directed cycle $\varphi$ 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(\mathcal{N})$ must be a directed cycle $\varphi$.

Remark 2. Consider the network (2) with a mass-flow digraph $M(\mathcal{N})$. Then, the flow of the set of materials $\mathcal{B}=\{{\beta }_1,\dots ,{\beta }_q,\dots ,{\beta }_{n_{\beta }}\}$ is thermodynamically circular if $M(\mathcal{N})$ is a directed cycle $\varphi$.

The mass-flow matrix $\mathbf{\Gamma }(\mathcal{N})$ associated with the network (2) is given by

(3) $\begin{array}{cc}\mathbf{\Gamma }(\mathcal{N})& =\left[\begin{array}{ccc}{\gamma }_{1,1}& \dots & {\gamma }_{1,n_v}\\ ⋮& \ddots & ⋮\\ {\gamma }_{n_v,1}& \dots & {\gamma }_{n_v,n_v}\end{array}\right]\\ \hspace{0.17em}& =\left[\begin{array}{cccc}{m}_1& {\dot{m}}_{1,2}& \dots & {\dot{m}}_{1,n_v}\\ {\dot{m}}_{2,1}& m_2& \dots & {\dot{m}}_{2,n_v}\\ ⋮& \ddots & \ddots & ⋮\\ {\dot{m}}_{n_v,1}& {\dot{m}}_{n_v,2}& \dots & m_{n_v}\end{array}\right],\end{array}$(3)

whose entries along the diagonal are the weights of the vertex-compartments $c_{i,j}^k\in \mathcal{R}$ (i.e. mass stocks) and whose off diagonal entries are the weights of the arc-compartments $c_{i,j}^k\in \mathcal{T}$ (i.e. mass flow rates). Hence, $\mathbf{\Gamma }(\mathcal{N})$ is a square matrix of size $n_v\times n_v$ with nonnegative real entries, i.e. $\mathbf{\Gamma }\in {\bar{\mathbb{R}}}_{+}^{n_v\times n_v}$.

The mass conservation principle (Kaminski and Jensen 2017) establishes the relationship between the entries of $\mathbf{\Gamma }(\mathcal{N})$, namely,

(4) $\frac{\mathrm{d}}{\mathrm{d}t}{m}_k=\sum _{i=1}^{n_v}{\dot{m}}_{i,k}-\sum _{j=1}^{n_v}{\dot{m}}_{k,j},$(4)

which can be further written in terms of the entries of $\mathbf{\Gamma }(\mathcal{N})$ as

(5) $\frac{\mathrm{d}}{\mathrm{d}t}{\gamma }_{k,k}=\sum _{i=1,i\ne k}^{n_v}{\gamma }_{i,k}-\sum _{j=1,j\ne k}^{n_v}{\gamma }_{k,j},$(5)

or, equivalently, in vector form as

(6) $\begin{array}{cc}\frac{\mathrm{d}}{\mathrm{d}t}m& =\frac{\mathrm{d}}{\mathrm{d}t}\left[\begin{array}{c}{m}_1\\ m_2\\ ⋮\\ m_{n_v}\end{array}\right]=\frac{\mathrm{d}}{\mathrm{d}t}\left[\begin{array}{c}{\gamma }_{1,1}\\ {\gamma }_{2,2}\\ ⋮\\ {\gamma }_{n_v,n_v}\end{array}\right]\\ \hspace{0.17em}& =\left[\begin{array}{c}\sum _{\begin{array}{c}\underset{i\ne 1}{i=1}\end{array}}^{n_v}{\gamma }_{i,1}-\sum _{\begin{array}{c}\underset{j\ne 1}{j=1}\end{array}}^{n_v}{\gamma }_{1,j}\\ \sum _{\begin{array}{c}\underset{i\ne 2}{i=1}\end{array}}^{n_v}{\gamma }_{i,2}-\sum _{\begin{array}{c}\underset{j\ne 2}{j=1}\end{array}}^{n_v}{\gamma }_{2,j}\\ ⋮\\ \sum _{\underset{i\ne n_v}{i=1}}^{n_v}{\gamma }_{i,n_v}-\sum _{\begin{array}{c}\underset{j\ne n_v}{j=1}\end{array}}^{n_v}{\gamma }_{n_v,j}\end{array}\right].\end{array}$(6)

Therefore,

(7) $\begin{array}{cc}m& =\left[\begin{array}{c}{\gamma }_{1,1}\\ {\gamma }_{2,2}\\ ⋮\\ {\gamma }_{n_v,n_v}\end{array}\right]\\ \hspace{0.17em}& =\left[\begin{array}{c}\underset{t_0}{\overset{t}{\int }}\left(\sum _{\underset{i\ne 1}{i=1}}^{n_v}{\gamma }_{i,1}-\sum _{\begin{array}{c}\underset{j\ne 1}{j=1}\end{array}}^{n_v}{\gamma }_{1,j}\right)\mathrm{d}\tau \\ \underset{t_0}{\overset{t}{\int }}\left(\sum _{\begin{array}{c}\underset{i\ne 2}{i=1}\end{array}}^{n_v}{\gamma }_{i,2}-\sum _{\begin{array}{c}\underset{j\ne 2}{j=1}\end{array}}^{n_v}{\gamma }_{2,j}\right)\mathrm{d}\tau \\ ⋮\\ \underset{t_0}{\overset{t}{\int }}\left(\sum _{\begin{array}{c}\underset{i\ne n_v}{i=1}\end{array}}^{n_v}{\gamma }_{i,n_v}-\sum _{\begin{array}{c}\underset{j\ne n_v}{j=1}\end{array}}^{n_v}{\gamma }_{n_v,j}\right)\mathrm{d}\tau \end{array}\right],t\ge t_0.\end{array}$(7)

Remark 3. In line with the standard nomenclature adopted in thermal engineering, the mass accumulation or depletion in vertex-compartments is denoted as $\frac{\mathrm{d}}{\mathrm{d}t}m$ and not as $\dot{m}$ (Kaminski and Jensen 2017). 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 $\lambda$.

Definition 11

(Graph-based circularity indicator). The graph-based circularity indicator of the network (2) associated with the mass-flow matrix (3) is

(8) $\lambda \left(\mathbf{\Gamma }\right)=\frac{{\sum }_{k=1}^{n_{\varphi }}\mathrm{C}\mathrm{M}\left({\varphi }_k\right)}{{\sum }_{k=1}^{n_{\varphi }}\mathrm{C}\mathrm{M}\left({\varphi }_k\right)+{\sum }_{{\gamma }_{i,j}\in \mathcal{Q}}{\gamma }_{i,j}},$(8)

where $n_{\varphi }$ is the number of directed cycles in $M$,

(9) $\mathcal{Q}=\{{\gamma }_{i,j}|{\dot{m}}_{i,j}\mathrm{d}\mathrm{o}\mathrm{e}\mathrm{s}\mathrm{n}\mathrm{o}\mathrm{t}\mathrm{b}\mathrm{e}\mathrm{l}\mathrm{o}\mathrm{n}\mathrm{g}\mathrm{t}\mathrm{o}\mathrm{a}\mathrm{n}\mathrm{y}\mathrm{d}\mathrm{i}\mathrm{r}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{e}\mathrm{d}\mathrm{c}\mathrm{y}\mathrm{c}\mathrm{l}\mathrm{e}\},$(9)

and

(10) $\mathrm{C}\mathrm{M}\left(\varphi \right)=\frac{1}{l}\sum _{{\gamma }_{i,j}\in \mathcal{Y}}{\gamma }_{i,j}$(10)

is the cycle mean of $\varphi$, with

(11) $\mathcal{Y}=\{{\gamma }_{i,j}|{\dot{m}}_{i,j}\in \varphi \}.$(11)

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 Figure 4 and involves three main steps; its output is a TMN, its goal is designing the flow of the material set $\mathcal{B}=\{{\beta }_1,\dots ,{\beta }_q,\dots ,{\beta }_{n_{\beta }}\}$, and it has two optional extensions at Steps 2 and 3 indicated by the black arrows.

Figure 4. Proposed methodology to design a circular flow of $\mathcal{B}$.

Summary of the proposed methodology, which consists of three steps and whose output is a TMN for modelling, planning, and control of circular material flows.

Specifically, the first step is the choice of the material set $\mathcal{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 Figure 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_{i,j}^k\in \mathcal{N}$ the dynamical form of the compartmental energy balance, that is,

(12) $\frac{\mathrm{d}{E}_k}{\mathrm{d}t}=\frac{\mathrm{d}}{\mathrm{d}t}(K_k+U_k+P_k)={\dot{Q}}_k-{\dot{W}}_k,$(12)

and/or the compartmental mass balance (4). Equation (12)(12) $\frac{\mathrm{d}{E}_k}{\mathrm{d}t}=\frac{\mathrm{d}}{\mathrm{d}t}(K_k+U_k+P_k)={\dot{Q}}_k-{\dot{W}}_k,$(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 ${\dot{Q}}_k$ and ${\dot{W}}_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 (Figure 4) if the power balance (12) and the mass balance (4) are written in state-space form

(13) $x_k^{\cdot }=F\left(x_k,u_k\right),$(13)

where $x_k\in {\mathbb{R}}^{n_k}$ is the state vector of the $k$-th compartment, $u_k\in {\mathbb{R}}^{z_k}$ is the control input of the $k$-th compartment, and $F:{\mathbb{R}}^{n_k}\times {\mathbb{R}}^{z_k}\to {\mathbb{R}}^{n_k}$ is a continuous nonlinear function.

V. Illustrative examples

A. Circularity Calculation

This example demonstrates the calculation of the circularity indicator (8) considering the net $\mathcal{N}$ depicted in Figure 5, where $\alpha \in [0,1]$.

Figure 5. Compartmental digraph in (a) and corresponding mass-flow digraph in (b) considered in Example 5.1 with $\alpha \in [0,1]$.

Compartmental digraph on the left side and mass-flow digraph on the right side corresponding to the first example.

Assume $m_i=0{|}_{i=2,3,4}$. Hence, the mass balance (7) is respected in $c_{2,2}^2$ and $c_{4,4}^4$, while in $c_{3,3}^3$ it yields

(14) $(1-\alpha ){\dot{m}}_{\alpha }={\dot{m}}_{\beta }+{\dot{m}}_{\alpha }\Rightarrow {\dot{m}}_{\beta }=-\alpha {\dot{m}}_{\alpha }<0,$(14)

which is non-physical. Therefore, the direction of $c_{3,5}^{10}$ must be inverted as $c_{5,3}^{10}$ to get ${\dot{m}}_{\beta }=\alpha {\dot{m}}_{\alpha }$. The remaining stocks $m_1$ and $m_5$ follow from the mass balance as $m_1(t)={\int }_0^t\alpha {\dot{m}}_{\alpha }\mathrm{d}\tau$ and $m_5(t)=-{\int }_0^t\alpha {\dot{m}}_{\alpha }\mathrm{d}\tau$, with the latter requiring that $\alpha {\dot{m}}_{\alpha }t\le m_{5,0}$, where $m_{5,0}$ is the initial stock in $c_{5,5}^5$. With this, all the stocks and flows different from zero are functions of ${\dot{m}}_{\alpha }$.

To calculate the circularity $\lambda$ given by (8), we observe that for $\alpha =1$ the network has no cycles, and hence $\lambda =0$. In contrast, if $\alpha =0$, then $\mathcal{Q}=\text{O}$, and hence $\lambda =1$. Finally, for $0<\alpha <1$, the net has a cycle with

(15) $\mathrm{C}\mathrm{M}\left(\varphi \right)=\frac{{\dot{m}}_{\alpha }\left(3-2\alpha \right)}{3},$(15)

and hence,

(16) $\lambda (\alpha )=\frac{3-2\alpha }{3+4\alpha }.$(16)

In summary,

(17) $\lambda (\alpha )=\left\{\begin{array}{cc}0,& \alpha =1\\ \frac{3-2\alpha }{3+4\alpha },& 0\le \alpha <1.\end{array}\right.$(17)

Figure 6 shows $\lambda (\alpha )$ as a function of $\alpha$.

Figure 6. Circularity indicator (8) as a function of $\alpha$ for the network in Figure 5.

Circularity indicator given in equation (8) as a function of the parameter α for the first example.

Observations:

1. The computation of $\lambda$ requires the preliminary steps to be completed as shown in Fig. 4: Step 1 first, then Step 2, then define $M(\mathcal{N})$ as indicated by the black arrow on the left, and finally compute $\lambda$. In this introductory example, the network $\mathcal{N}$ is given, hence Steps 1 and 2 are skipped. In general, for a particular case study, Steps 1 and 2 are required to compute $\lambda$.

2. A higher value of $\lambda$ indicates a higher circularity. An intuitive example with $\lambda =1$ (maximum value) is a network with a digraph configured as a closed planar figure, e.g. a triangle or a rectangle, because in this case all the flows belong to a directed cycle, and hence, $\mathcal{Q}=\text{O}$. If any of the edges of the planar figure (i.e. any of the arcs of the mass-flow digraph) are removed, then $\lambda =0$ (minimum value).

3. The mass balance underpinning MFA is used to calculate the weights of the digraph in Fig. 5, while the power balance has not been used. Hence, this example merges the mass balance of MFA with graph theory to calculate the circularity $\lambda$ for different material flows.

B. Subsystem of Biomethane 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 Figure 7): the biomass hub, the truck to transport the biomass, and the anaerobic digestion plant to convert 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{]}^{\mathrm{T}}$ (i.e. the material motion is along $x$ only), the position of the hub is $x_h$, $O\equiv 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$.

Figure 7. The subsystem under study in Example 5.2: (a) physical representation and (b) compartmental diagram.

Physical representation at the top and compartmental diagram at the bottom of the biomethane subsystem in the second example.

The first step of the methodology requires to choose the set of materials of interest, i.e. $\mathcal{B}$; in this case, $\mathcal{B}=\{{\beta }_1\}$, where ${\beta }_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, $\mathcal{N}=\{c_{1,1}^1,c_{2,2}^2,c_{1,2}^3\}$, $n_c=3$, $n_v=2$, and $n_a=1$. The thermodynamic modelling 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 modelled 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

(18) $m_1(n+1)-m_1(n)=-{\delta }_{n_l}(n)m_l,n\in {\bar{\mathbb{Z}}}_{+},$(18)

where ${\bar{\mathbb{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 ${\delta }_{n_l}(n)$ is the Kronecker delta, that is,

(19) ${\delta }_{n_l}(n)=\left\{\begin{array}{cc}0,& n\ne n_l\\ 1,& n=n_l.\end{array}\right.$(19)

2) $c_{1,2}^3$ (Truck): For $n_l\le n\le n_u$, the biomass with mass $m_l$ is in 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-wheeled vehicle shown in Figure 8 (Gabiccini 2011).

Figure 8. A three-wheeled vehicle model of the truck, that is, compartment $c_{1,2}^3$ (modified from (Gabiccini 2011)).

Three-wheeled vehicle model of the truck of the second example with geometrical dimensions and Lagrangian coordinates.

The equations of motion of the vehicle can be derived using Lagrange’s equations of motion (Siciliano et al. 2010)

(20) $\frac{\mathrm{d}}{\mathrm{d}t}{\left(\frac{\mathrm{\partial }L}{\mathrm{\partial }\stackrel{\cdot }{q}}\right)}^{\mathrm{T}}=\xi +{\left(\frac{\mathrm{\partial }L}{\mathrm{\partial }q}\right)}^{\mathrm{T}},$(20)

where

(21) $L=K-P$(21)

is the Lagrangian function of the mechanical system, $q\in {\mathbb{R}}^d$ is the vector of generalised coordinates, $\stackrel{\cdot }{q}\in {\mathbb{R}}^d$ is the vector of generalised velocities, and $\xi \in {\mathbb{R}}^d$ is the vector of generalised forces associated with the generalised coordinates $q$. Before we use the Lagrangian formulation (20) to define the dynamical equations of the three-wheeled 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,

(22) $\frac{\mathrm{d}}{\mathrm{d}t}(K+U+P)=\dot{Q}-\dot{W}.$(22)

Since the heat flow $\dot{Q}$ and the internal energy $U$ are neglected in solid mechanics, (22) reduces to

(23) $\frac{\mathrm{d}}{\mathrm{d}t}(K+P)=-\dot{W}.$(23)

Note that in solid mechanics the potential energy $P$ corresponds to the conservative work done by the gravitational force, and hence (Siciliano et al. 2010),

(24) $\frac{\mathrm{\partial }P}{\mathrm{\partial }\stackrel{\cdot }{q}}=0.$(24)

Taking the partial derivative on both sides of (23) with respect to $\stackrel{\cdot }{q}$ and transposing the resulting equation yields

(25) $\frac{\mathrm{d}}{\mathrm{d}t}{\left[\frac{\mathrm{\partial }\left(K+P\right)}{\mathrm{\partial }\stackrel{\cdot }{q}}\right]}^{\mathrm{T}}=-{\left(\frac{\mathrm{\partial }\dot{W}}{\mathrm{\partial }\stackrel{\cdot }{q}}\right)}^{\mathrm{T}},$(25)

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

(26) $\frac{\mathrm{d}}{\mathrm{d}t}{\left(\frac{\mathrm{\partial }L}{\mathrm{\partial }\stackrel{\cdot }{q}}\right)}^{\mathrm{T}}=-{\left(\frac{\mathrm{\partial }\dot{W}}{\mathrm{\partial }\stackrel{\cdot }{q}}\right)}^{\mathrm{T}},$(26)

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,

(27) $-{\left(\frac{\mathrm{\partial }\dot{W}}{\mathrm{\partial }\stackrel{\cdot }{q}}\right)}^{\mathrm{T}}=\xi +{\left(\frac{\mathrm{\partial }L}{\mathrm{\partial }q}\right)}^{\mathrm{T}},$(27)

Therefore, (26) can be written as (20).□

For the three-wheeled 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 (Gabiccini 2011)

(28) $\mathit{B}\left[\begin{array}{c}{\ddot{\theta }}_1\\ {\ddot{\theta }}_2\end{array}\right]=\left[\begin{array}{c}{\tau }_1\\ {\tau }_2\end{array}\right],$(28)

where ${\tau }_1$ and ${\tau }_2$ are the control torques of the engines $M_1$ and $M_2$, respectively, and

(29) $\mathit{B}=m_{vl}\left[\begin{array}{cc}\left(\frac{a^2{r}^2}{l^2}+\frac{r^2}{4}\right)+\epsilon & \left(-\frac{a^2{r}^2}{l^2}+\frac{r^2}{4}\right)-\epsilon \\ \left(-\frac{a^2{r}^2}{l^2}+\frac{r^2}{4}\right)-\epsilon & \left(\frac{a^2{r}^2}{l^2}+\frac{r^2}{4}\right)+\epsilon \end{array}\right],$(29)

where $m_{vl}=m_v+m_l$, $m_v$ is the vehicle mass, $\epsilon =I_z\frac{r^2}{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 ${\dot{\theta }}_1$ and ${\dot{\theta }}_2$ are determined by integration of (28), $\dot{\psi }$ and ${\dot{\theta }}_3$ are given by

(30) $\left[\begin{array}{c}{\dot{\theta }}_3\\ \dot{\psi }\end{array}\right]=F\left[\begin{array}{c}{\dot{\theta }}_1\\ {\dot{\theta }}_2\end{array}\right],$(30)

where

(31) $\mathit{F}=\left[\begin{array}{cc}\frac{\mathrm{c}\mathrm{o}\mathrm{s}\psi }{2}-\alpha \mathrm{s}\mathrm{i}\mathrm{n}\psi & \frac{\mathrm{c}\mathrm{o}\mathrm{s}\psi }{2}+\alpha \mathrm{s}\mathrm{i}\mathrm{n}\psi \\ \rho \left(-{\delta }_{-}-\alpha \mathrm{c}\mathrm{o}\mathrm{s}\psi \right)& \rho \left(-{\delta }_{+}+\alpha \mathrm{c}\mathrm{o}\mathrm{s}\psi \right)\end{array}\right],$(31)

with ${\delta }_{\pm }=\frac{\mathrm{s}\mathrm{i}\mathrm{n}\psi }{2}\pm \frac{d}{l}$, $\alpha =\frac{a+b}{l}$, and $\rho =\frac{r}{d}$. As mentioned above, here we assume that the material motion is along the $x$-axis only, and hence, $\psi (t)=0$, ${\theta }_1(t)={\theta }_2(t)=\theta (t)$, and ${\tau }_1(t)={\tau }_2(t)=\tau (t)$. Moreover, $x_G(t)=x_{O_2}(t)+a$, where $x_{O_2}(t)=r\theta (t)$, and hence, ${\ddot{x}}_G(t)={\ddot{x}}_{O_2}(t)$ and ${\ddot{x}}_G(t)=r\ddot{\theta }(t)$. Thus, the two equations of the system (28) are identical and give the dynamics of $\mathit{G}$

(32) $m_{vl}r{\ddot{x}}_G=2\tau ,$(32)

which is independent of $I_z$ as the motion is purely translational. Moreover, (30) gives ${\dot{\theta }}_3(t)=\dot{\theta }(t)$ and $\dot{\psi }(t)=0$.

3) $c_{2,2}^2$ (Chemical plant): As shown in Figure 7, the anaerobic digestion plant is divided into $P_r$ and $P_d$. The mass balance of $P_r$ in discrete-time yields

(33) $m_r(n+1)-m_r(n)={\delta }_{n_u}(n)m_l,n<n_d,n\in {\bar{\mathbb{Z}}}_{+},$(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\ge 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 modelled 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₄ and carbon dioxide CO₂ (Campos-Rodríguez et al. 2019). The set of four ordinary differential equations of the anaerobic digestion resulting from the mass balance are given by (Campos-Rodríguez et al. 2019)

${\dot{X}}_1(t)=\left[{\mu }_1(S_1(t))-\alpha D_1(t)\right]X_1(t),$

(34) $X_1(0)=X_{1,0},t\ge 0,$(34)

${\dot{X}}_2(t)=\left[{\mu }_2(S_2(t))-\alpha D_2(t)\right]X_2(t),$

(35) $X_2(0)=X_{2,0},$(35)

${\dot{S}}_1(t)=D_3(t)\left(S_{1\mathrm{i}\mathrm{n}}-S_1(t)\right)-k_1{\mu }_1(S_1(t))X_1(t),$

(36) $S_1(0)=S_{1,0},$(36)

${\dot{S}}_2(t)=D_4(t)\left(S_{2\mathrm{i}\mathrm{n}}-S_2(t)\right)$

$+k_2{\mu }_1(S_1(t))X_1(t)-k_3{\mu }_2(S_2(t))X_2(t),$

(37) $S_2(0)=S_{2,0},$(37)

where

(38) ${\mu }_1(S_1(t))={\mu }_{1\mathrm{m}\mathrm{a}\mathrm{x}}\frac{S_1(t)}{S_1(t)+K_{\mathrm{S}1}},$(38)

$D_j(t){|}_{j=1,2,3,4}$ is the dilution rate for the $j$-th state,

(39) ${\mu }_2(S_2(t))={\mu }_{2\mathrm{m}\mathrm{a}\mathrm{x}}\frac{S_2(t)}{S_2(t)+K_{\mathrm{S}2}+{(S_2(t)/K_{\mathrm{I}2})}^2},$(39)

and $\alpha ,S_{1\mathrm{i}\mathrm{n}}$, $S_{2\mathrm{i}\mathrm{n}}$, $k_1$, $k_2$, $k_3$, ${\mu }_{1\mathrm{m}\mathrm{a}\mathrm{x}}$, ${\mu }_{2\mathrm{m}\mathrm{a}\mathrm{x}}$, $K_{\mathrm{S}1}$, $K_{\mathrm{S}2}$, and $K_{\mathrm{I}2}$ are specific constants detailed in Bernard et al. (2001) and Campos-Rodríguez et al. (2019). The input material flow to $P_d$ supplied by $P_r$ is

(40) ${\dot{m}}_{d,\mathrm{i}\mathrm{n}}\left(u\left(t\right)\right)={\rho }_{b}u\left(t\right)V_d,$(40)

where ${\rho }_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{\mu }_2(S_2(t))X_2(t)$ (Bernard et al. 2001). Hence, for $n\ge n_d$, the mass balance of $P_r$ yields

(41) $\frac{\mathrm{d}{m}_r(t)}{\mathrm{d}t}=-{\rho }_b{V}_d\sum _1^4{D}_j(t),$(41)

4) Compartmental control: Here we introduce two compartmental controllers following the black arrow in Step 3 of our methodology (Figure 4): one for $c_{1,2}^3$ 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 ${\dot{x}}_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

(42) $\tau \left(t\right)=u\left(t\right)=\{\begin{array}{cc}\bar{\tau },& t_l\le t<\frac{t_u}{2}\\ -\bar{\tau },& \frac{t_u}{2}\le t<t_u,\end{array}$(42)

where

(43) $\bar{\tau }=\frac{2(m_v+m_l)rH}{t_u^2}.$(43)

To regulate the digester to the desired working point, that is, the equilibrium ‘SS6’ in Campos-Rodríguez et al. (2019), a nonlinear control system was designed and implemented. Specifically, we implemented the control law of Haddad and L’Afflitto (2016) because, for dynamical systems of the form

(44) $\stackrel{\cdot }{x}\left(t\right)=f\left(x\left(t\right)\right)+G\left(x\left(t\right)\right)u\left(t\right),x\left(0\right)=x_0,$(44)

where $\mathit{f}:{\mathbb{R}}^{n_k}\to {\mathbb{R}}^{n_k}$ and $\mathit{G}:{\mathbb{R}}^{n_k}\to {\mathbb{R}}^{n_k\times z_k}$ are continuous functions with $\mathit{f}(\mathbf{0})=\mathbf{0}$, it guarantees the stabilisation of the zero solution $\mathit{x}(t)\equiv \mathbf{0}$ of (44) in a finite-time, where $\mathbf{0}\in {\mathbb{R}}^{n_k}$ is a vector of zeros. In our case, $\mathit{x}(t)=[X_1(t),S_1(t),X_2(t),S_2(t){]}^{\mathrm{T}}$ and $u(t)={[D_1(t),D_2(t),D_3(t),D_4(t)]}^{\text{T}}$. We designed the controller by choosing the functions $L_2\left(\tilde{x}\right)$, $R_2^{-1}\left(\tilde{x}\right)$, and $V\left(\tilde{x}\right)$ (in the notation of Haddad and L’Afflitto (2016)) as

(45) $L_2\left(\tilde{x}\right)=2\left[f^{\mathrm{T}}\left(\tilde{x}\right)G\left(\tilde{x}\right)\right],$(45)

(46) $R_2^{-1}\left(\tilde{x}\right)=G^{-1}\left(\tilde{x}\right)[G^{\mathrm{T}}\left(\tilde{x}\right){]}^{-1},$(46)

and

(47) $V\left(\tilde{x}\right)=p^{\frac{2}{3}}{\left({\tilde{x}}^{\mathrm{T}}\tilde{x}\right)}^{\frac{2}{3}},$(47)

where $\tilde{x}\left(t\right)$ is the state vector translated to have the desired working point corresponding to the zero solution. This choice requires that $G{\left(\tilde{x}\left(t\right)\right)}^{-1}$ exists, but it has the following benefits: first, the condition (33) of Haddad and L’Afflitto (2016) simplifies significantly so that it can be easily checked for systems with complex expressions of $f\left(\tilde{x}\left(t\right)\right)$ and $G\left(\tilde{x}\left(t\right)\right)$ since now it only depends on $V\left(\tilde{x}\right)$; second, since $\mathit{f}(\mathbf{0})=\mathbf{0}$, the condition (34) of Haddad and L’Afflitto (2016) is satisfied regardless of the expressions of $f\left(\tilde{x}\left(t\right)\right)$ and $G\left(\tilde{x}\left(t\right)\right)$; and third, the dynamics of the closed-loop reduce to the gradient system

(48) $\stackrel{\cdot }{\tilde{x}}\left(t\right)=-\frac{1}{2}{\left[V{\hspace{0.17em}}_{\hspace{0.17em}}^{\prime }\left(\tilde{x}\left(t\right)\right)\right]}^{\mathrm{T}},\tilde{x}\left(0\right)={\tilde{x}}_0,$(48)

which is easy to interpret: the state dynamics $\stackrel{\cdot }{\tilde{x}}\left(t\right)$ is proportional to the negative gradient $V{\hspace{0.17em}}_{\hspace{0.17em}}^{\prime }\left(\tilde{x}\right)=\frac{\mathrm{\partial }V\left(\tilde{x}\right)}{\mathrm{\partial }\tilde{x}}$.

The continuous-time mass-flow matrix (3) for this network is

(49) $\begin{array}{cc}\mathbf{\Gamma }(t)& =\left[\begin{array}{cc}{m}_1(t)& {\dot{m}}_{1,2}(t)\\ {\dot{m}}_{2,1}(t)& m_2(t)\end{array}\right]\\ \hspace{0.17em}& =\left[\begin{array}{cc}\sum _{n=0}^{n_f}{m}_1(n)\chi (t-nT)& \sum _{n=0}^{n_f}{\dot{m}}_{1,2}(n)\chi (t-nT)\\ 0& \sum _{n=0}^{n_f}{m}_p(n)\chi (t-nT)\end{array}\right],\end{array}$(49)

where

(50) $\chi (t-nT)=\mathrm{r}\mathrm{e}\mathrm{c}\mathrm{t}\left(\frac{t-\frac{T}{2}-nT}{T}\right),$(50)

$\mathrm{r}\mathrm{e}\mathrm{c}\mathrm{t}(\cdot )$ is the rectangular function defined as (Luise and Vitetta 2009)

(51) $\text{rect}(\sigma )=\{\begin{array}{cc}1& \left|\sigma \right|\text{<}1/2\\ \raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.& \left|\sigma \right|=1/2\\ 0& \text{otherwise},\end{array}$(51)

$T$ is the sample time, $n_f=\frac{t_f}{T}$ is the number of samples and $t_f$ is the final time, $m_p(n)$ is the mass inside the digestion plant updated as

(52) $m_p(n+1)-m_p(n)={\delta }_{n_u}(n)m_l,n\in {\bar{\mathbb{Z}}}_{+},$(52)

$n_u=t_u/T$, and ${\dot{m}}_{1,2}(n)$ is the discrete-time mass flow rate between the hub and the digestion plant (i.e. the mass flow rate of the truck) updated as

(53) ${\dot{m}}_{1,2}(n+1)-{\dot{m}}_{1,2}(n)=\left[{\delta }_{n_l}(n)-{\delta }_{n_u}(n)\right]{\widehat{\dot{m}}}_{1,2},n\in {\overline{\mathbb{Z}}}_{+}.$(53)

Since the transportation by truck is performed in batches, the mass flow rate ${\dot{m}}_{1,2}$ is not defined. Hence, we introduced the pseudo-mass-flow-rate ${\widehat{\dot{m}}}_{1,2}=\frac{m_l}{t_u-t_l}$, whose behaviour is similar to the continuous case. Indeed, an increase in the delivery time $t_u-t_l$ with fixed $m_l$ yields a lower flow, an increase of $m_l$ with fixed $t_u-t_l$ yields a higher flow, and the unit is kg/s.

The entries of the discrete-time form of the mass-flow matrix, that is, of

(54) $\mathbf{\Gamma }(n)=\left[\begin{array}{cc}{m}_1(n)& {\dot{m}}_{1,2}(n)\\ 0& m_p(n)\end{array}\right],n\in {\bar{\mathbb{Z}}}_{+},$(54)

are shown in Figure 9 considering the values in Table 1.

Figure 9. Entries of $\mathbf{\Gamma }\left(n\right)$ achieved using the values in Table 1.

Dynamics of the entries of Γ(n): γ1,1(n) on the top left, γ1,2(n) on the top right, γ2,1(n) on the bottom left, and γ2,2(n) on the bottom right.

Table 1. Values of the parameters used for the numerical study.

Compartment	Parameter	Value	Description
$c_{1,1}^1$	$m_l$	200 kg	Biomass truckload
	$m_{1,0}$	5000 kg	Initial hub stock
	$n_l$	7	Loading time
$c_{1,2}^3$	$m_l$	200 kg	Biomass truckload
	$m_v$	3500 kg	Mass of truck without load
	$I_z$	3000 kg${\mathrm{m}}^2$	Yaw moment of inertia
	$a$	2 m	Chassis posterior length
	$b$	3 m	Chassis anterior length
	$r$	0.4 m	Radius of wheels
	$l$	2 m	Width of chassis
	$d$	0.1 m	Interaxis
	$H$	8000 m	Hub-plant distance
	$t_u$	600 s	Delivery time
	${\dot{m}}_{{1,2}_0}$	0 $\frac{\mathrm{k}\mathrm{g}}{\mathrm{s}}$	Initial truck mass flow rate
	${\theta }_{1,0}$, ${\dot{\theta }}_{1,0}$, ${\theta }_{2,0}$,	0 rad, 0 $\frac{\mathrm{r}\mathrm{a}\mathrm{d}}{\mathrm{s}}$	Initial conditions
	${\dot{\theta }}_{2,0}$, ${\theta }_{3,0}$, ${\psi }_0$
$c_{2,2}^2$	$\alpha ,S_{1\mathrm{i}\mathrm{n}}$, $S_{2\mathrm{i}\mathrm{n}}$, $k_1$,	See Bernard et al. (2001) and	See Bernard et al. (2001) and
	$k_2$, $k_3$, ${\mu }_{1\mathrm{m}\mathrm{a}\mathrm{x}}$, ${\mu }_{2\mathrm{m}\mathrm{a}\mathrm{x}}$,	Campos-Rodríguez et al. (2019)	Campos-Rodríguez et al. (2019)
	$K_{\mathrm{S}1}$, $K_{\mathrm{S}2}$, $K_{\mathrm{I}2}$, $k_6$
	$x_0$	${\left[3.43\frac{\mathrm{g}}{\mathrm{L}},-0.13\frac{\mathrm{g}}{\mathrm{L}},7.19\frac{\mathrm{g}}{\mathrm{L}},3.74\frac{\mathrm{m}\mathrm{m}\mathrm{o}\mathrm{l}}{\mathrm{L}}\right]}^{\mathrm{T}}$	Initial conditions in open loop
	${\tilde{\mathit{x}}}_0$	${\left[2\frac{\mathrm{g}}{\mathrm{L}},-2\frac{\mathrm{g}}{\mathrm{L}},1\frac{\mathrm{g}}{\mathrm{L}},-1\frac{\mathrm{m}\mathrm{m}\mathrm{o}\mathrm{l}}{\mathrm{L}}\right]}^{\mathrm{T}}$	Initial conditions in closed loop
	$p$	1	See (Haddad and L’Afflitto 2016)
	$\bar{D}$	0.5 $\frac{1}{\mathrm{d}\mathrm{a}\mathrm{y}}$	Dilution rate at equilibrium “SS6”
	$m_{p,0}$	0 kg	Initial plant stock
All	$T$	2 s	Sample time
	$t_f$	80 s	Final time

Further simulation results are shown in Figure 10, which were achieved considering the values in Table 1.

Figure 10. Simulation results of the plant digester $P_d$ and the truck.

Digester biomethane flow vs. time in open loop on the top left, digester translated states vs. time in closed loop on the top right, digester biomethane flow vs. time in closed loop on the bottom left, and truck position vs. time on the bottom right.

Observations:

1) Since $-{\delta }_{n_l}(n)m_l\le 0$, it follows from (18) that the truckload $m_l$ is subtracted from the hub biomass stock $m_1(n)$.

2) It follows from (32) that, for a fixed vehicle (i.e. $r$ and $m_v$ are fixed) and fixed acceleration ${\ddot{x}}_G(t)$, an increase of the truckload $m_l$ causes an increase of the engine torque $\tau$. Moreover, the multiplying factor ‘2’ on the right-hand side is because each of the two motors supplies $\tau$.

3) In this example, we assumed a straight trajectory of the vehicle, i.e. $\psi (t)=0$; in general, the road between the hub and the plant may turn, and hence, $\psi (t)\ne 0$.

4) Since ${\dot{\theta }}_2(0)={\dot{x}}_G(0)={\dot{\theta }}_2(t_u)={\dot{x}}_G(t_u)=0$, the truck begins and ends its delivery with null speed (Figure 10d); moreover, $x_G(t_u)=H+a=8002$ m.

5) The gravity force does not appear in the vehicle dynamics (28) because we assumed a planar trajectory. In general, the road may have undulations, and hence, the gravity force contribution depends on the road inclination and the total mass $m_{vl}=m_v+m_l$.

6) It is visible in Figure 9 that the truckload $m_l=200\mathrm{k}\mathrm{g}$ leaves the hub to enter the digestion plant through the truck. Indeed, ${\dot{m}}_{1,2}(n)\ne 0$ if $n_l\le n\le n_u$. The non-existing flow ${\dot{m}}_{2,1}$ is from the plant to the hub and corresponds to the non-existing arc-compartment $c_{2,1}^4$.

7) The total mass of the network is constant since the network is a closed system, i.e.

(55) $m_1(n)+m_2(n)+(t_u-t_l){\dot{m}}_{2,1}(n)=m_{1,0},n\in {\bar{\mathbb{Z}}}_{+}.$(55)

8) The sample time $T$ can be seen as the sample time with which real material flow data are stored and subsequently analysed using MFA. Hence, if a TMN accurately depicts the material flow dynamics, MFA can be performed on the simulated data generated by a TMN. This application of TMNs addresses the issue that the data for MFA are often limited and uncertain (Wang et al. 2022). 9) The controlled system stabilises the biomethane flow of the digester $q_M(t)$ to the desired value in approximately 5 days (Figure 10c), whereas the desired value is reached in approximately 12 days in the uncontrolled system (Figure 10a).

10) Since the mass-flow digraph for this network is equivalent to the one in Figure 3c, it yields $\lambda =0$ (see (8)). Therefore, the flow of ${\beta }_1$ is not thermodynamically circular with this network.

11) The adjunct of an arc-compartment $c_{2,1}^4$ bringing the biomethane from the digester back to the hub yields $\lambda =1$ with respect to the material set $\mathcal{B}=\{{\beta }_1,{\beta }_2\}$, where ${\beta }_1$ is the solid biomass (from $c_{1,1}^1$ to $c_{2,2}^2$) and ${\beta }_2$ is the biomethane (from $c_{2,2}^2$ to $c_{1,1}^1$). Instead, the thermodynamic circularity ${\lambda }^{\prime }$ with respect to $\mathrm{C}{\mathrm{O}}_2$, which is the other material output of the anaerobic digestion, would remain ${\lambda }^{\prime }=0$.

12) We used the mass balance underpinning MFA in (18), in (33), and to formulate the set of four ordinary differential equations of the anaerobic digestion occurring inside $P_d$. How our methodology extends MFA can be seen in the fact that we added the dynamical form of the energy balance (12) to define the dynamics of the truck, and we added the two control systems (42) and (45)–(47) to regulate the behaviour of the truck and the digester, respectively.

13) Since this second example focuses on the dynamics of the network, the network architecture is simpler than in the first example, where the dynamics of each compartment was modelled minimally. Specifically, the first example considers a network with 5 vertices and 5 arcs (10 compartments in total), whereas this example considers a network with 2 vertices and 1 arc (3 compartments in total).

C. Managerial Insight

The first example shows how to measure the material circularity from a physics perspective, that is, it measures to what extent the material flow returns back to where it was previously. Hence, this indicator is not useful to measure the circularity from a life-cycle perspective, that is, to what extent the material is re-used. This is because the indicator proposed here does not look at the material stages such as extraction, manufacturing, use, recycling, landfill, etc., but it looks at its motion in space and time. We started with this physics-based indicator to provide more solid foundations of circularity, which are currently missing in the literature.

In the second example, we focused on the dynamics of a subsystem of a biomethane supply chain. Since our model is dynamic, it can provide details not visible without the use of ordinary differential equations and dynamical energy balances. For example, the vehicle model can take into account the load of biomass being transported and the presence of road undulations; these factors affect the engine torque, which in turn allows us to estimate the fuel consumption and the carbon emissions (if the vehicle is not electric).

The controller we implemented to regulate the anaerobic digestion stabilises the digester to the desired production point in approximately 7 days less than without the use of a controller. Therefore, the controller makes the digestion process both more reliable and quicker than without a controller. Note that the paper is titled ‘[…] and control of circular material flows’ because, as we showed in the second example, the control system can be added to the supply chain model to regulate the material flow.

Note also that the paper title includes the planning of material flows (along with the modelling and control). This is because the design of a TMN allows to simulate different scenarios by varying the model parameters and the tuning of the implemented control systems. Once the TMN behaves as desired, it will show the way to plan its construction, i.e. it will show what are the thermodynamic compartments to be built and how they should be connected in order to circulate the target material. Essentially, a TMN is a tool that is meant to guide the actual design of circular material flows as also shown in Figure 4 (see ‘GOAL’ at the top).

Furthermore, since the data used to perform MFA are often limited and uncertain (Wang et al. 2022), a TMN could be used to generate the data needed. Clearly, the quality of the simulated data depends on the TMN accuracy.

VI. Conclusion

The TMN methodology has been proposed in this paper to extend MFA by adding the power balance, graph theory, and control theory. Our first example showed that a digraph is an effective mathematical formalism to depict the architecture of a material flow network; the knowledge of the architecture, then, allowed us to quantify the material circularity, which is maximum (i.e. $\lambda =1$) when all the flows belong to at least one directed cycle, and it is minimum (i.e. $\lambda =0$) when the digraph has no directed cycles. The mass conservation principle is needed to calculate any unknown weight of the network digraph. The importance of the mass balance is in common with the MFA methodology.

The second example focuses on the dynamics of each compartment and considers a subsystem of a biomethane supply chain to introduce the methodology. The stock in the hub is reduced by the truckload $m_l$ at the loading time $n_l$. The transported biomass $m_l$ affects the torques required to move the three-wheeled truck; the torque could be used to estimate the fuel consumption. The vehicle reaches the delivery point with null speed as desired using a simple control law. If a non-planar and non-linear motion of the truck trajectory is assumed, the vehicle control law is more complex.

The controller implemented to regulate the anaerobic digestion stabilises the system to the desiredbiomethane production approximately 7 days quicker than without controller. The flow of the solid biomass is not thermodynamically circular, i.e. $\lambda =0$. To increase the thermodynamic circularity, further flows should be added, e.g. a new flow from the digester back to the biomass hub. However, in practice, such a new flow may not be useful since the biomethane is produced without being used.

In future work, the general idea is that any system approachable with MFA can be seen as a TMN. The advantages of a TMN approach are the higher model accuracy, higher repeatability (the model accuracy is given in explicitly written dynamical equations rather than being ‘hidden’ in the data) and that it is less data-intensive. As TMN models are scalable (by adding or removing thermodynamic compartments), environmental researchers and engineers could develop and share their networks or compartments to be reused by the research community. For this reason, we published the source code of this paper.

Other future work is the development of a circularity indicator still based on the graph formalism, but that measures the circularity considering the life-cycle stages rather than the loops of material flows as the indicator provided here does. Indeed, increasing material loops may not be always relevant in practice as discussed above for the biomethane supply chain: the adjoint of a flow from the digester to the hub gives maximum thermodynamic circularity, but, in practice, this is not useful.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Notes on contributors

Federico Zocco

Federico Zocco received the B.S. in Mechanical Engineering and the M.S. in Robotics and Automation Engineering from University of Pisa in 2013 and 2016, respectively, and the Ph.D. in applied Machine Learning from Queen’s University Belfast in 2021. He is currently a postdoctoral researcher with the School of Mechanical, Electrical and Manufacturing Engineering and with the Centre for Sustainable Manufacturing and Recycling Technologies at Loughborough University. Dr. Zocco’s research merges computer vision, compartmental dynamical thermodynamics, control theory, and network science with the holistic perspective of industrial ecology to define the theoretical foundations of material flow design. The main area of application is the circularity of critical materials such as the atmospheric carbon dioxide and the rare-earth metals.

Pantelis Sopasakis

Pantelis Sopasak is received the Diploma (M.Eng.) degree in chemical engineering and the M.Sc. Degree (Hons.) in applied mathematics from the National Technical University of Athens (NTU Athens), Athens, Greece, in 2007 and 2009, respectively, and the Ph.D. degree from the School of Chemical Engineering, NTU Athens, in 2012. He held post-doctoral positions at IMT Lucca, Lucca, Italy, from 2013 to 2016, ESAT, KU Leuven, Leuven, Belgium, from 2016 to 2018, and KIOS, University of Cyprus, Nicosia, Cyprus, in 2018. He is currently a Lecturer with the School of Electronics, Electrical Engineering and Computer Science (EEECS), Queen’s University Belfast, where he is a member of the Centre for Intelligent Autonomous Manufacturing Systems (i-AMS). His research interests include the development of model predictive control methodologies, numerical algorithms for large-scale stochastic systems, and massive parallelization on general-purpose graphics processing units (GP-GPUs). The outcomes of his research have several applications in automotive and aerospace, microgrids, water distribution networks, and optimal, safe drug administration.

Beatrice Smyth

Beatrice Smyth is a senior lecturer in the School of Mechanical and Aerospace Engineering, Queen’s University Belfast, and is part of the Research Centre in Sustainable Energy. Her research is focused on environmental impacts and circular economy thinking, and specific areas of interest include energy and carbon life cycle analyses, optimisation of energy pathways, resource quantification and mapping, and economic assessment. Ongoing research projects include nutrient management of digestate combined with energy recovery, the use of short rotation coppice willow to reduce agricultural run-off and improve water quality, sustainability indicators for food, hydrogen for heavy duty transport, microplastic pollution in drinking water, and the sustainable use of plastics in a circular economy. Prior to moving to Queen’s in 2013, Beatrice worked in both the public and private sectors, mainly in energy/carbon management and in geotechnical and environmental engineering. Beatrice has wide experience in working in multi-disciplinary teams, with colleagues from psychology, management, public health, politics and sociology.

Wassim M. Haddad

Wassim M. Haddad received the B.S., M.S., and Ph.D. degrees in mechanical engineering from Florida Tech in 1983, 1984, and 1987. From 1987 to 1994 he served as a consultant for the Structural Controls Group of the Government Aerospace Systems Division, Harris Corporation, Melbourne, FL. In 1988 he joined the faculty of the Mechanical and Aerospace Engineering Department at Florida Tech, where he founded and developed the Systems and Control Option within the graduate program. Since 1994 he has been with the School of Aerospace Engineering at Georgia Tech, where he holds the rank of Professor, the David Lewis Chair in Dynamical Systems and Control, and Chair of the Flight Mechanics and Control Discipline. He also holds a joint Professor appointment with the School of Electrical and Computer Engineering at Georgia Tech. He is the Co-Founder, Chairman of the Board, and Chief Scientific Advisor of Autonomous Healthcare, Inc. Dr. Haddad has made numerous contributions to the development of nonlinear control theory and its application to aerospace, electrical, and biomedical engineering. His transdisciplinary research in dynamical systems and control is documented in over 700 archival journal and conference publications, and 8 books in the areas of science, mathematics, medicine, and engineering. Dr. Haddad is an NSF Presidential Faculty Fellow; a member of the Academy of Nonlinear Sciences; an IEEE Fellow; an AAIA Fellow; and the recipient of the AIAA Pendray Aerospace Literature Award.

Notes

1. https://github.com/fedezocco/TMNbiometh-SciPy.