A Lumped-Parameter Model of Multiscale Dynamics in Steam Supply Systems

TL;DR

This paper develops a simplified lumped-parameter model for multiscale steam supply system dynamics, capturing slow and fast processes, and validates it through theoretical analysis and numerical comparison with detailed simulations.

Contribution

It introduces a physically-derived lumped-parameter model that preserves multiscale dynamics and provides a theoretical basis for separating short-term and long-term system operations.

Findings

The model accurately captures multiscale dynamics of steam systems.

Existence of a normally hyperbolic invariant manifold confirms the model's validity.

Quantitative comparison shows the model's predictions align with detailed simulations.

Abstract

This paper focuses on multiscale dynamics occurring in steam supply systems. The dynamics of interest are originally described by a distributed-parameter model for fast steam flows over a pipe network coupled with a lumped-parameter model for slow internal dynamics of boilers. We derive a lumped-parameter model for the dynamics through physically-relevant approximations. The derived model is then analyzed theoretically and numerically in terms of existence of normally hyperbolic invariant manifold in the phase space of the model. The existence of the manifold is a dynamical evidence that the derived model preserves the slow-fast dynamics, and suggests a separation principle of short-term and long-term operations of steam supply systems, which is analogue to electric power systems. We also quantitatively verify the correctness of the derived model by comparison with brute-force simulation of the original model.

Tables

Table 1. Table 1: List of variables and parameters in the derived model ( 14 ). The values used for numerical simulations in Sec.  5 are also presented.

Meaning	Symbol	(Nominal) Value	Scaled value
Pressure of steam	$p$	$800\mathrm{kPa}$	1.0 (base)
Density of saturated steam	${\rho }_s$	$4.16\mathrm{kg}/m^3$	1.0 (base)
Density of saturated water	${\rho }_w$	$897\mathrm{kg}/m^3$
Specific enthalpy of saturated steam	$h_s$	$2768\mathrm{kJ}/\mathrm{kg}$	14.3
Specific enthalpy of saturated water	$h_w$	$721\mathrm{kJ}/\mathrm{kg}$	3.74
Specific enthalpy of the feed water	$h_f$
Mass flow rate of saturated steam from drum	$m_s^{\prime }$
Mass flow rate of feed water to drum	$m_f^{\prime }$
Coefficient of pressure variation given in (2)	$e$	$3073J/\mathrm{Pa}$	$1.8$
Temperature of saturated steam	$T_s$	$443K$
Total mass of evaporator and drum of a boiler	$m_t$	$50,000\mathrm{kg}$
Specific heat of the metal of boiler	$C_p$	$0.4\mathrm{kJ}/\left(K\cdot \mathrm{kg}\right)$
Total volume of steam	$V_s$	$10.0m^3$
Total volume of water	$V_w$	$10.0m^3$
Velocity of steam in the pipe	$u$	$30m/s$	1.0 (base)
Length of the steam pipe	$L$	$200m$	1.0 (base)
Diameter of the steam pipe	$d$	$0.2m$	1.0 (base)
Friction coefficient of the steam pipe	$\lambda$	0.016	16
Input rate of heat of a boiler	$Q^{\prime }$
Consumption rate of heat at a load	$Q_L^{\prime }$	$5.0\mathrm{MJ}/s$	5.2

Equations

$$e_v(p_v)\frac{\mathrm{d}{p}_v}{\mathrm{d}t}={ϵ}_1\left\{Q_v^{\prime }-m_{\mathrm{s}v}^{\prime }{h}_{\mathrm{c}}(p_v)\right\},$$

$$e_v(p_v)=$$

$$\frac{\partial {\rho }_l}{\partial t}+\frac{\partial }{\partial x}({\rho }_l{u}_l)=0,$$

$$\frac{\partial }{\partial t}({\rho }_l{u}_l)+\frac{\partial }{\partial x}({\rho }_l{u}_l^2)+\frac{1}{{ϵ}_2}\frac{\partial p_l}{\partial x}+{\lambda }_l\frac{{\rho }_l{u}_l|u_l|}{2d_l}=0,$$

$$\frac{\partial }{\partial t}({\rho }_l{h}_l)+\frac{\partial }{\partial x}({\rho }_l{h}_l{u}_l)=\frac{\partial p_l}{\partial t}+{ϵ}_3{Q}_{\mathrm{w}l},$$

$${\rho }_l\frac{\partial u_l}{\partial t}+\frac{1}{{ϵ}_2}\frac{\partial p_l}{\partial x}+\frac{{\lambda }_l{\rho }_l{u}_l|u_l|}{2d_l}=0.$$

$$L_l{\rho }_{\mathrm{av}l}\frac{\mathrm{d}{u}_l}{\mathrm{d}t}=\frac{p_{{\partial }^{+}l}-p_{{\partial }^{-}l}}{{ϵ}_2}-\frac{{\lambda }_l{\rho }_{\mathrm{av}l}{L}_l{u}_l|u_l|}{2d_l},$$

$${\rho }_{\mathrm{av}l}(t):=\frac{1}{L_l}{\int }_0^{L_l}{\rho }_l(x,t)\mathrm{d}x.$$

$$m_{\mathrm{s}v}^{\prime }(t)=m_{\mathrm{L}v}^{\prime }(t)$$

$$m_{\mathrm{s}v}^{\prime }(t)h_{\mathrm{s}}(p_v(t))=m_{\mathrm{L}v}^{\prime }(t)h_{\mathrm{L}v}(t)$$

$$-\sum _{l\in {\delta }^{-}v}\frac{\pi d_l^2}{4}{\rho }_l(L_l,t)u_l(t)h_l(L_l,t),$$

$$m_{\mathrm{L}v}^{\prime }{h}_{\mathrm{L}v}=Q_{\mathrm{L}v}^{\prime }+m_{\mathrm{L}v}^{\prime }{h}_{\mathrm{w}}(p_v).$$

$$m_{\mathrm{s}v}^{\prime }{h}_{\mathrm{c}}(p_v)=Q_{\mathrm{L}v}^{\prime }$$

$$-\sum _{l\in {\delta }^{-}v}\frac{\pi d_l^2}{4}{\rho }_l(L_l,t)(h_l(L_l,t)-h_{\mathrm{w}}(p_v))u_l.$$

$$h_l(0,t):=h_{\mathrm{s}}(p_{{\partial }^{+}l}(t)),h_l(L_l,t):=h_{\mathrm{s}}(p_{{\partial }^{-}l}(t)),$$

$${\rho }_l(0,t):={\rho }_{\mathrm{s}}(p_{{\partial }^{+}l}(t)),{\rho }_l(L_l,t):={\rho }_{\mathrm{s}}(p_{{\partial }^{-}l}(t)),$$

