Structure‐preserving methods for a coupled port‐Hamiltonian system of compressible non‐isothermal fluid flow

The port‐Hamiltonian (pH) formulation of partial‐differential equations and their numerical treatment have been elaborately studied lately. One advantage of pH‐systems is that fundamental physical properties, like energy dissipation and mass conservation, are encoded in the system structure. Therefore, structure‐preservation is most important during all stages of approximation and system coupling. In this context we consider the non‐isothermal flow of a compressible fluid through a network of pipes. Based on a pH‐formulation of Euler‐type equations on one pipe, we introduce coupling conditions, through which we can realize energy, mass and entropy conservation at the coupling nodes and thus, preserve the pH‐structure. We implement them through an input‐output‐coupling using the flow and effort variables of the boundary port. Thus, we can make use of the structure‐preserving model and complexity reduction techniques for the single pipe. This procedure becomes even more important for network simulations, as here, we deal with high dimensional and highly non‐linear dynamical systems. We explain the extension from a single pipe to a network and numerical examples are shown to support our findings.

This system of equations describes fluid flow through a pipe of length  and diameter  with friction factor  and heat transmission coefficient   .The ambient temperature of the pipe is denoted by  ∞ .The unknowns are the mass density , the velocity  and the internal energy density .The temperature  and the pressure  have to be given by state equations.We make use of the following three assumptions.First, the fluid flow is governed by the ideal gas law, such that we have for the pressure , the temperature , the specific entropy  and the total specific enthalpy ℎ, respectively.Here,  is the gas constant,   and   are the heat capacities and the adiabatic index is given by  =   ∕  .Secondly, the considered flow is subsonic, that is, || <  = √ ∕, where  denotes the speed of sound.Therefore, two boundary conditions need to be stated at the inflow and one at the outflow of the pipe, see [1].Lastly, we only consider networks without a change in flow direction.The total energy of (1), that is, the Hamiltonian, is given by (, , ) = ∫    2

PORT-HAMILTONIAN FORMULATION ON A SINGLE PIPE
In [2] we introduced an infinite dimensional port-Hamiltonian (pH)-formulation for (1) motivated by the finite dimensional pH-system proposed in [3].To make the pH-formulation favorable for structure-preserving approximation, we use the variable transformation  = ∕ and choose the mass flow  to be a state variable instead of the velocity . 1]  and the input  =  ∞ .The system operators are, The operator   is defined for functions  1 ,  2 ,  3 by 3 is self-adjoint and semi-elliptic, as   , ,  ≥ 0 for all  ≥ 0, and ()[⋅] ∶  1 () 3 ↦  1 () 3 is formally skewadjoint in the  2 inner-product, when the boundary terms are assumed to be zero.For ,  ∈  1 () 3 the latter can be proved using partial integration and by showing that (, ()) = −(, ()), that is, Since the boundary terms created by partial integration are usually not zero in applications, this gives rise to including them in a boundary port with boundary flow and effort variables at the inflow and outflow of the pipe, that is, Using the ansatz introduced in [4] these variables can be deduced from the skew-adjoint system operator (), which yields that is, the boundary flows are the mass flow at the inflow and the negative of the mass flow at the outflow, whereas the boundary efforts are given by the total specific enthalpy at both pipe ends.Let us additionally define the following boundary operators.Let  ∈  1 () and let [0] and [] be the evaluations of  at the pipe ends.Then we introduce and The structure-preserving incorporation of boundary conditions is possible through the boundary port variables and the boundary operators.Although there is some freedom in choosing the boundary conditions, we demand that one of the conditions at the inflow is stated in terms of the internal energy density.Thus,  e only acts at the inflow of the pipe.This becomes important for our coupling.For space-discretization with finite elements the weak formulation of System 2.1 is given as follows.
for  ∈  2 (), ,  ∈  1 () and  ∈ [0,   ],  = (0, ) with the boundary efforts The variational formulation can be deduced as follows.Testing the equations from System 2.1 with  ∈  2 () and ,  ∈  1 (), integrating over the pipe  and using partial integration on certain terms yields, For the incorporation of the boundary conditions into the formulation, we first use the equations for the boundary flows  B to state boundary conditions for the mass flow .Recall (3), such that we have We define the associated Lagrangian multiplier as ] .Thus, the boundary term in (5a) can be reformulated to As we want to preserve the pH-structure, we need to make the equations stating the boundary conditions dependent on the respective effort variables, that is, ẽ2 () =  and ẽ3 () = 1.As seen above, this dependency comes somewhat natural for (6).In the case of the boundary condition for the internal energy density we have to work a little bit more.Using the assumption that  B [0] ≠ 0, we can state the boundary condition for the internal energy as In this way, we can express it through the boundary flow.This might seem like unnecessarily blowing the system up, but by doing this the boundary effort variables act as output in the pH-sense.Finally, we also have to add this equation to our system via a Lagrange multiplier.We define  e = [0] Note that, energy dissipation and mass conservation can be proved analogously to [5].

