Topological solvability and DAE-index conditions for mass flow controlled pumps in liquid flow networks

(1)

www.oeaw.ac.at

www.ricam.oeaw.ac.at

Topological solvability and DAE-index conditions for mass flow controlled pumps

in liquid flow networks

A-K. Baum, M. Kolmbauer, G. Offner

RICAM-Report 2016-33

(2)

Topological solvability and DAE-index conditions for mass flow controlled pumps in liquid flow networks

Ann-Kristin Baum^∗† Michael Kolmbauer^† G¨unter Offner^‡ Thursday 29^th September, 2016

Abstract

This work is devoted to the analysis of a model for the thermal management in liquid flow networks consisting of pipes and pumps. The underlying model equation for the liquid flow is not restricted to the equation of motion and the continuity equation, describing the mass transfer through the pipes, but also includes thermodynamic effects in order to cover cooling and heating processes. The resulting model gives rise to a differential-algebraic equation (DAE), for which a proof of unique solvability and an index analysis is presented.

For the index analysis, the concepts of the Strangeness Index is pursued. Exploring the network structure of the liquid flow network via graph theoretical approaches allows to develop network topological criteria for the existence of solutions and the DAE-index.

The topological criteria are explained by various examples.

Keywords: Differential-algebraic equations, topological index criteria, hydraulic network.

AMS(MOS) subject classification: 65L80, 94C15

Introduction

Increasingly demanding emissions legislation specifies the performance requirements for the next generation of products from vehicle manufacturers. Conversely, the increasingly stringent emissions legislation is coupled with the trend in increased power, drivability and safety expectations from the consumer market. Promising approaches to meet these requirements are downsizing the internal combustion engines (ICE), the application of turbochargers, variable valve timing, advanced combustion systems or comprehensive exhaust aftertreatment but also different variants of combinations of the ICE with an electrical engine in terms of hybridization or even a purely electric propulsion. The challenges in the development of future powertrains do not only lie in the design of individual components but in the assessment of the powertrain as a whole. On a system engineering level it is required to optimize individual components globally and to balance the interaction of different sub-systems. A typical system engineering model comprises several sub-systems. For instance in case of a hybrid propulsion these can be the vehicle chassis, the drive line, the air path of the ICE including combustion and exhaust

∗Radon Institute for Computational and Applied Mathematics (RICAM), Austrian Academy of Sciences, Altenberger Str. 69, 4040 Linz, Austria.

†MathConsult GmbH, Altenberger Str. 69, 4040 Linz, Austria.

‡AVL List GmbH, Hans-List-Platz 1, 8020 Graz, Austria.

(3)

aftertreatment, the cooling system of the ICE, the electrical propulsion system including the engine and a battery pack and finally according control systems. Similar to the ICE, the battery pack requires a cooling system as well. Both cooling systems are typically represented by an according hydraulic network in the overall model. The simulation and optimization of hydraulic networks has been studied in various works, including [13, 12, 23, 4, 7] and the references therein. The considered models are motivated by drinking water supply systems, where the main target is to circulate an amount of water at any time, assuring a certain pressure at extraction points. The aim of this work is to consider and analyze hydraulic networks used for thermal management systems. Examples in automotive applications are the above mentioned cooling systems.

In contrast to water transportation networks, the primary interest is not the pressure distribution across the whole system, but the temperature distribution. Consequently, the models have to be equipped with energy balance laws in order to model the thermodynamic effects.

The intention of this work is to extend the results, which are already available for water transportation networks [13, 12] to cooling and heating systems used for thermal management and to networks including mass flow controlled pumps. Thermal flow networks consisting solely of pipes have been analyzed to the work in [2]. The extension to networks of pipes and pumps is not straight forward, since additional to Kirchhoff’s first law also Kirchhoff’s second law has to be considered. Especially Kirchhoff’s second law restricts the allowed pump constellations for a valid liquid flow models.

The model under consideration is a quasi-stationary pipe network, cp. [13], equipped with energy balance laws. This model is suited to describe circuits, which are filled with incompressible fluids (e.g. water). Here incompressible means, that density changes with respect to temperature changes or pressure changes are neglected.

While general networks consist of various types of elements (pipes, pumps, valves) [23], the model here is restricted to pipes and pumps only. Despite this simplification, the demanding issues are caused by the arbitrary network structure of the underlying model. Since valves can change the topology of the underlying network due to there discrete nature, they have to be treated separately.

State-of-the-art modeling and simulation packages such as Dymola¹, Matlab/Simulink², Flow- master³, Amesim⁴, SimulationX⁵or Cruise M⁶offer many excellent concepts for the automatic generation of dynamic system models, including hydraulic networks. Modeling is done in a modularized way, based on a network of subsystems which again consists of simple standardized sub-components. The network structure (topology) carries the core information of the network properties and therefore is predestinated to be exploited for the analysis and numerical simulation of those. In many applications the network describing equations are differential-algebraic equations (DAEs). Hence the analysis of existence and uniqueness of solutions, as well as rank considerations are a delicate issue.