$${\rho }_{\mathrm{av}l}(t):=\frac{{\rho }_{\mathrm{s}}(p_{{\partial }^{+}l}(t))+{\rho }_{\mathrm{s}}(p_{{\partial }^{-}l}(t))}{2}.$$

$$\frac{\mathrm{d}{p}_v}{\mathrm{d}t}=\frac{ϵ}{e_v(p_v)}\left\{Q_v^{\prime }-Q_{\mathrm{L}v}^{\prime }-\sum _{l\in {\delta }^{+}v}\frac{\pi d_l^2}{4}{h}_{\mathrm{c}}(p_v){\rho }_{\mathrm{s}}(p_v)u_l+\sum _{l\in {\delta }^{-}v}\frac{\pi d_l^2}{4}{h}_{\mathrm{c}}(p_v){\rho }_{\mathrm{s}}(p_v)u_l\right\},$$

$$\frac{\mathrm{d}{u}_l}{\mathrm{d}t}=\frac{1}{ϵ}\frac{2(p_{{\partial }^{+}l}-p_{{\partial }^{-}l})}{L_l\{{\rho }_{\mathrm{s}}(p_{{\partial }^{+}l})+{\rho }_{\mathrm{s}}(p_{{\partial }^{-}l})\}}-\frac{{\lambda }_l}{2d_l}{u}_l|u_l|,$$

$$p_v(t,ϵ)=p_v^{(0)}(t)+ϵp_v^{(1)}(t)+O({ϵ}^2),$$

$$u_l(t,ϵ)=u_l^{(0)}(t)+ϵu_l^{(1)}(t)+O({ϵ}^2).$$

$$\frac{\mathrm{d}{p}_v^{(0)}}{\mathrm{d}t}=0,p_{{\partial }^{+}l}^{(0)}=p_{{\partial }^{-}l}^{(0)}.$$

$$\mathbf{\psi }$$

