A Concise and Geometrically Exact Planar Beam Model for Arbitrarily Large Elastic Deformation Dynamics

The potential of large elastic deformations in control applications, e.g., robotic manipulation, is not yet fully exploited, especially in dynamic contexts. Mainly because essential geometrically exact continuum models are necessary to express these arbitrarily large deformation dynamics, they typically result in a set of nonlinear, coupled, partial differential equations that are unsuited for control applications. Due to this lack of appropriate models, current approaches that try to exploit elastic properties are limited to either small deflection assumptions or quasistatic considerations only. To promote further exploration of this exciting research field of large elastic deflection control, we propose a geometrically exact, but yet concise a beam model for a planar, shear-, and torsion-free case without elongation. The model is derived by reducing the general geometrically exact the 3D Simo–Reissner beam model to this special case, where the assumption of inextensibility allows expressing the couple of planar Cartesian parameters in terms of the curve tangent angle of the beam center line alone. We further elaborate on how the necessary coupling between position-related boundary conditions (i.e., clamped and hinged ends) and the tangent angle parametrization of the beam model can be incorporated in a finite element method formulation and verify all derived expressions by comparison to analytic initial value solutions and an energy analysis of a dynamic simulation result. The presented beam model opens the possibility of designing online feedback control structures for accessing the full potential that elasticity in planar beam dynamics has to offer.


INTRODUCTION
Dynamic robotic manipulation of highly deformable objects is still a rarely considered field in literature due to a lack of appropriate dynamic models. Elastic objects in robotic manipulation are usually either considered only quasistatically (Bretl and McCarthy, 2014) or dynamically under the assumption of small deflections which in return results in applications where elasticity is often treated as an undesired property that needs to be avoided or compensated (Tavasoli, 2015). There are some approaches that suggest exploiting elasticity instead, e.g., in terms of safety (Bicchi and Tonietti, 2004;Haddadin et al., 2009), payload estimation (Malzahn et al., 2015), or energy storage (Tantanawat and Kota, 2007;Haddadin et al., 2011). Some authors exploit elasticity also for quasistatic manipulation (McCarragher, 2000) or even to design new strategies for dynamic manipulation, e.g., Pekarovskiy et al. (2014) or Pham and Pham (2017). Nonetheless, we propose that elastic dynamics have even more potential that can be exploited when large deflections are considered.
Models that allow taking advantage of such arbitrarily large deformations need to fulfill two main criteria. They have to 1. describe the dynamics geometrically exact, i.e., independently of the magnitude of deformation 2. be concise and simple enough to allow the design of online feedback control However, such models are not yet available in literature. In this study, we start filling this gap for the case of a planar elastic beam undergoing pure bending, by proposing a single-dimensional dynamic continuum model.

