Geometrically exact planar Euler-Bernoulli beam and time integration procedure for multibody dynamics

A new formulation of geometrically exact planar Euler-Bernoulli beam in multi-body dynamics is proposed. For many applications, the use of the Euler-Bernoulli model is sufficient and has the advantage of being a nodal displacement-only formulation avoiding the integration of rotational degrees of freedom. In this paper, an energy momentum method is proposed for the nonlinear in-plane dynamics of flexible multi-body systems, including the effects of revolute joints with or without torsional springs. Large rotational angles of the joints are accurately calculated. Several numerical examples demonstrate the accuracy and the capabilities of the new formulation.


Introduction
Flexible multi-body dynamics has been an intensive topic of research for the last decades due to its many applications in different areas of engineering. Examples of multi-body dynamical systems can be found in aerospace, mechanical, manufacturing and transportation engineering. Manipulators, robots and space structures [1][2][3][4][5][6][7][8] are examples of multi-body dynamical systems. Also bio-dynamical systems e.g human bodies and animals are best understood as multi-body dynamical systems [9][10][11].
It's clear that a rigid multi-body formulation is no longer suitable for today's engineering applications and the flexibility of the members need to be considered [12,13]. For systems experiencing large deformations, special formulations should be used. The classical geometrically exact beam theory, first introduced by Reissner [14] in an intrinsic form and implemented by Simo and Vu-Quoc [15], introduces a rotation vector as a degree of freedom of the system. Different implementations can be found in the literature depending on the type of updating schemes adopted for the rotational degrees of freedom including co-rotational frames for planar beams [16,17] and inertial frames [18].
However, the integration of rotational parameters cannot be achieved using standard solvers for ordinary differential equations on vector space. In multi-body dynamics, the rotation tensor is also discontinuous at the joints. Problems associated with integration of rotational degrees of freedom are discussed in [17,19] to which the reader is referred. Alternatively, the Absolute Nodal Coordinate Formulation (ANCF), first proposed by Shabana [2,16,[20][21][22][23][24] for the analysis of large deformation problem in multi-body dynamics, has gained some popularity in the multi-body dynamics community. The beam element in ANCF is described by using a global slope coordinates that define the element orientation.
In most multi-body dynamics applications, the members are slender beams. Hence, an Euler-Bernoulli model could be completely sufficient with the great advantage of delivering a nodal displacement-only formulation. Nevertheless, for some applications, the shear flexibility of the beam cannot be neglected. The standard Timoshenko-type beams finite element formulations do suffer from locking phenomena [21,[25][26][27][28][29][30]. However, it should be mentioned that more complex formulations addressing locking are available among which those based on the helicoidal approach [31][32][33][34][35].
Indeed, besides the kinematics description, the time integration scheme is also at the heart of multi-body dynamics. Especially when dealing with high non-linearity, classical time integrations cannot be employed because of numerical instabilities [36][37][38][39][40][41][42]. The use of long-term stable time integration schemes is then crucial.
In an attempt to avoid rotational degrees of freedom while still accounting for full geometric non-linearity, the present authors resorted to the classical Bernoulli beam and developed a suitable energy-momentum method which provides both high levels of accuracy and numerical stability. In this paper, we extend our approach to account for multi-body dynamics, where beams are connected via revolute joints with or without torsional springs. At the heart of our approach is to consider the strain fields obtained by integrating the strain rates. Not only the displacements and strains, but also the rotation angles of the joints are integrated. It will be shown that the finite rotation angles of the joints can be accurately obtained. Furthermore, the conservation of energy will be proved formally and numerically. We restrict ourselves to planar motion and show that the discontinuity at the joint rotation can be removed using a displacement-only formulation.
The paper is organized as follows. A review of the equations describing the deformation of a planar rod is given in "Rod kinematics and strain measures" section. In "Revolute joint with torsional spring" section, we develop a new method to circumvent the rotation discontinuity for two planar rods connected by a revolute joint with a torsional spring. "Hamilton's principle and the dynamic equations" section is dedicated to Hamilton's principle and the corresponding dynamic equations. The finite element method and the energy-momentum time integration are detailed in "Finite element formulation and time integration scheme" section. In particular, formal proof of energy and angular momentum conservation are given. In "Numerical examples" section, we provide a range of numerical examples demonstrating the strength and applicability of the new approach.

Rod kinematics and strain measures
In this section, we briefly introduce the beam's strain measures. We consider a general Cartesian co-ordinate basis e i , i = 1, 2, 3 and a curve X 0 (s), with s being the arc-length. We understand X 0 (s) as describing the centre line of the rod cross section. We restrict ourselves to in-plane deformations taking place in the e 1 − e 2 plane and introduce the tangent vector G 0 = dX 0 (s) ds . Perpendicular to this vector, we define N to be the normal vector with z as the corresponding coordinate in the direction of N .
In addition to the Cartesian system, we define a suitable convected curvilinear coordinate system given by the triple (s, z, x 3 ). We define the vector X(s, z) = X 0 (s) + z N (s) as the position vector of points in the direction of N at the reference configuration. In a general cross section, the boundary is defined by z being a function of x 3 . Note that the third direction is not explicitly included in this equation, though it is implicitly understood that it is in the direction of e 3 and that the deformation is independent of that direction. Consequently, a local basis in the reference configuration is defined by the triple (G, N , e 3 ), We have also the following rela- where | • | denotes the absolute value of a vector, × and (·) denotes the cross and scalar product of vectors, respectively. The corresponding contra-variant basis vectors are then given by (G * , N , e 3 ), with G * = G |G| 2 . A deformation is given as x = ϕ(X) which defines the actual configuration. For our geometrically exact beam theory, we make use of the Bernoulli model where the cross sections are assumed to be rigid and remain perpendicular to the center line. Thus, the corresponding tangent vectors at the deformed configuration are defined as (g, n, e 3 ) with g 0 = g| z=0 , and the following relations hold ( Fig. 1) The displacement of material points in the direction of n can be characterized by the following relation where u(s) is the displacement vector of the centre line. From this we obtain immediately g = X 0 ,s + u ,s + z n ,s . The deformation gradient is written down in the curvilinear basis system as F = g ⊗G * +n⊗N +e 3 ⊗e 3 . The Green strain tensor is defined as E = 1 2 (C−I), where C = F T F. It has only one non-trivial component E 11 which is where the term in z 2 has been neglected since the thickness of the beam is small compared to its length. E 11 is split into two components as E 11 = ε 11 + z κ, the definition of which are given by The latter is the classical change of curvature.

Revolute joint with torsional spring
Let's consider first the frictionless planar revolute joint (see Fig. 2). This type of joints is very common in mechanisms and allows only a relative rotation between two rods. We call the two cross-sections with a hinge, A and B. The orientations of the cross sections are determined by their respective normal vectors. Denoting the relative rotation at the joint ω, one has ω = γ − , where is defined to be the angle between the vectors N B and N A . Similarly, the same holds for γ at the deformed configuration. We recall that and γ can be defined as the absolute value of the rotation vectors parallel to e 3 with the corresponding rotations from N B to N A and from n B to n A , respectively. These rotations are to be defined via the kinematic quantities N A , N B , n A and n B . From Fig. 2, one can derive the following relations: from which it follows with k being some integer. The expression for ω now reads Since the hinge can freely rotate, the only kinematic condition is that the displacement u defined in (3) is the same at sections A and B: Note that there is no condition about the continuity of the first derivative of u between the bodies. When an elastic torsional spring is added to the revolute joint, it generates an internal moment which is a function of the relative rotation between the two bodies of the system, see Fig. 3. The elastic energy stored in the torsional spring when it deforms is written as with C as the spring constant.

Hamilton's principle and the dynamic equations
We assume our multi-body system to be conservative. Hamilton's principle states that the action functional is stationary. That is where t refers to time, t 1 and t 2 are boundaries of the time interval, L is defined as the Lagrangian given by

Fig. 3 Revolute pair with torsional spring
Here T is the kinetic energy of the system and int , ext are the internal and the external potentials, respectively. The quantities in (17) are defined as where E is Young's modulus, L is the total length of all members, K is the number of torsional springs, q(s) is the distributed external force, P i , i = 1, ..., N and u i are the concentrated forces and the corresponding displacements, respectively, M 3j , j = 1, ..., M and θ j are concentrated external moments and the corresponding rotation angles, respectively. From (16), (17), together with an integration over the cross section, and due to the fact that the variations vanish at the boundaries, standard arguments of the calculus of variation deliver with I and A being the moment of inertia and the cross-section area, respectively. Since u is the only degree of freedom of the system, we have to relate the relative rotations to it and to its derivatives. Calculation of δθ starts from the fact that θ is the angle between N and n, we have Accordingly, the functional in (21) can be written as follows where u A ,s and u B ,s define the vectors n A and n B , respectively. In general, the above equation takes the form where m, k, f are the generalized mass vector, the generalized internal force vector and the generalized external force vector, respectively. The total energy of the system, which in most cases coincides with the Hamiltonian, is defined by Equation (26) is the starting point for a numerical treatment which includes also the time integration scheme. For more details about the mass matrix, internal and external force vector, please refer to the Appendix A.

Finite element discretisation
As a result of the Bernoulli hypothesis, functional (25) contains second derivatives. Hence, within a finite element method, the first derivative must exhibit continuity. We resort to the standard approach of utilizing Hermite polynomial within a finite element context. Each node comprises four degrees of freedom which are the displacement and their derivatives: u 1 , u 2 , u 1 , u 2 . For the displacement field the interpolations read where and with ζ = [−1, 1] and l is the length of the element.

An energy-momentum time integration scheme
Since the achievement of Simo and Tarnow [38], it has been realized that time integration schemes should reproduce conservation properties of the analytical model in the discrete case. Especially, energy conservation or energy control is considered key to long-term stable integration. In what follows, we will extend an energy-momentum method previously developed by [43][44][45] to cover the case of multibody dynamics. Assuming that at time step n, all kinematical fields and velocities are known, we need to find these quantities for time step n+1. Consider ξ to be a scalar which defines any position within the time interval t, with 0 ≤ ξ ≤ 1. We start with the following expressions: The first relation defines a convex set, the following two are true for some values of ξ . The midpoint rule corresponds to the choice ξ = 0.5. Consequently, to integrate u one has, where u = u n+1 − u n . The key step in our energy-momentum method is to integrate strain velocity fields to define the strain fields instead of using Eq. (5), (6). The same is done for the highly non-linear kinematic relations as in Eqs. (2) and (13), where velocities are defined for these too. The equations are used to merely define the velocity fields themselves. One haṡ Given the strain fields and the kinematic quantities n and ω at time t n , the strain fields and the same kinematic quantities at step n + 1 are then defined as follows: With these so defined quantities, the midpoint rule provides energy conservation. Using (48), the value of ω could be calculated exactly without limitation (for example 2π) because only the difference of ω is of importance.

The field equations
To obtain our field equations we start from functional (25) Using the following identities: Because the variations are arbitrary, one Euler-Lagrange equation with two boundary conditions follow from the above functional: EI κ ∂κ

Proof of conservation of angular momentum
Let's recall that x 0 = X 0 + u. To prove the conservation of angular momentum we start by taking the cross product of Eqs.
Integrating Eq. (58) gives Using the following identities: together with subsequent application of Gauss's theorem, Eq. (61) presents us with the following statement Using (59) Now we are going to prove that In component form, the normal vector n defined in (2) reads Here 3jk is the permutation symbol defined by 3jk = ⎧ ⎪ ⎨ ⎪ ⎩ 1 for even permutation 3,j,k -1 for even permutation 3,j,k 0 for two equal values of 3,j,k where δ ij is the Kronecker delta defined as The term N · n can be rewritten as follows which, together with (71), provides With this relation at hand, Eq. (69) gives Finally we get Hence, (68) is proven. and Consider next, Using similar procedure to get the results from (70) to (75) by replacing N by n B , we proved that Similarly for x 0,s n+ To further simplify this expression, we first are going to prove that Let's consider the second term in (81). After some lengthy algebraic manipulations, Eq. (69) deliver where use of x 0 3 ,s = 0 has been made. Consider further the expressionn × x 0,s · ∂n ∂x 0,s , which is equivalent to From (85) and (86), we conclude: Hence (82) is proven. Consider next From (87) and (88) we can write: Consider further the third term in (81) and recall the expression for ε defined in (6), one has ∂ε ∂u ,s = x 0,s . Hence, Finally, we are going to prove the following statement Consider the last two integrals in (81). The expression for κ, found in (6), can be rewritten as follows: which gives The term x 0,s × κ ∂κ ∂x 0,s results in Further, for x 0,ss × κ ∂κ ∂x 0,ss , it follows Hence, Eqs. (95) and (95) mean nothing but Using (96) which is equivalent to The expression is equivalent to the concise one which indicates that in the case of vanishing moments, we immediately infer This means nothing but the conservation of the total angular momentum.

Proof of conservation of the total energy
In the following, we will give a concise proof of energy conservation which is different than the proof presented earlier in [45]. Starting from Eq.
which, with the use of the midpoint rule, is equivalent to The last equation leads to Assuming that the external loading is conservative, from Eq. (104) we finally obtain which proves nothing but the conservation of the total energy. We refer the reader to Appendix for the proof of conservation of linear momentum. It is worthwhile to note that in some cases high frequency damping is a desirable quality of the integration scheme. This can be achieved by considering material-type viscosity in addition to the elastic response. However, this is out of the scope of this paper. Alternatively, slight modification of the present integration scheme can result in numerical damping, though this approach is not pursued here any further.

Numerical examples
To demonstrate the performance of the above method in geometrically exact multi-body dynamics, we investigate a range of examples, each of which will be run for one million time steps. The joint is modeled by revolute pair with or without a torsional spring.

Nonlinear dynamics of a flexible two-body system
The first example deals with a classical multi-body dynamics system consisting of two flexible rod-like bodies connected through a revolute joint without torsional spring. The two-body system is attached to a pinned support and subjected to a conservative horizontal force applied at one third of the beam length from the joint (see Fig. 4). The load increases linearly to a peak and decreases at the same rate to zero as shown in Fig. 5.  Figure 6 shows the number of iteration necessary in order for the Newton-Raphson method to converge. For t = 10 −2 s no more that 5 iterations are are needed to achieve convergence whereas for t = 10 −3 s 3 iterations are sufficient to reach dynamic equilibrium.

Fig. 6 Number of iterations until convergence
Further, we investigate the influence of the time step size on the results. In Fig. 7, the vertical displacement of the free edge (end) of the two-body system versus time is plotted for different time steps: t ≤ 10 −3 , t = 10 −2 , t = 0.1 and t = 0.2. For t ≤ 10 −3 and t = 1E−2, the results are similar and only a slight deviation can be observed for t = 0.1 and t = 0.2. However, in term of energy conservation, Fig. 8 shows that numerically the energy conservation is achieved with high degree of accuracy for t ≤ 10 −3 and t = 10 −2 whereas for larger time steps ( t ≥ 0.1), a numerical increase of energy can be observed. It should be noted however that for flexible beams with very high stiffness as considered in the present example, large time steps are physically not meaningful. The results of our energy-momentum method are now compared to solutions using classical midpoint rules. Figure 9 shows that the classical midpoint rule already suffers from blow up at t = 7 s with t = 10 −2 s. The energy blow-up can be seen in Fig. 10.

Nonlinear dynamics of a deformable multi-body system with flexible joint
To investigate the performance of the method for multi-body systems with torsional springs, we consider the previous example and equip the revolute joint with a flexible torsional spring (see Fig. 11). The system is subjected to a conservative horizontal force applied at free end of the multi-body system. The loading history is depicted in Fig. 12. The calculation was performed considering one million time steps with t = 10 −3 s.     Figures 13 and 14 show the horizontal and vertical displacements at the free end of the two-body system. The conservation of energy is depicted in Fig. 15. The time history of the rotation angle of the hinge is shown in Fig. 16. The rotation angle increases up to about 8 rad. It should be stressed that this smooth behaviour of the rotation angle is a direct outcome of the time integration method where the spring rotation angle is computed by integrating the rotational velocity via Eq. (44) and (48). Some snap shots of the two-body motion are shown in Fig. 17 together with the corresponding time and the value of the hinge rotation angle.
Finally, to demonstrate not only the stability and conservation properties of the proposed scheme but also its ability to compute finite spring rotation, we examine the response of this multi-body system using the standard mid-point time integration scheme. Figure 18 shows the displacement of the load application point versus time. An instability can be observed at about t = 95 s. The time history of the spring rotation angle is plotted in Fig. 19. It can be seen that ω has a correct value for ω = [−π, π]. When the spring rotation angle reaches ±π, we observe a sharp jump in its value, which is clearly incorrect. Hence, the displacements are also not inappropriate. This behaviour is a consequence of the use of Eq. (13) which provides a straightforward means to calculate the spring rotation. This equation is well defined only for angle amplitudes less than |π|. Again, a smooth behaviour is possible if the spring rotation angle is computed using the proposed

Free moving multi-body system
To further asses the performance of the formulation, we consider a free moving multibody system. More specifically, we investigate the motion of a two-elements mechanism connected to each other by a torsional spring (see Fig. 21). The mechanism is subjected to two external bending moments applied at each end. The time evolution of the loading is depicted in Fig. 22. The calculation has been performed considering one million time steps. In this example, not only the conservation of energy (see Fig. 23), but also the conservation of linear momentum (see Fig. 24), and the angular momentum (Fig. 25) are captured. Similarly, the time history of the rotation angle of the hinge and some snap-shots of the motion are illustrated in Figs. 26 and 27, respectively.

Slider-crank system
The last example deals with the classical slider-crank mechanism composed of two flexible beams connected by a revolute joint (see 28). On one side, the mechanism is attached to a pinned support. On the other side, the mechanical system is connected to a weightless slider using a revolute joint. The slider can move only in the horizontal direction and is   subjected to an externally applied moment acting at the beam end cross-section connected to the pinned support. As a consequence of the loading, the beam rotates about the outof-plane axis and the mechanism converts a rotary motion to a straight-line motion. The flexible beam element parameters are given below and the loading history is shown in Fig. 29.       Another excellent performance of the formulation is shown with this classical example of multibody mechanism. The energy is perfectly conserved for the whole time (see Fig. 30). The displacement time history of the slider is depicted in Fig. 31, here again no sign of instability is observed. Snap-shots of the motion are presented in Fig. 32.     where is the global degrees of freedom. T 1 and T 2 are the localised matrix of the two elements connected by the torsional spring. The global residual vector related to torsional springs is written as follows

Appendix B. Proof of conservation of linear momentum
To prove the conservation of linear momentum, we start by integrating the field Eq. Hence, the conservation of linear momentum is now proven.