$$R_{ij}:=\{\begin{array}{ll}1& \mathrm{if}{v}_i={\partial }^{+}{l}_j\ne {\partial }^{-}{l}_j,\\ -1& \mathrm{else}\mathrm{if}{v}_i={\partial }^{-}{l}_j\ne {\partial }^{+}{l}_j,\\ 0& \mathrm{otherwise}.\end{array}$$

$$\left[\begin{array}{cc}\mathsf{G}& 0\\ 0& \mathsf{H}\end{array}\right]\frac{\mathrm{d}}{\mathrm{d}t}\left[\begin{array}{c}\mathbf{\psi }\\ \mathbf{q}\end{array}\right]=\left[\begin{array}{cc}0& -\mathsf{R}\\ {\mathsf{R}}^{\top }& 0\end{array}\right]\left[\begin{array}{c}\mathbf{\psi }\\ \mathbf{q}\end{array}\right]+\left[\begin{array}{c}\mathbf{s}\\ -\mathbf{f}(\mathbf{q})\end{array}\right],$$

$$\mathsf{G}$$

$$\mathsf{H}$$

$$\mathbf{s}$$

$$\mathbf{f}(\mathbf{q})$$

$$\mathsf{R}{\mathbf{q}}^{\ast }=\mathbf{s},$$

$${\mathsf{R}}^{\top }{\mathbf{\psi }}^{\ast }=\mathbf{f}({\mathbf{q}}^{\ast }).$$

$$\mathbf{s}\in \mathrm{Im}(\mathsf{R})=(\mathrm{Ker}({\mathsf{R}}^{\top }){)}^{\perp }.$$

$$Q_1^{\prime }+\cdots +Q_n^{\prime }=Q_{\mathrm{L}1}^{\prime }+\cdots +Q_{\mathrm{L}n}^{\prime }.$$

$${\mathbf{\psi }}^{\prime }=\mathbf{\psi }-\left(\frac{s_1+\cdots +s_n}{G_{11}+\cdots +G_{nn}}t\right)1.$$

$$Q^{\ast}({\bm{s}})=\big{\{}\bm{q_{0}}+{\bm{q}}~{}\big{|}~{}{\bm{q}}\in{\rm Ker}({\sf R})\big{\}},$$

Full text

1 Introduction

In this paper we study a problem of mathematical modeling for dynamics occurring in steam supply systems that consist of distributed plants producing, consuming, and interchanging steam. This type of steam supply is crucial to realization of energy systems integration [1, 2], where multiple types of energy such as electricity, heat, and natural gas are managed consistently in order to satisfy specifications of stability, reliability, and energy efficiency. In this integration, the dynamics of energy transfer and conversion governed by different physical laws occur on a wide range of spatio-temporal scales [1, 2]. Modeling such multi-scale dynamics is of basic importance for establishing the control principle of the integrated energy systems. The problem which we study in this paper originates from the interaction between electricity and heat supply systems. This paper is a substantially-enhanced version of our conference papers [3, 4]. 111In this paper, we newly present the unified description of multiscale dynamics of steam supply systems by proper scaling of governing equations in Sec. 3, derivation of inner-limit of the derived model in Sec. 4, and numerical simulations for phase-space analysis and for verification of the model in Sec. 5, all of which are not reported in [3, 4].

The interaction between electricity and heat appears in the spatial deployment of Combined Heat and Power (CHP) [5] plants. The plants enable a local but spatially-distributed coupling between different energy systems because a CHP plant utilizes waste heat as a by-product of the conversion of fuel (natural gas, hydrogen, etc.) into electricity. It is stated in [1, 2] that this novel coupling would make it possible to design a new system architecture satisfying the above specifications. The overview of the target systems of this paper is shown in Fig. 1 with an example of three-site system. In this figure, multiple boilers (including heat recovery boiler within CHP) are connected via a steam pipe network, and the produced heat is exchanged between different sites. This type of the steam supply systems has appeared in practice, e.g. district heating systems [6]. In such conventional steam-based systems, the primary objective of the CHP operation is to supply a desirable fixed amount of steam, and hence the operation does not necessarily contribute to the electricity supply [6]. However, in the view point of energy systems integration, it is possible to consider a novel operation of the CHP plants that contributes to both the steam and electricity supply by utilizing their ability of rapid electric response. In fact, it is proposed in [5, 7, 8] that a CHP plant is used for compensating a variable output of renewable energy resources. This imposes a new problem on mathematical modeling for dynamics of steam supply systems against a large change of operating condition.

The contributions of this paper are twofold. The first one is to derive a lumped-parameter model that captures stability and multiscale properties of steam supply. For internal dynamics in a single boiler or single plant, much work has been reported on lumped-parameter modeling [9, 10, 11, 12], whereas for dynamics of steam flows in pipes, partial differential equations or distributed-parameter models [13, 14, 15, 16] are normally used. Although such models are crucial to plant design, detailed simulation, and commissioning, their simple coupling is too complicated to reveal the system-wide dynamics of interest. They are originally described by the model for fast steam flows over a pipe network coupled with the model for slow internal dynamics of boilers. Namely, the lumped-parameter model of boilers is regarded as a slowly time-varying boundary condition of the distributed-parameter model of steam pipes. In this paper, through physically-relevant approximations, we newly derive a lumped-parameter model that contains multiscale dynamic characteristics of the steam pipes and boilers as well as a graph-theoretic property of the pipe network.

The second contribution is to provide theoretical and numerical analyses of the derived model in terms of multiscale property of steam supply. The theoretical analysis is conducted with dynamical systems and graph theoretic methods [17, 18, 19, 20, 21]. We obtain the inner limit of the derived model using the standard regular expansion method [17] and locate a set of non-isolated equilibrium points of the inner-limit model. The set is thus proved to form a Normally Hyperbolic Invariant Manifold (NHIM) [18, 19] under mild technical conditions. The normal hyperbolicity characterizes the slow-fast vectorfield near the set, which is a dynamical evidence that the derived model preserves the slow-fast dynamics in the original model. Also, we conduct numerical simulations of the dynamics for an example of two-site system under a practical set of parameters. The correctness of the derived model is quantitatively verified by comparison with brute-force simulation of the original model, and the slow-fast dynamics near the NHIM are visualized. The existence of NHIM suggests a separation principle of short-term and long-term operations of steam supply systems, which is analogue to electricity supply operation [22].

The rest of this paper is organized as follows. In Sec. 2 we review the basic physical processes in steam supply systems. In Sec. 3, based on physical assumptions, we derive the lumped-parameter model. In Sec. 4 we theoretically analyze the derived model and prove the existence of NHIM. In Sec. 5 we perform numerical simulations of the two-site system for providing a technological implication and for verifying the derived model. Sec. 6 concludes this paper with a summary and future work.

2 Physical processes in steam supply systems

This section describes the basic physical processes of steam supply systems that are related to (i) producing, (ii) transporting, and (iii) consuming steam. Fig. 2 illustrates the three processes for the two-site system. The blocks in the figure represent (a) components of the two-site system and (b) components of a boiler. Each arrow describes the flow of steam, water, or heating gas. For (i), Fig. 2 (b) shows the production of steam from combustion gas and water in a boiler. Mixed steam and water are produced in the evaporator due to the boiling of water by high-temperature combustion gas. The mixed steam and water are supplied to the drum and stored at the phase equilibrium condition, while feedwater is supplied to the drum. Thus, the saturated steam in the drum is brought to the outside of boiler. For (ii), the transport of steam to a load in Fig. 2 (a) is realized due to the self-pressure of boiler. The transport of steam between the multiple sites is also realized by controlling the difference of pressures between the two sites connected to a pipe. Finally, for (iii), the transported steam is consumed at each load. The latent heat of the steam is extracted with a heat exchanger, and the resulting condensed water is returned to the boiler.

3 Derivation of lumped-parameter model

This section is devoted to the derivation of lumped-parameter model for dynamics of steam supply systems with the underlying physical processes in Sec. 2. The modeling procedure is exaplined in a dimensionless form in order to clearly describe the multiscale property of steam supply. The physical quantities with dimension are denoted by superscript $\ast$, and the reference quantities for scaling by subscript ${\rm r}$. The dimensionless time $t$ is scaled according to the time scale defined by the steam velocity $u_{\rm r}^{\ast}$ and the length scale $L_{\rm r}^{\ast}$ of the pipes: $t_{\rm r}^{\ast}:=L_{\rm r}^{\ast}/u_{\rm r}^{\ast}$. The detailed procedure of scaling and the reference quantities are presented in Appendix A.

3.1 Network description

In this and next sections, we will consider a general steam supply system with arbitrary number of sites based on the graph theory [20, 21]. The topology of a steam supply system is described by a directed graph $\mathcal{G}=(\mathcal{V},\,\mathcal{L})$, where $\mathcal{V}$ stands for a finite set of vertices representing sites, and $\mathcal{L}$ for a finite set of links representing steam transporting pipes. For a link $l\in\mathcal{L}$, the tail (or head) vertex is denoted by $\partial^{+}l$ (or $\partial^{-}l$). For a vertex $v\in\mathcal{V}$, the set of outgoing (or incoming) links is denoted by $\delta^{+}v$ (or $\delta^{-}v$). Below, the physical variables of boilers and pipes at each vertex and link are denoted by subscripts $v$ and $l$, respectively. The graph $\mathcal{G}$ is assumed to contain no self-loop and to be connected. The assumption of no self-loop is relevant because a steam transporting pipe normally connects different two sites. Under this assumption, the incidence relation of a graph is completely represented with the incident matrix [20], and hence the matrix will be used in our modeling and analysis. The assumption of connected graph is intended for simplifying the presentation in this paper and does not lose generality of the modeling and analysis here.

3.2 Physical assumptions

As described in Sec. 2, the pressure and flow rate of steam are important physical quantities. In order to simply describe essential characteristics of the system-wide dynamics, we make the following assumptions:

1. (A1)

Temperatures of drum, evaporator and their wall are equal to temperature of mixed steam and water, i.e. the saturation temperature [9, 10]. 2. (A2)

Change of the volumes of steam and water is negligible when the water level of a drum is well regulated [9]. 3. (A3)

Feedwater to a drum is at the condition of saturated liquid [10]. 4. (A4)

Pressure drop in a pipe is evaluated by the Darcy-Weisbach equation for steady flow [13]. 5. (A5)

No dominant effect of compressibility of steam appears on its velocity profile. This is relevant when the steam velocity is sufficiently smaller than the sound speed [23]. 6. (A6)

No dominant effect of heat loss appears on the pressure drop and volumetric flow. This is relevant when the mass fraction of vapor in the fluid, i.e. the quality of steam is sufficiently close to one [24, 25]. 7. (A7)

Pressure drop in the site’s components such as pressure regulators and valves are negligible.

The validity of (A1) to (A3) and (A4) have been tested in [9, 10] and [13, 14], respectively. Thus, we mainly discuss the assumptions from (A5) to (A7) in the rest of this paper.

3.3 Steam boiler

The dynamical model of a boiler is based on [9, 10]. In the model, $V$ represents volume, $\rho$ density, $h$ specific enthalpy, $T$ temperature, and $m^{\prime}$ mass flow rate. Furthermore, the three subscripts ${\rm s}$, ${\rm w}$, and ${\rm m}$ represent saturated steam, saturated water, and metal, respectively. The total mass of metals of the drum and the evaporator is represented by $m_{\rm t}$, and the specific heat of the metals by $C_{\rm p}$. It is stated in [9] that the dynamics of pressure are well captured by global mass and energy balance. This is because the internal energy is rapidly released or absorbed due to the uniform boiling and condensation inside the drum and evaporator. Thus, under the assumptions from (A1) to (A3), the dynamics of pressure $p_{v}$ at vertex $v\in\mathcal{V}$ are formulated as

[TABLE]

where $Q^{\prime}_{v}$ stands for the heat flow rate to the evaporator, and $h_{\rm c}:=h_{\rm s}-h_{\rm w}$ corresponds to the enthalpy of condensation. The small parameter $\epsilon_{1}:={d_{\rm r}^{\ast}}^{2}L_{\rm r}^{\ast}/e_{\rm r}^{\ast}$ describes the slowness of the pressure dynamics (1) in terms of the time scale $t_{\rm r}^{\ast}$ of steam flow. From [9], the coefficient $e_{v}(p_{v})$ represents the rate of change of internal energy stored in the boiler against a change of pressure, given by

[TABLE]

In this paper, according to [9], the thermodynamic properties $h_{\rm s}$, $h_{\rm w}$, $\rho_{\rm s}$, $\rho_{\rm w}$, and $T_{\rm s}$ are evaluated from the steam table [26] and are represented as functions of pressure $p_{v}$, for example, $h_{{\rm s}v}=h_{\rm s}(p_{v})$.

3.4 Steam pipe

The transient steam flow in a pipe is described by the one-dimensional continuity equations of mass, momentum, and energy [13, 14, 15]. For each link $l\in\mathcal{L}$, the mass balance is given by

[TABLE]

the momentum balance by

[TABLE]

and the energy balance by

[TABLE]

where $u_{l}$ stands for the velocity of steam in a pipe $l$, and $x$ for the displacement variable along the pipe. The parameter $\lambda_{l}$ stands for the friction coefficient of the Darcy-Weisbach equation [13, 25, 26] under (A4). The parameters $d_{l}$ and $Q_{{\rm w}l}$ stand for the diameter of the pipe and the heat flow through walls, respectively. The parameters $\epsilon_{2}:=\rho_{\rm sr}^{\ast}{u_{\rm r}^{\ast}}^{2}/p_{\rm r}^{\ast}$ and $\epsilon_{3}:={d_{\rm r}^{\ast}}^{2}{L_{\rm r}^{\ast}}Q_{\rm wr}^{\ast}/Q_{\rm r}^{\prime\ast}$ in (4) and (5) are small, and they reflect the assumptions (A5) and (A6), respectively. Under the two assumptions, the above original equations are simplified through the incompressibility condition $\partial u/\partial x=0$. It is widely accepted that the low Mach number (described by $\epsilon_{2}$) implies an incompressible model [23]. Further precise discussions are presented in [27, 28] for Navier-Stokes equations and in [29] for one-dimensional flow equations. The literature [27, 28, 29] shows that the simplification is relevant if both the constants $\epsilon_{2}$ and $\epsilon_{3}$ are sufficiently small. Hence, the energy equation (5) is decoupled from the other equations, and the dynamics of steam flows are described by the momentum equation (4) with the condition $\partial u/\partial x=0$:

[TABLE]

Therefore, by integrating (6) with respect to $x$ from $x=0$ to $x=L_{l}$ (see Fig. 3), the following ordinary differential equation is derived:

[TABLE]

where $\rho_{{\rm av}l}$ is given by

[TABLE]

3.5 Site

At each site, equations (1) and (7) are combined via the continuity equations of mass and energy. Under the assumption (A7) we have the following equations [13]: for each $v\in\mathcal{V}$ and all $t\in\mathbb{R}$, the mass balance is given by

[TABLE]

and the energy balance by

[TABLE]

where $m^{\prime}_{{\rm L}v}$ and $h_{{\rm L}v}$ stand for the mass flow rate and specific enthalpy consumed by the load at vertex $v$. They are related to the consumption rate $Q^{\prime}_{{\rm L}v}$ of heat as follows:

[TABLE]

By multiplying both the sides of (9) by $h_{\rm w}(p_{v})$ and using (10), we obtain

[TABLE]

The equation of $m^{\prime}_{{\rm s}v}h_{\rm c}(p_{v})$ is used in (1) and thereby combines the two equations (1) and (7).

3.6 Lumped-parameter model

Finally, in order to derive a model in a self-consistent manner, it is necessary to determine the thermodynamic quantities in the pipes. While they are calculated in the original model by the equations (3) for mass and (5) for energy, we give their values by functions of pressure $p_{v}$ based on the assumptions (A6) and (A7):

[TABLE]

The relevance of these approximations will be discussed in Sec. 5.3.

Consequently, the following model is derived for representing the steam supply dynamics: for each $v\in\mathcal{V}$ and $l\in\mathcal{L}$,

[TABLE]

where $\epsilon_{1}$ and $\epsilon_{2}$ were reset as $\epsilon$ by choosing the reference quantities as $e_{\rm r}^{\ast}:={d_{\rm r}^{\ast}}^{2}L_{\rm r}^{\ast}/\epsilon_{2}$. This resetting operation is relevant for a practical setting of parameters shown in Sec. 5.

Here, we discuss the multiscale property of the derived model: see Sec. 4 for its detailed analysis. The model (14) includes a single small parameter $\epsilon$, and the parameters and functions on the right-hand sides are order of $1$. From (14a), the pressure $p_{v}$ changes slowly in time due to the presence of the small parameter. In (14b), the first and second terms on the right-hand side should be order of $1$. That is, the pressure difference $p_{\partial^{+}l}-p_{\partial^{-}l}$ should be kept small (order of $\epsilon$). If this is not the case, for example if the pressure difference is $O(1)$, then the left-hand side of (14b) becomes $O(1/\epsilon)$, implying the steam velocity $u_{l}$ becomes large. This indicates the violation of (A5), and hence the derivation of the lumped-parameter model (14) loses its validity. Thus, the smallness of $p_{\partial^{+}l}-p_{\partial^{-}l}$ is necessary and will be assumed in the rest of this paper. Consequently, the short-term dynamics in $t\in[0,T]$ for $T=O(1)$ are described by the changes of $u_{l}$ in $O(1)$ and $p_{v}$ in $O(\epsilon)$. On the other hand, the long-term dynamics are related to the change of $p_{v}$ of $O(1)$, implying the changes of thermodynamic quantities such as $\rho_{\rm s}$.

4 Theoretical analysis

This section shows that the derived model (14) preserves the slow-fast dynamics in the original model. In phase-space geometric concepts, the presence of slow-fast dynamics can be described by the notion of normal hyperbolicity of invariant manifolds [18, 19]: slow dynamics along an invariant manifold and fast dynamics transversal to it. In order to theoretically analyze (14), here we simplify the model through the standard regular expansion method [17]. As will be shown later, it corresponds to the inner-limit of the derived model. Thus, we prove the existence of NHIM for the inner-limit model (19). This indicates that the slow-fast dynamics are involved in the model (14), because a NHIM persists under a perturbation of vectorfield [18, 19]. Indeed, in Sec. 5, we will numerically confirm the existence of NHIM for the model (14).

4.1 Notation

The notation frequently used in the rest of this paper is summarized below. The symbol $\mathsf{T}$ stands for the transpose operation of a vector or matrix. For a matrix $\mathsf{A}$, $\mathrm{Im}(\mathsf{A})$ represents the image space of linear mapping represented by $\mathsf{A}$, and $\mathrm{Ker}(\mathsf{A})$ the kernel space of the linear mapping $\mathsf{A}$. The symbol $\mathrm{diag}(\bm{v})$ stands for the diagonal matrix made from a vector $\bm{v}$. For a vector-valued function ${\bm{f}}=(f_{1},\ldots,f_{m})^{\top}$, its Jacobian is denoted by $D{\bm{f}}$. The constant vector ${\bm{1}}$ stands for all-one vector. For an Euclidean space $E$, $E^{\bot}$ stands for the orthogonal space of $E$.

4.2 Derivation of inner-limit model

First, we apply the regular expansion method to the derived model (14). The regular expansion starts with assuming a solution of (14) in the following form:

[TABLE]

By substituting them into (14) and equating the coefficient of each power of $\epsilon$, a series of differential equations is obtained. From the leading-order terms, we obtain the following conditions: for all $v\in\mathcal{V}$ and $l\in\mathcal{L}$,

[TABLE]

Thus, $p_{v}^{(0)}$ does not depend on both time and site, and is henceforth denoted by $p_{0}$. For the next-order of the series, we obtain the governing equations that contain the graph-theoretic property of the target steam supply system. For the graph $\mathcal{G}=(\mathcal{V},\,\mathcal{L})$, by labeling the vertices and links as $\mathcal{V}:=\{v_{1},\,\dots,\,v_{n}\}$ and $\mathcal{L}:=\{l_{1},\,\dots,\,l_{m}\}$, the first-order term of pressure and the zeroth-order term of volumetric flow of steam are described by

[TABLE]

The incidence matrix ${\sf R}=(R_{ij})\in\mathbb{R}^{n\times m}$ of the graph is defined by

[TABLE]

With the notation, the dynamics of $(\bm{\psi},\bm{q})$ are described as follows:

[TABLE]

with

[TABLE]

Note that (19) possibly has a unbounded solution, and thus the expansion (4.2) is not uniformly valid for all time $t\in\mathbb{R}$. Namely, there exists a finite $T$ such that the expansion is valid at all $t\in[0,T]$, and thus the limitation $\epsilon\to 0$ implies the inner limit [17]. In this sense, we call (19) as the inner limit model of (14).

4.3 Characterization of invariant manifold

In this subsection, we locate the invariant manifold in the inner-limit model (19) using graph theory: see Appendix B for its summary. As will be shown below, the manifold is located as a set of non-isolated equilibrium points. An equilibrium point $({\bm{\psi}}^{\ast},\,{\bm{q}}^{\ast})$ of (19) satisfies the following condition: for given $\bm{s}$, $\sf R$, and $\bm{f}$,

[TABLE]

The above conditions can be regarded as the combination of two linear equations defined by ${\sf R}$ and ${\sf R}^{\top}$. Thus, they are analyzed using the image and kernel spaces of ${\sf R}$ and ${\sf R}^{\top}$. First, (24a) has a solution if $\bm{s}$ belongs to the image of ${\sf R}$:

[TABLE]

Since (25) implies that $\bm{s}$ is orthogonal to ${\bm{1}}$ (see (54) in Appendix B), the sum of all the elements $s_{1}+\cdots+s_{n}$ should be zero. By the definition in (22), this condition is equivalent to

[TABLE]

This clearly indicates that the sum of all generation and consumption of steam at vertices is zero. Even if the above condition does not hold, the following analysis in this section is still applicable by introducing a new state variable $(\bm{\psi}^{\prime},\bm{q})$ with the following time-varying transformation:

[TABLE]

Thus, under the condition (25), the set of all solutions of (24a), denoted by $Q^{\ast}({\bm{s}})$, is described as follows:

[TABLE]

where \bm{q_{0}}standsforoneofthesolutionsof(\ref{eq:AEep-a}).Second,(\ref{eq:AEep-b})hasasolutionif$**f**(**q**)$belongstotheimageof$R^⊤$:\begin{equation}{\bm{f}}({\bm{q}})\in{\rm Im}({\sf R}^{\top}).\end{equation}Thisimpliesthat$**f**(**q**)$belongstothesubspaceofcutsets,whichisknownastheKirchhoff^{\prime}slawfortensions\cite[cite]{[@@bibref{}{iri69}{}{}]}.Underthecondition(\ref{eq:y-condition})andfor$**q**^∗∈Q^∗(**s**)$,theset$Ψ^∗(**q**^∗)$ofallsolutionsof(\ref{eq:AEep-b})isdescribedasfollows:\begin{equation}\Psi^{\ast}({\bm{q}^{\ast}})=\big{{}\bm{\psi_{0}}+{\bm{\psi}}{}\big{|}{}{\bm{\psi}}\in{\rm Ker({\sf R}^{\top})}=c,{\bm{1}}{},c\in\mathbb{R}\big{}}.\end{equation}where$**ψ_0$standsforoneofthesolutionsof(\ref{eq:AEep-b}).Consequently,underthecondition\begin{equation}{\bm{s}}\in({\rm Ker}({\sf R}^{\top}))^{\bot},\quad{\rm Im}({\sf R}^{\top})\cap{\bm{f}}(Q^{\ast}({\bm{s}}))\neq\emptyset,\end{equation}thesetofequilibriumpoints$**(**ψ**^∗, **q**^∗)$isdescribedasfollows:\begin{aligned} \Big{{}({\bm{\psi}}^{\ast},{\bm{q}}^{\ast})\in\mathbb{R}^{n+m}\Big{|}{}&{\bm{\psi}}^{\ast}={\bm{\psi}}_{0}+c,{\bm{1}},{}c\in\mathbb{R},{}{\bm{q}}^{\ast}\in{\bm{f}}^{-1}({\rm Im({\sf R}^{\top}))}\cap Q^{\ast}({\bm{s}})\Big{}},\end{aligned}where$**f**^-1 (Im(R^⊤)∩**f**(Q^∗(**s**)) = **f**^-1 (Im(R^⊤) ) ∩Q^∗(**s**)$holdsbecause$**f$isbijectivefromthedefinitionin\eqref{eq:definition_f}.Thisfactalsoindicatestheuniquenessof$q**^∗$ifitexists\cite[cite]{[@@bibref{}{iri69}{}{}]}.Asaresult,theset\eqref{eq:equilibrium_set}ofequilibriumpointsbecomesaone-dimensionalcurve(line)inthestatespaceoftheinner-limitmodel\eqref{eq:graph-model}.Thiscurveformsaninvariantmanifoldconsistingofallthenon-isolatedequilibriumpoints,whichwewilldenoteby$I$.\par InordertodiscussthetechnologicalimplicationinSec.,\ref{sec:numerical_example},thedynamicsofthemodel\eqref{eq:network_model}near$I$aredescribedbelow.Sincetheinvariantmanifoldhereinisone-dimensional,theslowdynamicswillbecharacterizedbyone-dimensionalreducedsystemonit.Aperturbationtosolutionsof\eqref{eq:graph-model}along$I$representstheuniformchangeofthepressures,andthustheslowdynamicscorrespondtotheresponsesofpressureleveloftheentiresystem.Ontheotherhand,thedynamicstransverseto$I$representthechangeofsteamvelocitiesandpressurefluctuationsnearthemanifold,andthefastdynamicscorrespondtothetransportofsteambetweenthedifferentsites.\par\par$

4.4 Proof of normal hyperbolicity of the invariant manifold

In this subsection, we prove that the located invariant manifold $\mathcal{I}$ in (4.3) is normally hyperbolic under certain conditions. An invariant manifold is called normally hyperbolic if the expansion or contraction rate of vectors transverse to the manifold dominates that of vectors tangent to the manifold. For a precise formulation of normal hyperbolicity, see [18, 30]. For the present discussion, since $\mathcal{I}$ consists of non-isolated equilibrium points, it is characterized by eigenvalues associated with the linearization of the model (19) at each equilibrium point. By substituting ${\bm{\psi}}={\bm{\psi}}^{\ast}+\bm{\Delta\psi}$ and ${\bm{q}}={\bm{q}}^{\ast}+\bm{\Delta q}$ into (19), the linearized system around $({\bm{\psi}}^{\ast},{\bm{q}}^{\ast})$ is obtained as follows:

[TABLE]

with

[TABLE]

In below, we will show that the center subspace of the linearized system (33) is one-dimensional and is tangent to $\mathcal{I}$. To do this, we analyze the eigenvector associated with zero eigenvalue under the following two assumptions. The first one is that the matrix $D{\bm{f}}(\bm{q}^{\ast})$ is non-singular. This is relevant if $q_{i}^{\ast}\neq 0$ for all $i=1,\dots,m$. Under the assumption, the eigenvector $(\bm{\Delta\psi_{0}},\bm{\Delta q_{0}})$ associated with zero eigenvalue satisfies the following equations:

[TABLE]

The second assumption states non-existence of pure imaginary eigenvalues of ${\sf A}$. As stated in [31], non-singularity of $D{\bm{f}}(\bm{q}^{\ast})$ is closely related to the non-oscillating condition, that is, non-existence of pure imaginary eigenvalues. However, in this paper, we will simply make both the assumptions. Under the two assumptions, one can verify that the center subspace is spanned by eigenvector associated to zero eigenvalue and is tangent to $\mathcal{I}$. By substituting (34b) into (34a), the following condition holds for $\bm{\Delta\psi_{0}}$:

[TABLE]

Since ${\sf\Sigma}$ is diagonal by definition, ${\sf R}{\sf\Sigma}{\sf R}^{\top}$ corresponds to the so-called Kirchhoff matrix [21]. From the assumption of connected graph, the kernel of ${\sf R}{\sf\Sigma}{\sf R}^{\top}$ is represented by

[TABLE]

From (34b), $\bm{\Delta q_{0}}={\bm{0}}$ holds because this space is also the kernel of ${\sf R}^{\top}$, that is, ${\sf R}^{\top}\bm{\Delta\psi_{0}}={\bm{0}}$. Note that, by discussion similar to above, one can verify that there is no generalized eigenvector associated with zero eigenvalue other than $(\bm{\Delta\psi_{0}},\bm{\Delta q_{0}})$. Thus, the center subspace of the linearized system (33) is explicitly represented as follows:

[TABLE]

This clearly indicates that the center subspace is one-dimensional and is tangent to the invariant manifold $\mathcal{I}$. From the definition, if the two assumptions—non-singularity of $D{\bm{f}}(\bm{q}^{\ast})$ and non-existence of pure imaginary eigenvalues of ${\sf A}$—are satisfied, then $\mathcal{I}$ is normally hyperbolic. The proof is thus completed.

5 Numerical simulations for two-site system

This section demonstrates the slow-fast dynamics near the NHIM and verifies the correctness of the derived model (14) by numerical simulations for the two-site system. For the system, in [32], we studied dynamics of electricity supply based on the classical formulation of power system swing equations [22] by assuming that the steam supply system is ideally operated. In this paper, following [32], we discuss the technological implication for the operation of steam supply systems by comparing to the electricity supply operation [22]. For the minimal two-site system, the dimension of the derived model is three, and hence it is possible to perfectly visualize the state space of the model including a NHIM and to apply phase-space geometric concepts [18] to it.

In [9] transient responses of a single boiler are examined experimentally as well as numerically under a setting of fuel profiles for a real plant. Following this, in this paper we provide responses of physical quantities under a step change and periodic change of the heat flow rates $Q^{\prime}_{i}$ to boilers. Although the abrupt change of $Q^{\prime}_{i}$ may not be possible in realistic operation, it provides basic information on the multiscale dynamics. The values of parameters used for the current simulations are shown in Tab. 1 and based on district heating systems [12, 25]. According to [11], the thermodynamic properties are calculated by using Xsteam package [33]. Numerical values of $e_{i}(p_{i})$ and $h_{\rm c}(p_{i})\rho_{\rm s}(p_{i})$ in (14) are shown in Fig. 4 for a practical range of pressure [12]: $0.03\,{\rm MPa}\leq p_{i}\leq 2\,{\rm MPa}$.

5.1 Time-response analysis

First, we simulate the short-term dynamics related to transport of steam between the two sites. Time-responses of physical quantities are provided under the following step change:

[TABLE]

This implies an abrupt change of operating conditions of the boilers, and thereby $1\,{\rm MJ/s}$ surplus (or deficit) of heat is caused in site $\#1$ (or site $\#2$). Note that in both cases the condition (26) is satisfied. Fig. 5 shows step responses of the state variables $(p_{1},\,p_{2},\,u)$ of (14) and heat output rates $Q^{\prime}_{{\rm o}i}$ from the boilers. The system is initially at a steady operating condition with no transport of steam between the sites, and then $Q^{\prime}_{i}$ changes at $t=10\,{\rm s}$. The pressures $p_{1}$ and $p_{2}$ and velocity $u$ move to a new operating condition after transients in a few tens of seconds. At the new operating condition, the values of $p_{1}-p_{2}$ and $u$ become positive. This clearly shows that the positive pressure drop $p_{1}-p_{2}$ induces the transport of steam from site $\#1$ to site $\#2$. The heat output rates $Q^{\prime}_{\rm o1}$ and $Q^{\prime}_{\rm o2}$ change symmetrically, and thus the surplus and deficit of heat at the two sites are compensated by the transport of steam.

5.2 Phase-space analysis

Second, we analyze the dynamics described by the model (14) from the viewpoint of phase space. Especially, the long-term dynamics described by (14) are considered in this subsection. The model (14) for the two-site system has the three independent variables. Based on the analysis in Sec. 4.3, we introduce the following variable transformation:

[TABLE]

Fig. 6 shows trajectories of the model (14) under the parameter setting as $(Q^{\prime}_{1},\,Q^{\prime}_{2})=(6\,{\rm MJ/s},4\,{\rm MJ/s})$. The red trajectory corresponds to the time response presented in Fig. 5 with subsequent long-term response from $t=0\,{\rm s}$ to $10,000\,{\rm s}$. This trajectory shows typical slow-fast dynamics as mentioned in Sec. 4. The mean pressure value $(p_{1}+p_{2})/2$ does not change dominantly while the pressure difference $p_{1}-p_{2}$ and the velocity $u$ exhibit fast oscillations shown in Fig. 5. After the fast oscillation is settled, the mean pressure $(p_{1}+p_{2})/2$ begins to decrease slowly. The blue trajectory in the figure shows a one-dimensional invariant manifold located with direct numerical integration of (14). It is confirmed that the trajectories converge to the located invariant manifold while exhibiting the slow and fast dynamics as mentioned above. This result becomes a numerical evidence that the located invariant manifold possesses the normal hyperbolicity.

Also, we analyze the long-term dynamics related to the boilers’ operation to offer several technological implications of the phase space analysis. Here, we consider the following periodic change of $Q^{\prime}_{1}$: for $n=1,\,2,\dots$,

[TABLE]

and $Q^{\prime}_{2}=4\,{\rm MJ/s}$. This periodic change is intended to the novel electricity-oriented operation of CHP plants mentioned in Sec. 1, and the similar profiles of the fuel flow rate are shown in [9] as experiment data of a real plant. Under the above setting, the condition (26) is not satisfied. Fig. 7 shows (a) the time responses of the state variables and (b) the corresponding trajectory in the phase space. The appearing short-term and long-term dynamics are characterized by the NHIM. For the short-term regime, the fast motion towards to the NHIM guarantees that the pressure and velocity oscillation are settled after a change of operating condition. This ensures that the amount of transport of steam between the two sites becomes bounded in the short-term regime. On the other hand, for the long-term regime, the mean pressure value changes due to the slow motion along the NHIM. Thus, the success of the novel electricity-oriented operation can be clarified as the existence of NHIM near which the separation of fast and slow motions holds.

It should be here noted that this type of slow-fast separation plays an important role in power system operation [22]. Conventionally, the frequency dynamics in power systems can be classified into three stages with different time-scales and are regulated separately with different mechanisms [22]. In this sense, we now identify a similar time-scale separation in the steam supply system, which has not been reported yet in literature. The time-scale separation enables independent operations of transport of steam between the sites in short-term regime and long-term supply-demand balancing of steam in order to maintain the pressure level in the system. The finding of the time-scale separation is thus expected to become a dynamical principle for operational design of steam supply.

5.3 Comparison with the original model

Lastly, in order to verify the correctness of the derived model (14), we present a comparison with the brute-force simulation of the original model represented by equations (1) to (5) and (3.5). The equations were implemented using the COMSOL Multiphysics$\textregistered$ Software. According to [15], the thermodynamic properties were given by the first order approximation around nominal values. From the procedure of lumped-parameter modeling, the approximations needed for deriving (14) are the incompressibility condition and the evaluation of thermodynamic quantities in (3.6). Since the relevance of these approximations is related to the parameter $\epsilon_{3}:=d_{\rm r}^{2}L_{\rm r}Q_{\rm wr}/Q^{\prime}_{\rm r}$, we present in Fig. 8 the simulation results for various setting of the heat loss, represented by $Q_{\rm w}$ in (5). The lines in the figure show the results of brute-force simulation of the original model with $-Q_{\rm w}\cdot\pi d^{2}/4=0\,{\rm W/m}$, $100\,{\rm W/m}$, and $200\,{\rm W/m}$, respectively. The sequence of points, denoted by $\times$, represents a sample plot of time responses of (14) presented in Fig. 5. The simulation result clearly shows that the lumped-parameter model (14) well approximates the dynamics of the original model when the heat loss is sufficiently small, i.e. when the steam pipes are well insulated.

For the cases with large heat loss, while the pressures $p_{1}$ and $p_{2}$ in the original model slowly decrease, the responses of $u_{0}$ and $u_{L}$ are well captured by the derived model (14). Also, the responses of pressures and velocities in $[10\,{\rm s},20\,{\rm s}]$ are correctly produced by the derived model. The result shows that the derived model describes the transport of steam in the shot-term regime even if the heat loss becomes large. This model is thus useful for designing the coordinated operation of heat and electricity supply because the typical time scale of the swing dynamics in power systems is also in the same time regime [22].

6 Conclusions

In this paper, we studied a problem on mathematical modeling for the dynamics in steam supply systems. The dynamics of interest were originally described by a distributed-parameter model for fast steam flows over a pipe network coupled with a lumped-parameter model for slow internal dynamics of boilers. Through physically-relevant approximations, we newly derived a lumped-parameter model that captured stability and multiscale properties of the dynamics. In order to describe the slow-fast dynamics, we used the notion of Normally Hyperbolic Invariant Manifold (NHIM). By theoretically analyzing the inner limit of the derived model, we located a set of non-isolated equilibrium points that formed a NHIM. Also, the numerical simulations under practical settings of parameters demonstrated the slow-fast dynamics near the NHIM. The existence of NHIM clearly suggests the so-called separation principle for operational design of steam supply, which is analogue to power system operation.

Future directions of this work are as follows. One is to verify the correctness and application limit of the derived model via experimental measurements. This is inevitable for practical use of the derived model. Another one is operational design of multiscale dynamics of steam supply based on the NHIM. The characterization is expected to be utilized for separation of different scales or model-order reduction.

Appendix A Derivation of dimensionless governing equations

This first appendix derives the dimensionless equations from (1) to (5). The physical quantities with dimension are denoted by superscript $\ast$, and the reference quantities by subscript ${\rm r}$: A physical quantity $z^{\ast}$ is associated by $z^{\ast}=z_{\rm r}^{\ast}z$ to a dimensionless quantity $z$. Similarly, a function $f^{\ast}$ is related by $f_{\rm r}^{\ast}f(z^{\ast}/z_{\rm r}^{\ast})=f^{\ast}(z^{\ast})$ to a dimensionless function $f$. The reference quantities are related by the following equations:

[TABLE]

A.1 Boiler model

It is well-known in [9, 10] that the dynamic behavior of boiler’s pressure is well captured by global mass and energy balances. The global mass balance is given by

[TABLE]

and the global mass balance by

[TABLE]

where the term $h^{\ast}-p^{\ast}/\rho^{\ast}$ corresponds to internal energy. Under (A2), by multiplying (42) by $h_{\rm w}^{\ast}$ and subtracting the result from (43) we have

[TABLE]

with

[TABLE]

In addition, (A1) implies $T_{\rm m}^{\ast}=T_{\rm s}^{\ast}$, and (A3) does $h_{\rm f}^{\ast}=h_{\rm w}^{\ast}$ in (44) and (45). As a result, by using the relation (41), the pressure dynamics of boiler are formulated as follows:

[TABLE]

Thus, (1) is derived by defining the small parameter $\epsilon_{1}:={d_{\rm r}^{\ast}}^{2}L_{\rm r}^{\ast}/e_{\rm r}^{\ast}$.

A.2 Steam pipe model

The continuity equations of mass, momentum, and energy with dimension are given as follows [13, 14, 15]:

[TABLE]

where $F_{\rm w}$ stands for the shear force acting on a steam element and is approximated by the Darcy-Weisbach equation [13, 26, 25] as follows:

[TABLE]

The coefficient $\lambda$ depends on the Reynolds number $Re$, pipe diameter $d$, and roughness of the inner surface of the pipe and therefore varies according to the steam velocity $u$. In the case of laminar flow under a low Reynolds number, the value of $\lambda$ is derived in [26] theoretically as $\lambda=64/Re$. On the other hand, in the case of turbulent flow under high Reynolds numbers, $\lambda$ is approximated by a constant that is determined by $d$ and the roughness [25, 26]. In this paper, the friction coefficient $\lambda$ is considered as a constant because the steam flow used in standard steam supply is turbulent [25, 26]. By using the relation (41), the above equations become

[TABLE]

Thus, (3) to (5) are derived by defining the parameters $\epsilon_{2}:=\rho_{\rm sr}^{\ast}{u_{\rm r}^{\ast}}^{2}/p_{\rm r}^{\ast}$ and $\epsilon_{3}:={d_{\rm r}^{\ast}}^{2}L_{\rm r}^{\ast}Q_{\rm wr}^{\ast}/Q_{\rm r}^{\prime\ast}$.

Appendix B Summarized Graph Theory

The second appendix provides a summarized theory of graph from [20]. Along the notation introduced in Sec. 4.1, consider a directed graph $\mathcal{G}$ with $n$ vertices and $m$ links. Assume $\mathcal{G}$ is connected and is represented by the incidence matrix ${\sf R}$. The matrix ${\sf R}$ (or ${\sf R}^{\top}$) is regarded as a linear map ${\sf R}:\mathbb{R}^{m}\to\mathbb{R}^{n}$ (or ${\sf R}^{\top}:\mathbb{R}^{n}\to\mathbb{R}^{m}$), and its image and kernel are related to the graph’s topology. ${\rm Ker}({\sf R})$ and ${\rm Im}({\sf R}^{\top})$ are subspaces of $\mathbb{R}^{m}$, and their dimensions coincide with the number of independent loops and cutsets, respectively. Since $\mathcal{G}$ is connected, we have ${\rm dim}({\rm Ker}({\sf R}))=m-n+1$ and ${\rm dim(Im}({\sf R}^{\top}))=n-1$ (see [20]), where ${\rm dim}(X)$ stands for the dimension of $X$. Also, ${\rm Ker}({\sf R}^{\top})$ is a subspace of $\mathbb{R}^{n}$ given as

[TABLE]

This result is derived from ${\rm dim(Im}({\sf R}^{\top}))=n-1$ and the fact that every link connects exactly two vertices. Finally, ${\rm Im}({\sf R)}$ is a subspace of $\mathbb{R}^{n}$ orthogonal to ${\rm Ker}({\sf R}^{\top})$, and its dimension is equal to $n-1$.

List of Tables

1. 1 List of variables and parameters in the derived model (14). The values used for numerical simulations in Sec. 5 are also presented.

List of Figures

1. 1 Overview of steam supply systems with multiple CHP plants that we consider in this paper. It shows an example of three-site system.
2. 2 Schematic diagram of the two-site steam supply. Each block in the figure shows (a) components of the two-site system and (b) components of a boiler. Red, orange, and black arrows describe the flow of steam, water, and heating gas, respectively.
3. 3 Schematic diagram of the steam transporting pipe $l$.
4. 4 Numerical values of (a) $e_{i}$ and (b) $h_{\rm c}\rho_{\rm s}$
5. 5 Responses of state variables $(p_{1},\,p_{2},\,u)$ and heat output rates $Q^{\prime}_{{\rm o}i}:=m^{\prime}_{{\rm s}i}h_{\rm c}(p_{i})$ from the boilers. The solid lines show the responses of (14), and the points ($\times$) show the sample plot of responses of (19). They are initiated by a step change of the parameters from $(Q^{\prime}_{1},\,Q^{\prime}_{2})=(5\,{\rm MJ/s},5\,{\rm MJ/s})$ to $(Q^{\prime}_{1},\,Q^{\prime}_{2})=(6\,{\rm MJ/s},4\,{\rm MJ/s})$.
6. 6 Trajectories of the derived model (14) with $(Q^{\prime}_{1},Q^{\prime}_{2})=(6\,{\rm MJ/s},4\,{\rm MJ/s})$. The red trajectory corresponds to the time response presented in Fig. 5 with subsequent long-term response. The blue trajectory shows a invariant manifold located with direct numerical integration of (14).
7. 7 Long-term dynamics for the periodic change of the parameter $Q^{\prime}_{1}$. Numerical simulation of (a) time responses of the state variables and (b) trajectory in the phase space are shown.
8. 8 Responses of $p_{1}$, $p_{2}$, $u_{0}:=u(t,0)$, and $u_{L}:=u(t,L)$ by the original coupled equations. The points ($\times$) show the sequence of time responses of the derived lumped-parameter model (14).