A Scalable port-Hamiltonian Model for Incompressible Fluids in Irregular Geometries

3rd IFAC/IEEE CSS Workshop on Control of Systems Governed by 3rd IFAC/IEEE CSS Workshop on Control of Systems Governed by Partial 3rd 3rd IFAC/IEEE IFAC/IEEE CSS CSS Workshop Workshop on on Control Control of of Systems Systems Governed Governed by by Partial Available online at www.sciencedirect.com Differential Equation, and XI Workshop Control of Distributed Partial 3rd IFAC/IEEE CSS Workshop on Control of Systems Governed by Partial Differential Equation, and XI Workshop Control of Distributed Parameter DifferentialSystems Equation, and and XI XI Workshop Workshop Control Control of of Distributed Distributed Partial Differential Equation, Parameter Systems Oaxaca, Mexico, May 20-24, Parameter Systems Differential Equation, and XI 2019 Workshop Control of Distributed Parameter Systems Oaxaca, Mexico, May 20-24, 2019 Oaxaca, Mexico, Mexico, May 20-24, 20-24, 2019 2019 Parameter Systems Oaxaca, May Oaxaca, Mexico, May 20-24, 2019 IFAC PapersOnLine 52-2 (2019) 102–107

ScienceDirect

A A A A

Scalable port-Hamiltonian Model for Scalable Scalable port-Hamiltonian port-Hamiltonian Model Model for for Incompressible Fluids Scalable port-Hamiltonian Model for Incompressible Fluids in in Irregular Irregular Incompressible Fluids in Irregular Geometries Incompressible Fluids in Irregular Geometries Geometries Luis A. Mora ∗∗ Geometries H´ ector Ram´ırez ∗∗ Juan I. Yuz ∗∗

