Port-Hamiltonian formulation and Symplectic discretization of plate models Part II : Kirchhoff model for thin plate

# Port-Hamiltonian formulation and Symplectic discretization of plate models Part II : Kirchhoff model for thin plate

Andrea Brugnoli Daniel Alazard Valérie Pommier-Budinger Denis Matignon ISAE-SUPAERO, Université de Toulouse, France.
10 Avenue Edouard Belin, BP-54032, 31055 Toulouse Cedex 4.
###### Abstract

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.

## Introduction

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

 ρ(x)\diffp[2]wt(x,t)+\diffp[2]x(EI(x)\diffp[2]wx)=0,x∈(0,L),t≥0 (1)

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

 αw =ρ(x)\diffpwt(x,t),Linear Momentum, (2) ακ =\diffp[2]wx(x,t),Curvature.

Those variables are collected in the vector , so that the Hamiltonian can be written as a quadratic functional in the energy variables

 H=12∫L0αTQαdx,whereQ=[1ρ(x)00EI(x)]. (3)

The co-energy variables are found by computing the variational derivative of the Hamiltonian

 ew :=\diffdHαw=\diffpwt(x,t),Vertical % velocity, (4) eκ :=\diffdHακ=EI(x)\diffp[2]wx(x,t),%Flexuralmomentum.

Those variables are again collected in vector , so that the underlying interconnection structure is then found to be

 \diffpαt=Je,whereJ=[0−\diffp[2]x\diffp[2]x0]. (5)

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

 f∂=⎛⎜ ⎜ ⎜ ⎜⎝ew(0)\diffpewx(0)\diffpeκx(L)eκ(L)⎞⎟ ⎟ ⎟ ⎟⎠,e∂=⎛⎜ ⎜ ⎜ ⎜⎝\diffpeκx(0)−eκ(0)−ew(L)\diffpewx(L)⎞⎟ ⎟ ⎟ ⎟⎠. (6)

The power flow is then easily evaluated as

 ddtH(t)=∫L0\diffpαt⋅edx=⟨e∂,f∂⟩IR4. (7)

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

 B:={(f,f∂,e,e∂)∈F×E}. (8)

The duality pairing between elements of is then defined as follows

 ≪((fa,fa∂),(ea,ea∂)),((fb,fb∂),(eb,eb∂))≫:=∫L0{(fa)Teb+(fb)Tea}dx+(fa∂)Teb∂+(fb∂)Tfa∂. (9)

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

 D={(f,f∂,e,e∂)∈F×E|f=−Je,f∂=(ew(0)\diffpewx(0)\diffpeκx(L)eκ(L))Te∂=(\diffpeκx(0)−eκ(0)−ew(L)\diffpewx(L))T}. (10)

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

 u(x,y,z)=−z\diffpwx,v(x,y,z)=−z\diffpwy,w(x,y,z)=w(x,y). (11)

and for the strains

 ϵ=⎛⎜⎝ϵxxϵyyγxy⎞⎟⎠=⎛⎜⎝\diffpx00\diffpy\diffpy\diffpx⎞⎟⎠(−z\diffpwx−z\diffpwy)=−z⎛⎜⎝\diffp[2]wx\diffp[2]wy2\diffpwx,y⎞⎟⎠. (12)

The curvature vector is defined as

 κ=⎛⎜⎝κxxκyyκxy⎞⎟⎠=⎛⎜⎝\diffp[2]wx\diffp[2]wy2\diffpwx,y⎞⎟⎠. (13)

The Hooke constitutive law for isotropic material is considered for the constitutive relation

 (14)

where the Poisson ratio and the Young modulus. The generalized momenta are found by integrating the stresses along the fiber

 M=⎛⎜⎝MxxMyyMxy⎞⎟⎠=∫h2−h2−zσdz=(∫h2−h2Ez2dz)κ

If it is assumed that the mechanical properties remain constant along the axis, then momenta are expressed vectorially as follows

 M=⎛⎜⎝MxxMyyMxy⎞⎟⎠=⎛⎜ ⎜⎝D(\diffp[2]wx+ν\diffp[2]wy)D(\diffp[2]wy+ν\diffp[2]wx)D(1−ν)\diffpwx,y⎞⎟ ⎟⎠. (15)

where is the bending rigidity and, in the most general case, . The relation between momenta and curvatures is expressed by matrix :

 M=Dκ,D:=∫h2−h2Ez2dz=D⎡⎢ ⎢⎣1ν0ν100012(1−ν)⎤⎥ ⎥⎦. (16)

The Kirchhoff-Love model arises as the Euler-Lagrange equation stemming from the following extrema problem (Hamilton principle):

 w0={w∈V;δF(w)|w0=0}

and

 F(w)=∫T0∫Ω{K−U+W}dΩdt=∫T0∫Ω{12μ(\diffpwt)2−12D[(\diffp[2]wx)2++(\diffp[2]wy)2+2ν(\diffp[2]wx)(\diffp[2]wy)+2(1−ν)(\diffpwx,y)2]+pw}dΩdt,

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

 W =pw, K =12μ(\diffpwt)2, U =12D[(\diffp[2]wx)2+(\diffp[2]wy)2+2ν(\diffp[2]wx)(\diffp[2]wy) +2(1−ν)(\diffpwx,y)2], =12κTDκ.

The total energy density is split into kinetic and potential energy

 H=K+U, (17)

and the corresponding total energies given by the following relations

 H=∫ΩH dΩ,K=∫ΩK dΩ,U=∫ΩU dΩ. (18)

 μ\diffp[2]wt+\diffp[2]Mxxx+2\diffpMxyx,y+\diffp[2]Myyy=p. (19)

The spatial derivatives of the acceleration have been neglected. If the coefficients are constant then the ruling EDP becomes

 μ\diffp[2]wt+DΔ2w=p, (20)

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

 α:=(μv, κxx, κyy, κxy)T, (21)

where . The Hamiltonian density is given by the following expression

 H=12αT[1μ00D]α,H=∫ΩH dΩ. (22)

So its variational derivative provides as co-energy variables

 e:=\diffdHα=(v, Mxx, Myy, Mxy)T, (23)

The port-Hamiltonian system and the formally skew-adjoint operator relating energy and co-energy variables are found to be

 \diffpαt=JeandJ:=⎡⎢ ⎢ ⎢ ⎢⎣0−\diffp[2]x−\diffp[2]y−(\diffpx,y+\diffpy,x)\diffp[2]x000\diffp[2]y000\diffpx,y+\diffpy,x000⎤⎥ ⎥ ⎥ ⎥⎦. (24)
###### Remark 1.

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

 ˙H =∫Ω\diffdHα⋅\diffpαtdΩ =∫Ω{v(−\diffp[2]Mxxx−\diffp[2]Mxxy−2\diffpMxyx,y)+Mxx\diffp[2]vx+Myy\diffp[2]vy+2Mxy\diffpvx,y}dΩ =∫Ω{e1(−\diffp[2]e2x−\diffp[2]e3y−2\diffpe4x,y)+e2\diffp[2]e1x+e3\diffp[2]e1y+2e4\diffpe1x,y}dΩ

In Figure 2 the notations for the different reference frames are introduced. By applying Green theorem, considering the split mixed derivative ()

You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters