Port-Hamiltonian formulation and
Symplectic discretization of plate models
Part II : Kirchhoff model for thin plate
The mechanical model of a thin plate with boundary control and
observation is presented as a port-Hamiltonian system (pHs111pHs stands for port-Hamiltonian systems.), both in vectorial and tensorial forms: the Kirchhoff-Love model of a plate is described by using a Stokes-Dirac structure and this represents a novelty with respect to the existing literature. This formulation is carried out both in vectorial and tensorial forms. Thanks to tensorial calculus, this model is found to mimic the interconnection structure of its one-dimensional counterpart, i.e. the Euler-Bernoulli beam.
The Partitioned Finite Element Method (PFEM222PFEM stands for partitioned finite element method.) is then extended to obtain a suitable, i.e. structure-preserving, weak form. The discretization procedure, performed on the vectorial formulation, leads to a finite-dimensional port-Hamiltonian system. This part II of the companion paper extends part I, dedicated to the Mindlin model for thick plates. The thin plate model comes along with additional difficulties, because of the higher order of the differential operator under consideration.
keywords:port-Hamiltonian systems, Kirchhoff plate, Partitioned Finite Element Method, geometric spatial discretization, boundary control.
As presented in part I of this companion paper, the port-Hamiltonian (PH) formalism BookZwart; bookPHs; Villegas allows the structured modeling and discretization of multi-physics applications involving interconnected finite- and infinite-dimensional systems (Cervera2007; ShaftIntInfinite; vanderShaftintFinInf). Preserving the port-Hamiltonian structure in the discretization process is a keypoint to take benefit of this powerful formalism. This issue was fisrt addressed in Golo, with a mixed finite element spatial discretization method, and in moulla:hal-01625008, with pseudo-spectral methods relying on higher-order global polynomial approximations. All those methods are difficult to implement, especially for those system the spatial dimension of which is bigger than one. Very recently weak formulations which lead to Galerkin numerical approximations began to be explored: in WeakForm_Kot, a structure preserving finite element method was introduced for the wave equation in a two-dimensional domain; this method exhibits good results, both in the spectral analysis and simulation part, though requiring of a primal and a dual mesh on the geometry of the problem. Another approach is the partitioned finite element method (PFEM) proposed in LHMNLCFlavio2018, already largely explored in part I of this companion paper; the advantages of this latter methodology are its simplicity of implementation and its potential to carry over to a wide set of examples, no matter the spatial dimension of the problem; the possible use of open source software like FreeFEM++ or FEniCS is also an appealing feature of this latter method
In part II of this companion paper, the modeling and discretization of thin plates described by the Kirchhoff-Love plate model is carried out within the PH framework, allowing for boundary control and observation. Compared to the existing literature on either the symplectic Hamiltonian formulation of LI2016984; LI2018310, or the port-Hamiltonian framework throgh jet theory in SchobertMindlin; SCHOBERL2015610, the main contribution of this paper concerns the PH formalism based on a Stokes-Dirac structure. This formalism is presented both in vectorial and tensorial forms. Moreover, the tensorial formalism (Grinfield, Chapter 16) highlights that this model mimics the interconnection structure of its one-dimensional counterpart, i.e. the Euler-Bernoulli beam. Compared to part I dedicated to thick plate Mindlin model in which first-order differential operators are explored in dimension two, and compared to LeGorrec2005 in which second- or higher-order differential operators were explored in dimension one only, the contribution of this paper is the PH formalism of systems of dimension two described with second-order differential operators, such as the Kirchhoff-Love model. The model, once written in a tensorial form, highlights new insights on second-order differential operators: especially the double divergence and the Hessian are proved to be adjoint operators one of another, which represents another important contribution of this paper. Finally, the extension of the PFEM method to the structure-preserving discretization of the Kirchhoff model is also a novelty of the paper. It allows simple implementation of numerical schemes compared to the jet theory formalism, while preserving the structure of PHS at the discrete level.
1 Second-order distributed PH systems: Euler-Bernoulli beam
The Euler-Bernoulli beam is the one-dimensional equivalent of the Kirchhoff-Love plate. The reader can refer to Gorrec2017 for a stability proof of a beam connected to non linear springs and dampers or to aoues:hal-01738092 for an illustration of rotating spacecraft with flexible appendages model as PH Bernoulli beams. This model consists of one PDE, describing the vertical displacement along the beam length
where is the transverse displacement of the beam. The coefficients and are the mass per unit length, Young’s modulus of elasticity and the moment of inertia of a cross section. The energy variables are then chosen as follows
Those variables are collected in the vector , so that the Hamiltonian can be written as a quadratic functional in the energy variables
The co-energy variables are found by computing the variational derivative of the Hamiltonian
Those variables are again collected in vector , so that the underlying interconnection structure is then found to be
For an infinite-dimensional system, boundary variables have to be defined as well. Those can be found by evaluating the energy rate flow across the boundary. One possible choice among others (see articleFlavio for a more exhaustive explanation) for this model is the following
The power flow is then easily evaluated as
The flow variables can now be defined as , so that the flow space is given by the tuples . Equivalently the effort space is given by . The bond space is therefore the Cartesian product of the two
The duality pairing between elements of is then defined as follows
The Stokes-Dirac structure for the Euler-Bernoulli beam is therefore the following
Theorem 1 (From LeGorrec2005, Stokes-Dirac structure for the Bernoulli beam).
Consider the space of power variables = and the bilinear form (+pairing operator) given by (9). Define the following linear subspace
Then , where is understood in the sense of orthogonality wit respect to the bilinear product i.e is a Dirac structure.
2 Kirchhoff-Love theory for thin plates
In this section the classical variational approach (Hamilton’s principle) to derive the equation of motions is first detailed. The physical quantities involved and the different energies, of utmost importance for the PH formalism, are reminded.
2.1 Model and associated variational formulation
The Kirchhoff-Love plate formulation rests on the hypothesis of small thickness compared to the in plane dimensions. The notations and symbols are borrowed form FEM_Cook and Oslo. The displacement field and the strains are defined by assuming that fibers orthogonal to the middle plane remain orthogonal (see figure 1). This leads to the following relations for the displacement field
and for the strains
The curvature vector is defined as
The Hooke constitutive law for isotropic material is considered for the constitutive relation
where the Poisson ratio and the Young modulus. The generalized momenta are found by integrating the stresses along the fiber
If it is assumed that the mechanical properties remain constant along the axis, then momenta are expressed vectorially as follows
where is the bending rigidity and, in the most general case, . The relation between momenta and curvatures is expressed by matrix :
The Kirchhoff-Love model arises as the Euler-Lagrange equation stemming from the following extrema problem (Hamilton principle):
where the surface density, the vertical load density. For example if gravity along the axis is considered . The domain is a connected region of . Furthermore, is a suitable functional space, which contains the essential boundary conditions, i.e. Dirichlet conditions; as usual, Neumann boundary conditions will be specified by means of the variational problem itself. This formulation is not complete since boundary terms may be added to take into account non homogeneous boundary conditions. Moreover, the external work , and , the kinetic and potential energy density per unit area, are respectively given by
The total energy density is split into kinetic and potential energy
and the corresponding total energies given by the following relations
The Euler-Lagrange equation reads
The spatial derivatives of the acceleration have been neglected. If the coefficients are constant then the ruling EDP becomes
where is the bilaplacian. In the sequel it is assumed . The specification of boundary conditions allows defining the function space .
3 PH formulation of the Kirchhoff plate
In this section the port-Hamiltonian formulation of the Kirchhoff plate is presented first in vectorial form in 3.1 and then in tensorial form in LABEL:sec:PH_ten_Kir.
3.1 PH vectorial formulation of the Kirchhoff plate
To obtain a port-Hamiltonian system (pHs) the energy variables as well as the underlying Stokes-Dirac structure, associated with the skew-adjoint operator , have to be properly defined. Since the Kirchhoff plate represents the 2D extension of the Euler-Bernoulli beam, it is natural to select as energy variables the linear momentum, together with curvatures. The energy variables are collected in the vector
where . The Hamiltonian density is given by the following expression
So its variational derivative provides as co-energy variables
The port-Hamiltonian system and the formally skew-adjoint operator relating energy and co-energy variables are found to be
From the Schwarz theorem for functions the mixed derivative could be be expressed as , instead of . However, in this way the symmetry intrinsically present in would be lost. The mixed derivative is here split to reestablish the symmetric nature of curvatures and momenta (that are of tensorial nature as explained in section LABEL:sec:PH_ten_Kir).
The boundary variables are obtained by evaluating the time derivative of the Hamiltonian
In Figure 2 the notations for the different reference frames are introduced. By applying Green theorem, considering the split mixed derivative ()