An Object-Oriented Library for Heat Transfer Modelling and Simulation in Open Cell Foams

This w ork has b een submitted to IF A C for p ossible publication An Ob ject-Orien ted Library for Heat T ransfer Mo delling and Sim ulation in Op en Cell F oams ? T obias M. Sc heuermann ∗ P aul Kot yczk a ∗ Christian Martens ∗ Haithem Louati ∗∗ Bernhard Masc hk e ∗∗ Marie-Line Zanota ∗∗∗ Isab elle Pitault ∗∗ ∗ T e chnic al University of Munich, Dep artment of Me chanic al Engine ering, Chair of Automatic Contr ol, Boltzmannstr ae 15, 85748 Gar ching, Germany. ∗∗ Univ. Lyon, Universit Claude Bernar d Lyon 1, CNRS, LAGEPP UMR 5007, 43 b oulevar d du 11 novembr e 1918, F-69100 Vil leurb anne, F r anc e. ∗∗∗ Univ. Lyon, CNRS, CPE Lyon, UCBL, LGPC UMR 5285, 43 b oulevar d du 11 novembr e 1918, F-69100 Vil leurb anne, F r anc e. Abstract: Metallic op en cell foams hav e m ultiple applications in industry , e. g. as catalyst supp orts in c hemical pro cesses. Their regular or heterogeneous microscopic structure determines the macroscopic thermo dynamic and c hemical prop erties. W e presen t an ob ject-oriented p ython library that generates state space models for sim ulation and control from the microscopic foam data, whic h can be imported from the image processing tool iMorph. The foam topology and the 3D geometric data are the basis for discrete mo deling of the balance la ws using the cell metho d. While the material structure imp oses a primal chain complex to define discrete thermo dynamic driving forces, the in ternal energy balance is ev aluated on a second c hain complex, which is constructed by top ological duality . The heat exchange b etw een the solid and the fluid phase is describ ed based on the a v ailable surface data. W e illustrate in detail the construction of the dual c hain complexes, and we show ho w the structured discrete mo del directly maps to the softw are ob jects of the python co de. As a test case, w e present simulation results for a foam with a Kelvin cell structure, and compare them to a surrogate finite elemen t mo del with homogeneous parameters. Keywor ds: P ort-Hamiltonian systems, metallic foam, cell metho d, distributed parameter systems, discrete mo deling, geometric discretization, pro cess systems, simulation 1. INTRODUCTION Metallic foams are a type of material that is used in mul- tiple w a ys for industrial purp oses. Two classes of metallic foams are distinguished: closed and open cell foams. In closed cell foams, the fluid phase is encapsulated in closed ca vities inside the foam. Op en cell foams ha v e connected p orous cells so that the fluid can flow through the material. In this pap er, we will concentrate on the latter. Due to their high surface to volume ratio, open cell metallic foams are used in catalytic reactors, see e. g. F rey et al. (2016). In order to design and control the chemical processes in a reactor, numerical mo dels of the thermodynamic b eha viour are needed. Existing approac hes use effective prop erties, e. g. from v olume a v eraging o v er Cartesian unit cells (Quin tard et al., 1997). ? This work was supp orted by Deutsche F orsch ungsgemeinschaft (pro ject num ber 317092854) and Agence Nationale de la Rec herc he (ID ANR-16-CE92-0028), pro ject DF G-ANR INFIDHEM. With the use of tomograph y , precise 3D vo xel data of a given foam sample can be generated and topological as w ell as geometric data can b e extracted using image pro cessing softw are like iMorph (Brun et al., 2008). W e will show an approac h to set up a numerical mo del for the heat transfer on op en cell foams that is directly based on the p ossibly heterogeneous foam top ology . Microscopic material parameters and the exact geometry complete the mo del in the discrete constitutiv e equations. The separation of a (Dirac) interconnection structure to describ e the structural exchange of pow er (or the time deriv ative of another appropriate p otential) via pairs of conjugated port v ariables from material-dep endent con- stitutiv e equations and energy storage, is at the heart of the Port-Hamiltonian (PH) framework, see e. g. Duindam et al. (2009) for an ov erview. In Seslija et al. (2014), and later for non-uniform boundary conditions in Kotyczk a and Maschk e (2017), the discrete mo delling of conserv a- tion la ws on dual chain complexes w as presented. The preliminary w ork (Sc heuermann et al., 2019) illustrates the discrete mo delling of heat transfer and exchange on op en cell foams. In this pap er, we adopt this paradigm for the computer- based modelling and simulation of heat transfer on op en cell foams. W e presen t the necessary extensions for the classification of topological ob jects from the regular 2D case as presented in Kotyczk a and Maschk e (2017) to irregular 3D meshes in Section 2. The structured repre- sen tation of the coupled heat equation on dual complexes is presented in Section 3, while we show ho w this mo del directly maps to the ob ject oriented p ython co de in Section 4. A n umerical example is given in Section 5, and the pap er closes with final remarks and an outlo ok in Section 6. 2. IMAGE PR OCESSING The input data for model generation and simulation is obtained from the image processing tool iMorph 1 . iMorph can extract the structure of the foam from 3D tomography pictures. A typical example is shown b elow. Fig. 1a shows the image of an open cell foam sample, while Fig. 1b displa ys the extracted solid graph. (a) Surface (b) Graph Fig. 1. Metallic op en cell foam (Kelvin cells) Besides the solid no des (Fig. 2a) and struts (Fig. 2b), whic h are represented by the edges of the solid graph, iMorph identifies cells (Fig. 2c) in the fluid phase. These cells are connected b y so-called “windo ws” (Fig. 2d). 3. DISCRETE HEA T EQUA TION ON DUAL COMPLEXES The mo delling is based on the cell metho d, see Alotto et al. (2013) for an in troduction to this numerical sc heme with references to the original works (T onti, 2001) and appli- cations. The PH framework explicitly considers op en sys- tems, i. e. systems with b oundary energy flo w, see Seslija et al. (2014) for the discrete modelling of conserv ation laws and v an der Schaft and Masc hk e (2013) for PH systems on graphs. W e follow the regular 2D approac h describ ed in Kot yczk a and Masc hk e (2017). The heterogeneous 3D case considered here requires some adaptations and additions, whic h are illustrated b elo w. 3.1 Cel ls, Chains and Chain Complex The top ology and geometry of the foam is describ ed in a structured w a y using j -dimensional cells, or in short 1 http://imorph.sourceforge.net/ (a) No de (b) Strut (c) Cell (d) Window Fig. 2. Structures in op en cell foams defined b y iMorph “ j -cells” 2 , see Arnold (1989), Section 35.D or Flanders (1989), Section 5.5. A j -cell is a geometric ob ject that consists of a conv ex p olyhedron D ⊂ R j , a differentiable f : D → M on the n -dimensional manifold M and an orien tation. A formal sum of j -cells is called j -chain. The linear vector space of j -chains on a tessellation K is denoted C j ( K, R ). The b oundary of each j -cell consists of a j − 1-chain and is found b y applying the b oundary op erator ∂ j . Applying the b oundary op erator t wice to a j - c hain results in an empty set, which is the central property of a chain complex, see e. g. J¨ anic h (2001), Section 7.6. The spaces of j -chains, j = n, . . . , 0, which, connected via the b oundary op erators, form an n -complex, can b e represen ted in a sequence diagram: C n ( K, R ) ∂ n − → C n − 1 ( K, R ) ∂ n − 1 − → . . . ∂ 1 − → C 0 ( K, R ) (1) In the follo wing, w e call a n -c hain with the collection of all j -cells, j = 0 ...n , app earing in the sequence ab o v e, an n - complex. The sym b ol ∂ will be used for both the boundary op erator and its matrix representation, i. e. an incidence matrix. F or our application, only the case with n = 3 is relev ant, so we will restrict ourselves to this case. 3.2 Definition of the Primal 3-Complex The primal 3-complex is initially given by the structure of the solid phase. Since an n -complex can b e seen as a gen- eralized dir e cte d graph, orientations hav e to b e assigned to all j -cells. The no des (0-cells) and edges (1-cells) of the primal 3-complex can b e tak en directly from the graph generated with iMorph. F aces (2-cells) corresp ond to the iMorph windows. The windo ws that enclose a fluid cell define a v olume (3-cells). The following classification of inner and b order j -cells is necessary for the direct imp osition of b oundary conditions 2 The term “cell” is used in tw o contexts, that should not b e confused with each other: It is used in iMorph to describ e a cavity in the foam or a j -dimensional geometric ob ject. Therefore, the latter is alwa ys denoted as j -cell. Fig. 3. Primal complex in the numerical model. T o realize Neumann b oundary conditions (NBCs), i. e. heat flux b oundary conditions on the appropriate dual ob jects, see Subsection 3.3, additional b order no des must b e defined, which lead to additional edges, faces and v olumes on a thin, artificial boundary la y er. Fig. 3 shows a minimal example for a 3-complex with the orien tation of the j -cells. The different categories of j -cells are describ ed b elo w. Inner no des: Solid no des inside the domain are called inner no des and are denoted b y n i ∈ N i . Bor der no des: Solid nodes on the boundary are called b order nodes and are denoted b y n b ∈ N b . A t these no des, a Diric hlet b oundary condition (DBC) is imp osed. A dditional b or der no des: These no des, denoted b y n b ∈ N B , are not a represen tation of a solid no de, but an in tersection of a strut with the b oundary (border edge, see b elow). Through the dual face to this edge, see next subsection, a NBC is imp osed. (a) Inner (b) Border (c) Additional b order Fig. 4. Primal no des Inner e dges: Inner edges e i ∈ E i are connections of the inner no des and b order no des. They represen t struts that are en tirely inside the domain or on its boundary ( n i with n i , n b with n b and n i with n b ). Bor der e dges: Border edges e b ∈ E b connect inner no des to additional b order no des an represen t struts that cross the system b oundary . A dditional b or der e dges: These edges e B ∈ E B ha v e no represen tation in the solid graph, but they are necessery to fill the en tire domain with v olumes. (a) Inner (b) Border (c) Additional b order Fig. 5. Primal edges Inner fac es: F aces b elonging to windows that are entirely inside the domain, are called inner faces f i ∈ F i . Bor der fac es: Border faces f b ∈ F b b elong to windows that are not completely inside the domain. A dditional b or der fac es: These faces f B ∈ F B do not b elong to a window, but they are necessery to fill the entire domain with v olumes. (a) Inner (b) Border (c) Additional b order Fig. 6. Primal faces Inner volumes V olumes that lie inside the domain or on the b oundary with a DBC are inner v olumes v i ∈ V i . Bor der volumes Only volumes at the b oundary with a NBC are b order v olumes v b ∈ V b . (a) Inner (b) Border Fig. 7. Primal v olumes R emark 1. On first sigh t it ma y seem, that some no des in Fig. 3, esp ecially at the corners, are missing. How ever, they were left out inten tionally . Similar to the additional b oundary edges in the 2D case in (Kot yczk a and Masc hk e, 2017, Fig. 5), that hav e no nodes at the corner of the face, the volumes in the 3D case can also hav e corners without no des. In 3D, there can even b e kinks in the faces without ha ving a “real” edge at that p osition. These kinks are dra wn with dotted lines and lie inside a face and hav e therefore no effect on the result of the b oundary op erator applied to the face. R emark 2. The categorization differs from Kotyczk a and Masc hk e (2017), because the physical v ariables are as- signed to the geometric ob jects in another wa y . This is b ecause the energy balance is ev aluated on the dual v ol- umes instead of the primal faces. Subsequen tly , the driving force is ev aluated on the primal instead of the dual edges. 3.3 Construction of the Dual 3-Complex The dual 3-complex is defined b y construction. F or b etter visibilit y , only one dual j -cell is drawn in Fig. 8. The same pro cedure is rep eated for all other primal j -cells. A barycentric dual is used, as in Alotto et al. (2013). This means, that the dual no de is located at the barycen tre of the primal volume (Fig. 8a). Accordingly , a dual edge in tersects with its primal face at the barycen tre of the face (Fig. 8b) and the dual face intersects with the primal edge also at the barycen tre of the edge (Fig. 8c). The dual complex is completed with the dual volumes around the primal no des (Fig. 8d). (a) No de (b) Edge (c) F ace (d) V olume Fig. 8. Primal and asso ciated dual cells 3.4 Discr ete PH R espr esentation F or a structured discrete mo del of the heat transfer on the foam, we start with the well-kno wn heat equation with distributed parameters on a single phase, x ∈ Ω ⊂ R 3 , t ∈ R + 0 , c ˙ T ( x , t ) = λ ∆ T ( x , t ) . (2) T ( x , t ) denotes the temp erature, the heat capacit y c and the thermal conductivity λ are assumed to be constant. W e rewrite (2) in p ort-Hamiltonian form (neglecting for the momen t the b oundary conditions) using the inner energy densit y u ( x , t ) as state and T ( x , t ) as co-state/effort 3 ,  ˙ u f  =  0 − div − grad 0   T φ  . (3) φ ( x , t ) and f ( x , t ) denote the vectors of heat flux and the temp erature gradien t as the thermodynamic driving force. The mo del is completed with the constitutiv e la ws φ = λ f u = cT . (4) The discrete model is found by integrating the equations o v er the appropriate j -chains of the dual and the primal complex, resp ectively , as indicated in T able 3 with sup er- script s or f referring to the solid or the fluid phase. 3 Which is the conjugate quantit y w. r. t. the artificial p otential R Ω 1 c u 2 ( x , t ) d x . T able 1. j -chains and asso ciated quantities j -chain (Integral) physical quantit y Primal no de n k T emp erature T s / f k Primal edge e k Driving force (temp erature difference) F s / f k Dual face ˆ f k Heat flow rate ˆ Φ s / f k Dual volume ˆ v k Energy ˆ U s / f k If k is the index for a dual con trol volume, and the set I ( k ) con tains the indices of the b oundary faces, the discrete energy balance on such a control volume can b e written for b oth the solid and the fluid phase as ∂ ∂ t ˆ U s k = − X l ∈I ( k ) ˆ Φ s k,l − ˆ Φ sf k (5a) ∂ ∂ t ˆ U f k = − X l ∈I ( k ) ˆ Φ f k,l + ˆ Φ sf k . (5b) The heat flow ˆ Φ sf k represen ts the heat transfer betw een b oth phases. The temperature differences along a strut (index k , 1 and 2 refer to the start and end node) for b oth phase, as w ell as b et w een b oth phases are F s k = − ( T s k, 2 − T s k, 1 ) , F f k = − ( T f k, 2 − T f k, 1 ) (6a) F sf k = T s k − T f k . (6b) Finally , the discrete appro ximations of the constitutiv e equations (4) for b oth phases, together with the heat transfer mo del b et w een b oth phases are ˆ Φ s k = λA s k F s k | r k, 2 − r k, 1 | , ˆ Φ f k = λA f k F f k | r k, 2 − r k, 1 | (7a) ˆ Φ sf k = αA sf k F sf k (7b) ˆ U s k = V s k c s T s k , ˆ U f k = V f k c f T f k . (7c) The discrete geometry parameters (note that (5) and (6) con tain only top ological information) are given in T able 2. T able 2. Geometry parameters Parameter Definition r k Position vector of no de n k A s / f k Solid / fluid part of the area of f k A sf k Contact area of the phases in v k V s / f k Solid / fluid part of the v olume of v k T o obtain a numerical model of the heat transfer in the complete foam, we collect the whole set of v ariables ˆ U k , F k , T k and ˆ Φ k in the vectors ˆ U i / b , F i / b , T i / b and ˆ Φ i / b , whic h represen t inner / b order co-c hains as algebraically dual ob jects to the j -c hains of the primal and the dual complex 4 . The result is the follo wing system of equations, where ˆ d 3 ii / bi = − ( d 1 ii / bi ) T and d 1 ii / ib denote the co-incidence matrices (i. e. the transp osed b oundary matrices) b etw een 4 F or a given j − 1-co-chain c j − 1 , which contains the integral v alues of a quantit y ov er j − 1-chains, the duality pairing, see Seslija et al. (2014), h c j − 1 , ∂ j c j i = h d j c j − 1 , c j i (8) defines the co-b oundary op erator d j . The sequence of spaces of co- chains and co-boundary op erators defines a co-c hain complex C 0 ( K, R ) d 1 − → C 1 ( K, R ) d 2 − → . . . d n − → C n ( K, R ) . (9) faces and volumes on the dual complex and no des and edges on the primal complex, resp ectiv ely 5 .        ˙ ˆ U s i ˙ ˆ U f i F s i F f i F sf i        =       0 0 ( − d 1 ii ) T 0 I 0 0 0 ( − d 1 ii ) T − I d 1 ii 0 0 0 0 0 d 1 ii 0 0 0 − I I 0 0 0              T s i T f i ˆ Φ s i ˆ Φ f i ˆ Φ sf i        +       0 0 ( − d 1 ib ) T 0 0 0 0 ( − d 1 ib ) T d 1 ib 0 0 0 0 d 1 ib 0 0 0 0 0 0           T s b T f b ˆ Φ s b ˆ Φ f b     (10) The subscripts i and b denote the lo cations (in the in terior or at the b oundary) of the j -c hains, on which the discrete quan tities are defined as presen ted in the previous sub- sections. Note that the skew-symmetry of the first matrix mimics the formal skew-adjoin tness of the matrix op erator in (3). The mo del is again completed b y the constitutiv e la ws ˆ U s / f = C s / f T s / f (11a) ˆ Φ s / f = Λ s / f F s / f (11b) with the diagonal matrices C s / f = diag( V s / f k c s / f ) (12a) Λ s / f = diag λA s / f k | r k, 2 − r k, 1 | ! (12b) 4. IMPLEMENT A TION F or the implementation of the 3-complexes and their j - cells, we chose an ob jected oriented approac h using the programming language Python . The goal of this implemen- tation is to represent the relations b etw een j -cells in the co de. The general structure of the core classes is sho wn as a UML diagram in Fig. 9. F or b etter treatment of j - cells with rev erse orien tation, the implemen tation includes some more classes than sho wn, but they follo w the same arc hitecture. T o a v oid redundancy , the classes of all j -cells inherit from a Cell class where common prop erties like num bering or lab eling are implemen ted. Node , Edge , Face and Volumes classes must b e instan tiated from top to bottom, since ev ery class needs an aggregation of its predecessor. This approac h relates to the application of the co-b oundary op erator as in (9). All ob jects of j -cells are collected in an instance of the PrimalComplex class, where the classification is imple- men ted and the incidence matrices are calculated. The DualComplex automatically generates all dual j -cells. 5. NUMERICAL EXAMPLE The presen ted approach is applied to a grid based on Kelvin cells as shown in Fig. 10, whic h has 848 degrees 5 F or the relations of co-incidence matrices between the dual com- plexes, see Kotyczk a and Masc hke (2017) or Seslija et al. (2014). 2..* 2 2..* 2 2..* 3..* 2..* 3..* 1..2 4..* 1..2 4..* 1 1..* 1..* 1..* 1..* 1 1..* 1..* 1..* 1..* Cell DualCell No de DualNo de Edge DualEdge F ace DualF ace V olume DualV olume PrimalComplex DualComplex Fig. 9. Simplified UML diagram of freedom (DOFs). F or b etter replicability , the grid is constructed, so that we can test the n umerical method without dep ending on user settings in iMorph or the need to com pensate p ossibly o ccuring defaults in the iMorph result. On the top and b ottom boundary , a DBC is applied (Fig. 11a). The other b oundaries hav e a NBC (Fig. 11b), in our case the heat transfer is set to 0, meaning it is p erfectly isolated at theses b oundaries. Fig. 10. Geometry of the example foam (a) Dirichlet b ound- ary condition (DBC) (b) Neumann b ound- ary condition (NBC) Fig. 11. Primal v olumes The material parameters used in the simulation are given in T able 3. T able 3. Material parameters Dimensions l × w × h 40 × 40 × 40 mm Mass m 16 . 463 g Density of aluminium ρ s 2 . 7 × 10 − 3 g mm 3 Density of air ρ f 1 . 204 × 10 − 6 g mm 3 Heat capacity of Al c s 0 . 897 J g K Heat capacity of air c f 1 . 005 J g K Thermal conductivity of Al λ s 0 . 2 W mm K Thermal conductivity of air λ f 2 . 6 × 10 − 5 W mm K Heat transfer co efficient α 1 . 0 × 10 − 4 W mm 2 K Surrogate thermal diffusivity a eff 1 . 85 mm 2 s Fig. 12 shows the transient b ehaviour of the temp erature on 4 selected no des. T 0 is the constant temp erature at the b ottom b oundary , while T 3 is increased at the top. T 1 and T 2 are the temp eratures of tw o no des at different heights close to the fron t b oundary . F or comparison, a Finite Element (FE) simulation with 18 081 DOFs was p erformed with a surrogate parameter for the diffusivity a eff = λ eff ρ eff c eff using FEniCS (Alnæs et al., 2015). The results are shown with markers and the sup erscript c . Fig. 12. T ransien t b eaviour of the foam The p erfect matc hing of b oth our simulation based on the structured mo del with the surrogate FE simulation is due to tw o facts: (a) the surrogate diffusivity has b een determined by curve fitting and (b) did we only consider the “harmless” case of pure heat conduction without the consideration of con v ectiv e transp ort. 6. CONCLUSION AND OUTLOOK W e show ed a structured approach to obtain a numerical mo del of heat transfer through metallic op en cell foams, in which the separation of top ology (expressed in terms of co-incidence matrices) on the one side and geometry and material parameters (constitutive equations) on the other side mimics the PH structure of the lo cal PDE model. The mo del allo ws to iden tify macroscopic foam parameters, and can b e used for design optimization and (after possible mo del reduction) for con trol. The model structure directly maps to the ob jects and dep endencies of the ob ject oriented python library , which can read top ology and geometry data ov er an interface to the iMorph image processing softw are. W e presented the sim ulation of a realistic foam mo del and its comparison to a FE sim ulation with surrogate effectiv e parameters. A t the moment, w e w ork in sev eral directions: (a) the sim- ulation of real foam data and comparison with the exp er- imen tal data obtained at LGPC Lyon, (b) the integration of con v ection in the mo del and (c) impro ving robustness of our mo del generation with resp ect to artefacts like not fully connected graphs from image pro cessing. A CKNO WLEDGEMENTS The authors cordially thank Jerme Vicente from Univer- sit y Aix-Marseille for the help with iMorph and in par- ticular for implementing mo difications in the new iMorph v ersion that allow us to directly access all necessary ob jects and parameters. REFERENCES Alnæs, M.S., Blech ta, J., Hake, J., Johansson, A., Kehlet, B., Logg, A., Ric hardson, C., Ring, J., Rognes, M.E., and W ells, G.N. (2015). The fenics pro ject version 1.5. Ar chive of Numeric al Softwar e , 3(100). doi:10.11588/ans.2015.100.20553. Alotto, P ., F reschi, F., Repetto, M., and Rosso, C. (2013). The Cel l Metho d for Ele ctric al Engine ering and Multiphysics Pr oblems . Springer Berlin Heidelb erg. doi:10.1007/978-3-642-36101-2. Arnold, V.I. (1989). Mathematic al Metho ds of Classical Me chanics . Springer New Y ork. doi:10.1007/978-1-4757-2063-1. Brun, E., Vicente, J., T opin, F., and Occelli, R. (2008). Imorph: A 3d morphological to ol to fully analyse all kind of cellular materials. Cel lular Metals for Structural and F unctional Applications . Duindam, V., Macchelli, A., Stramigioli, S., and Bruyninc kx, H. (2009). Mo deling and Contr ol of Complex Physic al Systems . Springer-V erlag GmbH. Flanders, H. (1989). Differ ential F orms with Applic ations to the Physic al Scienc es . Dov er Publications Inc. F rey , M., Romero, T., Roger, A.C., and Edouard, D. (2016). Open cell foam catalysts for co2 methanation: Presentation of coating pro cedures and in situ exothermicit y reaction study by infrared thermography . Catalysis T o day , 273. doi: 10.1016/j.cattod.2016.03.016. J¨ anich, K. (2001). V e ctor Analysis . Springer. Kotyczk a, P . and Maschk e, B. (2017). Discrete p ort-hamiltonian formulation and numerical approximation for systems of t w o conserv ation laws. at-Automatisierungste chnik , 308–322. Quintard, M., Kavian y , M., and Whitaker, S. (1997). Two-medium treatment of heat transfer in p orous media numerical results for effective prop erties. A dvanc es in Water R esour c es, . Scheuermann, T.M., Kotyczk a, P ., Zanota, M.L., Pitault, I., Louati, H., and Maschk e, B. (2019). Numerical appro ximation of heat transfer on heterogenous media. P AMM , 19(1). doi: 10.1002/pamm.201900372. Seslija, M., Sc herpen, J.M.A., and v an der Schaft, A.J. (2014). Explicit simplicial discretization of di stributed-parameter p ort- Hamiltonian systems. Automatic a , 50(2), 369–377. doi: 10.1016/j.automatica.2013.11.020. T onti, E. (2001). A direct discrete form ulation of field la ws: the cell method. Computer Modeling in Engine ering and Scienc es , 2(2), 237–258. doi:10.3970/cmes.2001.002.237. v an der Schaft, A.J. and Maschk e, B. (2013). Port-hamiltonian systems on graphs. SIAM Journal of Contr ol and Optimization , 51(2), 906–937.