Luis A. Mora ∗∗ H´ ector Ram´ırez∗∗∗∗ Juan I. Yuz ∗∗ Luis H´ e Ram´ ırez Le Gorrec Luis A. A. Mora Mora Yann H´ ector ctor ırez∗∗ Juan Juan I. I. Yuz Yuz Le Ram´ Gorrec ∗∗ ∗Yann ∗∗∗ Juan I. Yuz ∗ Le Gorrec Luis A. Mora Yann H´ e ctor Ram´ ırez Yann Le Gorrec ∗∗ ∗ Yann Le Gorrec Engineering, Universidad T´eecnica ∗ Department of Electronic cnica Federico Federico ∗ Department of Electronic Engineering, Universidad T´ ∗ Department of Electronic Electronic Engineering, Universidad T´ee(e-mail: cnica Federico Federico Santa Mar´ ıa, Av. Espa˜ n a 1680 Valpara´ ıso, Chile Department of Engineering, Universidad T´ cnica Santa Mar´ ıa, Av. Espa˜ n a 1680 Valpara´ ıso, Chile (e-mail: ∗ Santa Mar´ ıa, Av. Av. Espa˜ Espa˜ na a 1680 1680 [email protected], Valpara´ ıso, Chile Chile (e-mail: Department of Electronic Engineering, Universidad T´e(e-mail: cnica Federico [email protected], Santa Mar´ ıa, n Valpara´ ıso, [email protected], [email protected], [email protected], Santa Mar´ıa, Av. Espa˜ na 1680 [email protected], Valpara´ıso, Chile (e-mail: [email protected]). [email protected], [email protected], [email protected]). ∗∗ [email protected]). [email protected], [email protected], Univ. Bourgogne Franche-Comt´ ee,, CNRS, ∗∗ FEMTO-ST, [email protected]). CNRS, 24 24 rue rue ∗∗ FEMTO-ST, Univ. Bourgogne Franche-Comt´ ∗∗ FEMTO-ST, Univ. Bourgogne Franche-Comt´ e,, CNRS, CNRS, 24 24 rue rue [email protected]). Savary, F-2500 Besan¸ c on, France (e-mail: [email protected]) FEMTO-ST, Univ. Bourgogne Franche-Comt´ e Savary, F-2500 Besan¸ c on, France (e-mail: [email protected]) ∗∗ Savary, F-2500 Besan¸ Besan¸ on, France (e-mail: (e-mail: [email protected]) FEMTO-ST, Univ. Bourgogne Franche-Comt´ e, CNRS, 24 rue Savary, F-2500 ccon, France [email protected]) Savary, F-2500 Besan¸con, France (e-mail: [email protected]) Abstract: Abstract: The The behavior behavior of of a a fluid fluid in in pipes pipes with with irregular irregular geometries geometries is is studied. studied. Departing Departing Abstract: The behavior of a fluid in pipes with irregular geometries is studied. from the partial differential equations that describe mass and momentum balances aa scalable Abstract: The behavior of a fluid in pipes with irregular geometries is studied. Departing Departing from the partial differential equations that describe mass and momentum balances scalable from the partial differential equations that describe mass and momentum balances aa scalable Abstract: The behavior of a fluid in pipes with irregular geometries is studied. Departing lumped-parameter model is proposed. To this end the framework of port-Hamiltonian systems from the partial differential equations that describe mass and momentum balances scalable lumped-parameter model is proposed. To this end the framework of port-Hamiltonian systems lumped-parameter modelais ismodular proposed. To this endupon themass framework of port-Hamiltonian port-Hamiltonian systems from the partial differential equations that describe and momentum balances a scalable is instrumental to derive system which interconnection describes segments with lumped-parameter model proposed. To this end the framework of systems is instrumental to derive a modular system which upon interconnection describes segments with is instrumental instrumental to derive derive aismodular modular system which upon interconnection describes segments segments with lumped-parameter model proposed. To this end the framework of port-Hamiltonian systems different cross sections and dissipation effects. In order to perform the interconnection between is to a system which upon interconnection describes with different cross sections and dissipation effects. In order to perform the interconnection between different cross sections and dissipation effects. In order to perform the interconnection between is instrumental to derive a modular system which upon interconnection describes segments with segments the incompressibility hypothesis is relaxed in some infinitesimal section to different cross sections and dissipation effects. In order to perform the interconnection between different segments the incompressibility hypothesis is relaxed in some infinitesimal section to segments the incompressibility hypothesis is relaxed in some infinitesimal section to different cross sections and dissipation effects. In order to perform the interconnection between admit density variations and energy transference between segments. Numerical simulations segments the incompressibility hypothesis is relaxed in some infinitesimal sectionare to admit density variations and energy transference between segments. Numerical simulations are admit density density variations and energy energy transference between segments. Numerical simulations are different segments the incompressibility hypothesis is relaxed in some infinitesimal section to performed in order to illustrate the model. admit variations and transference between segments. Numerical simulations are performed in order to illustrate the model. performed in order order to illustrate illustrate the model. model. admit density variations and energy transference between segments. Numerical simulations are performed in to the © 2019, IFAC Federation Automatic Control) Hosting by Elsevier Ltd. All rights reserved. performed in (International order to illustrate the of model. Keywords: Keywords: Port-Hamiltonian Port-Hamiltonian systems, systems, PDE, PDE, approximation approximation of of PDEs, PDEs, computational computational methods methods Keywords: Keywords: Port-Hamiltonian Port-Hamiltonian systems, systems, PDE, PDE, approximation approximation of of PDEs, PDEs, computational computational methods methods Keywords: Port-Hamiltonian systems, PDE, approximation of PDEs, computational methods 1. INTRODUCTION INTRODUCTION niques, niques, such such as as finite-volumes finite-volumes and and LS-STAG, LS-STAG, among among othoth1. 1. INTRODUCTION niques, such as finite-volumes and LS-STAG, among others, that presents strong computational demands to ob1. INTRODUCTION niques, as finite-volumes and LS-STAG, among ers, thatsuch presents strong computational demands to othobers, that that presents strong computational demands to2016; ob1. INTRODUCTION niques, such as finite-volumes and LS-STAG, among othtain a detailed flow description (Bourantas et al., ers, presents strong computational demands to obA fluid is considered incompressible when its density is tain a detailed flow description (Bourantas et al., 2016; A fluid is considered incompressible when its density is tain that a detailed detailed flow description (Bourantas et al., al.,to2016; 2016; ers, presents strong computational demands obSharatchandra and Rhode, 1994). A fluid is considered incompressible when its density is tain a flow description (Bourantas et constant. In practice, to simplify the analysis of complex A fluid is considered incompressible when its density is Sharatchandra and Rhode, 1994). constant. In practice, to simplify the analysis of complex Sharatchandra 1994). a detailed and flowRhode, description (Bourantas et al., 2016; constant. practice, to simplify the analysis complex A fluid isitIn considered incompressible when itsof density is tain Sharatchandra and Rhode, 1994). systems, is common to consider as incompressible any constant. In practice, to simplify the analysis of complex To complexity of systems, it is common to consider as incompressible any To reduce reduce the the and complexity of fluid fluid systems, systems, aa common common Rhode, 1994). systems, it common consider as incompressible any constant. practice, toto the analysis of that complex fluid whose density variation is very very small, i.e., can Sharatchandra systems, itInis isdensity common tosimplify consider as small, incompressible any To reduce reduce the the complexityis of of fluid systems, systems, common engineering simplification to consider the fluid a fluid whose variation is i.e., that can To complexity fluid aa as common engineering simplification is to consider the fluid as a oneonefluid whose density variation is very small, i.e., that can systems, it is common to consider as incompressible any be considered as negligible. The incompressibility approxfluid whose density variation is incompressibility very small, i.e., that can To engineering simplification is to consider the fluid as a reduce the complexity of fluid systems, a common be considered as negligible. The approxdimensional flow described by a set of partial differential engineering simplification is to consider the fluid as a oneonedimensional flow described by a set of partial differential be considered as negligible. The approxfluid whose density variation is incompressibility verytosmall, i.e., that can imation depends on the conditions which the fluid is be considered as negligible. The incompressibility approxdimensional flow described aa set of partial differential engineering simplification is by to consider the fluid as a oneequations (PDEs). Additionally, to perform numerical imation depends on the conditions to which the fluid is dimensional flow described by set of partial differential equations (PDEs). Additionally, to perform numerical imation depends on the to fluid be considered The incompressibility approxsubjected Inasan annegligible. adiabatic process thewhich flow the of either either imation depends on the conditions conditions to which the fluid is is equations (PDEs). (PDEs). Additionally, toofperform perform numerical flow by aapproximated setto partial by differential computations the PDEs are sets subjected ... In adiabatic process the flow of aaa dimensional equations Additionally, numerical thedescribed PDEs are approximated by sets of of subjected In an adiabatic process the flow of either imation depends on the conditions to which the fluid is gas or a liquid can be considered incompressible when subjected . In an adiabatic process the flow of either a computations computations the PDEs PDEs are approximated approximated by sets to of equations (PDEs). Additionally, to perform numerical ordinary differential equations (ODEs). This allows gas or a liquid can be considered incompressible when computations the are by sets of ordinary differential equations (ODEs). This allows to gas or a liquid can be considered incompressible when subjected . In an adiabatic process the flow of either a its Mach number, the ratio between the fluid and the gas Mach or a liquid canthe be ratio considered incompressible ordinary differential equations (ODEs). This allows to the PDEs are approximated by sets of its number, between the fluid andwhen the computations describe the fluid in different pipe sections. ordinary differential equations (ODEs). This allows to the fluid in different pipe sections. its number, the ratio between the the gas orvelocities, a liquid can be small, considered incompressible sound is very very i.e., the the compressibility of a describe its Mach Mach number, the ratio between the fluid fluid and andwhen the describe the fluid in different pipe sections. ordinary differential equations (ODEs). This allows to sound velocities, is small, i.e., compressibility of a the fluid in different pipe sections. sound is very small, i.e., aa describe its Mach number, the speed ratio regimes. between the Mach fluid numbers and of the On hand, port-Hamiltonian (PH) flow is velocities, associated with For sound velocities, iswith very small, i.e., the the compressibility compressibility of On the the other other hand,inthe the port-Hamiltonian (PH) framework framework describe the fluid different pipe sections. flow is associated speed regimes. For Mach numbers On the other hand, the port-Hamiltonian (PH) flow is associated with speed regimes. For Mach numbers sound velocities, is very small, i.e., the compressibility of a is a useful mathematical tool to model systems (van less than 0.3 density variations are very small and the On the other hand, the port-Hamiltonian (PH) framework framework flow is associated with speed regimes. For Mach numbers is a useful mathematical tool to model systems (van der der less than 0.3 density variations are very small and the is useful mathematical tool to model systems (van der other hand, the port-Hamiltonian (PH) framework less 0.3 are very small and flow is associated with variations speed For Mach Schaft and Jeltsema, 2014). In aa PH model the incompressible assumption is very approximation is aa the useful mathematical tool model systems (van and der less than than 0.3 density density variations areuseful very small numbers and the the On Schaft and Jeltsema, 2014). In to PH model the inputs inputs and incompressible assumption is a aregimes. very useful approximation Schaft and Jeltsema, 2014). In to a PH PH model the inputs and a useful mathematical tool model systems (van der incompressible assumption is aa very useful approximation less than 0.3 density variations are very small and the is outputs are conjugated in the power sense, i.e., the prod(Johnson, 2016). This assumption for the airflow is usual Schaft and Jeltsema, 2014). In a model the inputs and incompressible assumption is very useful approximation outputs are conjugated in the power sense, i.e., the prod(Johnson, 2016). This assumption for the airflow is usual outputs are conjugated in the power sense, i.e., the prodSchaft and Jeltsema, 2014). In a PH model the inputs and (Johnson, 2016). This assumption for the airflow is usual incompressible assumption is and a very useful approximation uct between inputs and outputs represents the supplied in aerospace engineering (Wu Cao, 2015), and bioengioutputs are conjugated in the power sense, i.e., the prod(Johnson, 2016). This assumption for the airflow is usual between inputs and outputs represents the supplied in aerospace engineering (Wu and Cao, 2015), and bioengi- uct uct between inputs represents the supplied are conjugated inoutputs the of power sense, the prodin engineering (Wu and Cao, 2015), and (Johnson, 2016). This assumption for the airflow is usual outputs power. Additionally, dynamic the system is neering applications (Cal et al., al., 2017). uct between inputs and and outputs represents the supplied in aerospace aerospace engineering (Wu and2017). Cao, 2015), and bioengibioengipower. Additionally, dynamic of the systemi.e., is charactercharacterneering applications (Cal et power. Additionally, dynamic of the system is characteruct between inputs and outputs represents the supplied neering applications (Cal et al., 2017). in aerospace engineering (Wu and Cao, 2015), and bioengiized by a non-negative function that represents the power. Additionally, dynamic of the system is characterneering et treatment al., 2017).(Hager, 2010) or ar- ized by a non-negative function that represents the total total In areas applications such as as waste waste(Cal water ized by by Additionally, non-negative function that represents the total power. dynamic of thewhose system is characterneering applications (Cal et treatment al., 2017).(Hager, 2010) or ar- ized storage in the system and rate of change In areas such water aaenergy non-negative function that represents the total storage energy in the system and whose rate of change In areas such as waste water treatment (Hager, 2010) or arterial blood flow analysis (Guidoboni et(Hager, al., 2009) 2009) theorpipe pipe In areas suchflow as waste water treatmentet 2010) ar- storage in the system and whose rate of change ized by aenergy non-negative function that represents the total is bounded by the power supplied. A port-Hamiltonian terial blood analysis (Guidoboni al., the storage energy in the system and whose rate of change bounded by the power supplied. A port-Hamiltonian terial blood flow analysis (Guidoboni et al., 2009) the In areas such as waste water treatment 2010) orpipe ar- is presents geometric variations that the terial blood flow analysis (Guidoboni et(Hager, al., changes 2009) thein pipe is bounded by the power supplied. A port-Hamiltonian storage energy in the system and whose rate of change presents geometric variations that produce produce changes in the model is focused in describe the energy flux in aa system and is bounded by the power supplied. A port-Hamiltonian model is focused in describe the energy flux in system and presents geometric variations that produce the terialvelocities blood flow analysis et al., changes 2009) thein flow and energy(Guidoboni losses associated with the the pipe presents geometric variations thatassociated produce changes inpipe the is model is of focused infeatures describe the energyAflux flux in aa system system and bounded by the power supplied. port-Hamiltonian provide useful as stability in Lyapunov sense, flow velocities and energy losses with model is focused in describe the energy in and provide of useful features as stability in Lyapunov sense, flow energy with presents geometric variations thatassociated produce changes inpipe the model expansion and and contraction. The fluid behavior inthe irreguflow velocities velocities and energy losses losses associated within the pipe provide of usefulnorm features asthe stability in Lyapunov sense, is focused in describe energy flux in a system the state space is given by the system energy, the expansion and contraction. The fluid behavior irreguprovide of useful features as stability in Lyapunov sense, state space norm is given by the system energy, and the expansion and contraction. The fluid behavior in irreguflow velocities energy losses associated with pipe the lar geometries is commonly studied through distributedexpansion and and contraction. The fluidthrough behavior inthe irreguthe state space norm is given by the system energy, the provide of useful features as stability in Lyapunov sense, system is described through basic interconnectable energy lar geometries is commonly studied distributedthe state space norm is given by the system energy, the system is described through basic interconnectable energy lar geometries is commonly studied through distributedexpansion and contraction. The fluid behavior in irreguparameters models. These models use numerical techlar geometries is commonly studied use through distributedsystem is described through basic interconnectable energy the state space norm is given by the system energy, the storage elements, between others. These features provide parameters models. These models numerical techsystem is described through basic interconnectable energy storage elements, between others. These features provide parameters models. These numerical techlar geometries is commonly studied use through distributedparameters models. These models models use numerical tech- system storage elements, between others. These features(van provide is described through basic interconnectable energy nice properties such as passivity and scalability der storage elements, between others. These features provide  This work was supported in part by CONICYT through properties such as passivity and scalability (van der  parameters These inmodels useCONICYT numericalthrough tech- nice This work models. was supported part by nice properties properties such as et passivity and scalability (van der storage elements, between others. These features(van provide Schaft, 2017; Duindam al., 2009).  nice such as passivity and scalability der  grands Bec. Doc. Nac. by 2017-21170472, This CONICYT-PFCHA work was was supported supported in part part by CONICYT FONDEthrough Schaft, 2017; Duindam et al., 2009). This work in CONICYT through grands CONICYT-PFCHA Bec. Doc. Nac. 2017-21170472, FONDESchaft, 2017; Duindam et al., 2009). nice properties such as passivity and scalability (van der  Schaft, 2017; Duindam et al., 2009). CYT 1181090, MEC 80170066 and BASAL FB0008, also by the grands CONICYT-PFCHA Bec. Doc. Nac. 2017-21170472, FONDEThis work was supported in part by CONICYT through grands1181090, CONICYT-PFCHA Bec. Doc. Nac. 2017-21170472, The behavior of fluids using infinite-dimensional portCYT MEC 80170066 and BASAL FB0008, alsoFONDEby the The behavior of fluids using infinite-dimensional portSchaft, 2017; Duindam et al., 2009). Agence Nationale de la RechercheDeutsche Forschungs GemeinCYT 1181090, 1181090, MEC MEC 80170066 80170066 and BASAL BASAL FB0008, also also by by the the grands Bec. Doc. Nac. 2017-21170472, CYT and FB0008, Agence CONICYT-PFCHA Nationale de la RechercheDeutsche Forschungs FONDEGemeinThe behavior behavior of of fluids fluids using using infinite-dimensional portHamiltonian has studied as The infinite-dimensional portHamiltonian formulations formulations has been been studied in in works works as schat (ANR-DFG) project INFIDHEM, ID ANR-16-CE92-0028 Agence Nationale de la RechercheForschungs CYT 1181090, MEC 80170066 and Deutsche BASAL FB0008, alsoGemeinby and the Agence Nationale de la RechercheDeutsche Forschungs Gemeinschat (ANR-DFG) project INFIDHEM, ID ANR-16-CE92-0028 and Hamiltonian formulations has been studied in works as The behavior of fluids using infinite-dimensional portvan der Schaft and Maschke (2002); Altmann and Schulze Hamiltonian formulations has been studied in works as the Agence Nationale de la Recherche project EIPHI-BFC ANR-17schat (ANR-DFG) project INFIDHEM, ID ANR-16-CE92-0028 and van der Schaft and Maschke (2002); Altmann and Schulze Agence Nationale de la la RechercheForschungs ANR-17Gemeinschat (ANR-DFG) project INFIDHEM, ID ANR-16-CE92-0028 and the Agence Nationale de RechercheDeutsche project EIPHI-BFC van der Schaft and Maschke (2002); Altmann and Schulze Hamiltonian formulations has been studied in works as (2017). However to simulate these model it is necessary van der However Schaft and Maschke (2002); Altmann and Schulze EURE-0002 the Agence Agence Nationale de la la INFIDHEM, Recherche project project EIPHI-BFC ANR-17ANR-17schat (ANR-DFG) project ID ANR-16-CE92-0028 and (2017). to simulate these model it is necessary the Nationale de Recherche EIPHI-BFC EURE-0002 (2017). However to simulate these model it is necessary van der Schaft and Maschke (2002); Altmann and Schulze (2017). However to simulate these model it is necessary EURE-0002 the Agence Nationale de la Recherche project EIPHI-BFC ANR-17EURE-0002 (2017). However to simulate these model it is necessary EURE-0002 2405-8963 © 2019, IFAC (International Federation of Automatic Control) Hosting by Elsevier Ltd. All rights reserved.

Copyright © 2019 IFAC 102 Copyright 2019 IFAC 102 Control. Peer review© responsibility of International Federation of Automatic Copyright © under 2019 IFAC IFAC 102 Copyright © 2019 102 10.1016/j.ifacol.2019.08.018 Copyright © 2019 IFAC 102

2019 IFAC CPDE-CDPS Oaxaca, Mexico, May 20-24, 2019

Luis A. Mora et al. / IFAC PapersOnLine 52-2 (2019) 102–107

to apply spatial discretization methods that preserve the PH structure, as they proposed in Kotyczka and Maschke (2017) and Trenchant et al. (2018). In an incompressible flow, under adiabatic conditions and neglecting the gravitational effects, the fluid only stores kinetic energy, i.e., the incompressibility hypothesis neglects the potential energy, due to the change of density, and hence makes the interconnection between sections non-causal. In this paper we propose a scalable lumped port-Hamiltonian model to describe the behavior of incompressible fluids in irregular geometries. To this end we relax the incompressible hypothesis on some infinitesimal section to obtain a causal interconnection between different fluid sections. This paper is organized as follows. In Section 2 we present the basics on PHS, the fluid behavior assumptions and develop the PHS model describing the fluid dynamic in a pipe with two sections. In Section 3 we extend the model to describe a pipe with n internal sections and dissipation associated with geometry variations in the pipes. In Section 4 we present two numerical examples to illustrate the behavior of the proposed model. Finally, in Section 5 we present the conclusion and discussions on future work. 2. A SIMPLE MODEL FOR THE INCOMPRESSIBLE FLUID 2.1 Port-Hamiltonian systems (PHS) In this section we give a brief description on PHS. For a detailed overview see for instance van der Schaft and Jeltsema (2014). PHS are defined as follows Definition 1. Consider x ∈ Rn as the vector of state variables and H : Rn → R as a scalar energy function, also called Hamiltonian. Consider u ∈ Rm and y ∈ Rm as the input and output respectively, such that y and u are conjugated in power sense. Thus, a PHS is given by following differential equation x˙ = [J(x) − R(x)]∂x H + g(x)u (1) y = g(x)T ∂x H T

103

fluid sections k+1

k k−1

k

k+1

Nodes Fig. 1. Coupling incompressible fluid sections using a compressible behavior at the nodes and (5) represent the mass and momentum conservation laws for fluids, receptively (Bird et al., 2014). In a fluid the pressure is associated with the potential energy. Neglecting the gravitational effects and considering an isothermal fluid, the potential energy can be associated with the fluid work. However, under the incompressibility assumption the fluid work is zero and an instantaneous transmission of pressure from the boundary is assumed. We shall assume that the variations of the fluid density are negligible, i.e., ∆ρ ≪ ρ. Thus, we consider that the flow in one section behaves as an incompressible fluid and that the compressibility effects occur in some infinitesimal volume, that shall be denoted by node and which is located in the coupling zone between two adjacent sections, as shown in Figure 1. To relax the incompressible hypothesis we consider the following assumption. Assumption 1. During the expansion or contraction of a node the mass variations are much less than the volume variations. Hence, a uniform density distribution and a constant mass in each node is assumed, i.e., the following relationship is satisfied ρj Vj = κ

(6)

where ρj and Vj are the density and volume in node j respectively, and κ is the corresponding node mass. Then from this assumption, we develop a fluid model that allows a causal interconnection between adjacent fluid sections, as it’s shown below.

(2)

T

where J(x) = −J (x), R(x) = R (x) ≥ 0 are n × n real matrices called interconnection and dissipation matrix respectively, and g(x) ∈ Rn×m is the input matrix. The Hamiltonian H is a non-negative function whose rate of change is bounded as follow H˙ ≤ y T u (3) where the product y T u represents the instantaneous power supplied to the system. 2.2 The fluid and the incompressibility assumption The fluid dynamic is defined by a set of partial differential equations (PDEs) whose one-dimensional forms are given by ∂t ρ + ∂z (ρv) = 0 (4) ρ 2 ρ∂t v + ∂z v + ∂z p = 0 (5) 2 where ρ and v are, respectively, the fluid density and velocity, and p is the static pressure. The PDEs in (4) 103

2.3 Modeling a fluid section The balance equation (5) represents the momentum conservation law without losses. However, when the pipe presents variations in its geometry, as contraction, expansion and changes in direction, it is known that energy losses occurs in the flow (Brodkey and Hershey, 2003). These losses are given by local turbulences that contribute significantly to the pressure drop. Thus, to include these effects and considering the incompressibility assumption, we rewrite the momentum balance as ρ0 ∂t v = −∂z P − fd (z)

(7)

where ρ0 is the reference fluid density, P = 12 ρ0 v 2 + p is the total pressure field and fd (z) is a force associated with mechanical energy losses. Integrating (7) over the corresponding volume of section i, using the Leibnitz integral rule, the Gauss divergence theorem and considering an uniform cross-sectional area in the section, we obtain a finite-dimensional description of the fluid as follows,

2019 IFAC CPDE-CDPS 104 Oaxaca, Mexico, May 20-24, 2019



ρ0 ∂t vdVi = −



∂z P dVi −

Luis A. Mora et al. / IFAC PapersOnLine 52-2 (2019) 102–107



fd (z)dVi  ρ0 Vi v˙ i = Ai P1i − Ai P2i − Ai fd (z)dz

(8)

where v˙ i is the time derivative of the average velocity in section i, Ai is the cross-sectional area, and P1i and P2i are the average pressures of the inlet and outlet boundary surfaces S1i and S2i respectively. According to Mulley (2004) the loss of mechanical energy in pipe contractions and expansions has been empirically associated with the dynamical pressure of the fluid. Thus, in a volume Vi we use the following approximation  1 fd (z)dz ≈ λi ρ0 vi2 (9) 2 where the term (ρ0 /2)λi vi2 is the average pressure drop associated with the mechanical energy losses and λi is a dimensionless loss factor. Substituting (9) in (8), we can describe the dynamic behavior of the fluid in a section as ρ0 ρ0 Vi v˙ i = −Ai λi vi2 + Ai P1i − Ai P2i (10) 2 Defining πi = ρ0 Vi vi as the fluid momentum in section i and Ki = 12 Cxi πi2 as the corresponding kinetic energy, where Cxi = 1/(ρ0 Vi ) , the port-Hamiltonian formulation of (10) is given by   P (11a) π˙ i = −ri ∂πi Ki + [Ai −Ai ] 1i P2i     Q1j Ai (11b) = ∂ K −Ai πi i −Q2j

λ i πi is term associated with the dissipawhere ri = Ai ρ02V i tion in section i, and Q1i and Q2i are the inlet and outlet volumetric flows. The definition of the dissipation factor λi depends on geometric changes in the pipe and will be discussed below.

hypothesis of Assumption 1, we use the bulk modulus definition (Murdock, 1993), βS = −Vj (dpj /dVj ), and (6) to write the differential of the pressure in a node as βS dρj (14) dpj = ρj Considering that pj = 0 when ρj = ρ0 and solving (14), the pressure behavior in the j-th node is given by   ρj (15) pj (ρj ) = βS ln ρ0 Considering a very small mass in a node, the kinetic energy can be neglected. Thus, the energy in an arbitrary node j, is given by the potential energy associated with the corresponding density variations. From the first law of thermodynamics and considering an adiabatic process, the potential energy in node j, Ej , is given by the fluid work, dEj = −pj dVj . Using (6), the differential of the potential energy can be written as κ dEj = pj 2 dρj (16) ρj where pj is the node pressure. From (16), the pressure in node j can be described as pj = (ρ2j /κ)∂ρj Ej . The potential energy is then given by the following non-negative function ρj − ρ0 (1 + ln(ρj /ρ0 )) Ej (ρj ) = κβS (17) ρj ρ0 Finally, using (17) we can rewrite (13b) as the PHS    2  Q1j 2 ρ /κ /κ −ρ (18a) ρ˙j = 0∂ρj Ej + j j Q2j    2  ρj /κ p1j (18b) = ∂ E −p2j −ρ2j /κ ρj j 2.5 The complete model

2.4 The interconnection nodes The model proposed in (11) represents the fluid behavior in an arbitrary section i. However, note that due to the causality of the model the input and output ports are not compatible for interconnection width adjacent sections. To solve this obstacle, we define an infinitesimal coupling zone between fluid sections, that shall be denoted as node. The conservation mass law in (4) is then rewritten as ∂t ρj = −ρj ∂z vj (12) Integrating (12) in the volume of node j and using the Gauss divergence theorem, the model that describe the fluid behavior  in a node can be deduced as follow  (13a) ∂t ρj dVj = −ρj ∂z vj dVj

ρj ρ˙j = (Q1j − Q2j ) (13b) Vj where the Q1j and Q2j are the average flows in the inlet and outlet cross-sectional surfaces, S1j and S2j , of the corresponding node, respectively. Note that the inputs ports in (13b) are compatible with the outputs ports of adjacent fluid sections described by (11). To complete the coupling it is necessary to obtain an expression that relates the node pressure changes with the corresponding density variations. Then, using given the constant mass 104

Consider a pipe composed in two sections, i.e., i ∈ {1, 2} and j = 1. From (11), the sections are described by   P (19a) π˙ 1 = −r1 ∂π1 K1 + [A1 −A1 ] 11 P21     Q11 A1 (19b) = ∂ K −Q21 −A1 πi 1   P π˙ 2 = −r2 ∂π2 K2 + [A2 −A2 ] 12 (20a) P22     A2 Q12 (20b) = ∂ K −Q22 −A2 π2 2

where K1 = 12 Cx1 π12 and K2 = 12 Cx2 π22 are the kinetic energies in corresponding sections. On the other hand, from (18), the node is given by    Q11  (21a) ρ˙1 = 0∂ρ1 E1 + ρ21 /κ −ρ21 /κ Q21    2  ρ1 /κ p11 ∂ E (21b) = −p21 −ρ21 /κ ρ1 1 where E1 is given by (17). Note that the input flows in (21) are compatible with the outputs Q21 in (19) and Q12 in (20). Similarly, the outputs in (21) are compatible width

2019 IFAC CPDE-CDPS Oaxaca, Mexico, May 20-24, 2019

Luis A. Mora et al. / IFAC PapersOnLine 52-2 (2019) 102–107

the inputs P21 in (19) and P12 in (20). Thus, we obtain the following PHS   ρ21 0 −A −r 1 1    κ  ∂ π 1 H   π˙ 1 2    ρ π˙ 2 =  0 ∂π H −r2 A2 1   ∂ 2H  κ ρ˙1  ρ1  ρ2 ρ2 0 A1 1 −A2 1 κ κ   A1 0   P11 + 0 −A2 (22a) P22 0 0     ∂ π H  1 Q11 A1 0 0 ∂ π2 H (22b) = 0 −A2 0 −Q22 ∂ρ1 H where Q11 and Q22 are the inlet and outlet volumetric flows of the pipe, P11 and P22 the corresponding inlet and outlet pressures, and the total energy is given by 2 ρ1 − ρ0 (1 + ln(ρ1 /ρ0 ))  1 H = κβS Cxi πi2 + (23) ρ1 ρ0 2 i=1 3. A SCALABLE MODEL The model proposed in (22) describes the behavior of two fluid sections. However, it is straightforward to extend the model to n fluid sections. First, for n fluid sections we T define  π = [π1 , π2 , · · · , πn ] as the momenta vector and n K = i=1 Ki as the total kinetic energy, from (11) we describe the fluid behavior in n pipe sections as     u P (24a) π˙ = −R∂π K + [gπ1 −gπ2 ] π1 + g 1 uπ2 P2    T    gπ1 yπ1 Q1 (24b) K, = g T ∂π K = ∂ π T yπ2 −Q2 −gπ2

where P1 and P2 are the external pressures in the inlet and outlet boundary sections, respectively, Q1 and Q2 the T corresponding volumetric flows, uπ1 = [P12 , · · · , P1n ] is T  the inlet pressure vector, uπ2 = P21 , · · · , P2(n−1) is the T outlet pressure vector, yπ1 = [Q12 , · · · , Q1n ] and yπ2 =  T Q21 , · · · , Q2(n−1) are the internal inlet and outlet volumetric flow vectors, R ∈ Rn×n is a dissipation matrix associated with energy losses due to geometry changes, and matrices g ∈ Rn×2 and {gπ1 , gπ2 } ∈ Rn×(n−1) are given by     01×(n−1) diag ([A1 , . . . , An−1 ]) gπ1 = gπ1 = 01×(n−1) diag ([A2 , . . . , An ])   A1 0 g= 0 0 0 −An where 0 is a zero matrix of suitable dimensions. To obtain a description of the n − 1 nodes we define ρ = n−1 T [ρ1 , ρ2 , · · · , ρn−1 ] as the density vector and E = j=1 Ej as the fluid potential energy. Thus, from (18) the dynamic in the nodes is   uρ1 −g (25a) ρ˙ = 0∂ρ E + [gρ ρ] uρ2    T  gρ yρ1 (25b) = ∂ E yρ2 −gρT ρ 105

105

where uρ1 is the inlet flow vector of the nodes, uρ2 is the outlet flow vector, yρ1 and yρ2 are the inlet and outlet node pressures, respectively, and gρ ∈ R(n−1)×(n−1) is given by   2 ρ2n−1 ρ1 ,..., gρ = diag κ κ Notice that the input and output ports of (25) are compatible with the corresponding ports (24). Thus, for the interconnection of the fluid sections and the nodes we have      uπ1 0 0 0 −I yπ1 uπ2  0 0 I 0  yπ2  (26)  u  = 0 −I 0 0   y  ρ1 ρ1 I 0 0 0 uρ2 yρ2 Finally, using (24), (25) and (26), the complete model is         −R J ∂π H π˙ g P1 = (27a) + T ρ˙ 0 P2 −J 0 ∂ρ H      ∂π H  Q1 = gT 0 (27b) ∂ρ H −Q2 where P1 and P2 are the boundary pressures of pipe, Q1 and Q2 are the corresponding inlet and outlet volumetric flows, J = (gπ1 − gπ2 ) gρT and the total energy is given by H=

n n−1   1 πi2 ρj − ρ0 (1 + ln(ρj /ρ0 )) + κβS 2 ρ V ρj ρ0 0 i i=1 j=1

(28)

The definition of the loss factor λi depends on the geometry variation of the pipe and the class of bifurcation junctions (Mulley, 2004; P´erez-Garc´ıa et al., 2010). In this work we shall only consider the loss factors associated with sudden expansions and contractions. In a sudden expansion, when the fluid enters a section with enlarged cross-sectional area, a jet is formed as the fluid separates from the wall of the smaller pipe section. This jet expands until it fills the entire area and some fluid break away and circulates in the corner of the expanded section (Mulley, 2004). In this case the loss factor is given by 2  Ai e λi = 1 − (29) Ai+1 In a sudden contraction given the reduction of the pipe, the fluid accelerates as it enters the smaller section. The loss factor is given by (Mulley, 2004)   Ai 1 c λi = 1− (30) 2 Ai−1 According to Brodkey and Hershey (2003), the loss factor in the inlet depends of the entrance pipe geometry and is usually less than 0.78, while the loss factor associated with the outlet is equal to 1. From these considerations the dissipation matrix in (27) is   ρ 0 λ 1 π1 ρ0 πn R = diag A1 , . . . , An 2V1 2Vn where each λi , i ∈ [1, n−1] are define by λei or λci depending of the internal geometry and the loss factor in the last section is λn = 1. 4. NUMERICAL EXAMPLES In this section we present two examples to illustrate the model. First a pipe with uniform cross sectional area with four fluid sections and uniform geometry is considered. Secondly a pipe with two cross sectional area changes,

2019 IFAC CPDE-CDPS 106 Oaxaca, Mexico, May 20-24, 2019

4

0.6

π1 π2 π3 π4

2 0

0

5

10

15

∆ρ1 ∆ρ2 ∆ρ3

0.4 0.2 0

Normalized Energy

·10−4

∆ρi (%)

πi (kg × m/sec)

6

Luis A. Mora et al. / IFAC PapersOnLine 52-2 (2019) 102–107

0

time (msec)

(a)

5

10

15

1 Kin. Ener. Pot. Ener.

0.8 0.6 0.4 0.2 0

0

5

10

15

time (µsec)

time (msec)

(b)

(c)

Fig. 2. Simulation results for example 1. (a) Momentum behavior in each fluid section, (b) Density variation in each node, (c) Normalized kinetic (red line) and potential (blue dashed line) energies one contraction and one expansion, and six fluid sections is considered. The ODE15s solver of MATLAB is used in both examples. The fluid is air with density ρ0 = 1.1376kg/m3 at 35◦ C and bulk modulus βS = 142 × 103 P a. In practice a fluid is considered incompressible if the density variations are negligible. This implies that the potential energy of incompressible fluids is negligible in comparison with the corresponding kinetic energy. To evaluate that the model proposed in (27) satisfy these considerations, we analyze the density changes ∆ρi in each node using ρj − ρ0 (31) ∆ρi = 100 ρ0 and the normalized kinetic and potential energies in the ¯ and E, ¯ respectively, where system, K n 1 πi2 i=1 2 ρ0 V i ¯ = K (32) H n−1 ρj − ρ0 (1 + ln(ρj /ρ0 )) j=1 κβS ρj ρ0 ¯= (33) E H 4.1 Example: A pipe with equal cross-sectional areas Consider a pipe with transversal area A0 and length L, divided in 4 sections with volume V0 = A0 L/4. The PHS (27) is given by  2    −ρ1 0 0 A0 0 2 2  A0  0 0   ρ1 −ρ22 0 2  G= J=   0 0  0 ρ −ρ k 2 3 2 0 −A0 0 0 ρ3   ρ0 π4 R = diag 0, 0, 0, s4 A0 2V0 Setting A0 = 5×10−4 m2 , L = 0.1m, κ = 1×10−10 kg, P1 = 800P a and P2 = 0P a, the simulation results are shown in Figure 2. Notice that the momentum presents a fast convergence to a uniform behavior in all sections (Figure 2.a). The density variation is less than 0.6% (Figure 2.b), i.e., for practical considerations it can be neglected. The energy is shown in Figure 2.c. Here the potential energy is significant only in the first microseconds of the simulation, being the kinetic energy is dominant. All these results are in concordance with the conditions used traditionally to assume the incompressible hypothesis (Johnson, 2016). 106

4.2 Example. A pipe with different cross sectional areas

A0 , V0 1

A1 , V1 3

A2 , V2 4

2 5

6

Fluid sections Fig. 3. Pipe with irregular geometry Consider a pipe with sections that contracts and expands, as shown in Figure 3. The PHS is given by   −A0 ρ21 0 0 0 0  A0 ρ2 −A0 ρ2 0 0 0  1 2   2 2  1 0 A 0 0 ρ −A ρ 1 3 1 2   J=  0  0 A1 ρ23 −A1 ρ24 k 0   0 0 0 A2 ρ24 −A2 ρ25  0 0 0 0 A2 ρ25  T A0 0 0 0 0 0 G= 0 0 0 0 0 −A2   ρ0 λ 2 π2 λ 4 π4 π6 R = diag 0, s2 A0 , 0, s4 A1 , 0, s6 A2 2 V0 V1 V2 2

where λ2 = (1/2)(1 − A1 /A0 ) and λ4 = (1 − A1 /A2 ) are the loss factors associated with a sudden pipe contraction and expansion Hager (2010), respectively. Setting A0 = 5× 10−4 m2 , A1 = 0.5A0 , A2 = 2A0 , V0 = 1 × 10−5 m3 , V1 = 0.5V0 , V2 = 2V0 , κ = 1 × 10−10 kg, P1 = 800P a and P2 = 0P a, we obtain the simulation shown in Figure 4. The air momentum and velocity in each section are shown in Figure 4.a and Figure 4.c, respectively. Notice that the momentum presents a fast convergence to the same value for all sections, and the velocities depend on the section area. Figure 4.b shows the density changes in the nodes. Notice that the density variation is less than 1 percent, similarly to the previous example. Figure 4.d shows the potential and kinetic energies. 5. CONCLUSION In this paper, a simple scalable PHS model for incompressible fluids in irregular geometries has been proposed.

2019 IFAC CPDE-CDPS Oaxaca, Mexico, May 20-24, 2019

Luis A. Mora et al. / IFAC PapersOnLine 52-2 (2019) 102–107

0.4 0.2 0

0

10

20

time (ms)

(a)

v1 v2 v3 v4 v5 v6

40 30 20 10 30

0

Normalized Energy

50 ∆ρ1 ∆ρ2 ∆ρ3 ∆ρ4 ∆ρ5

vi (m/s)

∆ρi (%)

0.6

0

10

20 time (ms)

(b)

30

40

107

1 0.8

Kin. Ener. Pot. Ener.

0.6 0.4 0.2 0

0

5

10

15

time (µs)

(c)

Fig. 4. Simulation results for example 2. (a) Density variation in each node, (b) Velocities behavior in each fluid section, (c) Normalized kinetic (red line) and potential (blue dashed line) energies. Allowing small density variations in some infinitesimal pipe section we obtain a fluid description that allows the coupling of two adjacent incompressible fluid sections with a suitable energy transference. Two numerical examples have been carried out to show that the model preserve the features of incompressible fluids, density variations in coupling zones are very small, considered in practice as negligible, and the fluid potential energy is negligible compared with the kinetic energy. In future works, we consider to extend this method to pipes with time variable geometries, such as vocal folds and blood-vessels. REFERENCES Altmann, R. and Schulze, P. (2017). A port-Hamiltonian formulation of the Navier–Stokes equations for reactive flows. Systems & Control Letters, 100, 51–55. Bird, R.B., Stewart, W.E., Lightfoot, E.N., and Klingenberg, D.J. (2014). Introductory Transport Phenomena. Jhon Wiley & Sons, United States of America. Bourantas, G.C., Cheeseman, B.L., Ramaswamy, R., and Sbalzarini, I.F. (2016). Using DC PSE operator discretization in eulerian meshless collocation methods improves their robustness in complex geometries. Computers & Fluids, 136, 285–300. Brodkey, R. and Hershey, H. (2003). Transport Phenomena: A Unified Approach. Chemical engineering series. Brodkey Publishing. Cal, I.R., Cercos-Pita, J.L., and Duque, D. (2017). The incompressibility assumption in computational simulations of nasal airflow. Computer Methods in Biomechanics and Biomedical Engineering, 20(8), 853–868. Duindam, V., Macchelli, A., Stramigioli, S., and Bruyninckx, H. (2009). Modeling and Control of Complex Physical Systems. Springer Berlin Heidelberg, Berlin, Heidelberg. Guidoboni, G., Glowinski, R., Cavallini, N., Canic, S., and Lapin, S. (2009). A kinematically coupled timesplitting scheme for fluid–structure interaction in blood flow. Applied Mathematics Letters, 22(5), 684–688. Hager, W.H. (2010). Losses in Flow. In Wastewater Hydraulics, 17–54. Springer Berlin Heidelberg, Berlin, Heidelberg, second edition. Johnson, R. (2016). The Handbook of Fluid Dynamics. CRC Press, 2nd edition. Kotyczka, P. and Maschke, B. (2017). Discrete portHamiltonian formulation and numerical approximation for systems of two conservation laws. at - Automatisierungstechnik, 65(5), 308–322. 107

Mulley, R. (2004). Flow of Industrial Fluids: Theory and Equations. CRC Press. Murdock, J.W. (1993). Fundamental Fluid Mechanics for the Practining Enegineer, volume 82 of Mechanical Engineering. Marcel Dekker, Inc. P´erez-Garc´ıa, J., Sanmiguel-Rojas, E., and Viedma, A. (2010). New coefficient to characterize energy losses in compressible flow at T-junctions. Applied Mathematical Modelling, 34(12), 4289–4305. Sharatchandra, M.C. and Rhode, D.L. (1994). New, strongly conservative finite-volume formulation for fluid flows in irregular geometries using contravariant velocity components: Part 1. theory. Numerical Heat Transfer, Part B: Fundamentals, 26(1), 39–52. Trenchant, V., Ramirez, H., Le Gorrec, Y., and Kotyczka, P. (2018). Finite differences on staggered grids preserving the port-Hamiltonian structure with application to an acoustic duct. Journal of Computational Physics, 373, 673–697. van der Schaft, A.J. and Maschke, B.M. (2002). Hamiltonian formulation of distributed-parameter systems with boundary energy flow. Journal of Geometry and Physics, 42(1-2), 166–194. van der Schaft, A. (2017). L2-Gain and Passivity Techniques in Nonlinear Control. Communications and Control Engineering. Springer International Publishing, Cham, 3rd edition. van der Schaft, A. and Jeltsema, D. (2014). PortHamiltonian Systems Theory: An Introductory Overview. Foundations and Trends® in Systems and Control, 1(2-3), 173–378. Wu, Z. and Cao, Y. (2015). Numerical simulation of flow over an airfoil in heavy rain via a two-way coupled Eulerian–Lagrangian approach. International Journal of Multiphase Flow, 69, 81–92.