Network topology and coupling conditions
In the following we consider connected, finite and directed graphs  = ( , , ).A graph consists of a set of nodes  = { 1 , … ,   },  ∈ ℕ, and a set of pipes  = { 1 , … ,   } ⊂  ×  ,  ∈ ℕ.Every pipe  has a length   , 0 <   < ∞, which we collect in the set  = {  1 , … ,    }.The superscript  denotes the respective pipe .Furthermore, the structure of the graph  is given through an incidence mapping, which describes the direction of the edges.Further sets we work with are the set of edges adjacent to node , that is, () = { ∈  ∶  = (, ν) or  = ( ν, )}, the set of interior or coupling nodes  0 ⊂  , and the set of boundary nodes   =  ∖ 0 .Let  ∈  0 , we then have that is, the graph-topological directions of the pipes with respect to a node.From these sets and the flow direction, we can define,   () ∶= { ∈   () ∶   ≥ 0} ∪ { ∈   () ∶   ≤ 0},   () ∶= { ∈   () ∶   < 0} ∪ { ∈   () ∶   > 0}, which tell us the pipes  on which the flow travels into or out of the node , respectively.An example graph is given in Figure 1.Due to assuming that there is no change in flow direction, we are able to set   () =   () and   () =   () before the simulation.Otherwise, these sets would change during the simulation, which evidently has an impact on the computations and the implementation of the resulting algorithms.This is still a question of ongoing research.
As we need to state three conditions on each inner node, we choose the following three coupling conditions motivated by [6], (M) Conservation of mass: Since the system on each pipe is in pH-form we want to keep the pH-character also for the coupled network system.This is possible by using power-conserving input-output coupling.For the Euler equations ( 1), we already have that (M) and (H) imply energy conservation, that is, ∑ ∈()   [](  (  +   ))[] = 0, see [6].Therefore, the first two coupling conditions are sufficient for structure-preservation, whereas (S) is needed to close the system.

Structure-preserving incorporation of coupling conditions
The coupling conditions (M) and (H) are realizable using the boundary flow and effort variables, respectively.For mass conservation we have that For the enthalpy equality at each  ∈  0 we introduce a Lagrange multiplier   H , such that To realize (S) we use the ideal gas assumption and the equation for the specific entropy given in (2).Thus, we are able to rewrite the coupling conditions in terms of the internal energy density, that is, ) .
This yields a state equation for   [] at  ∈  0 which can be plugged into (4) for each pipe, that is, Therefore, the incorporation of the coupling conditions (M), (H) and (S) is done via an input-output-coupling and does not destroy the pH-structure.To write down the variational formulation on the network, we need to define the function spaces.Let  ∈  be identified with (0,   ), then the spatial domain of the network is given by Ω = { ∶  ∈ , for  ∈ }.Hence, we define  2 () = { ∶ Ω → ℝ with |  ∈  2 () for all  ∈ }.For ,  ∈  2 () the inner-product and the associated norm are given by (, )  = ∑ ∈ (|  , |  ) and ||||  = √ (, )  , respectively.We introduce the broken operator  ′  , which describes the edge-wise weak derivative, that is, ( ′  )|  =   |  for all  ∈ .Thus, The proof is analogous to Theorem 2.3 the power-conserving coupling.The structures of the variational formulations in 2.2 and 3.2 do not differ in general, except for the the function spaces and algebraic equations.Thus, the same structure-preserving numerical methods as introduced in the [7, 8] and [2] are applied to the concluding numerical example.

NUMERICAL EXAMPLE
Consider the network systems from Figure 2. Here, a coupling node is added in the middle of the single pipe, as we show that the coupling conditions do not change the simulation of the flow without or with model order reduction.We simulate the same flow on both networks.We work with two-atomic ideal gas, that is,   = 5∕2,  = 1.The system parameters are given by  = 40,   = 5∕4 and  ∞ = 1,  = 1.For space discretization  0 ,  1 ,  1 finite elements are used for the mass density, the mass flow and the internal energy density, respectively, with Δ = 0.01.As this leads to a pH differentialalgebraic equation we use an implicit Euler method with Δ = 0.1 for time integration.We use the same initial values and boundary conditions for both systems, that is,  0 = 3,  0 = 0.3,  0 = 9 and [ 1 ] = [ 3 ] = 0.3, [ 1 ] = 9.Consistent initial values for all Langrange multipliers have been computed.For the 1P system we choose   1 = 1 and   1 = 100 finite elements, which leads to a full order model (FOM) of dimension 305.For the 2P system we have   11 =   12 = 0.5 and   11 =   12 = 50, which yields a FOM of dimension 313, which is due to the algebraic coupling conditions.Applying model order reduction to the FOMs, yields the relative  2 -and projection errors shown in Figure 3. Reducing the 1P system to  = 20 finite elements yields a reduced order model (ROM) of dimension 65 and a relative error of (10 −5 ).For the 2P system there are two different cases for model reduction.Either we choose to reduce the network as a whole (net), that is,  = 20, or pipe-wise (pw), that is,   11 =   12 = 10, which yields in both cases a ROM of dimension 73, but different projection matrices.The net-approach yields an error in (10 −5 ) similar to the 1P scenario, whereas the pwapproach leads to an error of one order of magnitude higher, that is, (10 −4 ).Furthermore, the total energy and mass is the same for all full and ROMs, see Figures 4 and 5.The small errors in the mass conservation are due to rounding errors.Thus, adding the coupling node does not change the flow simulation as desired.Overall, we see that applying numerical methods designed for the single pipe system can be easily applied to network systems without bigger adjustments, as the coupling conditions where specifically chosen such that the energy and mass are conserved at the coupling nodes.

A C K N O W L E D G M E N T S
The support of the German BMBF, project ElAN, is acknowledged.Open access funding enabled and organized by Projekt DEAL.

F I G U R E 3 F I G U R E 4 F I G U R E 5
Relative  2 -and projection-errors.Total energy for all systems.Total mass for all systems.