Related Work
Large deflection dynamics are typically treated with one of two approaches-discretized approximations using a finite element method (FEM) or continuum models expressed as partial differential equations (PDEs). The FEM-based descriptions, for which recent literature such as the study of Duriez (2013) reaches real-time capable control, are limited to quasistatic deformable structures, as the computational cost of FEM descriptions for true dynamic large deformation online feedback control is still out of reach. Whereas, the geometrically exact continuum models result in multivariate and highly coupled nonlinear PDE systems, and control approaches are thus limited to oscillation damping (Ito, 2001;Hegarty and Taylor, 2012). Unlike existing control methods for linear beam models (Krstic et al., 2006a;Krstic et al., 2006b), the literature body on theory of nonlinear PDE systems is still too limited in its applicability for the complex expressions arising in these mechanical system models (Padhi and Ali, 2009). From a control point of view, however, it is well known that continuum models, unlike discretized approximations, do not face so called spillover phenomena that can lead to instabilities due to unmodeled high frequency dynamics (Meirovitch and Baruh, 1983). For these reasons, as a first contribution toward filling the discussed gap in literature, this work proposes a geometrically exact model of planar Euler-Bernoulli beam dynamics for arbitrarily large deflections that admit a surprisingly concise PDE formulation.
Although there are special purpose models for large deflection models, such as bullwhip dynamics by Bernstein et al. (1958) and McMillen and Goriely (2003) or dynamics of a fly fishing line by Spolek (1986), we are interested in more general beam dynamics that admit rope and cable dynamics as a special case of very low elastic stiffness. A vast literature body exists on 1D analytic beam theories alone. Therefore, because of the considerable computational advantage of 1D beam theories over 3D continuum mechanics theories, the latter will not be considered in this concise literature review. Starting with the first mathematical treatment of static elasticity by Galileo already in 1638, Hooke's treatise of linear elasticity is followed in 1678. The first precise definition of the elastica (a thin strip of elastic material) problem was carried out by Jakob Bernoulli (1691), and he published its first solution in Bernoulli (1694). His nephew Daniel Bernoulli (1742) did not himself solve the problem, but he suggested Euler the use of variational analysis, who delivers a closed-form solution of the elastica in the study by Euler (1744). A more detailed and insightful mathematical historical overview of the elastica can be found in the study by Levien (2008). While Euler's early work already predicts slender beam deformations with astounding precision, many authors built on this work to include further effects to account for more general conditions as well as geometries. It is said that Rayleigh (1877) added rotary inertia effects, and Timoshenko (1921) further enhanced the theory to account for shear effects. However, Elishakoff (2020) discusses original authors and naming of linear beam theories. He, e.g., mentioned that Bresse (1859) already included rotary inertia effects before Rayleigh, though the works were developed independently. Furthermore, a beam theory including shear effects was originally published in Timoshenko (1916), an earlier book in Russian language, where Timoshenko mentions to have developed the theory together with P. Ehrenfest. Elishakoff, therefore, suggests the historically justifiable name Bresse-Rayleigh-Timoshenko-Ehrenfest beam theory.
For large deflections, also referred to as finite strain, geometrically nonlinear models are necessary. Kirchhoff (1859) is a spatial generalization of the Euler-Bernoulli beam and allows modeling of 3D deformations through bending and torsion. The theory was later extended by Love (1892) to further account for axial tensions and is referred to as Kirchhoff-Love beam theory. Reissner (1972) added further measures to Kirchhoff's theory, accounting for shear deformations in planar curves and later for space curves in Reissner (1981). Simo (1985) enhanced Reissner's work in terms of approximations to what is nowadays known as Simo-Reissner beam theory or geometrically exact beam theory. To complete the overview, it is also worth mentioning that reduced versions of these two well-known theories have been proposed that neglect torsion modes, e.g., the study by Meier et al. (2015) for the Kirchhoff-Love and Romero et al. (2014) for the Simo-Reissner case. Meier et al. (2019) gave a more in-depth overview and analysis of these nonlinear beam theories. Table 1 shows how our proposed model fills the current literature gap of a concise model for geometrically exact descriptions.

Contribution
The main contribution of this work is twofold. We provide • the first single-dimensional, geometrically exact, dynamic beam model, and • a method for incorporating boundary conditions in a FEM formulation, for cases where the FEM model has only descriptive variables of higher-order derivatives than the boundary condition itself

Outline
In Section 2, we derive the model via step-by-step reduction, starting from the general Simo-Reissner beam theory, extracting first the Kirchhoff-Love beam theory, followed by special case assumptionsisotropic, torsion-free, inextensible, and planar bending. The translation into a FEM description is explained in Section 3, including our proposal for incorporating position-level boundary conditions into the tangent angle PDE system. This FEM description is applied for the simulation verification in Section 4. The work is concluded in the final Section 5.

MODELING
This section outlines the model derivation, starting from a general 3D theory. After reducing the model by gradually introducing assumptions, we couple the Cartesian coordinates to achieve an expression in a single PDE.

Nomenclature
This work follows the convention of using lowercase bold variables for vectors and uppercase bold variables for matrices. All nonbold variables are scalars. Furthermore, subscript annotations are reserved for index notation of multidimensional variables as well as expressing partial derivatives, whereas superscript annotations are part of the variable specification. Also, note that we omit explicit listing of function parameters whenever it is clear from the context to not unnecessarily clutter the notation. A list of the most frequently used variables at the end of this work is given.

