A port-Hamiltonian approach to modeling the structural dynamics of complex systems
Alexander Warsewa, Michael Böhm, Oliver Sawodny, Cristina Tarín
AA Port-Hamiltonian Approach to Modeling the Structural Dynamics ofComplex Systems (cid:73) , (cid:73)(cid:73) Alexander Warsewa a, ∗ , Michael Böhm a , Oliver Sawodny a , Cristina Tarín a a Institute for System Dynamics, University of Stuttgart, Waldburgstraße 17/19, 70563 Stuttgart, Germany
Abstract
With this contribution, we give a complete and comprehensive framework for modeling the dynamics of complexmechanical structures as port-Hamiltonian systems. This is motivated by research on the potential of lightweightconstruction using active load-bearing elements integrated into the structure. Such adaptive structures are ofhigh complexity and very heterogeneous in nature. Port-Hamiltonian systems theory provides a promisingapproach for their modeling and control. Subsystem dynamics can be formulated in a domain-independent wayand interconnected by means of power flows. The modular approach is also suitable for robust decentralizedcontrol schemes.Starting from a distributed-parameter port-Hamiltonian formulation of beam dynamics, we show the appli-cation of an existing structure-preserving mixed finite element method to arrive at finite-dimensional approxima-tions. In contrast to the modeling of single bodies with a single boundary, we consider complex structures com-posed of many simple elements interconnected at the boundary. This is analogous to the usual way of modelingcivil engineering structures which has not been transferred to port-Hamiltonian systems before. A block diagramrepresentation of the interconnected systems is used to generate coupling constraints which leads to differentialalgebraic equations of index one. After the elimination of algebraic constraints, systems in input-state-output(ISO) port-Hamiltonian form are obtained. Port-Hamiltonian system models for the considered class of systemscan also be constructed from the mass and stiffness matrices obtained via conventional finite element methods.We show how this relates to the presented approach and discuss the differences, promoting a better understand-ing across engineering disciplines. A Matlab framework is available on http://github.com/awarsewa/ph_fem/ to facilitate the application of the methods to different problems.
Keywords: port-Hamiltonian systems, modeling, finite element, structural dynamics, adaptive structures
1. Introduction
Lightweight structures have become a reality for mass-sensitive applications, such as large civil engineeringstructures and many more. For passive structures, these designs present in most cases a minimum in terms ofrequired mass under given safety limitations and user comfort constraints. However, it is possible to stay withinthese limits while further reducing the total embodied mass significantly by introducing adaptive structures ,which is referred to as ultralightweight design . Through their various actuators, these structures can react andadapt to external loads and disturbances – both static and dynamic – in order to minimize element stress andat the same time maximize lifetime expectancy [1].In light of worldwide population growth, demographic aging and ever increasing standard-of-living, construc-tion activities are at an all-time high and are expected to increase further drastically within the next 20 to 30years. This in mind, the technology of active and adaptive structures can help save millions of tons of concreteand steel and significantly reduce waste production and CO -emissions of the construction industry [2].Nevertheless, methods for analysis and control of adaptive structures need to be developed in order to makethis technology readily available. This requires accurate mathematical models of these structures. The currentstate of the art for modeling as well as static and dynamic analysis of civil engineering structures is finite elementmodeling (FEM). However, most often the elements are connected via their nodal displacements, not by theenergy flow from one to another. This technique impedes the integration of actuators into the dynamic modelof a structure. At the same time, it can be difficult for these structures to design decentralized control schemeswith guaranteed asymptotic stability. (cid:73) c (cid:13) http://creativecommons.org/licenses/by-nc-nd/4.0/ (cid:73)(cid:73) Published version: https://doi.org/10.1016/j.apm.2020.07.038 ∗ Corresponding author.
Email address: [email protected] (Alexander Warsewa)
Preprint submitted to Elsevier August 19, 2020 a r X i v : . [ phy s i c s . c o m p - ph ] A ug n this context, port-Hamiltonian systems have the potential to provide a framework for modeling bothpassive and active structures and for seamless and straightforward integration of actuators from different physicaldomains, such as electric motors or hydraulic cylinders. Furthermore, robust decentralized control strategiesbased on passivity based control can naturally be used for these structures.Truss structures and frames can be modeled as an interconnection of one-dimensional linear elastic beamelements and rods. Although our use case is adaptive structures, the structural dynamics of other complexsystems, such as in aerospace, can also be described this way. For simulation and control purposes, it isexpedient to describe the system dynamics in state space. In the case of port-Hamiltonian systems, a commonstate space representation is the input-state-output (ISO) form [3]˙ x = ( J − R ) d H d x + Gu , y = G T d H d x (1)with the state vector x ∈ R n which contains the energy variables of the physical system, the skew-symmetricinterconnection matrix J ∈ R n × n , the positive semi-definite symmetric dissipation matrix R ∈ R n × n and theoutput matrix G ∈ R n × n u . Inputs and outputs are collocated and n u is the number of both inputs and outputs.Evaluation of the Hamiltonian H ( x ) yields the energy stored in the system. For linear port-Hamiltonian systems,the Hamiltonian is quadratic and can be expressed as H ( x ) = x T Qx , where Q ∈ R n × n is a symmetric positivedefinite matrix. Flow and effort variables are defined as f := ˙ x and e := d H/ d x , respectively, with the powerproduct e T f ≤ M fe ¨ q + D fe ˙ q + K fe q = f ext , (2)where M fe ∈ R n DOF × n DOF is the mass matrix, D fe ∈ R n DOF × n DOF the damping matrix and K fe ∈ R n DOF × n DOF the stiffness matrix. The system’s global degrees of freedom (DOFs) are collected in the vector q ∈ R n DOF and the external forces acting on the structure in f ext ∈ R n DOF . Any such system can be written in ISO port-Hamiltonian form (1) as long as M fe and K fe are symmetric positive definite and D fe is symmetric positivesemi-definite, which is usually the case. By choosing the state vector as x = (cid:2) M fe ˙ q q (cid:3) , we obtain˙ x = (cid:20) − II (cid:21)| {z } J − (cid:20) D fe
00 0 (cid:21)| {z } R ! (cid:20) M − K fe (cid:21)| {z } Q x + (cid:20) I (cid:21)|{z} G f ext , y = G T Qx . (3)The concentrated parameter equations (2) are approximations of the dynamics of elastic mechanical systemsthat result from the spatial discretization of infinite-dimensional component models and their interconnection.In conventional FEM, this is done by applying the principle of virtual work and the Galerkin method. Bychoosing polynomial approximation bases for the DOFs, element mass and stiffness matrices directly resultfrom the integration of the potential and kinetic energy terms in the variational formulation of the systemHamiltonian. Transforming element DOFs to global coordinates and adding terms corresponding to the sameDOFs, the system mass and stiffness matrices and (2) or (3) are obtained.If we were only interested in a port-Hamiltonian formulation of linear elastic structures, there would be noneed to proceed any further. However, we are interested the interconnection of such systems with elementsfrom non-mechanical or hybrid domains. Suppose, that such elements are described by partial-differentialequations (PDEs) or include nonlinear relations. In those cases, it is not clear how to obtain a port-Hamiltonianformulation of the coupled system, given only the mechanical part (2). This is why we present the port-Hamiltonian way of modeling linear elastic structures using suitable FE approaches and a method for thecoupling of basic elements. In each step – from the infinite-dimensional equations via spatial discretization tothe interconnected concentrated parameter system – the port-Hamiltonian structure is preserved. From a controlengineer’s perspective, the passivity property of the port-Hamiltonian system representation can be exploited forcontrol design. Especially for nonlinear systems it eases stability proofs or finding suitable Lyapunov functions.Early lumping of infinite-dimensional systems with boundary controllers is also possible.The port-Hamiltonian approach does not directly produce element mass and stiffness matrices or equationsof the form (3). Even though the procedure is significantly different to established FE methods, we show that theresulting systems are identical and can be transformed to (3) when the same approximation bases are used. Tothe best of our knowledge, such a connection has not been established before. With a thorough understandingof the port-Hamiltonian way of modeling linear elastic mechanical systems and its relation to more establishedFE methods, extensions are achieved more easily and mutual understanding across disciplines is promoted. In2ddition to presenting a comprehensive method for systems composed of beams – which is the main focus ofthis contribution – coupling with the hydraulic domain is shown in an example. This emphasizes that it isstraightforward to construct models of complex multi-domain systems within the presented port-Hamiltonianframework. The dynamics of one-dimensional beam elements are by default described by PDEs. Physical systems ofthis type can be formulated as infinite-dimensional port-Hamiltonian systems without difficulties. For thecorresponding theory on distributed-parameter port-Hamiltonian systems, see e. g. [4, 5]. In order to obtainsystems of the form (1), a finite-dimensional approximation is necessary. A wide range of methods for thenumerical approximation of PDEs is available [6]. As stated above, the finite element method is a commontechnique in structural dynamics. For a comprehensive treatise of the approach, see e. g. [7] or [8]. Discretizationof infinite dimensional port-Hamiltonian systems, however, can not be done in the usual ways. Both the dualityand the passivity property need to be preserved in the finite-dimensional approximation. The former can beachieved by using mixed FEM approaches. First methods for the structure-preserving spatial discretizationof port-Hamiltonian systems were presented in [9] and [10]. They use differential forms for the numericalapproximation of one- and two-dimensional systems, where different approximation spaces are used for flowsand efforts. Their method was extended to allow for higher-order spatial derivatives by Bassi et. al [11] andfor its application to a class of irreversible thermodynamic systems by Baaiu et al. [12]. For an accurateapproximation, a relatively large number of finite elements is usually required with these approaches. Moullaet al. [13] adopted a collocation method to discretize 1D port-Hamiltonian systems. Using polynomial bases,this results in higher accuracy for lower-order approximations. In all of the approaches mentioned up to thispoint, the discretization happens in the strong form of the balance equations. This results in rather strictcompatibility conditions. Also, an input feed-through term is present in the finite-dimensional approximationwhich potentially complicates the application of control engineering methods.Cardoso-Ribeiro et al. [14] relax the compatibility conditions by using the method of [13] in the weakformulation in order to apply it to model the dynamics of a beam excited by a piezoelectric patch. Mixed finiteelement discretization in the weak form was introduced in a general way by Kotyczka et al. [15]. Their methodcan be applied to systems of any spatial dimension and resolves several issues that arise when approximatingdirectly in the strong form. It requires an additional projection to ensure non-degeneracy of the power productwhich introduces additional degrees of freedom that can be used to tune the discretized models. The boundaryports of resulting models are of alternating causality which is not always desired. In this contribution weuse the partitioned finite-element method (PFEM) of Cardoso-Ribeiro et al. [16]. It is closely related to themixed Galerkin method used in [15]. Integration by parts is, however, only applied to one of the dual balanceequations. This way, the port-Hamiltonian structure is obtained in a direct fashion and no further projection isrequired. Moreover, no feed-through term is produced in the discretization process. Causality of the boundaryports depends on which balance equation is selected for partial integration. It does not alternate for individualelements. The method was also shown to be compatible with existing finite element software. PFEM is appliedto the Mindlin plate in [17] and to the Kirchhoff-Love plate model in [18]. A symplectic Hamiltonian formulationof thin plates was also introduced by Li et al. in [19] to obtain analytic solutions for the free vibration problem.It has recently been extended to account for the buckling problem in [20] and thin cylindrical shell structures in[21]. While it constitutes a powerful method to obtain analytic solutions of boundary-value problems, it doesn’tesasily extend to systems interconnected in a complex way. All in all, PFEM is considered the most suitableapproach for the application described above.In available examples, the system usually consist of a single body and its boundary. Our use case differsin that we consider complex structures formed by an interconnection of many simple bodies at the boundary.This is a common approach in the modeling of civil engineering structures which we apply within the port-Hamiltonian framework in the following. To arrive at an automated way of generating models of arbitrarystructures composed of beam elements, we take a further look at the coupling of elements.
The structure of this contribution is as follows. In Sec. 2, distributed-parameter port-Hamiltonian modelsare presented for the load cases of a linear beam element, i. e. bending, axial loads and torsion. Regardingbending, both Euler-Bernoulli and the Timoshenko bending theory are considered. Depending on the natureof the problem examined, the formulation with higher validity can be chosen. Given the basic componentdynamics in PDE-form, the discretization process using PFEM is explained in Sec. 3. Instead of performingthe discretization process on each of the systems in Sec. 2 individually, it is stated on a more general systemclass level. The interconnection of subsystems is illustrated in Sec. 4. Based on a mechanical network diagramanalogy, an algorithm for the generation of coupling constraints is derived. After introducing and algebraicallymanipulating coupling constraints, the interconnected systems can be written as differential-algebraic equations3 bz ∂w∂z ϕw vω
AA A-A
Figure 1: A beam element defined on the interval Z = (cid:2) a b (cid:3) and its deflections (blue) due to bending ( w , ϕ ), torsional ( ω ) andaxial stress ( v ). (DAEs) of index one in ISO port-Hamiltonian form. Reduction to ordinary differential equations (ODEs) on aconstrained state space is outlined thereupon. Then, it is shown that the systems can be brought to the form(3) using simple transformations. In Sec. 5, a range of examples on the use of the approach is presented. First,a numerical analysis of the approximation error resulting from the application of PFEM on the beam equationsis carried out. This is followed by modeling the structural dynamics of a high-rise building composed of beamelements and an example of a multi-domain system with a hydraulic actuator. Concluding remarks are givenin Sec. 6. The models and methods of Secs. 2 - 4 were compiled into a Matlab framework which can be used forthe modeling of arbitrary truss structures and frames as port-Hamiltonian systems. It is available for downloadon GitHub .
2. Beam Element Models
Models of complex civil engineering structures are usually built by interconnecting basic element types withdifferent parameters. Usually, the number of elements forming a model is much higher than the number of basictypes it is composed of. In this article, we restrict ourselves to truss structures composed of rods and framescomposed of beam elements. Shell elements are not yet included in the framework presented in this work.Fig. 1 depicts a beam defined on the spatial domain Z := (cid:2) a b (cid:3) with spatial coordinate z which is subjectto bending, torsional and axial stress. When subject to bending, the beam experiences a deflection w ( z ) fromthe neutral axis and a rotation of the cross section by the angle ϕ ( z ). For torsional loads, the cross sectionrotates by the angle ω ( z ) about the neutral axis. Axial strain is determined by the deflection under axial loads, v ( z ). Since the load cases are considered linearly independent, their dynamics can be treated separately. Acomplete beam model thus consists of three basic components which are presented in the following. For bending,a port-Hamiltonian model is derived for both the Euler-Bernoulli and the Timoshenko beam equations. Thelatter are used when the classic assumption that the cross-section is always perpendicular to the neutral axis,i. e. ϕ = ∂w/∂z , does not hold (e. g. for thick beams). In Sec. 2.3, the model for a rod element deforming dueto axial loads is presented. Torsion is modeled according to the Saint-Venant equations in Sec. 2.4.All component models are given as distributed parameter port-Hamiltonian systems in this section. Con-centrated parameter models are obtained from the infinite dimensional systems via the procedure explained inSec. 3. In classical mechanics, the Euler-Bernoulli equation is widely used for describing the dynamics of a beamunder bending stress. Cordoso-Ribeiro et al. [14] derived a port-Hamiltonian model for a beam actuated bypiezoelectric patches using Euler-Bernoulli theory. We use the same formulation for the distributed parametersystem neglecting the actuator terms. The Euler-Bernoulli equation is given as µ ∂ w∂t = − ∂ ∂z (cid:18) EI ∂ w∂z (cid:19) , (4) http://github.com/awarsewa/ph_fem/ µ , Young’s modulus E and the second moment of area I . For the formulation ofthe system Hamiltonian, two conjugate energy variables are chosen as p A := µ ∂w∂t , q A := ∂ w∂z , (5)The letter A is used to distinguish between the variables in this system of equations and their respectivecounterparts in others. With x A = (cid:2) p A q A (cid:3) , the energy of system A can be stated in terms of it’s Hamiltonian H A ( x A ) = Z Z H A ( x A ) = 12 Z ba (cid:20) µ ( p A ) (cid:21) d z + 12 Z ba (cid:2) EI ( q A ) (cid:3) d z. (6)For all the mechanical systems modeled in this contribution, the Hamiltonian can be separated into a kineticenergy part and a potential energy part. Here, p A accounts for the former and q A for the latter form of energy.Taking the variational derivative of H A with respect to the energy variables yields the effort variables e A p := δH A δp A = ∂w∂t , e A q := δH A δq A = EI ∂ w∂z (7)where e A p ( z, t ) is the deflection velocity and e A q ( z, t ) the bending moment distributed over the interval Z . Definingthe flow variables as f A ( z, t ) := ˙ x A , system A can be written as a distributed parameter Hamiltonian systemof two conservation laws f A = (cid:20) f A p f A q (cid:21) = (cid:20) − ∂ z ∂ z (cid:21)| {z } J A (cid:20) e A p e A q (cid:21) . (8)The operator J A is formally skew symmetric according to [4] and therefore energy conserving. Substituting thedefinitions of flow and effort variables into (8) yields the Euler-Bernoulli equation (4) and an equality relation.The energy flow of this system is given by˙ H A = Z Z δH A δx A ∂x A ∂t d z = Z Z e A p ˙ x A p + e A q ˙ x A q d z = Z Z e A q ∂ z e A p − e A p ∂ z e A q d z = Z Z ∂ z e A p ∂ z e A q − ∂ z e A q ∂ z e A p d z + h e A q ∂ z e A p − e A p ∂ z e A q i ba . (9)From the second line of (9), we see that the energy flow depends only on the boundary values of the effortvariables and their derivatives. This motivates the definition of boundary inputs u A ∂ and outputs y A ∂ such thatthe Hamiltonian rate of change can be expressed as ˙ H A = ( u A ∂ ) T y A ∂ . A possible choice of the boundary variablesis given by u A ∂ = (cid:2) ∂ z e A q ( a ) − ∂ z e A q ( b ) − e A q ( a ) e A q ( b ) (cid:3) T , y A ∂ = (cid:2) e A p ( a ) e A p ( b ) ∂ z e A p ( a ) ∂ z e A p ( b ) (cid:3) T , (10)where the boundary inputs are forces and moments and the boundary outputs the collocated velocities andangular velocities. As stated in e. g. [22], other definitions of the boundary variables are also valid. However, aswe will show in Sec. 3, (10) is a “natural” choice as it has a normalizing property in the application of PFEM.The definition of u A ∂ and y A ∂ completes the formulation of (4) as distributed parameter port-Hamiltonian systemwith power exchange over the boundary. For thick beams, Eq. (4) does not accurately describe the deformation under bending. The rotation of thecross section with respect to the mid-surface normal needs to be taken into account. Including the rotationaldynamics, bending dynamics are described with a pair of equations in Timoshenko beam theory ρA ∂ w∂t = ∂∂z (cid:20) κAG s (cid:18) ∂w∂z − ϕ (cid:19)(cid:21) ρI ∂ ϕ∂t = ∂∂z (cid:18) EI ∂ϕ∂z (cid:19) + κAG s (cid:18) ∂w∂z − ϕ (cid:19) . (11)As in (4), E and I are Young’s modulus and the second moment of area, respectively. The beam’s cross sectionalarea is denoted by A , G s is the shear modulus and κ is the Timoshenko shear coefficient. A port-Hamiltonianformulation of the equations above can be found e. g. in [23]. We proceed accordingly by defining a set of energyvariables p B1 := ρI ∂ϕ∂t , p B2 := ρA ∂w∂t , q B1 := ∂ϕ∂z , q B2 := (cid:18) ∂w∂z − ϕ (cid:19) (12)5hich are collected in the state vector x B := (cid:2) p B1 p B2 q B1 q B2 (cid:3) . With the previous definitions, the systemHamiltonian can be written as H B ( x B ) = 12 Z ba (cid:20) ρI ( p B1 ) + 1 ρA ( p B2 ) (cid:21) d z + 12 Z ba (cid:2) EI ( q B1 ) + κAG ( q B2 ) (cid:3) d z. (13)Again, p B1 and p B2 are the kinetic energy variables of system B whereas q B1 and q B2 appear in potential energyterms. Parallel to the procedure in Sec. 2.1, we take the variational derivative of H B with respect to x B toobtain the effort variables for the Timoshenko beam e B p := ∂ϕ∂t , e B q := EI ∂ϕ∂ze B p := ∂w∂t , e B q := κAG s (cid:18) ∂w∂z − ϕ (cid:19) . (14)Here, e B p denotes the distributed angular velocity and e B p the deflection velocity whereas e B q is the distributedbending moment and e B q the shear force. An infinite-dimensional Hamiltonian system of two conservation lawsis obtained by defining the flow variables as f B := ˙ x B and establishing their relation to the effort variables as f B = f B p f B p f B q f B q = ∂ z
10 0 0 ∂ z ∂ z − ∂ z e B p e B p e B q e B p . (15)Compared to the Euler-Bernoulli beam, four equations instead of two are necessary to describe the dynamicsof the Timoshenko beam in Hamiltonian form, which results in higher-order concentrated parameter systems.However, the second-order spatial derivative is not present in (15). A convenient choice of the boundary variablesfor the Timoshenko beam is u B ∂ = (cid:2) − e B q ( a ) e B q ( b ) − e B q ( a ) e B q ( b ) (cid:3) T , y B ∂ = (cid:2) e B p ( a ) e B p ( b ) e B p ( a ) e B p ( b ) (cid:3) T . (16)It is easy to verify that the Hamiltonian rate of change can be expressed as ˙ H B = ( u B ∂ ) T y B ∂ . Modeling the dynamics of a linear beam is not complete without considering axial deformation. For a rodelement which only supports axial loads, the deformation along the neutral axis v ( z ) can be described by thefollowing dynamic equation µ ∂ v∂t = ∂∂z (cid:18) EA ∂v∂ z (cid:19) , (17)with the mass per unit length µ , Young’s modulus E and the beam’s cross sectional area A as defined in thepreceding sections. In accordance with the procedure for Euler-Bernoulli and Timoshenko beam, conjugateenergy variables are defined next p C := µ ∂v∂t , q C := ∂v∂z . (18)They separate the kinetic and potential energy domains and are used in the formulation of the system Hamil-tonian H C ( x C ) = 12 Z ba (cid:20) µ ( p C ) (cid:21) d z + 12 Z ba (cid:2) EA ( q C ) (cid:3) d z. (19)For the definition of the effort variables, we take the variational derivative of H C with respect to p C and q C e C p := ∂v∂t , e C q := EA ∂v∂z , (20)where e C p is the axial deflection velocity and e C q the axial force distributed over the interval Z . By defining theflow variables as f C := ˙ x C we arrive at a formulation of rod element as infinite-dimensional Hamiltonian systemof two conservation laws f C = (cid:20) f C p f C q (cid:21) = (cid:20) ∂ z ∂ z (cid:21) (cid:20) e C p e C q (cid:21) . (21)Again, we obtain a representation as distributed parameter port-Hamiltonian system by selecting boundaryinputs and outputs from the evaluation of the effort variables at the boundary u C ∂ = (cid:2) − e C q ( a ) e C q ( b ) (cid:3) T , y C ∂ = (cid:2) e C p ( a ) e C p ( b ) (cid:3) T , (22)which allows us to express the energy flow over the boundary as ˙ H C = ( u C ∂ ) T y C ∂ .6 bzˆ e p ˆ f p ˆ e p ˆ f p ˆ e p ˆ f p ˆ e p ˆ f p ˆ e p ˆ f p ˆ e p ˆ f p Figure 2: Division of a beam element into five intervals with N p = 6 supporting points with corresponding flows and efforts fornumerical approximation. Beam deformation due to torsion about the neutral axis is modeled using the Saint-Venant equation. Warpingtorsion is neglected here, assuming that cross sections are either warping-free or that warping is not constrained.Saint-Venant torsion is described by I p ∂ ω∂t = ∂∂z (cid:18) G s J t ∂ω∂ z (cid:19) , (23)where I p is the polar moment of inertia, G s the shear modulus, J t the torsion constant and ω ( z ) the torsionalangle as depicted in Fig. 1. We note that (23) and (17) are, except for the choice of parameters, identicalequations. Consequently, the port-Hamiltonian formulation will be identical to (21) when choosing p D := I p ∂ω/∂t and q D := ∂ω/∂z as the energy variables. The reader is referred to the previous section or to [24] forthe remaining steps to obtain (21).
3. Spatial Discretization
For the simulation of the dynamic response of beam elements, numerical approximation of the infinite dimen-sional port-Hamiltonian system equations derived in Sec.2 is necessary. When employing spatial discretizationmethods, care must be taken to preserve the port-Hamiltonian structure, which is why conventional approachesare not readily applicable. A spatial discretization method is structure preserving if the resulting concentratedparameter system can be written in ISO form (1). We use PFEM, as proposed by Cardoso-Ribeiro et al. [25],for reasons outlined in the introduction. In particular, due to its straightforward applicability and compatibilitywith standard FEM methods and software. In the following, the application of PFEM to systems of the form˙ x = f = J e, with J = N X i =0 (cid:20) A i ( − i +1 A T i (cid:21) ∂ i ∂z i (24)is presented. The matrices A i are quadratic with entries being either zero or one and N is the maximum orderof the spatial derivative. The system Hamiltonian can be expressed as H ( x ) = 12 Z Z x T L x d z, (25)where L can be a function of x , and the interval Z = (cid:2) a b (cid:3) as above. When observing the infinite dimensionalport-Hamiltonian systems in Sec. 2, we notice that they are all of this type with the order of the spatial derivative N ≤ L independent of x for isotropic beams. Taking for instance (15), the Timoshenko beam equationsin port-Hamiltonian form with N = 1, we get A = (cid:20) (cid:21) , A = (cid:20) (cid:21) . (26)By treating the spatial discretization process on a system class level, explicit formulation of the process for everysingle load case of the beam is avoided. At the same time, it is more evident that the method’s applicability isnot limited to truss structures and frames.In a first step, the energy variables x , flow variables f and effort variables e of a system of class (24) areapproximated using the method of weighted residuals x p ( z, t ) ≈ ˆ x T p ( t ) φ p ( z ) , x q ( z, t ) ≈ ˆ x T q ( t ) φ q ( z ) ,f p ( z, t ) ≈ ˆ f T p ( t ) φ p ( z ) , f q ( z, t ) ≈ ˆ f T q ( t ) φ q ( z ) ,e p ( z, t ) ≈ ˆ e T p ( t ) φ p ( z ) , e q ( z, t ) ≈ ˆ e T q ( t ) φ q ( z ) . (27)7ere, the basis functions φ p ( z ) and φ q ( z ) are chosen as the Lagrange polynomials of order N p − N q − p and q are used to distinguish between kinetic and potential energy related terms. Thevectors ˆ x p , ˆ f p , ˆ e p ∈ R N p and ˆ x q , ˆ f q , ˆ e q ∈ R N q are the energy variables, flows and efforts evaluated at uniformlydistributed points in Z = (cid:2) a b (cid:3) . For values of N p ≥ N q ≥
2, the boundary is always included. This isillustrated for the kinetic variables of a beam element in Fig. 2 where N p = 6 supporting points are used. Theapproximations of kinetic efforts ˆ e ip and flows ˆ f ip are included in the figure with the half arrow pointing in thepositive direction respectively. In order to be compatible with FEM-methods, PFEM proceeds via the weakform of (24) Z Z φ ( z ) f d z = Z Z φ ( z ) J e d z (28)using the approximation bases φ p ( z ) and φ q ( z ) as test functions. This transformation is similar to the Galerkinmethod, which is widely used for numerical approximation in FEM (see e. g. [8, 26]). Substituting the approx-imations (27) for energy variables, flows and efforts into (28) yields (cid:18)Z Z φ p φ T p d z (cid:19) ˆ f p = N X i =0 Z Z φ p A i ∂ i ∂z i φ T q d z ! ˆ e q (cid:18)Z Z φ q φ T q d z (cid:19) ˆ f q = N X i =0 ( − i +1 Z Z φ q A T i ∂ i ∂z i φ T p d z ! ˆ e p . (29)For systems of two conservation laws such as (24), structure-preserving mixed Galerkin approximation as pre-sented in [15] can be used. However, with PFEM, the finite-dimensional port-Hamiltonian system can beobtained in a direct fashion without further projection. This is achieved by applying integration by parts toonly one of the equations in (29). In principle, either equation can be integrated by parts. By choosing one orthe other, we decide on the causality of the boundary ports – i. e. which variables will be the boundary inputsand which the boundary outputs. When integrating the respective dual equation by parts instead, causalityof the boundary ports is reversed. For the systems defined in Sec. 2, integrating the upper equation by partsresults in forces and torques as boundary inputs and velocities and angular velocities as boundary outputs.Integrating the lower equation by parts instead, the role of boundary inputs and outputs is reversed. To beconsistent with the FEM approach, we choose to apply integration by parts N times to the equation containingthe kinetic flows. Intermediate steps of this calculation are shown in Appendix A. As a result we obtain (cid:18)Z Z φ p φ T p d z (cid:19) ˆ f p = N X i =0 ( − i Z Z ∂ i ∂z i φ p A i φ T q d z ! ˆ e q + N X i =1 i X j =1 ( − j − (cid:20) ∂ j − ∂z j − φ p A i ∂ i − j ∂z i − j φ T q ˆ e q (cid:21) ba . (30)Observing the term in round brackets on the right hand side, we note that it is the negative transpose of theterm in front of ˆ e p in the lower equation in (29). Applying Gaussian quadrature to evaluate the integrals, aconcentrated parameter approximation of the system dynamics is obtained (cid:20) M p M q (cid:21)| {z } M ˆ f = (cid:20) P − P T (cid:21)| {z } J ˆ e + (cid:20) G p (cid:21)| {z } G u ∂ (31)with the symmetric mass matrices M p ∈ R N p × N p and M q ∈ R N q × N q , the square skew-symmetric interconnectionmatrix J with P ∈ R N p × N q , the input matrix G and the boundary input variables u ∂ . The latter are the efforts φ T q e q and their spatial derivatives, appearing in the second term on the right hand side in (30), evaluated atthe boundary. According to the effort definitions in Sec. 2 in case of the beam dynamics, the inputs are forcesand bending moments. For one-dimensional systems such as the beam elements considered in this publication,the number of boundary inputs amounts to 2 N with u ∂ ∈ R N and G p ∈ R N p × N . With the boundary inputsstated for the individual systems in Sec. 2, it can be easily verified that G p is directly obtained by evaluation of φ p (and ∂ z φ p in case of the Euler-Bernoulli beam) at the boundary. This is why we regard them as “natural”choices for the application of PFEM. For the systems B, C and D, we thus get G p = I for N p = 2 N .Equation (31) is already a valid finite-dimensional approximation of the port-Hamiltonian system, but it isnot yet in the desired input-state-output form (1). To bring it to this form, we recall that f := ˙ x and transformthe energy variables ˜ x := M ˆ x . (32)8 z θ y θ z θ y θ z θ x θ z θ x z xy A CA CM MF F F F
Figure 3: Interconnection of two linear elastic elements with boundary mechanical ports.
For a concentrated-parameter port-Hamiltonian system, the effort variables are obtained by taking the gradientof the system Hamiltonian with respect to the energy variables. Thus, it remains to show thatˆ e = ∇ ˜ x H (˜ x ) . (33)First, the system Hamiltonian needs to be written in terms of ˜ x which is achieved by substitution of theapproximations (27) into (25) and application of the transformation (32) H (˜ x ) = 12 ˜ x T M − (cid:18)Z Z φ L φ T d z (cid:19) M − ˜ x . (34)From Sec. 2, we remember that e = ∂H/∂x . After substitution of the approximations (27) and transformationto the weak form, this expression becomes (cid:18)Z Z φφ T d z (cid:19) ˆ e = (cid:18)Z Z φ L φ T d z (cid:19) ˆ x . (35)Since the integral on the left hand side is M , it is easily seen that (33) is true. We define Q := M − Z Z φ L φ T d z (36)such that system (31) can be written as ˙˜ x = JQ ˜ x + Gu ∂ y ∂ = G T Q ˜ x (37)with the boundary output variables y ∂ being the input-collocated velocities and angular velocities for all systemsconsidered in this article.
4. Assembly of Complex Systems
With the systems defined in Sec. 2 and the numerical approximation method described in Sec. 3, the dynamicsof one-dimensional primitive elements can be obtained in ODE form (37) for all load cases of a linear beam. Tomodel more complex truss structures and frames, the interconnection and coupling of multiple basic elementsneeds to be taken into account. In the following, a method for the automatic assembly of arbitrary 3D structurescomposed of beams and rods is derived. Throughout this section we make use of Fig. 3 to illustrate the concepts.It depicts the interconnection of a Euler-Bernoulli beam and a rod element respectively labeled system “A” and“C”. The system equations of the Euler-Bernoulli beam are˜ f A = J A ˆ e A + G A u A ∂ , y A ∂ = ( G A ) T ˆ e A , (38)with u A∂ ∈ R the shear forces and bending moments at the boundary and y A∂ the corresponding collocatedvelocities and angular velocities. The dimension of ˜ f A depends on the order of the approximation polynomials φ p and φ q in (27). For the rod element we obtain˜ f C = J C ˆ e C + G C u C ∂ , y C ∂ = ( G C ) T ˆ e C , (39)9here u C ∂ ∈ R are the axial forces at the boundary and y C ∂ the collocated velocities.In the bottom part of Fig. 3 the interconnection of A and C is depicted in a mechanical network diagram.Specific to this block diagram representation is the appearance of a normalized orientation vector θ i at eachboundary port, where i denotes the port number. A port labeled “M” indicates a torque/angular velocitypair at the boundary and “F” marks ports with a force/velocity pair. Ports of the same type with a non-zeroscalar product of the orientation vectors can be coupled. Given the information in Fig. 3, it is possible towrite an algorithm for the automatic assembly of the coupled system. In the remainder of this section, thenecessary steps such an algorithm needs to perform are explained in more detail. First, a way to concatenatethe system matrices of individual elements is shown, followed by the generation of coupling constraints. Finally,the elimination of the generated constraints is outlined in Sec. 4.3. Without considering the coupling of the dynamic equations yet, there are several ways to concatenate thesystem matrices of systems A and C. When performing integration by parts on the equations involving thekinetic flows, the boundary inputs always act on f p and never on f q . Thus, we propose to maintain theseparation of kinetic and potential energy on concatenation. As we will later see, this also facilitates notation infurther steps. For both (38) and (39), the interconnection matrix and input matrix have the following structure J = (cid:20) J p J q (cid:21) , G = (cid:20) G p (cid:21) . (40)When defining the concatenated energy variables as ˜ x AC := (cid:2) ˜ x C p ˜ x A p ˜ x A q ˜ x C q (cid:3) T , the resulting joint systemmatrices are J AC = J C p J A p J A q J C q , G AC = G C p G A p
00 00 0 (41)and the structural separation of energy domains is maintained.
According to Newton’s second law of motion, forces and torques at each nodal point of a structure have tosum up to zero respectively for each degree of freedom considered. At the same time, velocities and angularvelocities with the same orientation have to be identical. Looking again at the example in Fig. 3, this implies v A b = v z , v C a = θ z v z (42)where v A b is the right boundary velocity in z -direction of system A, v C a the left boundary velocity in z -directionof the rod and v z the velocity in z -direction of the second node. The components of the orientation vector ofthe fifth port θ are depicted in Fig. 3. Since the Euler-Bernoulli beam cannot take loads in the x -direction, θ x does not appear in the formulation of constraints. For the boundary forces F A b and F C a we get F A b + θ z F C a + F z = 0 (43)with the additional contribution of an external force F z in z -direction.Since the boundary outputs of the concatenated system include all boundary velocities and angular velocities,we can write (42) in a more general fashion. For an arbitrary structure, we define its vector of global nodalvelocities ˙ q ∈ R n DOF with n DOF the total number of DOFs. It contains a velocity or angular velocity variablefor each DOF. This way, we can express each boundary output as a linear combination of entries of ˙ q y i∂ = ( c iv ) T ˙ q , i = 1 . . . n u (44)where n u is the total number of boundary inputs or outputs and c iv maps θ i to the corresponding nodal velocityor angular velocity. From the relation between ˙ q and y ∂ , algebraic constraints on the state variables x areobtained by eliminating ˙ q y ∂ = C v ˙ q ⇒ C ⊥ v y ∂ = 0 (45)using a left annihilator C ⊥ v of C v and the output equation of (1). Resulting algebraic constraints are of theform 0 = B T Q ˜ x . (46)For statically determinate structures, the row rank of C ⊥ v is n DOF and it follows that B ∈ R n × n c with n c = n u − n DOF . 10n a similar fashion, we can find a generalized expression for (43) and use it to reformulate the boundaryinputs u ∂ . For each global DOF, an external force or torque is introduced and they are collected in the vectorof external forces and torques f ext ∈ R n DOF . Consequently, each external force or torque can be equated with alinear combination of boundary inputs ( c ju ) T u ∂ = f j ext , j = 1 . . . n DOF C u u ∂ = f ext . (47)Usually, the number of DOFs is lower than the number of boundary inputs. With n DOF input constraints, thismeans that free boundary inputs u f ∈ R n c can be chosen. This enables a reformulation of the input term Gu ∂ as Gu ∂ = G f u f + Kf ext , (48)where G f is the resulting input matrix of the free inputs u f and K that of the external forces f ext . Assuming G f has full column rank, we can construct its left inverse G +f = ( G Tf G f ) − G Tf and rewrite u f as u f = G +f B λ (49)where λ ∈ R n c are the Lagrange multipliers corresponding to the algebraic constraints (46). Given (46), (48)and (49), the coupled structure can be written as a DAE system˙˜ x = JQ ˜ x + Kf ext + B λ ˙ q = K T Q ˜ x B T Q ˜ x (50)where the collocated outputs to the external force inputs are the global nodal velocities ˙ q . It can be shownthat the DAE system (50) is of index one in case the matrix B T ∂ H∂ ˜ x B has full rank [3], which is always thecase for the systems considered here. For now, we do not make further use of the DAE formulation. Instead,in the next section, we show how to eliminate the algebraic constraints to obtain an ODE formulation of thestructural dynamics. For a thorough treatise of differential algebraic port-Hamiltonian systems, see [27]. Systems of the form (50) can be reduced to the constrained state space X c with uniquely defined dynamics.As a consequence, the algebraic constraints are eliminated. A procedure for doing so is described in detail ine. g. [28] and [29]. We repeat the necessary steps here, to give a complete description of the methodology usedto obtain the dynamic equations for arbitrary structure models in ODE form.In a first step, the transformation V = (cid:20) B ⊥ B + (cid:21) (51)is introduced. It is composed of both the left annihilator B ⊥ and the left inverse B + of B . Left-multiplyingthe uppermost equation in (50) with V and defining a new state vector z := V ˜ x , yields˙ z = VJQ ˜ x + VKf ext + (cid:20) I (cid:21) λ (52)To retain the port-Hamiltonian form, the term Q ˜ x needs to be replaced with the gradient of H with respect to z ∇ z H = (cid:18) ∂ ˜ x ∂ z (cid:19) T ∇ ˜ x H = V − T QV − z = ˜ Qz (53)Defining ˜ J := VJV T and ˜ K := VK and separating z into z ∈ R n − n c and z ∈ R n c , the system can bereformulated as follows (cid:20) ˙ z ˙ z (cid:21) = (cid:20) ˜ J ˜ J ˜ J ˜ J (cid:21) (cid:20) ˜ Q ˜ Q ˜ Q ˜ Q (cid:21) (cid:20) z z (cid:21) + (cid:20) ˜ K ˜ K (cid:21) f ext + (cid:20) I (cid:21) λ ˙ q = (cid:2) ˜ K T1 ˜ K T2 (cid:3) (cid:20) ˜ Q ˜ Q ˜ Q ˜ Q (cid:21) (cid:20) z z (cid:21) (cid:2) I (cid:3) (cid:20) ˜ Q ˜ Q ˜ Q ˜ Q (cid:21) (cid:20) z z (cid:21) . (54)11he third equation implies ∇ z H = 0 and that z can be expressed as z = − ˜ Q − ˜ Q z . (55)Therefore, the system dynamics on the constrained state space X c are given as˙ z = ˜ J ( ˜ Q − ˜ Q ˜ Q − ˜ Q ) z + ˜ K f ext ˙ q = ˜ K T1 ( ˜ Q − ˜ Q ˜ Q − ˜ Q ) z . (56)Due to the separation of energy and co-energy related terms according to (41), the interconnection matrix ˜ J of the ODE system has the following structure˜ J = (cid:20) J p − ˜ J T p (cid:21) (57)with ˜ J p ∈ R ( n p − n c ) × ( n − n c ) , where n p is the number of potential energy variables in the DAE system (50)and ( n p − n c ) = n DOF . Since the external forces f ext and the constraint forces λ in (50) act on ˜ x p only, thetransformation (51) retains ˜ x q in the new state z . With ˜ x q retained in the state vector z , the system dynamics of the reduced order ODE system can beformulated as follows using z = (cid:2) z p ˜ x q (cid:3) T (cid:20) ˙ z p ˙˜ x q (cid:21) = (cid:20) J p − ˜ J T p (cid:21) (cid:20) ˜ Q p Q q (cid:21) (cid:20) z p ˜ x q (cid:21) + (cid:20) ˜ K p (cid:21) f ext , ˙ q = (cid:2) ˜ K T p (cid:3) (cid:20) ˜ Q p Q q (cid:21) (cid:20) z p ˜ x q (cid:21) . (58)In this representation, it becomes clear that the n q state derivatives ˙˜ x q are linearly dependent since they arecomputed using a reduced number of potential energy variables z p ∈ R n p − n c . For ease of notation, we assume n p = n q in the following, without loss of generality as long as n q > n p − n c . We can then proceed to constructan orthonormal basis T q of − ˜ J T p such that ˜ x q = T q z q . (59)The effort variables ˜ e q are then related to z q as follows˜ e q = ∇ ˜ x q H ( z ) = Q q ˜ x q = Q q T q z q . (60)This allows a formulation of algebraic constraints on ˜ e q using B ∗ = ker (cid:2) ( Q q T q ) T (cid:3) such that( B ∗ ) T Q q ˜ x q = 0 . (61)Given this relation, the process outlined in the previous section can be repeated, which results in the eliminationof linearly dependent efforts or rather state variables. Afterward, the system (58) is expressed in the newcoordinates ¯ z = (cid:2) z p z q (cid:3) T as ˙¯ z = ¯ J ¯ Q ¯ z + ¯ Kf ext , ˙ q = ¯ K T ¯ Q ¯ z (62)with ¯ z ∈ R n DOF . Note, that the reduced order system (62) is of the same order as the second order mechanicalsystem in ISO port-Hamiltonian form (3), given the global DOFs q are identical. In the next section it is shown,how to transform (62) to obtain (3) without dissipation. In the following, we assume that (3) – without the damping term – and (62) describe the same system. Inthis case, (3) with R = 0 can be obtained from (62) by means of two consecutive transformations. With thefirst transformation T u = (cid:20) ¯ K + ¯ K ⊥ (cid:21) (63)a normalization of (62) with respect to the input f ext is achieved T u ˙¯ z = T u ¯ JT T u T − T u ¯ QT − u T u ¯ z + T u ¯ Kf ext . (64)120 − − − − − − Number of elements E rr o r o f s t e i g e n f r e q u e n c y i n % (a) N p = N q = 2 for axial loads and N p = N q = 4 forbending. − − − − − − Number of elementsTimoshenkoBernoulliRodSaint-Venant (b) N p = N q = 4 for axial loads and N p = N q = 6 forbending. Figure 4: Analysis of the approximation error of PFEM for the beam models for a different number of elements and supportingpoints.
Since T u ¯ K = (cid:2) I (cid:3) T , it follows that T u ˙¯ z = (cid:2) M fe ¨ q ˙˜ z q (cid:3) T , which results in (cid:20) M fe ¨ q ˙˜ z q (cid:21) = (cid:20) P − ¯ P T (cid:21)| {z } T u ¯ JT T u (cid:20) M −
00 ¯ S (cid:21)| {z } T − T u ¯ QT − u (cid:20) M fe ˙ q ˜ z q (cid:21) + (cid:20) I (cid:21) f ext , ˙ q = (cid:2) I (cid:3) (cid:20) M −
00 ¯ S (cid:21) (cid:20) M fe ˙ q ˜ z q (cid:21) . (65)The final step to obtain (3) is the construction of a second transformation T p = (cid:20) I
00 ( − ¯ P T ) + (cid:21) , with T p T u ¯ JT T u T T p = (cid:20) − II (cid:21) . (66)When T p is applied to (65) in the same fashion as T u in (64), it follows that the transformed system must beidentical to (3) with R = 0. In case D fe in (2) is a Rayleigh or Caughey damping term, R can now be easilyconstructed from M fe and K fe . Other types of damping might involve adding dissipation terms on the PDElevel, but this is not considered here.
5. Numerical Examples
In this section, we provide numerical examples that illustrate the use of the presented methods. First, theapproximation error of PFEM for the beam equations in Sec. 2 is analyzed numerically. This is followed bya dynamic model of a high-rise building composed of many different elements which illustrates all of Sec. 4.Finally, the coupling of a mechanical system with hydraulic cylinders (i. e. a multi-domain system) is presented.
For each of the systems presented in Sec. 2, a numerical analysis of the approximation error is performedfor the case of a single beam with a quadratic cross section and
L/h = 50, where h is the cross sectionalheight. Different boundary conditions are chosen for bending and axial loads. In case of the Timoshenko andEuler-Bernoulli equations, the beam is simply supported and otherwise it is clamped at one end and free atthe other. A Poisson ratio of ν = 0 . κ = 5 /
6. Thebeam is then approximated using an increasing number of elements and the error in the first eigenfrequency iscomputed. Analytical solutions for the eigenfrequencies of each of the considered load cases can e. g. obtainedfrom [26]. The results are depicted in Fig. 4 for all of the beam equations with the number of elements rangingfrom N e = 1 . . . N p = N q = 2 supporting points were chosen in case of the rod and the Saint-Venant elementand N p = N q = 4 for both Euler-Bernoulli and Timoshenko beam. For all systems, the error starts at a value of ≈
10 % and decreases monotonically until N e ≈
40 where it starts to oscillate around a more or less stationary13 yz (a) High-rise structure
10 mm30 cm c m
12 mm 10 mm c m
20 cm 8 mm c m columns hor. diagonalsvert. diagonals bars (b) Beam and rod element cross sections Figure 5: 12 story high-rise building composed of beams (columns) and rod elements (all others) including a depiction of the beamand rod element cross sections. value for the bending elements. The results for the rod and Saint-Venant element are not discernible fromeach other. This was to be expected as systems C and D in Sec. 2 only differ in terms of parameters but areotherwise identical. A saturation of the error is not visible for the axial loads. The oscillating behavior aboutan error of ≈ × − % for the Euler-Bernoulli beam and ≈ × − % for the Timoshenko beam is explainedby the fact that the composite system becomes numerically ill-conditioned for a high number of elements. Thiseffect can be reduced by stabilizing the numerical methods involved in the system assembly and eigenfrequencycalculation.As the approximation error does not only depend on the number of elements but also on the order of theapproximation polynomials, an additional analysis was conducted with a higher number of supporting points.This time, N p = N q = 4 supporting points were chosen for the axial loads and N p = N q = 6 for bending. Theresults are shown in Fig. 4 b). With only a single element, the error is now about three orders of magnitudeslower than before for all systems. It also decreases more rapidly and is already below the minimum value visiblein Fig. 4 a) for N e = 5 elements (with the exception of the Timoshenko beam, for which a lower value resultsat N e = 74). However, for N e >
10 the systems start to become numerically ill-conditioned, which results ina saturation of the approximation error for the rod and Saint-Venant elements. Instead of a saturation, anincreasing trend is visible for the error of the simply supported Timoshenko and Euler-Bernoulli beams.The numerical analysis of the approximation error is limited by the accuracy of the employed numericalmethods and the numerical conditioning of the systems for which the error is calculated. A mathematicalanalysis of the convergence properties of PFEM for the systems considered in this publication would alleviatethose issues but is beyond the scope of this publication. With the results shown in Fig. 4, we conclude, thatreasonable errors can be achieved for a low number of approximating elements, especially when using more thanthe minimum number of supporting points required per element. However, numerical ill-conditioning becomesa problem for a high number of elements which, depending on the application, needs to properly addressed.Optimizing the algorithms for system assembly on this account is to be considered in further work.
With the model components and methods described in Secs. 2-4, a port-Hamiltonian system model forarbitrary linear elastic 3D truss structures and frames can be built. To facilitate the application of the proceduresdescribed in this article, they were compiled into a Matlab framework which is available on GitHub .In this section we illustrate the creation of a port-Hamiltonian model for a 12 story high-rise structure.The structure shown in Fig. 5 a) is entirely composed of steel bars, 36 m tall and has a 4 . × .
75 m footprint.Floor plates are represented by diagonal horizontal bars. The cross sections of columns, horizontal diagonals,horizontal bars and diagonal bracings are depicted in Fig. 5 b). Diagonals and horizontal bars are modeled as http://github.com/awarsewa/ph_fem/ − x - d i s p l a ce m e n t i n mm (a) x -displacement of node 52 H ( x ) i n k J (b) System Hamiltonian Figure 6: Structural response of the structure shown in Fig. 5 a) for an initial displacement in x -direction.Table 1: First six eigenfrequencies of the structure in Fig. 5 a). Eigenmode N p = N q = 2 and can only take axial loads. Columns are composed of two Euler-Bernoullielements – one for each bending axis – with N p = N q = 4, a rod element and a torsion element with N p = N q = 2each. We use Lagrange polynomials for Φ p and Φ q and Gauss-Legendre quadrature for evaluating the integralsin Sec. 3. The structure is composed of 4 identical modules. Each of them is composed of 12 columns of length L = 3 m, 12 horizontal bars, 8 diagonal bracings and 6 horizontal diagonals. This accounts for n = 392 states permodule and n = 1568 states in total with n DOF = 312 DOFs after concatenation of system matrices accordingto (41).Given the node and element table of the structure, the method described in Sec. 4.2 is used to automaticallyformulate coupling constraints. All translational DOFs of the lowermost 4 nodes in including the rotation aboutthe z -axis are locked. This is achieved by setting the corresponding velocity outputs to zero in (44). In total, n c = 488 constraints are formulated for the sample structure in the DAE form (50). After eliminating allalgebraic constraints using the procedure in Sec. 4.3, the system is reduced to n = 1080 ODEs. New systeminputs are the external forces f ext ∈ R acting on the structure’s DOFs with the collocated outputs being theglobal nodal velocities ˙ q . Proceeding by elimination of linearly dependent states according to the method inSec. 4.4, the number of ODEs is further reduced to n = 592, which equals twice the number of DOFs. Then, thetransformations from Sec. 4.5 are applied to the system equations to obtain the mass and stiffness matrices M fe and K fe . A Rayleigh-type damping term D fe = α M fe + α K fe is added as in (3), with the damping coefficientsset to α = 0 .
05 and α = 0 . x , which was calculated byassuming a stationary wind load acting in the positive x -direction. The wind load increases linearly with thestory number and amounts to f wind = 20 kN for the uppermost left nodes. The structural displacement of node52 (top of the structure) is shown in Fig. 6 a). Oscillation caused by the initial displacement is damped almostcompletely within five seconds. This is also reflected in the system energy shown in terms of the Hamiltonian H ( x ) in Fig. 6 b).For comparison, M fe and K fe were assembled according to a conventional FE method by adding up theelement mass and stiffness matrices after transforming them to global coordinates. Element mass and stiffnessmatrices for the Euler-Bernoulli beam can be found e. g. in [7]. The structural response of the system obtainedthe conventional way, when simulating it with the same initial condition, is almost identical to the one shownin Fig. 6. There are minor differences caused by the respective numerical methods involved in the calculationof M fe and K fe .A comparison between the first six eigenfrequencies of the structure obtained by the method presented inthis article and by the conventional FE approach is shown in Tab. 1. Except for the second eigenmode, theeigenmodes deviate slightly in the fourth decimal place. The first and second, as well as the fourth and fifthmodes, are bending modes in the x - and y -direction. Due to the symmetry of the structure, their eigenfrequenciesshould be identical, as observed in case of the conventional FE approach. The fact that they differ slightly forthe presented approach is again explained by the numerics involved. Mass and stiffness matrix are obtainedafter a number of transformations, while analytical formulations of the approximated element mass and stiffnessmatrices are used in conventional FEM. The only transformation involved in assembling the global mass andstiffness matrix in the latter case is the rotation of local coordinate frames which makes the method lesssusceptible to numerical errors. Thus, the higher number of transformations and numerical methods involved in15 S p T ms, w s A A p p x v ˙ q ˙ q L F ext (a) Double-acting hydraulic cylinder with valve xz (b) Steel frame with parallel hydraulic cylinder (red) Figure 7: Multi-Domain model of an adaptive structure that is actuated by a hydraulic cylinder. the presented framework can be viewed as a disadvantage. Since it is not strictly necessary, the transformationto global coordinates explained in Sec. 4.5 can be omitted to reduce numerical error. Especially when themechanical parts of the system are coupled with non-mechanical or hybrid components. An example of such amulti-domain system is presented in the following section.
In a recent publication of some of the authors, a method for the integration of actuators into FE modelsof adaptive structures is presented [30]. To allow for an easy integration of active elements, a quasi-stationaryassumption is made for the local dynamics of passive elements. While this is a valid simplification for mostlarge structures, the physical coupling between the dynamics of the actuators and those of the structure islost this way. If this is not desired, a significant increase in complexity results. In this section, we showthat it is straightforward to integrate hydraulic actuators into a mechanical structure within the presentedport-Hamiltonian modeling framework. Simplifications such as the ones in [30] need not be made.For the coupling of a mechanical structure with a hydraulic actuator, a port-Hamiltonian model of the latteris required. We use a slightly modified version of the model of a valve-controlled double-acting hydraulic cylinderpresented by Kugi and Kemmetmüller in [31] for which we will briefly present the necessary equations in thefollowing. For details, the reader is referred to [31]. A schematic representation of the piston actuator with athree-land-four-way valve is shown in Fig. 7 a). The volumetric flows into the cylinder chambers are given by˙ q = Γ x v = ( k v √ p S − p x v for x v ≥ ,k v √ p − p T x v for x v < , ˙ q = Γ x v = ( k v √ p − p T x v for x v ≥ ,k v √ p S − p x v for x v < x v and k v are the valve displacement and the valve coefficient, p S and p T are the supply and the tankpressure and p and p the respective chamber pressures. With the common simplification of a linearizedconstitutive law of the bulk modulus β , the system dynamics of the hydraulic cylinder can be formulated asfollows ˙ s = 1 m w s , ˙ w s = p A − p A − F ext , (68)˙ p = βA s ( − A w s + Γ x v ) , ˙ p = − βA ( L − s ) ( A w s − Γ x v ) . (69)Here, s is the piston position, m its mass and w s its impulse (in [31] the velocity). Multiplying the chamberpressures with the effective piston areas A and A and taking into account the external force F ext yields theforce balance equation in (68). For the system Hamiltonian we get H = U + 12 m w s , with U = A s (cid:18) β (cid:18) exp (cid:18) p β (cid:19) − (cid:19) − p (cid:19) + A ( L − s ) (cid:18) β (cid:18) exp (cid:18) p β (cid:19) − (cid:19) − p (cid:19) . (70)See [31] or [32] for the derivation of the internal energy U . Everything necessary for the formulation of theactuator dynamics as a nonlinear port-Hamiltonian system has now been defined. Choosing the state vector as x = (cid:2) s w s p p (cid:3) T and the input as u = (cid:2) x v F ext (cid:3) T , we get˙ x = J ( x ) ∂H∂ x ( x ) + g ( x ) u , y = g T ( x ) ∂H∂ x ( x ) , (71)16 0 . . . . d i s p l a ce m e n t i n mm n7 x n7 z n8 x n8 z (a) Displacements of the frame’s topmost nodes . . . . p i s t o nd i s p l a ce m e n t i n mm . . p r e ss u r e i nb a r ∆ sp p (b) Hydraulic cylinder dynamics Figure 8: Structural response of the system shown in Fig. 7 for sinusoidal excitation of the hydraulic cylinder where J ( x ) = − βs − βL − s − βs βL − s , g ( x ) = − βA s Γ βA ( L − s ) Γ . (72)The hydraulic cylinder is connected to the steel frame shown in Fig. 7 b) in parallel to the left vertical column.To simplify things a little, the frame is taken to be one side of the lowermost module (three stories) of the high-rise building shown in Fig. 5 a). It is composed of the same elements and the same parameters as in Sec. 5.2 areused for both approximation and damping. Coupling of the frame and the cylinder can be done just like beforeusing the procedure described in Sec. 4 until a system of the form (50) is obtained. Since the coupled systemis now nonlinear, we cannot proceed with the elimination of algebraic constraints in an automated fashion.However, it is still possible to obtain an explicit formulation by solving for the Lagrange multipliers (see e. g.[33] and [3]) λ ( x , u ) = − (cid:18) B T ( x ) ∂ H∂ x ( x ) B ( x ) (cid:19) − B ( x ) ∂ H∂ x ( x ) (cid:18) J ( x ) ∂H∂ x ( x ) + g ( x ) u (cid:19) . (73)A simulation of the coupled multi-domain system is shown in Fig. 8, where the valve displacement x v followsa sinusoidal excitation. For the cylinder, the length was set to L = 35 cm, the effective piston areas were chosenas A = 133 cm and A = 94 . and the piston mass as m = 10 kg. The system is driven with a supplypressure of p S = 200 bar and a tank pressure of p T = 0 bar. The initial pressures in the chambers were set to p = p = 1 bar and the initial piston position was set to s = 10 cm. A bulk modulus of β = 1 . x v = 2 l / min were assumed.In Fig. 8 a), the x - and z -displacements of the topmost nodes of the frame are shown, whereas Fig. 8 b)depicts the chamber pressures and the piston displacement ∆ s . The coupling between piston displacement andthe motion of the frame is immediately visible. Comparing the z -displacement of node 7 in Fig. 7 a) with ∆ s in Fig. 7 b), especially at the peaks at t = 0 . t = 1 .
6. Conclusion and Further Work
In this contribution, a selection of methods for the modeling of the dynamics of complex structures as linearISO port-Hamiltonian systems was presented. Starting from the distributed parameter formulations of axial,torsional and bending motion of a beam element, the corresponding infinite-dimensional port-Hamiltonian sys-tems were derived. In order to enable numerical simulations of these systems, finite-dimensional approximationsare necessary. Those were obtained via a structure preserving mixed finite element method. A numerical analy-sis of the approximation error was presented for the beam equations. Given the spatially discretized models forindividual load cases, complex structures can be assembled by interconnecting them via coupling constraints.This leads to index one DAEs that can be reduced to an ODE system on a constraint state space by means of atransformation. Afterward, further reduction of the system order is achieved by eliminating linearly dependent17ffort variables. The resulting system coordinates can be transformed to the space of global DOFs, which allowscomparison of the system’s mass and stiffness matrices to the ones obtained using conventional FE methods.The process was presented for a multi-story high-rise building.When a system is only composed of mechanical components, the mass and stiffness matrices from a con-ventional FE model can also be used to formulate a port-Hamiltonian system. In fact, this is advantageous,since the number of transformations and numerical routines involved in the presented approach is higher, whichresults in a higher numerical error. If the aim is, however, to model heterogeneous multi-domain systems asport-Hamiltonian systems, the merit of the presented method is, that non-mechanical elements can be includedin each step of the process. Since the port-Hamiltonian structure is preserved in each step, this makes it easy toadd such elements on both the PDE and ODE level. We demonstrated the coupling of the hydraulic with themechanical domain which is straightforward with the presented methods. It also shows that the methodologyis not limited to linear systems.A Matlab framework which encompasses all the models and methods presented in this contribution isavailable online. Not only can it be used to reproduce the results in Sec. 5, but also to build arbitrary structurescomposed of the elements introduced in Sec. 2. The framework is written in an object-oriented way, such thatthe user only needs to provide the geometry and material parameters of the elements including a node andelement table. Different examples are provided which illustrate the usage of the methods. Extension of the codeto account for e. g. different types of damping or model order reduction is a matter of ongoing work. Also anoptimization with respect to numerical stability and precision.Currently, coupling of simple elements is based on a block diagram analogy. In further work, we aim torepresent component interconnection in bond graphs which is more in line with port-Hamiltonian systems theoryand facilitates interfacing with non-mechanical components and comprehension of system characteristics. Tomodel a wider range of adaptive structures, further subsystems such as plate and shell elements or integratedfluidic actuators need to be included in the framework. Furthermore, we also want to develop control and stateestimation methods based on the models presented in this article.
Acknowledgment
The authors gratefully acknowledge the generous funding of this work by the German Research Foundation(DFG - Deutsche Forschungsgemeinschaft) as part of the Collaborative Research Center 1244 (SFB) "AdaptiveSkins and Structures for the Built Environment of Tomorrow"/projects A06 and B02.
Appendix A. Intermediate Steps in Discretization
Since it is not immediately clear how to obtain (30) from (29) using integration by parts, some intermediatesteps are shown here to enhance comprehensibility. We begin with the term including the first order spatialderivative for which we obtain Z Z φ p A ∂∂z φ T q e q d z = − Z Z ∂∂z φ p A φ T q e q d z + h φ p A φ T q e q i ba (A.1)and are finished in case our system has no higher-order spatial derivatives. For the second order spatialderivative, integration by parts needs to be applied twice. The term Z Z φ p A ∂ ∂z φ T q e q d z (A.2)is omitted on the left hand side in the following due to limited space and we proceed with . . . = − Z Z ∂∂z φ p A ∂∂z φ T q e q d z + (cid:20) φ p A ∂∂z φ T q e q (cid:21) ba = Z Z ∂ ∂z φ p A φ T q e q d z + (cid:20) φ p A ∂∂z φ T q e q (cid:21) ba − (cid:20) ∂∂z φ p A φ T q e q (cid:21) ba . (A.3)Accordingly, integration by parts needs to be repeated N times to Z Z φ p A N ∂ N ∂z N φ T q e q d z (A.4)18hich yields · · · = − Z Z ∂∂z φ p A N ∂ N − ∂z N − φ T q e q d z + (cid:20) φ p A N ∂ N − ∂z N − φ T q e q (cid:21) ba = . . . =( − N Z Z ∂ N ∂z N φ p A N φ T q e q d z + (cid:20) φ p A N ∂ N − ∂z N − φ T q e q (cid:21) ba − . . . + ( − N − (cid:20) ∂ N − ∂z N − φ p A N φ T q e q (cid:21) ba . (A.5)Merging (A.1), (A.3) and (A.5) and combining the bracketed integrals in sums respectively, it can be seen thatthe result of applying integration by parts N times is indeed given by (30). References [1] S. Korkmaz, A Review of Active Structural Control: Challenges for Engineering Informatics, Computers & Structures 89 (23)(2011) 2113–2132.[2] W. Sobek, Ultraleichtbau, Stahlbau 83 (11) (2014) 784–789. doi:10.1002/stab.201410211 .URL https://doi.org/10.1002/stab.201410211 [3] A. van der Schaft, D. Jeltsema, et al., Port-Hamiltonian systems theory: An introductory overview, Foundations and Trends R (cid:13) in Systems and Control 1 (2-3) (2014) 173–378.[4] H. Zwart, B. Jacob, Distributed-parameter port-Hamiltonian systems, CIMPA.[5] A. J. van der Schaft, B. M. Maschke, Hamiltonian formulation of distributed-parameter systems with boundary energy flow,Journal of Geometry and Physics 42 (1-2) (2002) 166–194.[6] A. Quarteroni, A. Valli, Numerical approximation of partial differential equations, Vol. 23, Springer Science & Business Media,2008.[7] R. Schwertassek, O. Wallrapp, Dynamik flexibler Mehrkörpersysteme: Methoden der Mechanik zum rechnergestützten Entwurfund zur Analyse mechatronischer Systeme, Springer-Verlag, 2017.[8] T. J. R. Hughes, The finite element method: linear static and dynamic finite element analysis, Courier Corporation, 2012.[9] V. Talasila, G. Golo, A. J. van der Schaft, The wave equation as a port-Hamiltonian system and a finite dimensional approxi-mation, in: Proceedings of 15th international symposium mathematical theory of networks and systems (MTNS), South Bend,2002.[10] G. Golo, V. Talasila, A. Van Der Schaft, B. Maschke, Hamiltonian discretization of boundary control systems, Automatica40 (5) (2004) 757–771.[11] L. Bassi, A. Macchelli, C. Melchiorri, An algorithm to discretize one-dimensional distributed port-Hamiltonian systems, in:Lagrangian and Hamiltonian methods for nonlinear control 2006, Springer, 2007, pp. 61–73.[12] A. Baaiu, F. Couenne, L. Lefevre, Y. Le Gorrec, M. Tayakout, Structure-preserving infinite dimensional model reduction:Application to adsorption processes, Journal of Process Control 19 (3) (2009) 394–404.[13] R. Moulla, L. Lefevre, B. Maschke, Pseudo-spectral methods for the spatial symplectic reduction of open systems of conserva-tion laws, Journal of computational Physics 231 (4) (2012) 1272–1292.[14] F. L. Cardoso-Ribeiro, D. Matignon, V. Pommier-Budinger, Piezoelectric beam with distributed control ports: a power-preserving discretization using weak formulation., IFAC-PapersOnLine 49 (8) (2016) 290–297.[15] P. Kotyczka, B. Maschke, L. Lefèvre, Weak form of Stokes–Dirac structures and geometric discretization of port-Hamiltoniansystems, Journal of Computational Physics 361 (2018) 442–476.[16] F. L. Cardoso-Ribeiro, D. Matignon, L. Lefevre, A structure-preserving Partitioned Finite Element Method for the 2D waveequation, IFAC-PapersOnLine 51 (3) (2018) 119–124.[17] A. Brugnoli, D. Alazard, V. Pommier-Budinger, D. Matignon, Port-Hamiltonian formulation and symplectic discretization ofplate models Part I: Mindlin model for thick plates, Applied Mathematical Modelling 75 (2019) 940–960.[18] A. Brugnoli, D. Alazard, V. Pommier-Budinger, D. Matignon, Port-Hamiltonian formulation and symplectic discretization ofplate models Part II: Kirchhoff model for thin plates, Applied Mathematical Modelling 75 (2019) 961–981.[19] R. Li, Y. Tian, P. Wang, Y. Shi, B. Wang, New analytic free vibration solutions of rectangular thin plates resting on multiplepoint supports, International Journal of Mechanical Sciences 110 (2016) 53–61.[20] R. Li, X. Zheng, H. Wang, S. Xiong, K. Yan, P. Li, New analytic buckling solutions of rectangular thin plates with all edgesfree, International Journal of Mechanical Sciences 144 (2018) 67–73.[21] R. Li, X. Zheng, Y. Yang, M. Huang, X. Huang, Hamiltonian system-based new analytic free vibration solutions of cylindricalshell panels, Applied Mathematical Modelling 76 (2019) 900–917.[22] Y. Le Gorrec, H. Zwart, B. Maschke, Dirac structures and boundary control systems associated with skew-symmetric differentialoperators, SIAM journal on control and optimization 44 (5) (2005) 1864–1892.[23] A. Macchelli, C. Melchiorri, Modeling and control of the Timoshenko beam. The distributed port-Hamiltonian approach, SIAMJournal on Control and Optimization 43 (2) (2004) 743–767.[24] F. L. Cardoso-Ribeiro, D. Matignon, V. Pommier-Budinger, Modeling of a fluid-structure coupled system using port-Hamiltonian formulation, IFAC-PapersOnLine 48 (13) (2015) 217–222.[25] F. L. Cardoso-Ribeiro, D. Matignon, L. Lefèvre, A Partitioned Finite Element Method for power-preserving discretization ofopen systems of conservation laws, arXiv preprint arXiv:1906.05965.[26] W. Schnell, D. Gross, W. Hauger, P. Wriggers, Technische Mechanik: Band 4: Hydromechanik, Elemente der HöherenMechanik, Numerische Methoden, Springer-Verlag, 2006.[27] A. J. Van der Schaft, Port-Hamiltonian differential-algebraic systems, in: Surveys in Differential-Algebraic Equations I,Springer, 2013, pp. 173–226.[28] F.-L. Cardoso Ribeiro, Port-Hamiltonian modeling and control of a fluid-structure system, Ph.D. thesis, University of Toulouse(12 2016).
29] Y. Wu, B. Hamroun, Y. Le Gorrec, B. Maschke, Port-Hamiltonian system in descriptor form for balanced reduction: Appli-cation to a nanotweezer, IFAC Proceedings Volumes 47 (3) (2014) 11404–11409.[30] M. Böhm, J. L. Wagner, S. Steffen, J. Gade, F. Geiger, W. Sober, M. Bischoff, O. Sawodny, Input modeling for active structuralelements – extending the established FE-workflow for modeling of adaptive structures, in: 2020 IEEE International Conferenceon Advanced Intelligent Mechatronics (AIM), 2020, (accepted for publication).[31] A. Kugi, W. Kemmetmüller, New energy-based nonlinear controller for hydraulic piston actuators, European journal of control10 (2) (2004) 163–173.[32] A. Kugi, W. Kemmetmüller, Energy based modelling of lumped-parameter hydraulic systems, in: 4th Mathmod, 2003, pp.825–834.[33] V. Duindam, A. Macchelli, S. Stramigioli, H. Bruyninckx, Modeling and control of complex physical systems: the port-Hamiltonian approach, Springer Science & Business Media, 2009.29] Y. Wu, B. Hamroun, Y. Le Gorrec, B. Maschke, Port-Hamiltonian system in descriptor form for balanced reduction: Appli-cation to a nanotweezer, IFAC Proceedings Volumes 47 (3) (2014) 11404–11409.[30] M. Böhm, J. L. Wagner, S. Steffen, J. Gade, F. Geiger, W. Sober, M. Bischoff, O. Sawodny, Input modeling for active structuralelements – extending the established FE-workflow for modeling of adaptive structures, in: 2020 IEEE International Conferenceon Advanced Intelligent Mechatronics (AIM), 2020, (accepted for publication).[31] A. Kugi, W. Kemmetmüller, New energy-based nonlinear controller for hydraulic piston actuators, European journal of control10 (2) (2004) 163–173.[32] A. Kugi, W. Kemmetmüller, Energy based modelling of lumped-parameter hydraulic systems, in: 4th Mathmod, 2003, pp.825–834.[33] V. Duindam, A. Macchelli, S. Stramigioli, H. Bruyninckx, Modeling and control of complex physical systems: the port-Hamiltonian approach, Springer Science & Business Media, 2009.