Topology based index analysis for networks connects the research fields ofAnalysis for DAEs [22] andGraph Theory [8] in order to provide the appropriate base to analyze DAEs stemming from automatic generated system models. So far it has been established for various types of

1http://www.dynasim.com

2http://www.mathworks.com

3http://www.mentor.com

4http://www.plm.automation.siemens.com

5http://www.iti.de

6http://www.avl.com

(4)

networks, including electric circuits [24], gas supply networks [10] and water supply networks [13, 12, 23]. Although all those networks share some similarities, an individual investigation is required due to their different physical nature. Recently, a unified modeling approach for different types of flow networks has been introduced in [14], aiming for a unified topology based index analysis for the different physical domains on an abstract level.

The structure of this work is the following. In Section 1, the main two concepts required for the basic ingredients of the analysis are introduced. First an introduction to graph theory is given, then the application to equations imposed on networks is described and the core tools for the following analysis are proven. The network model and arising DAE is formulated in Section 2. Furthermore some basic properties are derived, which lead to the full DAE analysis in Section 3. Beside existence and uniqueness results, DAE-index considerations are performed. Throughout the analysis, the sufficient algebraic conditions are linked to necessary conditions imposed on the network structure. Those topological conditions are explained in term of examples. A summary of the results with comments on their practical relevance in commercial simulation software concludes the paper in Section 4.

1 Graphs and their application in network dynamics

In this section, we introduce the notation and graph theoretical concepts that we need in our analysis and prove some additional results in Lemma 1.1 and Lemma 1.2.

For a detailed introduction to graph theory, we refer the reader to, e.g., [8] or [3]. Agraph G is a pair G= {V,E} of subsets V,E⊂N such that E⊂ V× V, i.e., each element ej ∈ E corresponds to a pair (v_i₁,v_i₂) ∈ V× V, cp. [8, p. 2]. If the pairs (v_i,v_k)∈ Eare ordered, then Gis called anoriented graph, cp. [8, p. 25]. If Gis oriented, then vi and vk are called originating and terminating vertex of ej = (vi,v_k), respectively, [8, p. 25]. If Gcontains no self-loops or parallel edges, then Gis called simple, cp. [8, p. 25].

Two vertices vi,vk ∈ Vare called adjacent if there exists an edge ej ∈ Esuch that ej = (vi,vi2) [8, p. 13]. The edge ej is called incident to vi and vi, respectively [8, p. 13]. Two edges e_j,e_l ∈ Eare called adjacent if they are incident with a common vertex v_i [8, p. 13].

For vi ∈ V, the incident edges are summarized in the set

Einc(vi) :={e_j ∈ E| ∃v` ∈ V: ej = (vi,v`)}.

If Einc(vi) =∅, then vi is isolated and if |E_inc(vi)|= 1, then vi is anend vertex [8, p. 2].

The connection structure of Gis described by theincidence matrix A, which, if Gis oriented, is defined as

A_ij =







1, if vi is the left vertex of ej,

−1, if v_i is the right vertex of e_j, 0, else.

A subset Gs={Vs,Es}with Vs⊂ Vs is a subgraph of Gif Es⊂ Vs×Vs [8, p. 3]. If Vs = V, then Gs spanns G[8, p. 3]. The incidence matrix of Gs is given byAs = [Aij]_(v_i_,e_j)∈Vs×E_s [8, p. 3].

In our analysis, we consider a simple, oriented graph Gwhose vertices Vand edges Eare composed from subsets V1, ...,V_v_ˆand E1, ...,E_ˆ_esuch that V=∪^ˆ^v_I=1V_I, E=∪^e_J=1^ˆ E_J. Accord- ingly, the incidence matrix A is composed from submatrices A_IJ describing the connection

(5)

structure of the subsets G_IJ:={V_I,E_J}. In general, a set G_IJ is not a proper subgraph as the edges EJ may be incident to vertices outside VI. Then, the connection matrix AIJ does not have the usual pattern of two non-zero entries per column. To characterize the fundamental subspaces ofA_IJ, we partition the edges into

EJ=E_IJînner∪E_IJ^loose∪E_IJîsolated, where E_Jînner contains the edges incident to vertices in VI, i.e.,

E_J^inner :={e_j ∈ E_J|e_j = (v_j₁,v_j₂) with v_j₁,v_j₂ ∈ V_I},

E_IJ^loose contains theloose edges incident to a vertex in V_I and a vertex outside V_I, i.e., E_IJ^loose :={ej ∈ EJ|ej = (v_j₁,vj2) with vj1 ∈ VI,vj2 ∈ V\ VI},

and E_IJ^isolated contains theisolated edges incident to vertices outside VI, i.e., E_IJ^isolated:={e_j ∈ E_J|ej = (vj1,vj2) with vj1,vj2 ∈ V\ V_I}.

For simplicity, we assume that there is at most one loose edge per vertex. Using an equivalence relation, the following results can be extend to the case of multiple loose edges per vertex.

Furthermore, we set

V_IJ^outer:= V_I∪ V_IJ^c,