Model Reduction
The presented geometrically exact model for a planar Euler-Bernoulli beam is found via reduction of the general 3D Simo-Reissner beam theory. We gradually introduce further assumptions to simplify the dynamic governing equations of the beam. Please note that we only define the individual components of the equations once necessary to keep the derivation clear and easy to follow. The resulting model forms a special case of a planar Kirchhoff-Love beam theory parametrized solely in the curve tangent angle, in an analog to Euler's elastica (Euler, 1744). A more in-detail analysis and discussion of the general Simo-Reissner and Kirchhoff-Love beam theories, including the special cases of isotropic cross-sections and torsion-free formulations, can be found in Meier et al. (2019).

General Simo-Reissner Beam Theory
The first theory that accounts for very general 3D beam deformations including spatial bending, torsion, axial tension, and shear deformation was published by Simo (1985). The consideration of shear effects makes it an adequate theory for thick rod dynamics. It is also been denoted by geometrically exact beam theory because the description is consistent at the deformed state regardless of the magnitude of displacements, rotations, and strains, cf. Crisfield and Jelenić (1999). Simo himself also used the term finite strain beam formulation. The strong form of the Simo-Reissner beam theory, cf. Simo (1985), is a system of six coupled PDEs and can be stated with the equilibrium equations where the internal force f internal and moment vector f internal results form internal stresses acting on the beam cross-section area at point r of the beam center line. The quantities f ext and m ext account for externally imposed forces, and f intertia and m intertia are the components due to inertia effects. The detailed constitutive equations that relate f internal and m internal to the first Piola-Kirchhoff stress tensor require an introduction into 3D continuum mechanics and is omitted in this work. Instead, we point the interested reader to related text books such as by Gurtin (1982) and define the expressions only after reduction to the specified special cases.
In the following, all objective deformation measures are chosen to be work conjugated to the material stress resultants in Eq. 1. We further assume a hyperelastic constitutive relation between these kinetic and kinematic quantities.

Assumption: Vanishing Shear Strains (Kirchhoff-Love Beam Theory)
Neglecting shear deformations and assuming that the crosssection is always perpendicular to the center line of the beam, the change in internal forces f internal l can be split up into a parallel component f internal l and ⊥ f internal l , a component perpendicular to the center line. Furthermore, the moment balance Eq. 1b reduces to the projection onto the center line tangential base vector g 1 , Figure 1 for an illustration. The Kirchhoff-Love beam equations are thus given with where Eq. 2b is now a scalar expression, and the beam model is thus reduced to 4 PDEs.

Assumption: Initially Straight and Isotropic
If now an initially straight beam with an isotropic cross-section is assumed for the hyperelastic beam, the components of Eq. 2 are given with with the Young modulus E, inertia I, density ρ, cross-section area A, and the curvature vector The components κ 1 , m ext 1 , and ω 1 in Eq. 3b relate to the curvature, externally imposed moment, and angular velocity along the tangential direction of the beam.

Assumption: Torsion-Free
Assuming pure bending and no torsional effects, the inertia moment and the moment balance Eq. 3b vanish completely. Moreover, only the perpendicular component of the external moment affects the force balance equation with the axial tension parameter ε : (r) l 2 − 1, which considers that the relation of the undeformed beam, does not in general hold for the deformed case zr zl 2 ≠ 1, due to possible elongations in the beam structure.

Assumption: Inextensible Beam
If it is assumed that the beam does not undergo axial elongations, the gradient Eq. 7 does always equals 1. Hence, the axial tension parameter ϵ evaluates to and thus, f internal l vanishes. The beam model Eq. 5 consequently further simplifies to However, unlike the previous assumptions that can be incorporated implicitly with an adequate choice of parametrization variables, it is in general difficult to find such a set of variables that fulfill the inextensibility constraint Eq. 8a by construction. A common practice to enforce the equality constraint Eq. 8a on the simulation result in a weak sense, i.e., in an integral form instead of point-wise, is by means of extending the weak form of the model Eq. 8b with a Lagrange multiplier potential, cf. Meier et al. (2019).
The following assumption of pure planar bending, however, does again permit a parametrization that fulfills this constraint directly in the strong sense, i.e., for every point along the beam.

