TL;DR
This paper introduces a structure-preserving finite element discretization method for hyperbolic conservation law systems, ensuring power preservation and applicable to complex geometries and higher-order systems.
Contribution
A novel partitioned finite element method for port-Hamiltonian systems that maintains power conservation across diverse geometries and extends to higher-order models.
Findings
Successfully applied to nonlinear 2D Shallow Water Equation
Preserves power and structure in discretization
Extensible to curvilinear coordinates and higher-order systems
Abstract
This paper presents a structure-preserving spatial discretization method for distributed parameter port-Hamiltonian systems. The class of considered systems are hyperbolic systems of two conservation laws in arbitrary spatial dimension and geometries. For these systems, a partioned finite element method is derived, based on the integration by parts of one of the two conservation laws written in weak form. The nonlinear 1D Shallow Water Equation (SWE) is first considered as a motivation example. Then the method is investigated on the example of the nonlinear 2D SWE. Complete derivation of the reduced finite-dimensional port-Hamiltonian system is provided and numerical experiments are performed. Extensions to curvilinear (polar) coordinate systems, space-varying coefficients and higher-order port-Hamiltonian systems (Euler-Bernoulli beam equation) are provided.
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Code & Models
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
A Partitioned Finite Element Method for power-preserving
discretization of open systems of conservation laws
Flávio Luiz Cardoso-Ribeiro
Instituto Tecnológico de Aeronáutica, Brazil
Denis Matignon
ISAE-SUPAERO, Université de Toulouse, France
Laurent Lefèvre
Univ. Grenoble Alpes, Grenoble INP, LCIS, F-26000, Valence
Abstract
This paper presents a structure-preserving spatial discretization method for distributed parameter port-Hamiltonian systems. The class of considered systems are hyperbolic systems of two conservation laws in arbitrary spatial dimension and geometries. For these systems, a partioned finite element method is derived, based on the integration by parts of one of the two conservation laws written in weak form. The nonlinear 1D Shallow Water Equation (SWE) is first considered as a motivation example. Then the method is investigated on the example of the nonlinear 2D SWE. Complete derivation of the reduced finite-dimensional port-Hamiltonian system is provided and numerical experiments are performed. Extensions to curvilinear (polar) coordinate systems, space-varying coefficients and higher-order port-Hamiltonian systems (Euler-Bernoulli beam equation) are provided.
keywords:
geometric spatial discretization , structure-preserving discretization , port-Hamiltonian systems , partitioned finite element method.
††journal: arXiv
1 Introduction
The port-Hamiltonian formalism has been proven to be a powerful tool for the modeling and control of complex multiphysics systems. In many cases, spatio-temporal dynamics must be considered and infinite-dimensional port-Hamiltonian models are needed. Classical academic examples such as the transmission line, the shallow water or the beam equations have been investigated in the port-Hamiltonian framework [1].
Besides 2D and 3D problems have been recently considered [2, 3, 4]. In many of these examples, e.g. those arising from mechanics, systems of two balance equations are considered such as mass and momentum or volume and momentum balance equations.
In order to simulate and design control laws, obtaining a finite-dimensional approximation which preserves the port-Hamiltonian structure of the original system can be advantageous. It may serve as a design guide such as in Control by Interconnection (CbI) or in Interconnection and Damping Assignment Passivity Based Control (IDA-PBC). Besides, preserving the underlying Dirac interconnection structure results in energy conservation properties and associated dynamical properties (e.g. stability, controllability, etc.).
Mixed finite element methods were introduced a long time ago to perform structure-preserving spatial discretization of the Maxwell field equations [5, 6]. An extension to open port-Hamiltonian systems was presented in [7] where the authors proposed a mixed finite element structure-preserving spatial discretization for 1D hyperbolic systems of conservation laws, making use of disctinct low-order Whitney bases functions to approximate respectively the energy and co-energy variables. These ideas were applied later for the discretization of a parabolic diffusion problem related to pressure-swing-adsorption columns [8], piezo-electric beams [9] and connected to finite volume and staggered grid finite difference methods for 1D problems [10]. They were also generalized for 2D systems in [2] where they are applied to a vibro-accoustic systems and in [4] where connection was made with 2D finite difference staggered grids schemes. Finally, these structure-preserving mixed finite element methods, applied to the spatial discretization of general port-Hamiltonian systems with boundary energy flows, were stated in a geometry independent form making use of discrete exterior calculus results [11].
In these previous works, the central idea was to define different discretization bases for the energy and co-energy variables such that the strong-form of the equations was exactly satisfied in the corresponding spanned finite-dimensional approximation spaces. This idea was extended to geometric pseudo-spectral methods using conjugated high-order polynomial bases ([12]) or Bessel ([13]) bases functions, globally defined on the whole spatial domain. Still the same idea has been considered in [14, 15, 16] for the 1D transmission line equation and for the Maxwell equations. In these latter works, one of the balance law is kept in strong-form (with exact spatial derivation) while the other one is considered in the weak sense only. As it was noticed in [17], defining these compatible spaces - with power-conjugated approximation bases for the energy and co-energy variables - is (relatively) straightforward for 1D systems, but seems to be cumbersome for higher spatial dimensions or higher order methods. Indeed, the kernel of the exterior derivatives in -D dimensional domains is not anymore trivial and the discretization of the trace operator on boundaries with non trivial (i.e. non rectangular) geometries often leads to dimensionality problems. As suggested in [18] the discretization of the weak formulation of the considered port-Hamiltonian system may be a practical solution to deal with these higher dimensional problems or more complex geometries. We propose in this paper to follow this approach but to perform integration by parts - which was used in [18] to get the weak formulation - on one of the two balance equations only, defining in this way a partitioned mixed finite element method. Doing so, the discretization in the chosen bases for the energy and co-energy variables (and the associated test functions) directly leads to a finite-dimensional Dirac interconnection structure and no further projection is required to get a finite-dimensional port-Hamiltonian system equations with reversible causality. Besides boundary conditions are naturally handled, even in the case of higher-order finite element bases. Finally, the use of this weak-form formulation enables the use of standard finite-element software to perform the proposed structure-preserving discretization and consequently paves the way for further applications to more involved higher-order problems with complex geometries.
This paper starts with a motivation example detailed in Section 2: a structure-preserving spatial discretization for the one-dimensional (1D) Shallow Water Equations (SWE). In Section 3, the approach is generalized in a general setting, where the initial model is stated in a coordinate-free form, independent of the specific geometry. The covariant formulation for systems of two conservation laws with boundary energy flows is presented, as well as its weak form. In Section 4 the proposed partitioned finite element method (PFEM) is applied to the 2D SWE example. Numerical experiments are presented in Section 5. Finally, extensions of the method to curvilinear (polar) coordinate systems, space-varying coefficients and higher-order port-Hamiltonian systems (Euler-Bernouilli beam equations) are provided in Section 6. The paper ends with conclusions and open questions which are discussed in Section 7.
2 An introductory example
The aim of this section is to present the general idea of this paper - partial integration by parts of the weak form for systems of two conservation laws and structure-preserving projections in finite-dimensional approximation spaces - applied on a simple one-dimensional example, namely the 1D SWE written in the port-Hamiltonian formulation. First, the port-Hamiltonian formulation for these equations is recalled (subsection 2.1). Then the partial integration by parts idea is performed on the weak form for this 1D SWE example and a structure-preserving finite element spatial discretization method is applied to obtain the finite-dimensional port-Hamiltonian system (pHs) (subsection 2.2). This general idea differs from previous works (as [7] and [12]) where the central idea was to define different discretization spaces for the energy and co-energy variables such that the strong form of the equations were exactly satisfied in these finite-dimensional spaces. Instead, we use a weak-form representation for the equations, where only one of the conservation laws is integrated by parts. This partitioned approach naturally leads to a skew-symmetric interconnection matrix between the energy and co-energy variables. Furthermore, the use of weak form enables the use of classical finite-element methods to perform the discretization.
2.1 Port-Hamiltonian strong formulation for the 1D SWE
The Shallow Water Equations are a set of partial differential equations that can be used to represent an incompressible fluid with free-surface motion. These equations are typically used to model fluid motion in water channels [19], wave propagations in oceans and lakes, and sloshing in fluid tanks [20, 21]. When one considers the frictionless flow in a horizontal channel with uniform rectangular cross-section, the one-dimensional mass and momentum balance equations may be written as:
[TABLE]
where is the fluid height, , the fluid average velocity in a cross-section, , the spatial coordinate, , the time and , the gravitational acceleration.
The total energy of the system inside the domain is given by the sum of kinetic and potential (gravitational) energy:
[TABLE]
where is the width of the water channel (or tank) rectangular cross-section and , the water density. Defining the energy-variables and , the system Hamiltonian (total energy) is given by:
[TABLE]
Using these newly defined variables, (1) can be rewritten as:
[TABLE]
where and are the co-energy variables (respectively, the total pressure and the water flow) which are defined as the variational derivatives of the Hamiltonian with respect to and :
[TABLE]
From the above definitions of energy and co-energy variables, using the SWE written in the canonical Hamiltonian form (4) and Stokes theorem, one obtains for the power balance equation:
[TABLE]
where boundary port input variables, , are defined as the values of the co-energy variables evaluated at the spatial domain boundary:
[TABLE]
while the power-conjugated boundary output variables are defined as:
[TABLE]
The power balance equation (6) defines a natural pairing or bilinear form
[TABLE]
where the Bond space is defined as the product of the effort real vector space
[TABLE]
and its dual flow real vector space
[TABLE]
with and denoting respectively the Sobolev space of functions with square integrable derivatives on and the usual Lebesgue space of square integrable functions on . Using the bilinear form (9), the power balance equation (6) simply reads
[TABLE]
with and (or the reverse). Besides, the pairing (9) may be symmetrized to obtain the associated inner product:
[TABLE]
It may be shown that the Hamiltonian formulation (4) for the SWE, together with the boundary conditions (7) and output (8) may be equivalently implicitely defined as
[TABLE]
where is the linear subspace which is maximally isotropic (i.e. ) with respect to the inner product (13) [12]. In that sense, the natural pairing (9) fully describes the geometric structure of the port-Hamiltonian system (4) with boundary values (7). Therefore, in this paper, structure-preserving (or symplectic) spatial discretization will be understood as approximations (projections) which preserve this power form (9). Symplecticity in that sense not only implies preservation of the power balance (6) or (12) (i.e. isotropy) but preservation of the whole geometric structure of the system (e.g. the Poisson structure in the example of closed systems or the Dirac structure in the case of open systems with time varying boundary conditions) [12, 18].
Note that the particular input and output port variables chosen here above in (equations (7) and (8)) is only one among other possible choices. A description of all the possible choices of input/output variables which lead to well-posed problems (in the linear case) is described in [22].
2.2 Partitioned weak form and structure-preserving discretization for the 1D SWE
We will now introduce a weak formulation for the 1D SWE and then perform integration by parts on the mass balance equation. Let and denote any arbitrary test functions, we may obtain from the strong formulation (4) the following weak form:
[TABLE]
Integrating by part the mass balance equation only, we get the following partitioned weak form:
[TABLE]
Remark 1**.**
In the specific case where and , we get:
[TABLE]
which shows that the two conservation laws for the total mass and the total momentum in the spatial domain are preserved in the weak formulation. When and are chosen, one gets:
[TABLE]
Therefore the power-balance equation (6) reads:
[TABLE]
which shows that the power balance is also preserved in the weak formulation.
We will now project the partitioned weak formulation (16) into finite-dimensional approximation spaces chosen in such a way as to preserve the total mass and momentum conservation laws, the power balance equation and the underlying Dirac structure of the original port-Hamiltonian model (4). Unlike in [12, 18] where different approximation bases are chosen for the energy and co-energy variables, we obtain the mass, momentum, power and structure-preservation by the selection of different approximation bases for the mass and momentum densities. This “partitioned” choice for the approximation bases lead us to square skew-symmetric interconnection matrices.
Let us approximate the energy variables and as
[TABLE]
where are the chosen approximation basis functions in , the chosen approximation basis functions in , while and are the approximation coordinates for and in the approximation bases and . The test functions and are approximated in the same bases as and , respectively. From the substitution of the approximated variables (20) in the weak form (16), the following finite-dimensional equations are obtained:
[TABLE]
where the effort functions have been approximated in the basis, while the effort functions have been approximated in the basis. Similarly, the flow functions and have been approximated respectively in the and bases. Since these equations should remain valid for any choices of test functions coordinates and , one gets:
[TABLE]
where and are square mass matrices (of size and , respectively) defined as
[TABLE]
Matrix is of size , and defined as
[TABLE]
and is an matrix. Using the input-output conjugated boundary port variables as defined in (7) and (8), since the boundary values of may be written
[TABLE]
the approximation (22) may be written using the following finite-dimensional Dirac structure representation:
[TABLE]
where and denote the vector coordinates for the flow approximations respectively in the in the and approximation bases, that is and .
We obtained the finite-dimensional Dirac structure representation (26) from the projection of the Stokes-Dirac structure on a the chosen approximation spaces. We will now derive the corresponding approximation of the port-Hamiltonian system dynamics (14) by restricting the Hamiltonian functional to the same approximation spaces. From the definition of co-energy variables as variational derivatives of the Hamiltonian with respect to and , we get:
[TABLE]
Using the approximations (20) for the energy and co-energy variables, this power balance may be approximated as:
[TABLE]
where new energy variables
[TABLE]
have been defined. Therefore, in order to write the power balance (28) as the total time derivative of the discrete Hamiltonian written as a function of the finite-dimensional vector coordinates for the energy variables, the following reduced co-energy variables must be defined:
[TABLE]
where
[TABLE]
We obtain the finite-dimensional pHs formulation for the proposed structure-preserving reduction scheme by combining equations (26) (for the linear finite-dimensional Dirac interconnection structure) and the nonlinear constitutive equations (30) and (31):
[TABLE]
The state space model (32) for the reduced dynamics exhibits the usual general pHs form for finite-dimensional systems [23]. It is a generic formulation which may be obtained for any 1D systems of conservation laws written in the port-Hamiltonian formulation. In the 1D SWE example, the Hamiltonian function is neither quadratic, nor separable. Nevertheless, an explicit form may be obtained for the constitutive equations (30). Since the Hamiltonian function restricted to the approximation spaces for and reads:
[TABLE]
One obtains for the reduced effort variables the expressions:
[TABLE]
Note that both constitutive equations exhibits nonlinear terms. In order to compute them, the following procedured was used. The first equation can be written as:
[TABLE]
where the components of can be computed as:
[TABLE]
Note that , for are matrices of dimension , that can be computed once and remains constant. Similarly, the second constitutive relationship can be written as a function of constant matrices.
Remark 2**.**
We may deduce from the pHs representation (32) that the power balance equation reads:
[TABLE]
Hence, the power balance (and the corresponding power product value) is preserved by the proposed partitioned spatial discretization scheme. In that sense, we call it a structure-preserving or symplectic scheme. When the bases functions satisfy
[TABLE]
the mass and momentum conservation laws are also satisfied in the finite-dimensional approximation spaces:
[TABLE]
In the beginning of this section, we motivated this work by the fact that previous work on structure-preserving spatial discretization that relies on exact satisfaction of the strong form of the equations usually lead to difficulties when generalizing to 2D or 3D systems. The following questions arise: can the proposed PFEM method be easily generalized to higher-dimensional problems (2D and 3D)? Does it work with different coordinate systems? What about convergence? We will answer these questions in the following sections.
3 A general setting
In this section we will define a port-Hamiltonian system of two conservation laws with time-varying port boundary variables in the so-called co-variant form, that is defining energy and co-energy variables as differential forms [24]. This will allow us to define, independently from the particular spatial dimension, geometry or coordinate system, the class of problems which can be solved by using the structure-preserving spatial discretization scheme proposed in this paper. In the next subsection (§ 3.1), to make the paper self-contained, we give a short introduction to exterior calculus with differential forms and their functional spaces. Then, in § 3.2, we generalize the definition of port-Hamiltonian systems, given in the previous section for the 1D SWE example, for general systems of two conservation laws, using these differential forms. In § 3.4, we define de Stokes-Dirac interconnection structures associated to these port-Hamiltonian systems and the corresponding power pairings which will be preserved in the discretization. We also detail the 2D SWE which will be used later in the numerical experiments. Finally, in § 3.5, we will give the general weak form which will be used for the structure-preserving discretization scheme presented in section 4. Readers could refer to [25, chapter 5], [26] or [27] for a general intuitive introduction to differential forms and their use in physical systems modelling. Numerical approximations using finite element spaces of differential forms are described in [28, 29].
3.1 Differential forms and exterior calculus
Let be an open, bounded and connected -dimensional spatial domain (manifold) with a -dimensional Lipschitz boundary .
We denote the space of smooth differential -forms (i.e. field of alternated -forms defined on the tangent subspace ) with smooth coefficients in .
Differential -forms are endowed with a product, called exterior product, denoted as , used for instance to define the power as the product of -forms flow variables (time derivatives of energy variables) and -forms co-energy variables (variational derivatives of the Hamiltonian with respect to the energy variables). The exterior product (also called wedge product or Grassman product) is a skew symmetric exterior product such that
[TABLE]
for all -form , -form and with .
Differential forms are also endowed with a derivation which is called exterior derivation, denoted by d, which generalize the gradient, divergence and curl operators respectively for [math]-forms (i.e. functions), -forms and -forms in . The exterior derivative (or co-boundary map) is a derivation of degree one
[TABLE]
such that
[TABLE]
for all and . The exterior derivative is defined such that the Stokes theorem applies, that is
[TABLE]
where the trace denotes the continuous extension of to the boundary . The Stokes theorem is the central result which allows us to define, from the natural power pairing, the associated Stokes-Dirac interconnection structure (see § 3.4). Symplectic discretization schemes will be defined in order to preserve either the power product or equivalently the Stokes theorem in finite-dimensional approximation spaces. In the sequel, we will make use of the duality product (natural power pairing) in between -forms and -forms :
[TABLE]
In this duality product (power pairing), one of the argument (differential form) will denote a flow variable while the other will be an effort variable. These two variables will be related together, on one side by balance equations (conservation laws) which are structural (metric independent) and on the other side by constitutive equations. The constitutive equations are related to the physical properties of the considered material domain (magnetic, visco-elastic, etc.). Therefore, they are metric dependent and will be formulated using the so-called Hodge star operator. For instance, in Section 3.3, the constitutive equations (60) or (61) are derived from the Hamiltonian density which describes, for the 2D SWE, the surfacic potential and kinetic energy densities. These constituve equations relates pressure tensors and water flows with mass and momentum densities. Therefore, metric dependent relations are needed and the Hodge star operator is introduced. The Hodge star operator induces an inner product on the space of differential -forms on the -dimensional manifold from the duality product:
[TABLE]
This chosen product is not necessarily the standard inner product for the norm. For instance, non uniform metric are required when the material domain properties are not uniform in space or isotropic. However, unless otherwise stated, we will consider in this paper the usual standard (metric). For instance, using Cartesian coordinates in some 2-dimensional domain , we will get , and .
The generalized Stokes’ theorem (44), together with the product rule (42), gives the following integration by parts formula for the power pairing:
[TABLE]
This formula will be helpful, in Section 3.5, to obtain the generalized weak port-Hamiltonian formulation for systems of two conservation laws with boundary energy flows.
When dealing with finite elements approximation spaces, one has to pay some attention to the regularity of differential forms if we want to apply the previous formula. Readers can refer to [29, section 4] for an introduction to differential forms whose coefficients are in Lebesgue spaces and Sobolev spaces . We will here point out only two comments on the weak exterior derivative and the trace theorem for differential forms.
Let us consider a -form , that is a -form with coefficients in . Then, its weak exterior derivative may be defined, using the integration by parts formula (47), as the unique -form (when it exists) such that:
[TABLE]
for all -form , that is a -form with coefficients in the space of continuous functions with compact support in . In the sequel of the paper, when dealing with the weak formulation of port-Hamiltonian system, exterior derivatives will be understood in this weak sense.
As far as we are dealing with systems of conservation laws with boundary energy flows, boundary values of the co-energy variables will be necessary to define input/output pairs of power conjugated boundary port variables (see section 3.2 hereafter). Hopefully, we can make use of the so-called trace theorem which has been extended to Sobolev spaces of differential forms [29, section 4]:
[TABLE]
The injection in equation (49) is continuous when Lebesgue space and the Sobolev space are equipped with the usual inner products. Throughout the paper, we will make use of the following compact notation:
[TABLE]
for all in .
3.2 The covariant port-Hamiltonian formulation for systems of two conservation laws with boundary energy flows
We will now extend the port-Hamiltonian formulation which has been presented in Section 2.1 only for the 1D SWE example. Hyperbolic systems of two conservation laws will be stated using differential forms, as introduced in the previous section, both for the energy and co-energy variables. Let us consider the two conserved quantities and with and denotes the dimension of the open, connected domain with Lipschitz boundary . These variables and are vector-valued distributed energy state variables defined for any ( is the position vector) and time .
Let the Hamiltonian functional be defined as:
[TABLE]
where denotes the Hamiltonian density -form which is assumed to be a smooth function. The variational derivatives of with respect to and are the unique differential and forms, denoted respectively and , such that:
[TABLE]
Therefore, from the Hamiltonian defined in (51), we may define the co-energy variables (efforts):
[TABLE]
with and . The Hamiltonian system of two canonically interacting conservation laws for and may be defined as:
[TABLE]
where the exponent ensures the formal skew symmetry of the matrix-valued differential operator (that is skew symmetry assuming zero boundary conditions for the arguments). According to (52), the time derivative of the energy functional (power balance) reads:
[TABLE]
According the state equations (54) and to the integration by part formula (47), this power balance may be written:
[TABLE]
This latter formula suggests of the following boundary port variables:
[TABLE]
defined in such a way that the power balance equation (55) and (56) may be written:
[TABLE]
where the flows variables and are defined as
[TABLE]
The power balance equation (58) simply states that the time derivative of the energy increase inside the domain equals the power supplied through the boundary . As it may be seen, we have extended the systems of two conservation laws (54) with a boundary power supply and related boundary port variables, obtaining an open system of two conservation laws with boundary energy flows. The explicit definition of an open port-Hamiltonian system of two canonically interacting conservation laws is given by the distributed state equation (54) together with the definition of the boundary port variables (57). It leads to the structural power balance equation (58) which is independent of the specific considered Hamiltonian function (i.e. from the effort constitutive equations (52)). Many 1D, 2D or 3D examples, either linear or nonlinear may be recast in that framework and satisfy this definition [23, 1]. Even parabolic systems may be formally represented with the skew-symmetric operator with a appropriate definition of the effort variables [3].
3.3 The irrotational 2D Shallow Water Equation example
We will consider as a running example for this paper the two-dimensional irrotational Shallow Water Equations (2D SWE) which describe the flow of an inviscid liquid where the horizontal components of the velocity field may be averaged on the water level and where the vertical velocity component may be omitted (low depth or shallow water assumption). Besides, we will consider a “non rotating” flow. It is known that the corresponding 2D SWE express then the mass and momentum balance equations. Therefore, we will choose for the energy state variables the mass density (which is a 2-form proportional to the water level ) and the momentum density (which is a 1-form). For instance, using Cartesian coordinates, one would choose (where denotes the water level) and where and denote the horizontal components of the fluid velocity while denotes the fluid mass density. According to the previous notations, is a 2-form and a 1-form, both defined in the 2-dimensional () horizontal spatial domain of the flow. Using these energy state variables, one gets for the total (kinetic and potential) energy111Since the Hamiltonian density and the resulting constitutive equations relate differential forms of different degrees and are metric dependent, the Hodge star is introduced. We will consider in this paper the usual uniform metric. inside the domain :
[TABLE]
Therefore, the co-energy variables are defined as:
[TABLE]
which are respectively the hydrodynamic pressure (a 0-form or intensive variable in the 2D domain ) and the volume flow (which is indeed a 1-form in the 2D spatial domain). For instance, using the same Cartesian coordinates as previously, one gets and . Using these co-energy variables, the Hamiltonian system of two canonically interacting conservation laws (54) reads:
[TABLE]
which are exactly the usual irrotational 2D SWE. For instance (62) reads in Cartesian coordinates using the usual vector calculus notations:
[TABLE]
3.4 The geometric structure of port-Hamiltonian systems
In the previous section, we proposed a port-Hamiltonian formulation for open systems of two canonically interacting conservation laws (distributed state space equations (54) and the boundary equations (57)). Since boundary energy flows are considered, boundary port variables are needed to derive the power balance equation (58). The port-Hamiltonian model (54, 57) may then be implicitly defined as a linear subspace in the Bond space of effort and flow variables which embeds boundary effort and flow variables. In turn, this linear subspace may be geometrically defined as the linear subspace which is maximally isotropic with respect to some inner product associated to the natural power product - or power form in the Bond space - between effort and flow variables. Therefore, in the sequel, we aim at structure-preserving discretization which will preserve this power form in the approximation spaces. We will speak about symplectic discretization in the sense that this power form is preserved. In this section, we will define the Bond space, the power symplectic form, the associated inner product and the associated Stokes-Dirac structure which implicitely defines the port Hamiltonian model (54, 57). Readers are referred to [24] for details about this representation.
Let the Bond space of extended flow and effort variables be with
[TABLE]
We may define on this Bond space the real power pairing or power form which maps any effort-flow vector to
[TABLE]
in such a way that every pair of extended effort and flow variables in the Bond space, satisfying the port Hamiltonian equations (54) and (57)), also satisfies the power balance equation . From the power pairing (65), we may define the following symmetric bilinear form:
[TABLE]
With the help of this symmetric bilinear form (inner product on ) we may define the Dirac structure associated to the power pairing (65) as the linear subspace which is maximally isotropic, that is such that where the orthogonality is defined with respect to the inner product . In particular, any satisfies , hence the power balance .
Dirac interconnection structure may used to define implicitely the dynamics of port Hamiltonian systems. In particular, the linear subspace in the Bond space , with and as in (64), which is defined by:
[TABLE]
is a Dirac structure associated to the natural power pairing (65). This is proved by using the generalized Stokes theorem [24]. Therefore, in this particular case, the interconnection structure is called a Stokes-Dirac structure. The dynamics (54), generated by the Hamiltonian function (see definition (51)), with boundary energy flow and port boundary variables (57), may be implicitly defined by:
[TABLE]
In that sense, we will say that (68) defines a boundary port Hamiltonian system of two canonically interconnected conservation laws. In order to define a well-posed Cauchy problem, boundary conditions still need be chosen either for or .
Remark 3**.**
The Stokes-Dirac structure (67), associated to the natural power bilinear form (65), may be used to represent hyperbolic systems, either linear (when the Hamiltonian density is quadratic) or nonlinear (in the other cases). It may even be used to represent parabolic systems when the effort differential forms do not derive from the same Hamiltonian operator [30]. The choice of Dirichlet boundary conditions for or will lead to a well-posed system. They are however many other possible choices for admissible boundary conditions which lead to well-posed systems. For linear systems (quadratic Hamiltonian), these admissible boundary conditions have been parameterized elegantly [22]. Many examples of physical systems have been represented using the port Hamiltonian formulation and its implicit representation using Stokes-Dirac structures, including Maxwell field equations and Navier-Stokes flow problems [24], beam and membrane equations [1], vibro-acoustic problems [31], shallow water flow problems [32, 33, 34], advection-diffusion or adsorption problems [30], tokamak plasma MHD problems [3],etc.
3.5 The partitioned weak form for port-Hamiltonian systems of two conservation laws
We will follow the approach presented in section 2.2 for the 1D SWE and generalize it to the general covariant formulation (54) for port Hamiltonian systems of two conservation laws. Let and denote any arbitrary test differential forms. We may obtain from the strong formulation (54) the following weak form:
[TABLE]
Integrating by part the first state equation only, and using the integration by part formula (47), we get the following partitioned weak form:
[TABLE]
Remark 4**.**
With the particular choice and , one gets:
[TABLE]
Therefore the power-balance equation (6) reads:
[TABLE]
which shows that the power balance is also preserved in the weak formulation.
Remark 5**.**
Instead of integrating the first equation by parts in (70), we could integrate the second equation, which would also lead to another skew-symmetric structure; the boundary inputs and outputs would not be the same either.
In the case of the irrotational 2D SWE, written in Cartesian coordinates, using usual vector calculus notations, the partitioned weak formulation (70) transforms into
[TABLE]
In this latter equation, notations have been introduced both for simplification and to distinguish between scalar coordinates (for functions and 2-forms) and vector coordinates (for 1-forms). For instance, has been used instead of and denotes the water level (proportional to the water surfacic density), while and directly stands for the corresponding functions (0-forms). The bold vector notations , and denote respectively the vectors of Cartesian coordinates for the momentum (), water flow () and test () which are all 1-forms. The notation denotes the usual 2D Nabla differential operator while is the vector of Cartesian coordinates for the length 1-form (perpendicular to the boundary ). Note that the term denotes the boundary port variable (perpendicular water flow at the boundary) which will be chosen as the boundary input .
In the following section, we discretize the partioned weak-form (73), and we show that the resulting system is a finite-dimensional port-Hamiltonian system preserving the power-balance of the original system.
4 The Partitioned Finite-Element Method (PFEM) for the 2D SWE
In this section, the partitioned weak-form representation for the Shallow Water Equations, presented in § 3.5, will be projected into finite-dimensional approximation spaces in such a way as to preserve the total mass and momentum conservation equations and the underlying Dirac structure of the original port-Hamiltonian model. This section is divided in two parts, firstly, the weak-form is discretized in § 4.1 and a finite-dimensional port-Hamiltonian system is obtained. Secondly, we show how to obtain the discrete constitutive relationships in § 4.2.
4.1 Structure-preserving finite element discretization
Let us approximate the energy variables using the following basis with elements:
[TABLE]
where , are the chosen approximation basis functions in , and are the approximation coordinates for . The test functions and the co-energy variables are also approximated using .
Similarly, the vectorial variable is approximated as:
[TABLE]
where represents a 2D-vectorial basis function and, consequently, is an matrix. The variables are the approximation coordinates for in the approximation space. The test functions and the co-energy variable are also approximated using .
Remark 6**.**
Note that a particular choice for is :
[TABLE]
such that we can decompose the variables with index in their Cartesian components as: .
With this particular choice, we recover the case where the basis functions components of the vectorial variables are decoupled. We studied this case in [35].
Finally, the boundary input can be discretized using any one-dimensional set of basis functions, say :
[TABLE]
Remark 7**.**
In the sequel, in the implementation of our finite element method, we conveniently chose as evaluated on the boundary. Other choices could be investigated.
From the substitution of the approximated variables, (74),(75) and (77), in the weak form (73), the following finite-dimensional equations are obtained:
[TABLE]
Since these equations should remain valid for any test functions coordinates and , one gets:
[TABLE]
where and are square mass matrices (of size and , respectively), defined as:
[TABLE]
The matrix is of size , defined as:
[TABLE]
and is an matrix:
[TABLE]
Defining , the output conjugated to the input as:
[TABLE]
the approximated system can be written using the following finite-dimensional Dirac structure representation:
[TABLE]
where and denote the vector coordinates for the flow approximations coordinates in the and approximation bases, that is and .
From the definition of the co-energy variables as the variational derivatives of the Hamiltonian with respect to and , the time-derivative of the continous Hamiltonian is given by:
[TABLE]
Using the approximations for the energy and co-energy variables, this power balance can be approximated as:
[TABLE]
The approximated equations (84), together with the power balance (86), provides a Dirac structure representation that is a projection of the Dirac-Stokes structure. As we did for the 1D SWE in § 2.2, we will now derive the corresponding approximation of the port-Hamiltonian dynamics by restricting the Hamiltonian functional to the same approximation spaces. New energy and co-energy variables will be defined, so that we can find an explicit port-Hamiltonian representation for the approximated system.
Since and are symmetric positive definite matrices, we can use Cholesky decomposition222The use of Cholesky is advantageous since it simplifies the numerical inversion of the mass matrices. Note that the procedure presented in this section slightly differs from the one presented in § 2.2, since different energy and co-energy variables are defined here.. There exist triangular matrices and , such that: and . The power balance becomes:
[TABLE]
where new energy and co-energy variables:
[TABLE]
have been defined. Therefore, in order to write the power-balance as the total time derivative of the discrete Hamiltonian written as a function of the finite-dimensional vector coordinates for the energy variables, the newly defined reduced co-energy variables must be equal to the gradient of the Hamiltonian with respect to the energy variables:
[TABLE]
where the approximated Hamiltonian is defined as:
[TABLE]
Rewriting the finite-dimensional equations (116), we get the following explicit finite-dimensional port-Hamiltonian system:
[TABLE]
where is the conjugated output of the discretized system.
From (86) and (91), the time derivative of the approximated Hamiltonian is given by:
[TABLE]
Note that is the discrete analog of the continuous power-balance equation (19). Furthermore, this power-balance is exactly preserved in the finite-dimensional approximation spaces. From the definition of the matrix (78), the definition of the approximated boundary input and approximated co-energy variable , we get:
[TABLE]
Remark 8**.**
Note that using classical finite-elements 1D discretization basis for the boundary input, the coordinates provide the values of the inflow () at the boundary nodes. For instance, in the case of shallow water equations, these are the values of volumetric influx into the system. The conjugated output is related with the curve integral of along the elements. The co-energy variable is the pressure, thus the discretized outputs coordinates are related to the forces per unit length applied along the external boundary.
Note that we can also define an output that is given by the point-wise values of the co-energy variables evaluated on the boundary. In this case, a convenient choice of basis function for the approximation would be the same as the input, i.e.:
[TABLE]
where the coordinates represents the values of pressure at the boundary nodes. The power-balance through the boundary is computed as:
[TABLE]
where
[TABLE]
is a symmetric positive-definite mass matrix.
Furthermore, since the power-balances (92) and (95) must coincide, the following relationship between these two output definitions must hold:
[TABLE]
Consequently, from (87) and (95), the following power product must hold (satisfying (84)):
[TABLE]
Alternatively, the power balance can be reduced to (satisfying (91)):
[TABLE]
where and .
With the help of a symmetric bilinear form as (66), using any of the previous power products, finite-dimensional Dirac interconnection structures can be defined.
4.2 Obtaining the non-linear constitutive relationships: discretization of the Hamiltonian
In the previous section, a finite-dimensional Dirac structure was obtained for the 2D Shallow Water Equations, relating the energy and co-energy variables as well as the boundary inputs and outputs. The next step is to obtain the constitutive relationships from the Hamiltonian.
The Hamiltonian of the 2D SWE (60) can be rewritten using the coordinate variables as:
[TABLE]
The energy variables are restricted to the approximation spaces for and . From (90), the discretized Hamiltonian reads:
[TABLE]
The reduced effort variables are obtained from the gradient of the discretized Hamiltonian:
[TABLE]
Both constitutive equations exhibits nonlinear terms. In order to compute them, the following procedured was used. The first equation can be written as:
[TABLE]
where the components of can be computed as:
[TABLE]
Note that , for are matrices of dimension , that can be computed once and remains constant. Similarly, the second constitutive relationship can be written as a function of constant matrices.
5 Numerical experiments
In this section, we present numerical experiments to test the Partitioned Finite Element Method. Firstly, results for the 1D SWE are presented in § 5.1. Then, the 2D case is presented in § 5.2.
5.1 One-dimensional Shallow Water Equations
The example presented in Section 2 was implemented using finite elements with polynomial basis functions.
Firstly, a spectral convergence analysis of the numerical method was done. The eigenvalues obtained from the linearized numerical model were compared to the exact eigenfrequencies of the linear wave equation with constant coefficients. The inputs of (32) were considered to be zero.
Recall that the variables , and must be discretized in with polynomials of order at least one (since they are derived once on (16)). Fig. 1 shows the relative error of the first modal frequency for four different choices of polynomial approximations. stands for first order polynomial for the variables related to , and order zero for the variables related to . uses first-order polynomial for both all variables. uses first-order polynomial for the variables, and variables. Finally, uses third-order polynomial for all variables.
Time-domain simulations for the Shallow Water Equations considering the following boundary conditions:
[TABLE]
such that these conditions represent a harmonic influx through both boundaries.
Two different amplitudes were used in the simulations. Snapshots of the simulations are presented in Figs. 2 and 3. The first figure shows snapshots for a small amplitude value for the input inflow. In the second figure, the amplitude is multiplied by 100. Nonlinear phenomena is observed in this case.
5.2 Two-dimensional Shallow Water Equations
The example presented in § 4 was implemented using quadrilateral finite elements with polynomial basis functions for a square domain.
As we did for the 1D SWE in the previous subsection, firstly, a spectral convergence analysis of the numerical method was done.
For the 2D SWE, recall that the variables with index (, , ) must be discretized with polynomials of order at least one (since they are derived once on (73)). Fig. 4 shows the relative error of the first modal frequency for three different choices of polynomial approximations. stands for ”i-th” order polynomial for the variables with index , and ”j-th” order for the two components of the vector variables of index (, and ).
Numerical simulations were performed using the discretized system under boundary-port excitation. The following boundary conditions were considered:
[TABLE]
and . The boundary is split in four sides: . These conditions impose a harmonic inflow on one side of the boundary and the opposite condition on the other side. Simulations for two different values of amplitude are presented. First, snapshots for small amplitudes are presented in Figs. 5 and 6. Secondly, snapshots for large amplitudes are presented in Figs. 7 and 8.
Figure 9 shows how the approximated Hamiltonian (101) and total volume of fluid () change with time. As expected, the Hamiltonian only changes during the first second of simulation, while the system is excited through the boundary ports. After that, since , the Hamiltonian is constant. Finally, the total volume is kept nearly constant and the changes are only due to numerical precision (of order ).
6 Extensions
In this final section, the extension of PFEM to more general configurations is addressed, in order to illustrate the flexibility of the proposed numerical method: in § 6.1 the specific expression in polar coordinates is investigated, in § 6.2 the case of heterogeneous medium with variable coefficients is presented, and in § 6.3 the method is applied to second-order differential operators, such as those involved in the port-Hamiltonian formulation of the Euler-Bernoulli beam in 1D.
6.1 Polar coordinates
The goal of this subsection is to prove the applicability of PFEM in 2D, when the chosen coordinate system is not Cartesian: polar coordinates are presented. In order to avoid unnecessary technicalities, the geometry of a disc has been chosen to illustrate this extension.
The 2D Shallow Water Equations in a disc, as a Port-Hamiltonian system in polar coordinates
Let us consider the disc of radius with boundary , the circle of radius . Polar coordinates and will be used. In vector calculus, the 2-form and the 1-form are represented by the scalar function , and the vector function respectively. The Hamiltonian reads
[TABLE]
The effort or co-energy variables can be computed as , and (which is associated to the 1-form . Finally, the 2-form translates into the scalar , and the 1-form translates into the vector . With these notations and definitions, we get the same system as (63) for the strong form of the pHs, namely:
[TABLE]
with boundary control and collocated boundary observation at the boundary . Let us conclude with the energy balance for this system:
[TABLE]
The Partitioned Finite Element Method directly applies to the port-Hamiltonian system in polar coordinates
Now, let us first rewrite the weak form (73) with test functions and .
[TABLE]
Let us approximate the scalar energy variables using the following basis with elements:
[TABLE]
The variables and are also approximated using .
Similarly, the vectorial energy variable is approximated as:
[TABLE]
where represents a 2D-vectorial basis function and, consequently, is an matrix. Remark 6 does apply here also. Furthermore, and are also approximated using .
Finally, the boundary input, localized on the circle of radius can be discretized using any one-dimensional set of basis functions, say , provided -periodicity is ensured (trigonometric polynomials are a fair trial approximation basis, see e.g. [36, Chapter 18]):
[TABLE]
Introducing the notation and for the matrices of partial derivatives of the functions , we define matrix
[TABLE]
where the apparent singularity at has been removed. Then, with classical mass matrices , , together with the control matrix , the finite-dimensional equations become:
[TABLE]
where and are square matrices (of size and , respectively). is an matrix and is an matrix.
Defining , the output conjugated to the input as:
[TABLE]
with boundary mass matrix , the approximated system can be written using the finite-dimensional Dirac structure representation given by (84), and as found in remark 8, the global energy balance reads
[TABLE]
which mimicks that at the continuous level, namely (110).
6.2 Heterogeneous case with variable coefficients
The goal of this subsection is to prove the applicability of PFEM when the coefficients are space varying. In order to avoid unnecessary technicalities, the choice has been made to tackle the 1D model, first derived in § 2 as introductory example. Another fully worked out example can be found in [37] on the anisotropic heterogeneous wave equation in 2D.
The variable-coefficient physical model as a port-Hamiltonian system
Let us consider the Shallow Water equation in a water channel with a space-varying cross section, i.e. with the width of the channel, it is easy to understand that the energy given by (2) remains unchanged, with function instead of coefficient , but then then the balance equations (1) must be modified as follows:
[TABLE]
With the appropriate choice of energy variables , and , the Hamiltonian (2) now reads:
[TABLE]
The co-energy variables are found to be , and ; thus system (119) becomes in compact form:
[TABLE]
to be compared with (4) in the uniform case. Hence the new interconnection operator reads:
[TABLE]
Since for smooth scalar functions and , a straightforward computation shows that for smooth vector-valued functions and vanishing at the ends of the interval, and with the standard scalar product in . Hence, the unbounded matrix-valued differential operator proves skew-symmetric in .
Remark 9**.**
The port-Hamiltonian formulation (121) makes use of a non canonical Stokes-Dirac structure with non uniform coefficents. Instead, one could use the intrinsic geometric formulation presented in subsection 3.2, together with constitutive equations in covariant form for the SWE in subsection 3.3 and make use of the same energy state variables as in the uniform case, that is:
[TABLE]
In this approach, information about the geometry - here the space-varying reach width - will be embedded in the constitutive equation, that is in the Hamiltonian functional, through the definition of the appropriate Hodge star operator. More precisely, using the non uniform metric derived from the duality product , one gets
[TABLE]
Therefore the Hamiltonian functional (60) reads
[TABLE]
and the corresponding co-energy variables defined in (61) read
[TABLE]
which are the usual hydrodynamic pressure and the volume flow conjugated variables in the 1D hydraulic domain. Using these co-energy variables, the 1D SWE with non uniform width reads:
[TABLE]
All the geometry-dependent parameters have been embedded in the constitutive equations and the structural interdomain coupling between potential and kinetic energies (i.e. between mass and momentum balance equations) in the SWE pops up again in the canonical 1D Stokes-Dirac structure in (122).
The Partitioned Finite Element Method directly applies to the port-Hamiltonian system with variable coefficients. Here, the same procedure as in § 2.2 is being followed. We begin with a weak formulation of (121), then two complementary choices can be made.
If we choose to integrate by parts the mass balance equation only, i.e. the first line of the obtained weak form, we get exactly the same finite-dimensional pHs as (22), but with a new matrix, defined by
[TABLE]
and a new control matrix . The boundary control remains , and the new collocated boundary observation reads .
If instead, we choose to integrate by parts the momentum balance equation only, i.e. the second line of the obtained weak form, we get the following finite-dimensional pHs:
[TABLE]
but with another matrix, defined by
[TABLE]
and a new control matrix . The boundary control is now defined by , and the new collocated boundary observation reads .
Finally, note that in both the above cases, the following power balance is met:
[TABLE]
6.3 Higher-order systems
In the previous sections, the PFEM was applied to first-order (1D and 2D) formally skew-symmetric differential operators. Indeed, the method seems to be much more general and can be applyied similarly to higher-order equations.
The Euler-Bernoulli beam equation can be written as a port-Hamiltonian system of second order
The equations are given by (see, e.g., [9]):
[TABLE]
where and are obtained from the variational derivative of the Hamiltonian:
[TABLE]
From the definition of the variational derivatives, the time-derivative of the Hamiltonian is computed as:
[TABLE]
Note that depends only on the boundary values of (vertical speed), (moment), (rotation speed) and (force). This motivates the definition of the boundary ports, which allows writing the infinite-dimensional equations as port-Hamiltonian systems. From (128), one possible definition for the boundary ports is as follows:
[TABLE]
The final power balance () can thus be written as:
[TABLE]
Weak-form representation of Euler-Bernoulli beam equation
Let us use arbitrary test functions and and develop a weak form of (126):
[TABLE]
Integrating the first equation by parts twice, we get the following partitioned weak form:
[TABLE]
Finite-dimensional port-Hamiltonian system
Similarly to the development in § 2.2, we chose finite-dimensional bases functions (20) as and , for the variables with index 1 and 2, respectively. From the substitution of the approximation functions in the weak form (132), we find:
[TABLE]
where and are square mass matrices (of size and , respectively), equivalent to (23). The matrix is of size :
[TABLE]
and is an matrix:
[TABLE]
Finally, the conjugated-output can also be written in terms of the previous matrix:
[TABLE]
Defining the flow variables as and , we find the following finite-dimensional Dirac structure representation:
[TABLE]
Following the same procedure presented in the previous sections for the 1D and 2D SWE, from the discretization of the Hamiltonian using the energy variables approximation spaces, one gets the underlying port-Hamiltonian dynamics for the approximated Euler-Bernoulli beam equations.
The analogue of the Euler-Bernoulli beam in 1D is the Kirchhoff plate in 2D, one can refer to [38] for the modelling as a port-Hamiltonian system using tensor calculus, and the application of PFEM to it, with various boundary controls ; note that the analogue of the Timoshenko beam in 1D is the Mindlin plate in 2D, and PFEM can also be applied to this model [39].
7 Conclusion and Open questions
The Partitioned Finite Element Method provides a full-rank structure-preserving representation of port-Hamiltonian systems in 2D and 3D: a general setting has been proposed here, written in the language of differential forms, and also translated into vector calculus for common PDE applications. This method can be easily implemented thanks to ready to use FEM software to compute the matrices of the representation, which are all sparse. It applies to complex geometries, works in any coordinate systems, and allows for space-varying coefficients; moreover higher order differential operators can also be tackled. Although PFEM has already been succesfully applied to linear PDEs with quadratic Hamiltonian functionals, e.g. vibrating membranes and plates, here the methodology carries over to a non-linear PDE with non-quadratic and non-separable Hamiltonian functional, the irrotational Shallow Water Equation in 2D.
Future work will deal with mixed boundary control (possibly leading to differential algebraic problems as pHDAEs), and mathematical convergence analysis (choice of the finite element bases and theoretical rate of convergence). Some worked-out 2D test cases will be studied on coupled systems, e.g. fluid-structure interaction (FSI), or thermal-structure coupling. Lastly, structure-preserving model reduction techniques will be tested on the high-fidelity finite-dimensional systems obtained by PFEM, see e.g. [40] for pHs or [41] for pHDAEs. The reduced order system will then be most useful to apply dedicated control laws for pHs, like IDA-PBC, which do take advantage of the specific structure of these dynamical systems with collocated inputs and outputs, see e.g. [42].
THANKS Part of this work has been performed in the framework of the Collaborative Research DFG and ANR project INFIDHEM n∘ ANR-16-CE92-0028 (http://websites.isae.fr/infidhem).
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1[1] V. Duindam, A. Macchelli, S. Stramigioli, H. Bruyninckx, Modeling and Control of Complex Physical Systems: The Port-Hamiltonian Approach, Springer Berlin Heidelberg, Berlin, Heidelberg, 2009.
- 2[2] Y. Wu, B. Hamroun, Y. Le Gorrec, B. Maschke, Power preserving model reduction of 2D vibro-acoustic system: A port Hamiltonian approach, Proceedings of the 5th IFAC Workshop on Lagrangian and Hamiltonian Methods for Nonlinear Control (LHMNC) 48 (13) (2015) 206–211.
- 3[3] N. M. T. Vu, L. Lefèvre, B. Maschke, A structured control model for the thermo-magneto-hydrodynamics of plasmas in tokamaks, Mathematical and Computer Modelling of Dynamical Systems 22 (3) (2016) 181–206.
- 4[4] V. Trenchant, H. Ramirez, Y. L. Gorrec, P. Kotyczka, Structure preserving spatial discretization of 2D hyperbolic systems using staggered grids finite difference, in: Proceedings of the 2017 American Control Conference, Seattle, USA, 2017.
- 5[5] A. Bossavit, Whitney forms: A class of finite elements for three-dimensional computations in electromagnetism, IEE Proceedings A (Physical Science, Measurement and Instrumentation, Management and Education, Reviews) 135 (8) (1988) 493–500.
- 6[6] A. Bossavit, Computational electromagnetism: variational formulations, complementarity, edge elements, Academic Press, 1998.
- 7[7] G. Golo, V. Talasila, A. J. van der Schaft, B. Maschke, Hamiltonian discretization of boundary control systems, Automatica 40 (5) (2004) 757–771.
- 8[8] 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.