where V_IJ^c contains the vertices outside V_I that are incident to edges in E_J, i.e., V_IJ^c :={vi ∈ V\ VI|Eadj(vi)∩EJ6=∅}.

With this notation, we set

G_IJôuter:={V_IJôuter,E_J}, G_IJînner :={V_I,E_IJînner},

where G_IJôuter is the minimal subgraph containing G_IJ and G_IJînner the maximal subgraph con- tained in GIJ. Using G_IJôuter, G_IJînner we can straightforward extend the standard definitions, cp., e.g., [6, 8], to the set G_IJ.

A subset P:= {V_P,E_P} ⊂ G_IJ,k is called a path in G_IJ if it is a path in G_IJ^outer, i.e., if the vertices in VPare pairwise distinct and there exists a numbering such thatvi,ej are adjacent to vi+1,ej+1 for (i, j)∈N^|V^P^|−1×N^|E^P^|−1, respectively. If Gis oriented, with respect to this numbering, we assign asign to every edge e_j ∈Pby

sgn_P(ej) =

(1, ej = (vj1, vj2),

−1, ej = (vj2, vj1), (1) and define the path matrix P = P

ej∈EPsgn_P(ej)ej, where e1, ..., e|E| ∈ R^|E| denotes the standard canonical basis.

If sgn_P(e_j) = sgn_P(e_l) for e_j,e_l ∈ E_P, then Pis calleddirected. If v_i₁,v_i_|

EP| ∈ V_IJ^c, then P is called acrossing path. If vi1,vi|EP| ∈ VI with vi1 =vi|EP|, then C:=Pis called acycle in G_IJ.

The set G_IJisconnected if E_IJîsolated =∅and if G_IJînner is connected, i.e., if every pair of vertices vi,vk ∈ VI can be connected by a path. If GIJ is not connected, then it is composed of connected components G_IJ,k ={V_IJ,k,E_IJ,k},k= 1, ..., K, containing the connected components G_IJ,kînner of G_IJînner, respectively, plus loose edges E_IJ,k^loose incident to vertices in G_IJ,kînner.

(6)

A subgraph T₁ ⊂ Gof a connected graph Gthat contains no cycles and spans Gis called a spanning tree [8, p. 13]. Every connected graph has at least one spanning tree [8, p. 14].

If G_IJ is connected, then a subset T₁ ⊂ G_IJ is called a spanning tree if T₁ = T₁înner∪ {e₀}, where T₁înner is a spanning tree of G_IJînner and e₀ ∈ E_IJ^loose is a reference loose edge. The complement T2 is called the chord set. For the incidence matrix, the associated edges are selected by the permutation Π = [Π₁,Π₂] with Π_i = [e_j]_e_j∈T_i, i = 1,2, where e₁, ..., e_|E_J_| ∈ R^|E^J^| denotes the standard canonical basis. Every interior chord ek ∈ ET2 ∩E_IJînner defines a uniquefundamental cycle C_k={V_k,E_k}with E_k\ {e_l} ⊂ T₁∩E_IJînner. Thefundamental cycle matrix is defined as C = [C₁, ..., C_c], where C_k is the cycle matrix of C_k. Similarly, every loose chord ek ∈ E_IJ^loose \ {e₀} defines a unique fundamental crossing path Pk starting and ending (or vice versa) in e_k and e₀, respectively. The fundamental crossing path matrix is defined asP_IJ = [P₁, ..., P_|Eloose

IJ |−1],whereP_k denotes the path matrix ofP_k.

If G_IJ is connected and E_IJ^loose =∅, then choosing a ground node v0 ∈ V, we set V2 :={v0} and denote the associated reduced vertex set by V₁ := V\ V₂. If E_IJ^loose 6= ∅, then V₂ = ∅ and the reduced vertex set is given by V1 = VI. For the incidence matrix, these vertices are selected by the permutation Γ = [Γ1,Γ2] with Γi= [e_k]_v_k∈Vi,i= 1,2, where e1, ..., e_|V|∈R^|V^I^| denotes the standard canonical basis.

For a subset Vs ⊂ VI, the vertex identification of Vs merges the elements of Vs into a new vertex ¯v := S

vi∈Vsvi. Removing all edges connecting vertices vi,v_k ∈ Vs, we obtain the contraction of G_IJ with respect to V_s, which is the graph ¯G_IJ := {V¯_I,¯E_J} with ¯V_I := (V_I\

Vs)S

{¯v}and ¯EJ:= EJ\{e_j|ej = (vj1,vj2)|vj1,vj2 ∈ Vs}. Note that ¯Gmight have multiple edges and self loops even if Gis simple. The associated identification matrix is given by1^TΠ, where 1 = [1, ...,1]∈ R^|V^s^| and Π∈ R^|V^s^|×|^V^s^| is a permutation such that [1^TΠ]_i = 1 if and only if vi∈ Vs.

In the following, we assume that GIJ is numbered such that

EIJ =EIJ,1∪. . .∪EIJ,K∪E_IJîsolated, VIJ = VIJ,1∪. . .∪ VIJ,K, (2) where GIJ,k = {EIJ,k,VIJ,k}, k = 1, ..., K, are the connected components of GIJ. These are ordered such that, fork= 1, ...,ˆk1, G_IJ,kcorresponds to proper subgraphs, fork= ˆk1+ 1, ...,ˆk to subsets with loose edges and for k = ˆk+ 1, ..., K to isolated vertices. Accordingly, we denote by G_IJ,kôuter, G_IJ,kînner the corresponding subgraphs of GIJ,k. Then, the connection matrix is given as

A_IJ=





 A_IJ,1

. ..

A_IJ,_k_ˆ 0







with A_IJ,k=

A^inner_IJ,k A^loose_IJ,k

, (3)

where Aînner_IJ,k is the incidence matrix of G_IJ,kînner and A^loose_IJ,k denotes the connection matrix of {V_IJ,k,E_IJ,k^loose}. The zero block row and column correspond to the isolated vertices and edges, respectively. For each G_IJ,k, we number G_IJ,kôuter such that V_IJ,kôuter = V_IJ,k ∪ V_IJ,k^c and the incidence matrix Aôuter_IJ,k is given by Aôuter_IJ,k = [A^T_IJ,k,(A^c_IJ,k)^T]^T, where A^c_IJ,k describes the connection structure of {V_IJ,k^c ,EIJ,k}.

Considering the substructures introduced above, we define the reduced vertex set V_IJ,k,1 :=

∪^K_k=1V_IJ,k,1 withground node V_IJ,k,2 :=∪^K_k=1V_IJ,k,2 with selection matrices

Γ_IJ= [Γ_IJ,1,Γ_IJ,2], (4a)

(7)

the spanning tree T_IJ,k,1 := ∪^K_k=1T_IJ,k,1 with chord set T_IJ,k,2 := ∪^K_k=1T_IJ,k,2 with selection matrices

ΠIJ= [ΠIJ,1,ΠIJ,2], (4b)

thefundamental cycles CIJ,k := ∪^K_k=1CIJ,k,l and the crossing paths PIJ,k :=∪^K_k=1PIJ,k,l with selection matrices

V₂= [C_IJ, P_IJ, L_IJ], (4c) whereV_IJ,k,i,T_IJ,k,i,i= 1,2, C_IJ,k ={C_IJ,k,1, ...,C_IJ,k,|C_IJ,k_|}andP_IJ,k ={P_IJ,k,1, ...,P_IJ,k,|P_IJ,k_|} denote the respective structure in each component GIJ,k with selection matrices Γ_IJ,k,i, Π_IJ,k,i, CIJ,k, PIJ,k, such that ΓIJ,i = diag (ΓIJ,k,i)_k=1,...,K, ΠIJ,i = diag (ΠIJ,k,i)_k=1,...,K, i= 1,2 and P_IJ= diag (P_IJ,k)_k=1,...,K,C_IJ= diag (C_IJ,k)_k=1,...,K andL_IJ= [0, I_|Eloose

IJ,k |]^T. For the connected components that are proper subgraphs, we further consider theidentification V¯IJ =Sk^ˆ

k=1v¯ik

with ¯vi_k :=vi_k,0∪(∪_v_i_∈_V_I,kvi) and identification matrix 1_IJ= diag

1_|V_IJ,k_|

k=1,...,ˆk, (4d)

Like for a proper graph and its incidence matrix, cp., e.g., [8, pp. 23], we can interpret the fundamental subspaces of a submatrix AIJ as substructures of the set GIJ.

Lemma 1.1. Let G={V,E} be a simple, oriented graph with V=∪_I∈I_VV_I, E=∪_J∈J_EE_J. Consider a subset G_IJ:={V_I,E_J}with connection matrixA_IJ. Then,rank(A_IJ) =P^ˆk

k=1|V_IJ,k|−

k, whereˆ GIJ,k, k= 1, ...,kˆ denotes the connected components of GIJ that itself are subgraphs.

For the matrices defined in (4), it holds thatker(AIJ) = span(V2),corange(AIJ) = span(ΠIJ,1) and coker(A_IJ) = span(1_IJ), range(A_IJ) = span(Γ_IJ,1).

The matrices U := [Γ_IJ,1,1_IJ] andV := [Π_IJ,1, V_IJ,2] are nonsingular with U⁻¹ =

U₂⁻ Γ^T_IJ,2

, V⁻¹ =







V₂⁻ 0 Π^T_IJ,2 0

0 I_|E^loose

IJ,k |





, where U₂⁻= Γ^T_IJ,1−1IJΓ^T_IJ,2 and V₂⁻= Π^T_IJ,1−Π^T_IJ,1[CIJ, PIJ]Π^T_IJ,2. Proof. From (3), we get that rank(A_IJ) =Pk^ˆ

k=1rank(A_IJ,k). Noting that G_IJ,kôuter is connected as G_IJ,k is connected and using that V_IJôuter= V_I∪V_IJ^c, we have that rank(Aôuter_IJ,k ) =|V_IJ,k|+

In conclusion, we have proven that rank(A_IJ) = P^ˆk1

k=1(|V_IJ,k| −1) +Pk^ˆ

k=ˆk1+1|V_IJ,k|, i.e., rank(A_IJ) =P^ˆk

k=1|V_IJ,k| −ˆk1.

Now, we consider a connected component GIJ,k. With the given numbering, the fundamental cycle and crossing path matrices are given by

C_IJ,k =



 C_IJ^inner

0 0



, P_IJ,k =





∗ 1|E_IJ,k|−1

I_|E_IJ,k_|−1



, (5)

(8)

whereC_IJ,kînner denotes the fundamental cycle matrix of G_IJ,kînner. From (3) and (5) it follows that AIJ,kCIJ,k=Aînner_IJ,k C_IJ,kînner= 0,

as the fundamental cycles of G_IJ^inner lie in ker(A_IJ). Similarly, we get from (3) and (5) that

A^outer_IJ,k P_IJ,k = 0

A^c_IJ,kP_IJ,k

=



 0 1|E_IJ,k|−1

I_|E_IJ,k_|−1





as the incidence matrix applied to a path matrix returns exactly the starting and end vertices of the path, cp. [6, p. 157]. Thus, A_IJ,kP_IJ,k = 0, implying that span([C_IJ,k, P_IJ,k]) ⊂ ker(AIJ,k). Considering (5), we find that rank([CIJ,k, PIJ,k]) = rank(CIJ,k) + rank(PIJ,k) with rank(CIJ,k) = |E_IJ,k^inner| − |VIJ,k|+ 1 and rank(PIJ,k) = |E_IJ,k^loose| −1. For k = 1, ...,ˆk1, we thus get that

rank([C_IJ,k, P_IJ,k]) =|E_IJ,k| − |V_IJ,k|+ 1 = dim(ker(A_IJ,k)), while fork= ˆk₁+ 1, ...,k, we get thatˆ

Hence, span([C_IJ,k, P_IJ,k]) = ker(A_IJ,k) and it follows that span(V₂) = ker(A_IJ).

For the left nullspace, we note that fork= 1, ...,k,ˆ GIJ,k is a proper subgraph, implying that every column ofAIJ,k contains exactly two nonzero entires 1,−1. Hence,1^T_|_V

If ker(AIJ) = span(V2), corange(AIJ) = span(ΠIJ,1) and coker(AIJ) = span(1IJ), range(AIJ) = span(Γ_IJ,1), then the matricesU := [Γ_IJ,1,1_IJ] andV := [Π_IJ,1, V₂] are nonsingular. To verify the representation of U⁻¹, V⁻¹, we check the properties of the inverse. For convenience, we drop the index IJ of the subset GIJ. First, we note that

Γ₁ 1

Γ^T₁ −1Γ^T₂ Γ^T₂

= Γ₁Γ^T₁ + (1−Γ₁1_IJ)Γ^T₂ = Γ₁Γ^T₁ + Γ₂1Γ^T₂.

Noting that [Γ₁1]_i = 1 for v_i ∈ V_IJ,1 and [Γ₁1]_i = 0 for v_i∈ V_IJ,2, we have that1−Γ₁1= Γ₂, such that

Γ₁ 1

Γ^T₁ −1Γ^T₂ Γ^T₂

= Γ₁Γ^T₁ + Γ₂Γ^T₂ =I_|V_IJ_|.

(9)

On the other hand, we have that Γ^T₁ −1Γ^T₂

Γ^T₂

Γ₁ 1

=

I_|V_red_|−1 (Γ^T₁ −1Γ^T₂)1

0 Γ^T₂1

.

In Γ^T₁ −1Γ^T₂, every row contains exactly two non-zero entries given by 1,−1, implying that (Γ^T₁ −1Γ^T₂)1= 0. With Γ^T₂1_IJ=I_|V_IJ,2_|, we thus get that [Γ1,1]⁻¹[Γ1,1] =I_|_V_IJ_|.

As E(T)∩E_IJ^loose=∅, where E(T) denotes the edges of a spanning tree T, we have that Π₁ V₂

=

"

Π˜₁ C P

0 0 0 I_|E^loose

IJ,k |

# . Thus, it suffices to show that

Π˜₁ [C, P]−1

=

Π^T₁ −Π^T₁[C, P]Π^T₂ Π^T₂

.

We show the assertion by checking the properties of the inverse. First, using that I|V_IJ|− Π₁Π^T₁ = Π₂Π^T₂, we get that

Π₁ [C, P]

= Π₁Π^T₁ −Π₁Π^T₁[C, P]Π^T₂ + [C, P]Π^T₂

= Π₁Π^T₁ + Π₂Π^T₂[C, P]Π^T₂.

The matrix Π₂Π^T₂ is a projection onto the edges of the chord set T_IJ,2. As the fundamental cycles and crossing paths Ck_l,Pkm contain exactly one edge ek_l,ekm ∈ TIJ,2, we have that Π₂Π^T₂[C, P] = Π₂. Hence,

Π1 C P

= Π1Π^T₁ + Π2Π^T₂ =I_|E_IJ_|. On the other hand, we have that

Π₁ [C, P]

=

I_|E_(T_IJ_)| Π^T₁[C, P](I_n_J−n_I+1−Π^T₂[C, P])

0 Π^T₂[C, P]

. Again, as every fundamental cycle and crossing path contains exactly one chord, we have that Π^T₂[C, P] =I_n_J−n_I+1. Then, it follows that [Π₁,[C, P]]⁻¹[Π₁,[C, P]] =I_|E_IJ_|.

Hence, the fundamental cycles, crossing paths and loose edges spann the right nullspace ker(A_IJ), while the edges in the spanning tree build a basis of corange(A_IJ). The identification of connected components with a ground node spanns the left nullspace coker(AIJ), while the vertices of the reduced vertex set build a basis of corange(AIJ).

Now, we equip the vertices and edges of Gwith potentials and flows, respectively. To each vertex vi ∈ V, we assign a potentialpi that are summarized asp= [pi]i=1,...,|V|. Similarly, to each edgeej ∈ G, we assign aflow qj and setq= [qj]_j=1,...,|E_|. The flow is directed: A flowqj

is calledpositive, i.e.,q_j >0, ifq_j agrees with the direction of its associated edge e_j. Ifq_j 6= 0 is opposed to the direction of edge ej, thenqJ is called negative, i.e.,qj <0. If G={V,E}

with V=∪_I∈I_VV_I, E=∪_J∈I_EE_J, we partition the flows and potential accordingly and write q_J,j ∈ E_J and p_I,i∈ V_I.

(10)

The flow and potential satisfy the following fundamental relations that generalize Kirchhoff’s circuit laws, i.e.,

Γ^T₁Aq= 0, (6a)

V₂^TA^Tp= 0. (6b)

The equations (6) allow to give a physical interpretation of Lemma 1.1. The fundamental cycles, crossing paths and loose edges correspond structures on which the potential difference vanishes, while the spanning tree selects a structure on which the potential difference is well- defined. The reduced vertex set denotes those vertices on which the potential is fixed in relation to the reference value given by the ground node. Thus, the identification of the reduced vertex set with its ground node summarizes all vertices on which the potential is not fixed, like in an isolated vertex or a subgraph without connection to a ground node.

We transform the flow and potential with respect to these substructures by setting

˜

q := [ΠIJ,1, VIJ,2]⁻¹q, p˜:= [ΓIJ,1,1IJ]⁻¹p, (7) such that

˜

q1 = (Π^T_IJ,1−Π^T_IJ,1[CIJ, PIJ]Π^T_IJ,2)q, p˜1= (Γ^T_IJ,1−1IJΓ^T_IJ,2)p,

˜

q₂ = Π^T_IJ,2q, p˜₁= Γ^T_IJ,2p.

The flows q₂ belong to edges on the chord set T_IJ,2, while the flowsq₁ denote the difference between a branch flowq_1,j ∈TIJ and the chord flows q_2,l ∈TIJ,2 associated with fundamental cycles and crossing paths q1,j is part of. Similarly, p2 denote the potentials in the ground nodes V_IJ,2 while p₁ correspond the difference between a potential p_1,j ∈ V_IJ,1 and to its associated ground nodep_2,j ∈ VIJ,2.

Now, we think of the flow as information running through the network. To describe the structure of the subset GIJ on this informational level, we partition the set of incident edges Einc(vi), vi∈ VI, into those along which vi receives andsends information, respectively, i.e., we set

E_inc,s(v_i) :={e_j ∈E_inc(v_i)|A_ijsgn(q_j(t))>0, i.e., q_j starts in v_i}, E_inc,e(v_i) :={e_j ∈E_inc(v_i)|A_ijsgn(q_j(t))<0, i.e., q_j ends in v_i}.

Defining theflow matrix

BIJ,i`=





 P

ej∈Einc,s(vi)∩EJ|q_j|, i=`,

−|q_j|, ej ∈ Einc,e(vi)∩Einc(v_`)∩E_J, i6=`,

0, else,

(8)

the information flow in G_IJis graphically represented by theflow graph G_IJ^flow := G(B_IJ^T),where the graph G(A) of a matrixA∈R^n×nis defined as G(A) =

{v₁, ...,v_n},{(v_i,v_j)|a_ij 6= 0} , i.e., whenever theij-th entry is nonzero, there is an edge from vertex vitovj, cp. [20, p. 528].

Hence,

G_IJ^flow={V_I,E_IJ^flow} (9) E_IJ^flow:={e_j := (v_`,v_i)} |E_inc,e(v_i)∩E_inc(v_`)6=∅,v_`6= v_i ∨ E_inc,s(v_i)6=∅,v_` =v_i}.

(11)

The graph G_IJ^flowbasically has the same connection structure as the set G_IJ. At vertices v_i∈ V_I sending a non-zero mass flow intoGIJ, however, G_IJ^flowhas self loops and edgesej ∈ EJequipped with a zero mass flowq_j = 0 are absent in E_IJ^flow. The orientation of G_IJ^flow is determined by the direction of the mass flows, where e_j ∈ E_IJ^flow is directed from v_` to v_i if v_i receives a mass flow from v`, cp. Figure 1.

Figure 1: A graph (left) and its flow graph (right).

vI,1

vI,2

vI,3

vK,2

vK,1

vI,6

vK,3

q1>0

q2>0 q3>0

q6= 0 q4>0

q5<0 q7>0

vI,1

vI,2

vI,3

vK,2

vK,1

vI,6

vK,3

q1

q2

q3

q4

q5

q7

The connectivity of G_IJ^flow on this informational level is described by the concept of strong connectivity, cp. [20, p. 528], i.e., we assume that G_IJ^flow is composed from strongly connected components Gf low,IJ,k:={Ef low,IJ,k,Vf low,IJ,k}, i.e., every pair of vertices vi,vk∈ Vf low,IJ,k is connected by a directed path from v_i tov_kand a directed path fromv_k tov_i, cp. [20, p. 528].

For each G_{f low,IJ,k}, we denote the interior subgraph by G_{f low,IJ,k}înner ={E_{f low,IJ,k}înner ,V_{f low,IJ,k}înner }for k= 1, ..., K.

The flow matrix B_IJ is nonsingular, if in every strongly connected component G_IJ,k of G_IJ there exists at least one vertex v_ˆ_i sending a nonzero flow into G_IJ\ G_IJ,k, i.e., outside the strongly connected component.

Lemma 1.2. Let G={V,E} be a simple, oriented graph with V=∪_I∈I_VV_I, E=∪_J∈J_EE_J. Consider a subset GIJ:={VI,EJ} with connection matrix AIJ. Then,

B_IJ= 1

2A_IJdiag (q_J,j)_j diag (sgn(q_J,j))_jA^T_IJ+|A^T_IJ| . For v_i ∈ V_I, if P

ej∈Einc,s(vi)∩EJ|q_j|>0 and X

ej∈Einc,s(vi)∩(EJ\E_{IJ,f low,k}^inner )

|q_j| ≥0, k= 1, ..., K, (10)

where G_{f low,IJ,k} :={E_{f low,IJ,k},V_{f low,IJ,k}} denotes the strongly connected components of G_IJ^flow with interior subgraphs G_{f low,IJ,k}înner = {E_{f low,IJ,k}înner ,V_{f low,IJ,k}înner }, and for every k = 1, ..., K there exists vî∈ Vf low,IJ,k, such that (10) is strictly satisfied, then BIJ is nonsingular.

Proof. To prove the representation ofB_IJ, we set ˜A_IJ:=A_IJdiag (q_J,j)_j diag (sgn(q_J,j))_jA^T_IJ+

|A^T_IJ|

and note that

B_IJ,i` := X

ej∈EJ

|q_j|A_IJ,ij A_IJ,`j+|A_IJ,`j|sgn(q_j) .

(12)

For v_i,v_`∈I_Vand j∈J_E, the entries of the incidence matrix A satisfy

A_ijA_`j =







1, ej ∈ Einc(v_i), i=`,

−1, ej ∈ Einc(vi)∩Einc(v_`), i6=`, 0, else,

A_ij|A_`j|=

(Aij, ej ∈ Einc(vi)∩Einc(v_`), 0, else

and together with the definition of Einc,s(vi),Einc,e(vi), we verify that 2BIJ= ˜AIJ.

IfG_{f low,IJ}is composed fromKstrongly connected components G_{f low,IJ,k} :={V_{f low,IJ,k},E_{f low,IJ,k}} it follows that BIJ is congruent to a block upper triangular matrix with irreducible diagonal blocks B_IJ,kk, k= 1, ..., K, i.e., there exists a permutation Π, such that Π^TBIJΠ = [B_IJ,kl]_kl withB_IJ,kl = 0,k > l and B_IJ,kk irreducible, cp. [20]. We show that under the given conditions, each B_IJ,kk is irreducibly diagonally dominant and hence nonsingular for k= 1, ..., K, cp., e.g., [11, p. 403] or [5, p. 67]. Then,BIJ is nonsingular.

Thei-th column ofB_IJ,kk is given by

[B_IJ,kke_i]_`=





 2P

ej∈Einc,s(vi)∩EJ|q_j|, i=l,

−2|q_j|, ej ∈ Einc,s(vi)∩Einc(v`)∩EJ, i6=`,

0, else,

`= 1, ...,|V_{f low,IJ,k}|,

and noting that [

v`∈V_{IJ,f low,k}\{v_i}

E_inc,s(v_i)∩E_inc(v_`)∩E_J

= E_inc,s(v_i)∩ E_{IJ,f low,k}^inner

for vi∈ VIJ,f low,k, the i-th column sum of BIJ,kk is given by X

v`∈VIJ,f low,k

|B_IJ,kk,i`|= X

ej∈Einc,s(vi)∩(E_J\E_{IJ,f low,k}^inner )

|q_j|.

Hence, if condition (10) is satisfied forvi ∈ Vf low,IJ,kandk= 1, ..., K, thenB_IJ,kkis diagonally dominant fork = 1, ..., K. For k = 1, ..., K, if there exists vˆi ∈ Vf low,IJ,k, such that (10) is strictly satisfied, then, together with the irreducibility, it follows that B_IJ,kk is irreducible diagonally dominant for k= 1, ..., K. Then,B_IJ is nonsingular.

Example 1.1. For the graph of Figure 1, the flow matrix is given by

B_IJ=







v_I,1 v_I,2 v_I,3 v_I,6 v_I,1 |q₁|+|q₄| 0 −|q₃| 0 vI,2 −|q₁| |q₂| 0 0 vI,3 0 −|q₂| |q₃| 0 vI,6 0 0 0 |q₇|.







The strongly connected components ofG_IJ^floware given byE_{f low,IJ,1}^inner ={{vI,1,vI,2,vI,3},{e₁,e2,e3}}

and E_{f low,IJ,2}^inner ={{v_I,4}}. The matrix B_IJ is irreducible as there exists no permutation that

(13)

would transform this matrix to an upper triangular matrix and we observe that B_IJ is nonsingular only if|q₄|,|q₇|>0, i.e., only if the flows |q₄|,|q₇|>0 start in vI,1, vI,6. This agrees with the conditions of Lemma 1.2, claiming that this is reflected in the solvability condition, claiming that P

ej∈Einc,s(vI,1)∩(EJ\E_{IJ,f low,k}^inner )|q_j|= |q₄| >0,, P

ej∈Einc,s(vI,6)∩(EJ\E_{IJ,f low,k}^inner )|q_j|=

|q₇|>0,whereas P

ej∈Einc,s(vI,i)∩(E_J\E_{IJ,f low,k}^inner )|q_j|= 0 for i= 1,2,3.

Besides these graph theoretical results, we frequently use the following identity for the rank of a block matrixA= [A_ij]_i,j=1,2, cp., e.g., [11, p. 25]. If A₁₁, A₂₂ are nonsingular, then

rank(A) = rank(A₁₁) +S_A₁₁(A) = rank(A₂₂) +S_A₂₂(A). (11) whereSA11(A) :=A22−A21A⁻¹₁₁A12,SA22(A) =A11−A12A⁻¹₂₂A21 denotes the Schur comple- ments.

2 A network model for incompressible flow networks

We consider a network

N={Pi,Pu,De,Jc,Re} (12) that is composed of pipes, pumps, demands, junctions and reservoirs and that is filled by an incompressible fluid, e.g. water. The pipes Pi := {Pi₁, ...,Pi_n_Pi} and pumps Pu :=

{Pu₁, ...,Pu_n_Pu} are connected by junctions Jc := {Jc₁, ...,Jc_n_Jc}, in which the mass flow of the fluid is split or merged. We distinguish between virtual connection points Jc₀ :=

{Jc_0,1, ...,Jc_0,n_Jc,0} and connection points Jc_V := {Jc_V,1, ...,Jc_V,n_Jc,V} posessing a volume V_i > 0. Those virtual connections point have a certain importance in the design of system simulation software, since they allow to connect standardized sub-components without introducing additional volumes (and as a consequence additional thermal inertia). The connection to the environment is modeled by reservoirs Re := {Re₁, ...,Re_n_Re} and demand branches De:= {De₁, ...,DenDe} that impose predefined pressures enthalpies as well as mass and enthalpy flows into the network. The number of each element in N is de- noted by n_Pi, n_Pu, n_Jc, where n_Jc =n_Jc,0 +n_Jc,V and n_De and n_Re, respectively, and we set n:=nPi+nPu+nJc+nRe+nDe.

Given boundary conditions ¯p_Re= [¯p_Re_i]i=1,...,n_Re, ¯h_Re= [¯h_Re_i]i=1,...,n_Reand ¯q_De = [¯q_De_i]i=1,...,n_De, H¯_De = [ ¯H_De_i]_i=1,...,n_De, the task is to compute the mass and enthalpy flowsq_Pi= [q_Pi_i]_i=1,...,n_Pi, qPu = [qPui]i=1,...,nPu, qDe = [qDei]i=1,...,nDe and HPi = [HPii]i=1,...,nPi, HPu = [HPui]i=1,...,nPu, HDe = [HDei]i=1,...,nDe in the pipes, pumps and demand branches as well as the pressures and specific enthalpies p_Jc = [p_Jc_i]_i=1,...,n_Jc, p_Re = [p_Re_i]_i=1,...,n_Re and h_Jc = [h_Jc_i]_i=1,...,n_Jc, hRe= [hRei]i=1,...,nRe in the junctions and reservoirs.

To set up the governing equations for the network, we consider the characteristic relation that every element imposes on the enthalpy flow and specific enthalpy as well as on the mass flow and pressure. In a pipe Pij, the enthalpy flow Hj agrees with the product of the mass flow qj and the specific enthalpyhi in the originating vertex vi. Hence, if the pipe Pij is directed from v_j₁ to v_j₂, we have that

H_j = q_Pi,j

2 ((sgn(q_j) + 1)h_j₁−(sgn(q_j)−1)h_j₂) =:f_Pi^∗(q_Pi,j, h_j₁, h_j₂). (13a)