Assumption: Pure Planar Bending
The last step in the model reduction is the restriction to pure planar bending. For the remainder of this section, we will switch to a component-wise notation in Cartesian coordinates. The beam model from Eq. 8b reduces to two coupled PDEs and is fully described by and the additional inextensibility constraint because the third row of the vector equation Eq. 8b vanishes for the purely planar problem. Thus, the remaining external inputs are forces f ext in the xy-plane, as well as the external moment ⊥ m ext z is perpendicular to this plane. With the curve tangent angle φ defined as the angle between the tangent vector of the deformed beam center line g 1 and its undeformed counterpart e 1 , the beam curvature Eq. 4 can be expressed as for the shear-free, planar, inextensible case. Furthermore, it allows stating the geometric identities The planar beam model Eq. 9 can thus be stated as a beam model in a mixed form containing Cartesian coordinates as well as the curve tangent angle as describing variables. In the remainder of this Section, Eq. 12 is the starting point to first develop the static case followed by the general dynamic case, both entirely expressed in the curve tangent angle.

Static Beam Model Expressed in the Curve Tangent Angle
Only considering solutions in a static equilibrium, the Cartesian acceleration terms in Eq. 12 vanish and remains. By computing all derivatives, and rotating the equations from their Cartesian xy coordinate system by φ around the z-axis, by premultiplying Eq. 14 with the rotation matrix which allows extracting the components perpendicular ⊥ and parallel to the beam center line The fact that no Cartesian xy parameter of the beam description remains but is rather described in the curve tangent angle cta alone means that the geometric identities Eq. 11 and the planar inextensibility constraint Eq. 9b are now implicitly fulfilled by construction. No further treatment such as Lagrangian multipliers are thus necessary, in contrast to the above beam models (Eq. 8b and Eq. 9). In case of an absent external force f ext , the static beam model Eq. 16 even admits a simple analytic solution. If a nontrivial curvature (φ) l ≠ 0 is assumed, Eq. 16 reduces to which can be integrated twice and yields a unique solution if boundary conditions are applied.

Dynamic Beam Model Expressed in the Curve Tangent Angle
To also fully state the dynamic planar beam model in terms of the curve tangent angle, the Cartesian xy acceleration terms in Eq. 12 remain to be expressed in terms of the curve tangent angle φ. This is achieved by differentiating the system of equations Eq. 12 w.r.t. the beam parameter l. Assuming no buckling of the object, x and y have continuous derivatives, and thus, Schwarz's theorem allows changing the order of the derivations. Applying the geometric identities (Eq. 11) now also to the acceleration terms, leads to with the auxiliary variable a PDE system entirely expressed in terms of the curve tangent angle φ. As for the static case (Eq. 16), the geometric identities fulfill the inextensibility constraint Eq. 9b by construction; thus, no special consideration is necessary. Expanding all partial derivatives and grouping the trigonometric terms, Eq. 18a yields Similar to the static case, premultiplying the entire system with the rotation matrix R z (φ) from Eq. 15 again extracts the components parallel and perpendicular to the beam center line. The only acceleration term (φ) tt , however, appears solely in the perpendicular direction which in the case of no external inputs admits the very concise strong form Frontiers in Robotics and AI | www.frontiersin.org May 2021 | Volume 7 | Article 609478 as a single PDE governing the beam dynamics in a single parameter φ. While this reduced model is relevant for PDE controller development, it is not directly applicable for use in simulations. We therefore present in the following section a respective approximation with a system of ordinary differential equations (ODEs), in terms of a FEM formulation.

FEM FORMULATION
In this section, we outline the development of a FEM simulation procedure, starting from the development of the weak form of the beam model (Eq. 21), without considering external forces. After transforming the integrodifferential weak form into a system of nonlinear ODEs of thee second order in time via a Bubnov-Galerkin approximation, a finite element discretization leads to a simulation procedure.

Weak Form of Large Deformation in the Curve Tangent Angle
The weak form of Eq. 21 is found by multiplying the equation with the test function δφ and integrate over the beam length l [0, L]: A sequence of integrations by parts will lead to the final weak form. In the first intermediate step, Using this result, the terms on the right-hand side of Eq. 22 result in that builds the basis for the following FEM formulation. The function φ as well as the variation δφ have to be the members of the Sobolev space H 2 , where with the function space of square integrable functions such that they are twice continuously differentiable in l, cf. Reddy (2013). A common method to choose candidates for φ and its variation δφ is given by the Bubnov-Galerkin approximation and eventually leads to a system of ODEs.

Bubnov-Galerkin Approximation
In the sense of the Bubnov-Galerkin method, the function φ as well as the test function δφ will be approximated by using the same set of n weighted orthogonal spatial basis functions ψ 1...n ∈ H 2 , together with n timedependent scaling coefficients a φ i for the approximation of the curve tangent angle and b φ j for the test function. The weak formulation (Eq. 24) thus reads As the coefficients b φ j from the variation Eq. 25b are arbitrary, the weak formulation result in the system of n equations for j [1, n]. Reorganizing the terms and using a matrix representation finally leads to with the component-wise definitions Note the index order of the matrix component definitions M ji and F ji that is important for the vector notation in Eq. 28a. In the context of FEM formulations, a particular choice of piecewise orthogonal basis functions ψ is used.

Finite Element Discretization
The core idea of the FEM is to discretize the continuous body into a finite number of N elements e ∈ [1, N], connected at the N + 1 node locations N 0...N , and choose a special set of piecewise function elements ψ 1...N that form the global spline approximation φ h (l) of Eq. 25a.
For the weak formulation (Eq. 24), the orthogonal basis functions need to be at least of class H 2 , such that the integral-containing second derivatives can be evaluated. In this work, we use, for demonstration purposes, the prominent choice of a Hermite cubic splines with cubic elements ψ 1...n , constructed by the Hermite local cubic polynomial basis functions: They have the special property that those coefficients that build the spline element, ψ e (λ) : ξ e 1 (λ) ξ e 2 (λ) ξ e 3 (λ) ξ e 4 (λ) directly correspond to the function value φ and its derivative (φ) l at the node locations N e−1 and N e . Figure 2 illustrates the interplay of the local element basis functions and the function values at the node locations to approximate the full curve tangent angle function. Although the full dynamics are encoded in the resulting system of ODEs, what are missing for a well-posed problem are the boundary conditions of a particular simulation case.

Boundary Conditions Expressed in the Curve Tangent Angle
Incorporate the boundary conditions on the curve tangent angle φ, and its derivatives follow standard procedures, e.g., Hughes (2012). Because treating position-based boundary conditions is not directly possible in the curve tangent angle beam model (Eq. 21), this chapter focuses on the development of a strategy to express boundary conditions in terms of higher order derivatives only. Without loss of generality, we consider for our beam model a fixed end, and, respectively, its dynamic counterpart Fixing an elastic beam in its position at one end introduces point-wise reaction forces from the mounting onto the beam in x and/or y directions. While these boundary conditions are straight forward to be incorporated in a FEM formulation for a model in the parameters x and y, e.g., Eq. 9a, the FEM description of the reduced model (Eq. 28a) directly acts on the tangent angle function φ and its derivatives, thus not offering any parameter to incorporate position boundary conditions. However, except for the a-priori known position at the boundary condition itself, the curve tangent angle dynamics Eq. 21 do not require any position parameters to govern the beam profile.
We thus propose to transform the position boundary condition at node N 0 , which cannot be incorporated directly, into a dynamic boundary condition for the neighboring node N 1 entirely expressed in curve tangent coefficients a φ . To define this substitutional boundary condition, we first derive another FEM formulation for the beam model parametrized in Cartesian coordinates (Eq. 9a). Not considering external forces for simplicity and recalling the geometric identities Eq. 11, the beam equation Eq. 9a reduces to Formulating the weak forms and applying another integration by parts reads and the Bubnov-Galerkin approximation with the same set of orthogonal functions ψ that leads to the systems of equations with the components Note that for both systems of equations in Eq. 36, the matrices M are the same as for the FEM in the curve tangent angle φ from Eq. 28a.
Again using Hermite cubic spline basis functions (Eq. 30), the equations of interest in the two systems of N equations (Eq. 36) are the ones relating to the position basis function of N 1 , i.e., j 3: As illustrated in Figure 2, this requires an integration from N 0 until N 2 to fully account for all associated local basis functions and thus involves the first six coefficients of the acceleration vectors, a x 1..6 tt : and a y 1..6 tt : relating to the function values and first derivatives at the first three nodes. While the right hand side of the FEM formulations (Eq. 36) is already defined in the curve tangent angle φ, what remains is to also rewrite coefficient vectors (a x 1..6 ) tt and (a y 1..6 ) tt in terms of φ instead of x and y. Starting again from the geometric identities (x) l ≡ cos(φ) and (y) l ≡ sin(φ), the coefficients in Eq. 39 can be expressed as Frontiers in Robotics and AI | www.frontiersin.org May 2021 | Volume 7 | Article 609478 y ltt (l, t) : cos φ(l, t) φ(s, t) tt ds − sin φ(l, t) φ(l, t) 2 t , (40c) containing Volterra integrals with an upper limit l. Recalling the spline approximation φ h from Eq. 30 and using it for the acceleration terms, allows to express the entire FEM balance equations (Eq. 38) in terms of φ h with the additional boundary values (x) tt (0, t) and (y) tt (0, t) which are, however, known a-priori from the position boundary condition we are incorporating. Eventually, Eq. 38 can be evaluated entirely in φ with the left hand sides and M 3,1..6 a y 3,1..6 tt where ○ , denotes an element-wise product, considering the piecewise function element definitions (ψ e ) tt from Eq. 30 and the local Hermite cubic base functions ξ e (λ, t) expressed in the global arc length coordinate ξ e (l, t) : ξ e (λ l/N, t). Note that only the first six columns of M affect the position value at N 1 , and all remaining entries of the third row of M are thus zero. The right hand side load value in x direction of N 1 expressed in φ reads and the load value in y direction reads These expressions are used to incorporate the position boundary conditions into the FEM formulation of the curve tangent angle beam equation.
The resulting system of nonlinear ODEs is used to verify the proposed beam model as well as the presented strategy for incorporating position boundary conditions.

SIMULATION VERIFICATION
The proposed beam model of Section 2 together with the method for incorporating boundary conditions on lower level derivatives than the descriptive variables of the FEM formulation from Section 3.2 is verified in simulation. First, the developed FEM strategy is tested in two different initial value problems and demonstrated for different magnitudes of external forces and momenta. The second part verifies a dynamic simulation in terms of energy consistency of the beam profile. For both cases, we consider the case of a clamped end at l 0 and a free end at l L. This is expressed in the following boundary conditions: φ λλ (l L, t) 0.
While the essential boundary condition Eq. 43c and natural boundary conditions Eq. 43d and Eq. 43e are directly considered in the FEM formulation with conventional techniques, conditions Eq. 43a and Eq. 43b are incorporated with the expressions developed in Section 3.4. All simulations have been conducted in MATLAB R2020a.

Initial Value Problem
The static solutions of the FEM formulation (Eq. 28a) with external nodal forces f ext , are verified in two scenarios. First, the analytic solution of Eq. 17 is replicated with an external moment at the free end, and in the second example, an external nodal force f ext is applied to the middle of the beam. The initial value problem for both cases is formulated as the nonlinear least-squares problem where the last two rows impose the boundary conditions Eq. 43a and Eq. 43b using the expressions from Eq. 42a and Eq. 42b. The results show that stem from MATLAB's nonlinear least-squares solver lsqnonlin() for a nominal beam of L 1 m and material parameter c 1, with a FEM discretization into 10 beam elements.

External Moment
Applying an external nodal moment at the free end of the beam results in and thus a constant curvature (φ) l along the entire beam, according to the analytic solution (Eq. 17). This is also consistent with analytical solutions known in literature Antman (2005). To verify the proposed formulation with this test case, various external momenta proportional to the material parameters ⊥ m z ext (l L) ∈ {0.5, 1, 2, 4}π/EI in N m are applied to the free end of the nominal beam. They are incorporated directly as boundary condition (φ) l (l L) ⊥ m z ext , substituting Eq. 43d. The results are shown in Figure 3 and depict segments of a perfect circle due to the constant curvature (φ) l . While (x) l and (y) l components show a FEM approximation of cos(·) and sin(·) functions, respectively, the according angle φ is a purely linear function in this special case and thus can be approximated with arbitrary accuracy, even for a low number of elements. Note that the analytical trigonometric solution is not visually distinguishable from the simulation result and thus is not shown in Figure 3. consisting of the potential energy E pot , the kinetic energy E kin , and the externally injected energy E ext . While the potential energy is directly proportional to the curvature of the beam center line, It can be directly evaluated from the simulation results. The kinetic energy however, needs to be rewritten in terms of φ first. The geometric identities (x) l ≡ cos(φ) and (y) l ≡ sin(φ) together with Schwarz's theorem on changing the order of derivatives again allow to rewrite the Cartesian components at an arc length l of the beam center line, with the Volterra integrals and thus, the kinetic energy of the entire beam can be evaluated with Note that the rotational component of the kinetic energy is not considered, due to the Euler-Bernoulli assumption of slender beams.

Simulation Results
For the dynamic case, an initially unloaded beam with a normalized beam material parameter c 1 and unit length L 1 m is perturbed by an external nodal force, and again at the beam center, with a maximal unit force f ext max 1 over a perturbation time of t max 0.5 s. The result of the dynamic simulation with 10 beam elements is shown in Figure 5. While in the first 0.5 s, energy is injected in the system via the external perturbation, the resulting beam movement results in an energy exchange between kinetic and potential along the beam while conserving the total amount of energy in the system. Note that no dissipative terms such as damping are considered in the simulated beam (Eq. 21). The right plot of Figure 5 shows the energy distribution within the entire beam at t 1 m. The accumulation of potential energy at the clamped site of the beam again demonstrates the effectiveness of the position boundary condition expressions Section 3.4.

CONCLUSION
A model for large planar deformation dynamics of Euler-Bernoulli beams was presented and put into context with well-known more general beam models. Literature does already offer various models that account for arbitrarily large deformations; however, they typically result in a system of coupled nonlinear PDEs expressions. Whereas, the presented approach admits a single-dimensional PDE in one variable, i.e., the curve tangent angle of the beam center line to describe planar beam dynamics under the common Euler-Bernoulli assumptions of shear-free constant cross-sections.
While boundary conditions on the beam profile derivatives-which is sufficient for sliding and/or free ends-can be directly encoded in a simulation algorithm of the curve tangent angle beam model, there is no descriptive variable available to directly incorporate boundary conditions on the beam position-needed for clamped and/or hinged ends. These cases are, however, of course of FIGURE 5 | Energy conservation during an FEM simulation of the dynamic curve tangent angle model in a clamped-free scenario, where an external torque is applied to the center of an initially straight beam for 0.5 s. The left plot shows the energy distribution. The right plot shows the energy along the beam at a snapshot taken at t 1 s.
Frontiers in Robotics and AI | www.frontiersin.org May 2021 | Volume 7 | Article 609478 highly practical relevance. To also address these cases, we additionally outlined a novel method that allows incorporating boundary conditions in FEM formulations, using solely descriptive variables of higher order derivatives. We apply this method to impose positionbased boundary conditions in the FEM formulation of the presented beam model, but it is not limited to solely use this case. The strategy is verified in initial value problems, where it replicates analytical solutions and a time-variant FEM simulation by evaluating energy conservation. The presented beam model does not only reduce the computational effort due to the dimensional reduction to a single parameter, but the beam profiles are at the same time less complex in this curve tangent angle parametrization and thus require fewer elements in the FEM description.
Although being nonlinear, the derived beam model provides a concise continuum model for future control theory applications of large deformations, where other more complex model descriptions are not appropriate for current model-based PDE controller development. The model reduction process in this work, however, also outlines more general beam models considering, e.g., shearing, axial torsion, elongation, and/or 3D spatial deformations. This work is thus also a good source for extended beam models including these additional effects, which might be relevant with the advancement of PDE controller development.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding authors.

AUTHOR CONTRIBUTIONS
GH, DW and MB initiated the study. GH formulated, implemented and conducted the research. All authors contributed to manuscript revision, read, and approved the submitted version.