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

Carlo Sansour, T.L.U. Nguyen, Mohammed Hjiaj, S. Chhang

► To cite this version:

Carlo Sansour, T.L.U. Nguyen, Mohammed Hjiaj, S. Chhang. Geometrically exact planar Euler-Bernoulli beam and time integration procedure for multibody dynamics. Advanced Modeling and Simulation in Engineering Sciences, 2020, 7 (1), pp.33. 10.1186/s40323-020-00166-1. hal-02932029

HAL Id: hal-02932029

https://hal.science/hal-02932029v1

Submitted on 7 Sep 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
<figure><img src="image_2.png" /></figure>

RESEARCH ARTICLE

Open Access

<figure><img src="image_3.png" /></figure>

Check for
updates

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

Carlo Sansour, Tien Long Nguyen, Mohammed Hjiaj\* and Sophy Chhang

\*Correspondence:
Mohammed.Hjiaj@insa-rennes.fr
INSA Rennes, LGCGM/Structural
Engineering Research Group, 20
avenue des Buttes de Coësmes,
CS 70839, 35708 Rennes Cedex 7,
France

**Abstract**

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.

**Keywords:** Multi-body dynamics, Geometrically exact beam theory, Euler-Bernoulli beam, Energy-momentum method

**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–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–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

<figure><img src="image_4.png" /></figure>
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 Sha-bana [2, 16, 20–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–30]. However, it should be mentioned that more complex formulations addressing locking are available among which those based on the helicoidal approach [31–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–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 = \frac{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) + zN(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)$, with $G = \frac{\partial X}{\partial s}$, $N = \frac{\partial X}{\partial z} = \frac{\partial X}{\partial z}|_{z=0}$ and $G_0 = G|_{z=0}$. We have also the following relations $G \cdot N = G_0 \cdot N = 0$, $|N| = 1$, $N = e_3 \times \frac{G}{|G|} = e_3 \times G_0$, where $|\cdot|$ denotes the absolute value of a vector, $\times$ and $(\cdot)$ 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^* = \frac{G}{|G|^2}$.

A deformation is given as $x = \varphi(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)

$$g = \frac{\partial x}{\partial s}, \qquad (1)$$

$$n = e_3 \times \frac{g}{|g|} = e_3 \times \frac{g_0}{|g_0|}. \qquad (2)$$

The displacement of material points in the direction of $n$ can be characterized by the following relation

$$x = X(s) - zN(s) + u(s) + zn(s) = X_0(s) + u(s) + zn(s). \qquad (3)$$

where $u(s)$ is the displacement vector of the centre line. From this we obtain immediately $g = X_{0,s} + u_s + zn_s$. The deformation gradient is written down in the curvilinear basis system as $F = g \otimes G^* + n \otimes N + e_3 \otimes e_3$. The Green strain tensor is defined as $E = \frac{1}{2}(C-I)$, where $C = F^T F$. It has only one non-trivial component $E_{11}$ which is

$$E_{11} = u_s \cdot X_{0,s} + \frac{1}{2} u_s \cdot u_s + z [n_s \cdot (X_{0,s} + u_s) - N_s \cdot X_{0,s}], \qquad (4)$$

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} = \varepsilon_{11} + z\kappa$, the definition of which are given by

$$\varepsilon = u_s \cdot X_{0,s} + \frac{1}{2} u_s \cdot u_s, \qquad (5)$$

$$\kappa = n_s \cdot (X_{0,s} + u_s) - N_s \cdot X_{0,s}. \qquad (6)$$

The latter is the classical change of curvature.
<figure><img src="image_5.png" /><figcaption>Fig. 1 Kinematic of beam element</figcaption></figure>

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 $\omega$, one has $\omega = \gamma - \Gamma$, where $\Gamma$ is defined to be the angle between the vectors $N_B$ and $N_A$. Similarly, the same holds for $\gamma$ at the deformed configuration. We recall that $\Gamma$ and $\gamma$ 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:

$$ (N_B \times N_A) \cdot e_3 = |N_B||N_A| \sin \Gamma = \sin \Gamma, \qquad (7) $$

$$ N_B \cdot N_A = |N_B||N_A| \cos \Gamma = \cos \Gamma, \qquad (8) $$

$$ (n_B \times n_A) \cdot e_3 = |n_B||n_A| \sin \gamma = \sin \gamma, \qquad (9) $$

$$ n_B \cdot n_A = |n_B||n_A| \cos \gamma = \cos \gamma, \qquad (10) $$

from which it follows

$$ \Gamma = \left( \arcsin \left( (N_B \times N_A) \cdot e_3 \right) - \frac{\pi}{2} \right) \frac{N_B \cdot N_A}{|N_B \cdot N_A|} + \frac{\pi}{2}, \qquad (11) $$

$$ \gamma = \left( \arcsin \left( (n_B \times n_A) \cdot e_3 \right) - \frac{\pi}{2} \right) \frac{n_B \cdot n_A}{|n_B \cdot n_A|} + \frac{\pi}{2} + 2k\pi, \qquad (12) $$
<figure><img src="image_6.png" /></figure>

with k being some integer. The expression for ω now reads

$$ \omega = \left( \arcsin \left( (\mathbf{n}_B \times \mathbf{n}_A) \cdot \mathbf{e}_3 \right) - \frac{\pi}{2} \right) \frac{\mathbf{n}_B \cdot \mathbf{n}_A}{|\mathbf{n}_B \cdot \mathbf{n}_A|} - \left( \arcsin \left( (N_B \times N_A) \cdot \mathbf{e}_3 \right) - \frac{\pi}{2} \right) \frac{N_B \cdot N_A}{|N_B \cdot N_A|} + 2k\pi. \quad (13) $$

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:

$$ \mathbf{u}_A = \mathbf{u}_B. \quad (14) $$

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

$$ U_s = \frac{1}{2} C \omega^2, \quad (15) $$

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

$$ \delta \int_{t_1}^{t_2} \mathcal{L} dt = 0, \quad (16) $$

where t refers to time, $t_1$ and $t_2$ are boundaries of the time interval, $\mathcal{L}$ is defined as the Lagrangian given by

$$ \mathcal{L} = T - \Psi_{\text{int}} - \Psi_{\text{ext}}. \quad (17) $$
<figure><img src="image_7.png" /><figcaption>Fig. 3 Revolute pair with torsional spring</figcaption></figure>

Here T is the kinetic energy of the system and $\Psi_{int}$, $\Psi_{ext}$ are the internal and the external potentials, respectively. The quantities in (17) are defined as

$$T = \frac{1}{2} \int_B \rho \dot{\mathbf{x}} \cdot \dot{\mathbf{x}} dV, \qquad (18)$$

$$\Psi_{int} = \frac{1}{2} \int_B (E E_{11} E_{11}) dV + \frac{1}{2} \sum_{k=1}^{K} C_k \omega_k^2, \qquad (19)$$

$$\Psi_{ext} = - \int_L \mathbf{q}(s) \cdot \mathbf{u}(s) ds - \sum_{i=1}^{N} \mathbf{P}_i \mathbf{u}_i - \sum_{j=1}^{M} M_{3j} \theta_j, \qquad (20)$$

where E is Young's modulus, L is the total length of all members, K is the number of torsional springs, $\mathbf{q}(s)$ is the distributed external force, $\mathbf{P}_i, i = 1, ..., N$ and $\mathbf{u}_i$ are the concentrated forces and the corresponding displacements, respectively, $M_{3j}, j = 1, ..., M$ and $\theta_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

$$\int_L \rho A \ddot{\mathbf{u}} \cdot \delta \mathbf{u} \, ds + \int_L \rho I \ddot{\mathbf{n}} \cdot \delta \mathbf{n} \, ds + \int_L (E A \varepsilon_{11} \delta \varepsilon_{11} + E I \kappa \delta \kappa) \, ds + C \omega \delta \omega \\ - \int_L \mathbf{q}(s) \cdot \delta \mathbf{u}(s) \, ds - \sum_{i=1}^{N} \mathbf{P}_i \cdot \delta \mathbf{u}_i - \sum_{j=1}^{M} M_{3j} \delta \theta_j = 0, \qquad (21)$$

with I and A being the moment of inertia and the cross-section area, respectively.

Since $\mathbf{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 $\delta\theta$ starts from the fact that $\theta$ is the
angle between $N$ and $n$, we have

$$ \sin \theta = (N \times n) \cdot e_3, \qquad (22) $$

$$ \cos \theta \delta\theta = \delta((N \times n) \cdot e_3), \qquad (23) $$

$$ \delta\theta = \frac{\partial((N \times n) \cdot e_3)}{\partial \mathbf{u}_{,s}} \frac{1}{N \cdot n} \cdot \delta\mathbf{u}_{,s}. \qquad (24) $$

Accordingly, the functional in (21) can be written as follows

$$ \int_L \rho A \ddot{u} \cdot \delta \mathbf{u} \, ds + \int_L \rho I \ddot{n} \cdot \frac{\partial \mathbf{n}}{\partial \mathbf{u}_{,s}} \delta \mathbf{u}_{,s} \, ds + \int_L EA \varepsilon_{11} \frac{\partial \varepsilon_{11}}{\partial \mathbf{u}_{,s}} \delta \mathbf{u}_{,s} \, ds \\ + \int_L EI \kappa \left( \frac{\partial \kappa}{\partial \mathbf{u}_{,s}} \delta \mathbf{u}_{,s} + \frac{\partial \kappa}{\partial \mathbf{u}_{,ss}} \delta \mathbf{u}_{,ss} \right) \, ds + \sum_{k=1}^{K} C_k \omega_k \left( \frac{\partial \omega}{\partial \mathbf{u}_{,s}^A} \delta \mathbf{u}_{,s}^A + \frac{\partial \omega}{\partial \mathbf{u}_{,s}^B} \delta \mathbf{u}_{,s}^B \right) \\ - \int_L \mathbf{q}(s) \cdot \delta \mathbf{u}(s) \, ds - \sum_{i=1}^{N} \mathbf{P}_i \cdot \delta \mathbf{u}_i - \sum_{j=1}^{M} M_{3j} \frac{\partial((N \times n) \cdot e_3)}{\partial \mathbf{u}_{,s}} \frac{1}{N \cdot n} \cdot \delta \mathbf{u}_{,s} = 0, \qquad (25) $$

where $\mathbf{u}_{,s}^A$ and $\mathbf{u}_{,s}^B$ define the vectors $n^A$ and $n^B$, respectively. In general, the above equation takes the form

$$ \rho \mathbf{m}(\ddot{u}, \dot{u}, \ddot{u}_{,s}, \dot{u}_{,s}) + \mathbf{k}(\mathbf{u}_{,s}, \mathbf{u}_{,ss}) - \mathbf{f}(\mathbf{u}, \mathbf{u}_{,s}) = \mathbf{0} \qquad (26) $$

where $\mathbf{m}, \mathbf{k}, \mathbf{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 $\mathcal{H} = T + \Psi_{\text{int}} + \Psi_{\text{ext}} = \text{constant}$.

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 formulation and time integration scheme

### 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

$$ \mathbf{u} = N \mathbf{U}, \qquad (27) $$

$$ \mathbf{u}_{,s} = N_{,s} \mathbf{U}, \qquad (28) $$

$$ \mathbf{u}_{,ss} = N_{,ss} \mathbf{U}, \qquad (29) $$

where

$$ N = \begin{pmatrix} N1 & 0 & N2 & 0 & N3 & 0 & N4 & 0 \\ 0 & N1 & 0 & N2 & 0 & N3 & 0 & N4 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \end{pmatrix}, \\ \mathbf{U}^T = (u_1^I \ u_2^I \ u'_1^I \ u'_2^I \ u_1^{II} \ u_2^{II} \ u'_1^{II} \ u'_2^{II}), \qquad (30) $$
and

$$N1 = \frac{\zeta^3}{4} - \frac{3\zeta}{4} + \frac{1}{2}, \qquad (31)$$

$$N2 = \frac{l\zeta^3}{8} - \frac{l\zeta^2}{8} - \frac{l\zeta}{8} + \frac{l}{8}, \qquad (32)$$

$$N3 = \frac{-\zeta^3}{4} + \frac{3\zeta}{4} + \frac{1}{2}, \qquad (33)$$

$$N4 = \frac{l\zeta^3}{8} + \frac{l\zeta^2}{8} - \frac{l\zeta}{8} - \frac{l}{8}, \qquad (34)$$

with $\zeta = [-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–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 $\xi$ to be a scalar which defines any position within the time interval $\Delta t$, with $0 \le \xi \le 1$. We start with the following expressions:

$$\mathbf{x}_{n+\xi} = \xi \mathbf{x}_{n+1} + (1-\xi) \mathbf{x}_{n}, \qquad (35)$$

$$\dot{\mathbf{x}}_{n+\xi} = \frac{\mathbf{x}_{n+1} - \mathbf{x}_{n}}{\Delta t}, \qquad (36)$$

$$\ddot{\mathbf{x}}_{n+\xi} = \frac{\dot{\mathbf{x}}_{n+1} - \dot{\mathbf{x}}_{n}}{\Delta t}. \qquad (37)$$

The first relation defines a convex set, the following two are true for some values of $\xi$. The midpoint rule corresponds to the choice $\xi = 0.5$. Consequently, to integrate $\mathbf{u}$ one has,

$$\dot{\mathbf{u}}_{n+\frac{1}{2}} = \frac{\Delta u}{\Delta t}, \qquad (38)$$

$$\dot{\mathbf{u}}_{n+1} = \frac{2\Delta u}{\Delta t} - \dot{\mathbf{u}}_n, \qquad (39)$$

$$\ddot{\mathbf{u}}_{n+\frac{1}{2}} = \frac{2\Delta u}{\Delta t^2} - \frac{2\dot{\mathbf{u}}_n}{\Delta t}, \qquad (40)$$

where $\Delta u = \mathbf{u}_{n+1} - \mathbf{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 has

$$\dot{\xi} = \frac{\Delta \mathbf{u}_s}{\Delta t} \cdot X_{0,s} + \mathbf{u}_s \cdot \frac{\Delta \mathbf{u}_s}{\Delta t}, \qquad (41)$$

$$\dot{\kappa} = \frac{\partial \kappa}{\partial \mathbf{u}_s} \cdot \frac{\Delta \mathbf{u}_s}{\Delta t} + \frac{\partial \kappa}{\partial \mathbf{u}_{ss}} \cdot \frac{\Delta \mathbf{u}_{ss}}{\Delta t}, \qquad (42)$$

$$\dot{\mathbf{n}} = \frac{\partial \mathbf{n}}{\partial \mathbf{u}_s} \cdot \frac{\Delta \mathbf{u}_s}{\Delta t}, \qquad (43)$$
$$ \dot{\omega} = \frac{\partial \omega}{\partial \mathbf{u}_s^A} \cdot \frac{\Delta \mathbf{u}_s^A}{\Delta t} + \frac{\partial \omega}{\partial \mathbf{u}_s^B} \cdot \frac{\Delta \mathbf{u}_s^B}{\Delta t}, \quad (44) $$

Given the strain fields and the kinematic quantities $n$ and $\omega$ at time $t_n$, the strain fields and the same kinematic quantities at step $n+1$ are then defined as follows:

$$ \varepsilon_{n+\frac{1}{2}} = \varepsilon_n + \frac{1}{2} \Delta t \dot{\varepsilon}_{n+\frac{1}{2}} \quad (45) $$

$$ \kappa_{n+\frac{1}{2}} = \kappa_n + \frac{1}{2} \Delta t \dot{\kappa}_{n+\frac{1}{2}} \quad (46) $$

$$ \ddot{\mathbf{u}}_{n+\frac{1}{2}} = \frac{\dot{\mathbf{u}}_{n+1} - \dot{\mathbf{u}}_n}{\Delta t} \quad (47) $$

$$ \omega_{n+\frac{1}{2}} = \omega_n + \frac{1}{2} \Delta t \dot{\omega}_{n+\frac{1}{2}} \quad (48) $$

With these so defined quantities, the midpoint rule provides energy conservation. Using (48), the value of $\omega$ could be calculated exactly without limitation (for example $2\pi$) because only the difference of $\omega$ is of importance.

### The field equations

To obtain our field equations we start from functional (25)

$$ \int_L \rho A \ddot{\mathbf{u}} \cdot \delta \mathbf{u} \, ds + \int_L \rho I \ddot{\mathbf{n}} \cdot \frac{\partial \mathbf{n}}{\partial \mathbf{u}_s} \delta \mathbf{u}_s \, ds + \int_L EA \varepsilon_{11} \frac{\partial \varepsilon_{11}}{\partial \mathbf{u}_s} \delta \mathbf{u}_s \, ds \\ + \int_L EI \kappa \left( \frac{\partial \kappa}{\partial \mathbf{u}_s} \delta \mathbf{u}_s + \frac{\partial \kappa}{\partial \mathbf{u}_{ss}} \delta \mathbf{u}_{ss} \right) \, ds + \sum_{k=1}^K C_k \omega_k \left( \frac{\partial \omega}{\partial \mathbf{u}_s^A} \delta \mathbf{u}_s^A + \frac{\partial \omega}{\partial \mathbf{u}_s^B} \delta \mathbf{u}_s^B \right) \\ - \int_L \mathbf{q}(s) \cdot \delta \mathbf{u}(s) \, ds - \sum_{i=1}^N \mathbf{P}_i \cdot \delta \mathbf{u}_i - \sum_{j=1}^M M_{3j} \frac{\partial ((N \times \mathbf{n}) \cdot \mathbf{e}_3)}{\partial \mathbf{u}_s} \frac{1}{N \cdot \mathbf{n}} \cdot \delta \mathbf{u}_s = 0. \quad (49) $$

Using the following identities:

$$ \ddot{\mathbf{n}} \cdot \frac{\partial \mathbf{n}}{\partial \mathbf{u}_s} \cdot \delta \mathbf{u}_s = \left( \ddot{\mathbf{n}} \cdot \frac{\partial \mathbf{n}}{\partial \mathbf{u}_s} \cdot \delta \mathbf{u} \right)_s - \left( \ddot{\mathbf{n}} \cdot \frac{\partial \mathbf{n}}{\partial \mathbf{u}_s} \right)_s \cdot \delta \mathbf{u}, \quad (50) $$

$$ \varepsilon \frac{\partial \varepsilon}{\partial \mathbf{u}_s} \cdot \delta \mathbf{u}_s = \left( \varepsilon \frac{\partial \varepsilon}{\partial \mathbf{u}_s} \cdot \delta \mathbf{u} \right)_s - \left( \varepsilon \frac{\partial \varepsilon}{\partial \mathbf{u}_s} \right)_s \cdot \delta \mathbf{u}, \quad (51) $$

$$ \kappa \frac{\partial \kappa}{\partial \mathbf{u}_s} \cdot \delta \mathbf{u}_s = \left( \kappa \frac{\partial \kappa}{\partial \mathbf{u}_s} \cdot \delta \mathbf{u} \right)_s - \left( \kappa \frac{\partial \kappa}{\partial \mathbf{u}_s} \right)_s \cdot \delta \mathbf{u}, \quad (52) $$

$$ \kappa \frac{\partial \kappa}{\partial \mathbf{u}_{ss}} \cdot \delta \mathbf{u}_{ss} = \left( \kappa \frac{\partial \kappa}{\partial \mathbf{u}_{ss}} \cdot \delta \mathbf{u}_s \right)_s - \left( \left( \kappa \frac{\partial \kappa}{\partial \mathbf{u}_{ss}} \right)_s \cdot \delta \mathbf{u} \right)_s + \left( \kappa \frac{\partial \kappa}{\partial \mathbf{u}_{ss}} \right)_{ss} \cdot \delta \mathbf{u}, \quad (53) $$

together with Gauss's theorem, the above functional can be rewritten to give

$$ \int_L \left( \rho A \ddot{\mathbf{u}} - \rho I \left( \ddot{\mathbf{n}} \cdot \frac{\partial \mathbf{n}}{\partial \mathbf{u}_s} \right)_s \right) \cdot \delta \mathbf{u} \, ds \\ - EA \left( \varepsilon \frac{\partial \varepsilon}{\partial \mathbf{u}_s} \right)_s - EI \left( \kappa \frac{\partial \kappa}{\partial \mathbf{u}_s} \right)_s + EI \left( \kappa \frac{\partial \kappa}{\partial \mathbf{u}_{ss}} \right)_{ss} - \mathbf{q}(s) \cdot \delta \mathbf{u} \, ds \\ + \left( \ddot{\mathbf{n}} \cdot \frac{\partial \mathbf{n}}{\partial \mathbf{u}_s} + EA \varepsilon \frac{\partial \varepsilon}{\partial \mathbf{u}_s} + EI \kappa \frac{\partial \kappa}{\partial \mathbf{u}_s} - EI \left( \kappa \frac{\partial \kappa}{\partial \mathbf{u}_{ss}} \right)_{ss} \right) \cdot \delta \mathbf{u} \Big|_0^L $$
$$ + EI \kappa \frac{\partial \kappa}{\partial \mathbf{u}_{,ss}} \cdot \delta \mathbf{u}_{,s} \Big|_{0}^{L} - \sum_{i=1}^{N} \mathbf{P}_{i} \cdot \delta \mathbf{u}_{i} + \sum_{k=1}^{K} C_{k} \omega_{k} \left( \frac{\partial \omega}{\partial \mathbf{u}_{,s}^{A}} \delta \mathbf{u}_{,s}^{A} + \frac{\partial \omega}{\partial \mathbf{u}_{,s}^{B}} \delta \mathbf{u}_{,s}^{B} \right) - \sum_{j=1}^{M} M_{3j} \frac{\partial ((N \times \mathbf{n}) \cdot \mathbf{e}_{3})}{\partial \mathbf{u}_{,s}} \frac{1}{N \cdot \mathbf{n}} \cdot \delta \mathbf{u}_{,sj} = 0. \quad (54) $$

Because the variations are arbitrary, one Euler-Lagrange equation with two boundary conditions follow from the above functional:

$$ \rho A \ddot{u} - \rho I \left( \ddot{n} \cdot \frac{\partial \mathbf{n}}{\partial \mathbf{u}_{,s}} \right)_{,s} - EA \left( \varepsilon \frac{\partial \varepsilon}{\partial \mathbf{u}_{,s}} \right)_{,s} - EI \left( \kappa \frac{\partial \kappa}{\partial \mathbf{u}_{,s}} \right)_{,s} + EI \left( \kappa \frac{\partial \kappa}{\partial \mathbf{u}_{,ss}} \right)_{,ss} = \mathbf{q}(s), \quad (55) $$

$$ \left( \rho I \ddot{u} \cdot \frac{\partial \mathbf{n}}{\partial \mathbf{u}_{,s}} + EA \varepsilon \frac{\partial \varepsilon}{\partial \mathbf{u}_{,s}} + EI \kappa \frac{\partial \kappa}{\partial \mathbf{u}_{,s}} - EI \left( \kappa \frac{\partial \kappa}{\partial \mathbf{u}_{,ss}} \right)_{,s} \right) \Big|_{0}^{L} = \mathbf{P}_{i} \Big|_{0}^{L}, \quad (56) $$

$$ EI \kappa \frac{\partial \kappa}{\partial \mathbf{u}_{,ss}} \Big|_{0}^{L} = - \sum_{k=1}^{K} C_{k} \omega_{k} \left( \frac{\partial \omega}{\partial \mathbf{u}_{,s}^{A}} + \frac{\partial \omega}{\partial \mathbf{u}_{,s}^{B}} \right) + M_{3j} \frac{\partial ((N \times \mathbf{n}) \cdot \mathbf{e}_{3})}{\partial \mathbf{u}_{,s}} \frac{1}{N \cdot \mathbf{n}} \Big|_{0}^{L}. \quad (57) $$

### Proof of conservation of angular momentum

Let's recall that $\mathbf{x}_0 = X_0 + \mathbf{u}$. To prove the conservation of angular momentum we start by taking the cross product of Eqs. (55), (56) with $\mathbf{x}_{0n+\frac{1}{2}}$ and the cross product of Eq. (57) with $\mathbf{x}_{0,sn+\frac{1}{2}}$:

$$ \begin{aligned} & \mathbf{x}_{0n+\frac{1}{2}} \times \rho A \ddot{u}_{n+\frac{1}{2}} - \mathbf{x}_{0n+\frac{1}{2}} \times \rho I \left( \ddot{n}_{n+\frac{1}{2}} \cdot \frac{\partial \mathbf{n}}{\partial \mathbf{u}_{,s}} \right)_{,s} - \mathbf{x}_{0n+\frac{1}{2}} \times EA \left( \varepsilon_{n+\frac{1}{2}} \frac{\partial \varepsilon}{\partial \mathbf{u}_{,s}} \right)_{,s} \\ & \quad - \mathbf{x}_{0n+\frac{1}{2}} \times EI \left( \kappa_{n+\frac{1}{2}} \frac{\partial \kappa}{\partial \mathbf{u}_{,s}} \right)_{,s} + \mathbf{x}_{0n+\frac{1}{2}} \times EI \left( \kappa_{n+\frac{1}{2}} \frac{\partial \kappa}{\partial \mathbf{u}_{,ss}} \right)_{,ss} = \mathbf{x}_{0n+\frac{1}{2}} \times \mathbf{q}(s), \end{aligned} \quad (58) $$

$$ \begin{aligned} & \mathbf{x}_{0n+\frac{1}{2}} \times \rho I \ddot{u}_{n+\frac{1}{2}} \cdot \frac{\partial \mathbf{n}}{\partial \mathbf{u}_{,s}} \Big|_{0}^{L} + \mathbf{x}_{0n+\frac{1}{2}} \times EA \varepsilon_{n+\frac{1}{2}} \frac{\partial \varepsilon}{\partial \mathbf{u}_{,s}} \Big|_{0}^{L} + \mathbf{x}_{0n+\frac{1}{2}} \times EI \kappa_{n+\frac{1}{2}} \frac{\partial \kappa}{\partial \mathbf{u}_{,s}} \Big|_{0}^{L} \\ & \quad - \mathbf{x}_{0n+\frac{1}{2}} \times EI \left( \kappa_{n+\frac{1}{2}} \frac{\partial \kappa}{\partial \mathbf{u}_{,ss}} \right)_{,s} \Big|_{0}^{L} = \mathbf{x}_{0n+\frac{1}{2}} \times \mathbf{P}_{i} \Big|_{0}^{L}, \end{aligned} \quad (59) $$

$$ \begin{aligned} & \mathbf{x}_{0,sn+\frac{1}{2}} \times EI \kappa_{n+\frac{1}{2}} \frac{\partial \kappa}{\partial \mathbf{u}_{,ss}} \Big|_{0}^{L} = - \sum_{k=1}^{K} C_{k} \omega_{k} \left( \mathbf{x}_{0,sn+\frac{1}{2}} \times \frac{\partial \omega}{\partial \mathbf{u}_{,s}^{A}} + \mathbf{x}_{0,sn+\frac{1}{2}} \times \frac{\partial \omega}{\partial \mathbf{u}_{,s}^{B}} \right) \\ & \mathbf{x}_{0,sn+\frac{1}{2}} \times M_{3j} \frac{\partial ((N \times \mathbf{n}) \cdot \mathbf{e}_{3})}{\partial \mathbf{u}_{,s}} \frac{1}{N \cdot \mathbf{n}} \Big|_{0}^{L} \cdot W \end{aligned} \quad (60) $$

Integrating Eq. (58) gives

$$ \int_{L} \mathbf{x}_{0n+\frac{1}{2}} \times \rho A \ddot{u}_{n+\frac{1}{2}} \, ds - \int_{L} \mathbf{x}_{0n+\frac{1}{2}} \times \rho I \left( \ddot{n}_{n+\frac{1}{2}} \cdot \frac{\partial \mathbf{n}}{\partial \mathbf{u}_{,s}} \right)_{,s} \, ds $$
$$ - \int_L \mathbf{x}_{0,n+\frac{1}{2}} \times EA \left( \varepsilon_{n+\frac{1}{2}} \frac{\partial \varepsilon}{\partial \mathbf{u},s} \right)_{,s} \, ds - \int_L \mathbf{x}_{0,n+\frac{1}{2}} \times EI \left( \kappa_{n+\frac{1}{2}} \frac{\partial \kappa}{\partial \mathbf{u},s} \right)_{,s} \, ds \\ + \int_L \mathbf{x}_{0,n+\frac{1}{2}} \times EI \left( \kappa_{n+\frac{1}{2}} \frac{\partial \kappa}{\partial \mathbf{u},ss} \right)_{,ss} \, ds = \int_L \mathbf{x}_{0,n+\frac{1}{2}} \times \mathbf{q}(s) \, ds. \quad (61) $$

Using the following identities:

$$ \mathbf{x}_{0,n+\frac{1}{2}} \times \left( \ddot{\mathbf{u}}_{n+\frac{1}{2}} \cdot \frac{\partial \mathbf{n}}{\partial \mathbf{u},s} \right)_{,s} = \left( \mathbf{x}_{0,n+\frac{1}{2}} \times \ddot{\mathbf{u}}_{n+\frac{1}{2}} \cdot \frac{\partial \mathbf{n}}{\partial \mathbf{u},s} \right)_{,s} - \left( \mathbf{x}_{0,s,n+\frac{1}{2}} \times \ddot{\mathbf{u}}_{n+\frac{1}{2}} \cdot \frac{\partial \mathbf{n}}{\partial \mathbf{u},s} \right)_{,s}, \quad (62) $$

$$ \mathbf{x}_{0,n+\frac{1}{2}} \times \left( \varepsilon_{n+\frac{1}{2}} \frac{\partial \varepsilon}{\partial \mathbf{u},s} \right)_{,s} = \left( \mathbf{x}_{0,n+\frac{1}{2}} \times \varepsilon_{n+\frac{1}{2}} \frac{\partial \varepsilon}{\partial \mathbf{u},s} \right)_{,s} - \left( \mathbf{x}_{0,s,n+\frac{1}{2}} \times \varepsilon_{n+\frac{1}{2}} \frac{\partial \varepsilon}{\partial \mathbf{u},s} \right)_{,s}, \quad (63) $$

$$ \mathbf{x}_{0,n+\frac{1}{2}} \times \left( \kappa_{n+\frac{1}{2}} \frac{\partial \kappa}{\partial \mathbf{u},s} \right)_{,s} = \left( \mathbf{x}_{0,n+\frac{1}{2}} \times \kappa_{n+\frac{1}{2}} \frac{\partial \kappa}{\partial \mathbf{u},s} \right)_{,s} - \left( \mathbf{x}_{0,s,n+\frac{1}{2}} \times \kappa_{n+\frac{1}{2}} \frac{\partial \kappa}{\partial \mathbf{u},s} \right)_{,s}, \quad (64) $$

$$ \begin{aligned} \mathbf{x}_{0,n+\frac{1}{2}} \times \left( \kappa_{n+\frac{1}{2}} \frac{\partial \kappa}{\partial \mathbf{u},ss} \right)_{,ss} = & \left( \mathbf{x}_{0,n+\frac{1}{2}} \times \left( \kappa_{n+\frac{1}{2}} \frac{\partial \kappa}{\partial \mathbf{u},ss} \right)_{,s} \right)_{,s} \\ & - \left( \mathbf{x}_{0,s,n+\frac{1}{2}} \times \left( \kappa_{n+\frac{1}{2}} \frac{\partial \kappa}{\partial \mathbf{u},ss} \right)_{,s} \right)_{,s} \\ & + \mathbf{x}_{0,ss,n+\frac{1}{2}} \times \left( \kappa_{n+\frac{1}{2}} \frac{\partial \kappa}{\partial \mathbf{u},ss} \right)_{,s}. \end{aligned} \quad (65) $$

together with subsequent application of Gauss's theorem, Eq. (61) presents us with the following statement

$$ \begin{aligned} & \int_L \mathbf{x}_{0,n+\frac{1}{2}} \times \rho A \ddot{\mathbf{u}}_{n+\frac{1}{2}} \, ds + \int_L \mathbf{x}_{0,s,n+\frac{1}{2}} \times \rho I \left( \ddot{\mathbf{u}}_{n+\frac{1}{2}} \cdot \frac{\partial \mathbf{n}}{\partial \mathbf{x}_{0,s,n+\frac{1}{2}}} \right) \, ds \\ & + \int_L \mathbf{x}_{0,s,n+\frac{1}{2}} \times EA \left( \varepsilon_{n+\frac{1}{2}} \frac{\partial \varepsilon}{\partial \mathbf{x}_{0,s,n+\frac{1}{2}}} \right) \, ds + \int_L \mathbf{x}_{0,s,n+\frac{1}{2}} \\ & \times EI \left( \kappa_{n+\frac{1}{2}} \frac{\partial \kappa}{\partial \mathbf{x}_{0,s,n+\frac{1}{2}}} \right) \, ds \\ & + \int_L \mathbf{x}_{0,ss,n+\frac{1}{2}} \times EI \left( \kappa_{n+\frac{1}{2}} \frac{\partial \kappa}{\partial \mathbf{x}_{0,ss,n+\frac{1}{2}}} \right) \, ds - \mathbf{x}_{0,n+\frac{1}{2}} \times \rho I \ddot{\mathbf{u}}_{n+\frac{1}{2}} \cdot \frac{\partial \mathbf{n}}{\partial \mathbf{u},s} \Big|_0^L \\ & - \mathbf{x}_{0,n+\frac{1}{2}} \times EA \varepsilon_{n+\frac{1}{2}} \frac{\partial \varepsilon}{\partial \mathbf{u},s} \Big|_0^L - \mathbf{x}_{0,n+\frac{1}{2}} \times EI \kappa_{n+\frac{1}{2}} \frac{\partial \kappa}{\partial \mathbf{u},s} \Big|_0^L \\ & + \mathbf{x}_{0,n+\frac{1}{2}} \times EI \left( \kappa_{n+\frac{1}{2}} \frac{\partial \kappa}{\partial \mathbf{u},ss} \right)_{,s} \Big|_0^L - \mathbf{x}_{0,s,n+\frac{1}{2}} \times EI \kappa_{n+\frac{1}{2}} \frac{\partial \kappa}{\partial \mathbf{u},ss} \Big|_0^L \\ & = \int_L \mathbf{x}_{0,n+\frac{1}{2}} \times \mathbf{q}(s) \, ds. \end{aligned} \quad (66) $$

Using (59) and (60), the last statement results in

$$ \begin{aligned} & \int_L \mathbf{x}_{0,n+\frac{1}{2}} \times \rho A \ddot{\mathbf{u}}_{n+\frac{1}{2}} \, ds + \int_L \mathbf{x}_{0,s,n+\frac{1}{2}} \times \rho I \left( \ddot{\mathbf{u}}_{n+\frac{1}{2}} \cdot \frac{\partial \mathbf{n}}{\partial \mathbf{x}_{0,s,n+\frac{1}{2}}} \right) \, ds \\ & + \int_L \mathbf{x}_{0,s,n+\frac{1}{2}} \times EA \left( \varepsilon_{n+\frac{1}{2}} \frac{\partial \varepsilon}{\partial \mathbf{x}_{0,s,n+\frac{1}{2}}} \right) \, ds \end{aligned} $$
$$ 
\begin{aligned}
&+ \int_L \mathbf{x}_{0,sn+\frac{1}{2}} \times EI \left( \kappa_{n+\frac{1}{2}} \frac{\partial \kappa}{\partial \mathbf{x}_{0,sn+\frac{1}{2}}} \right) ds \\
&+ \int_L \mathbf{x}_{0,ssn+\frac{1}{2}} \times EI \left( \kappa_{n+\frac{1}{2}} \frac{\partial \kappa}{\partial \mathbf{x}_{0,ss}} \right) ds = \int_L \mathbf{x}_{0,n+\frac{1}{2}} \times \mathbf{q}(s) ds + \sum_{i=1}^{N} \mathbf{x}_{0,n+\frac{1}{2}} \times \mathbf{P}_i \\
&- \sum_{k=1}^{K} C_k \omega_k \left( \mathbf{x}_{0,sn+\frac{1}{2}} \times \frac{\partial \omega}{\partial \mathbf{u}_s^A} + \mathbf{x}_{0,sn+\frac{1}{2}} \times \frac{\partial \omega}{\partial \mathbf{u}_s^B} \right) \\
&+ \sum_{j=1}^{M} \mathbf{x}_{0,sn+\frac{1}{2}} \times M_{3j} \frac{\partial ((N \times \mathbf{n}) \cdot \mathbf{e}_3)}{\partial \mathbf{u}_s} \frac{1}{N \cdot \mathbf{n}}.
\end{aligned} 
\quad (67) $$

Now we are going to prove that

$$ \mathbf{x}_{0,s} \times \left( \frac{\partial((N \times \mathbf{n}) \cdot \mathbf{e}_3)}{\partial \mathbf{u}_s} \frac{1}{N \cdot \mathbf{n}} \right) = \mathbf{e}_3. \quad (68) $$

In component form, the normal vector **n** defined in (2) reads

$$ n_j = \mathbf{x}_{0,k,s} \epsilon_{kj3} \frac{1}{|\mathbf{x}_{0,s}|}. \quad (69) $$

Here $\epsilon_{3jk}$ is the permutation symbol defined by $\epsilon_{3jk} = \begin{cases} 1 & \text{for even permutation 3,j,k} \\ -1 & \text{for even permutation 3,j,k} \\ 0 & \text{for two equal values of 3,j,k} \end{cases}$.

Let us introduce $V_0 = (N \times \mathbf{n}) \cdot \mathbf{e}_3$ :

$$ V_0 = N_i n_j \epsilon_{ij3} = N_i \mathbf{x}_{0,k,s} (\delta_{ik} \delta_{jj} - \delta_{ij} \delta_{kj}) \frac{1}{|\mathbf{x}_{0,s}|} = \frac{N_i \mathbf{x}_{0,i,s}}{|\mathbf{x}_{0,s}|}. \quad (70) $$

where $\delta_{ij}$ is the Kronecker delta defined as

$$ \delta_{ij} = \begin{cases} 1 & \text{if } i = j \\ 0 & \text{if } i \neq j \end{cases}. \quad \text{Hence, one has now} \\ \frac{\partial V_0}{\partial \mathbf{x}_{0,i,s}} = \frac{\mathbf{x}_{0,j,s} N_m \mathbf{x}_{0,l,s} \epsilon_{ml3} \epsilon_{ij3}}{|\mathbf{x}_{0,s}|^3}. \quad (71) $$

The term $N \cdot \mathbf{n}$ can be rewritten as follows

$$ N \cdot \mathbf{n} = N_m \mathbf{x}_{0,l,s} \epsilon_{lm3} \frac{1}{|\mathbf{x}_{0,s}|} \quad (72) $$

which, together with (71), provides

$$ \left( \frac{\partial((N \times \mathbf{n}) \cdot \mathbf{e}_3)}{\partial \mathbf{u}_s} \frac{1}{N \cdot \mathbf{n}} \right) \cdot \mathbf{e}_j = \mathbf{x}_{0,k,s} \epsilon_{kj3} \frac{1}{|\mathbf{x}_{0,s}|^2}. \quad (73) $$

With this relation at hand, Eq. (69) gives

$$ \frac{\partial((N \times \mathbf{n}) \cdot \mathbf{e}_3)}{\partial \mathbf{u}_s} \frac{1}{N \cdot \mathbf{n}} = \frac{\mathbf{n}}{|\mathbf{x}_{0,s}|} = \mathbf{e}_3 \times \frac{\mathbf{x}_{0,s}}{|\mathbf{x}_{0,s}|^2}. \quad (74) $$

Finally we get

$$ \mathbf{x}_{0,s} \times \left( \frac{\partial((N \times \mathbf{n}) \cdot \mathbf{e}_3)}{\partial \mathbf{u}_s} \frac{1}{N \cdot \mathbf{n}} \right) = \mathbf{x}_{0,s} \times \left( \mathbf{e}_3 \times \frac{\mathbf{x}_{0,s}}{|\mathbf{x}_{0,s}|^2} \right) = \mathbf{e}_3. \quad (75) $$

Hence, (68) is proven.
Then, we are going to prove that

$$ \sum_{k=1}^{K} C_k \omega_k \left( \mathbf{x}_{0,s_{n+\frac{1}{2}}} \times \frac{\partial \omega}{\partial \mathbf{u}_{s}^{A}} + \mathbf{x}_{0,s_{n+\frac{1}{2}}} \times \frac{\partial \omega}{\partial \mathbf{u}_{s}^{B}} \right) = 0. \quad (76) $$

Consider first $\mathbf{x}_{0,s_{n+\frac{1}{2}}} \times \frac{\partial \omega}{\partial \mathbf{u}_{s}^{A}}$ with $\omega$ defined in (13), one has

$$ \frac{\partial \omega}{\partial \mathbf{u}_{s}^{A}} = \frac{\partial}{\partial \mathbf{u}_{s}^{A}} \left( \arcsin \left( (\mathbf{n}_B \times \mathbf{n}_A) \cdot \mathbf{e}_3 \right) \frac{\mathbf{n}_B \cdot \mathbf{n}_A}{|\mathbf{n}_B \cdot \mathbf{n}_A|} \right) - \frac{\pi}{2} \frac{\partial}{\partial \mathbf{u}_{s}^{A}} \left( \frac{\mathbf{n}_B \cdot \mathbf{n}_A}{|\mathbf{n}_B \cdot \mathbf{n}_A|} \right) \quad (77) $$

and

$$ \frac{\partial}{\partial \mathbf{u}_{s}^{A}} \left( \frac{\mathbf{n}_B \cdot \mathbf{n}_A}{|\mathbf{n}_B \cdot \mathbf{n}_A|} \right) = \frac{\mathbf{n}_B}{|\mathbf{n}_B \cdot \mathbf{n}_A|} \frac{\partial \mathbf{n}_A}{\partial \mathbf{u}_{s}^{A}} - \frac{(\mathbf{n}_B \cdot \mathbf{n}_A)^2}{\sqrt{(|\mathbf{n}_B \cdot \mathbf{n}_A|)^6}} \mathbf{n}_B \frac{\partial \mathbf{n}_A}{\partial \mathbf{u}_{s}^{A}} = 0. \quad (78) $$

Consider next,

$$ \begin{aligned} & \mathbf{x}_{0,s_{n+\frac{1}{2}}} \times \frac{\partial}{\partial \mathbf{u}_{s}^{A}} \left( \arcsin \left( (\mathbf{n}_B \times \mathbf{n}_A) \cdot \mathbf{e}_3 \right) \frac{\mathbf{n}_B \cdot \mathbf{n}_A}{|\mathbf{n}_B \cdot \mathbf{n}_A|} \right) \\ &= \frac{\mathbf{x}_{0,s_{n+\frac{1}{2}}} \times \left( (\mathbf{n}_B \times \frac{\partial \mathbf{n}_A}{\partial \mathbf{u}_{s}^{A}}) \mathbf{e}_3 \right)}{\sqrt{1 - ((\mathbf{n}_B \times \mathbf{n}_A) \cdot \mathbf{e}_3)^2}} \frac{\mathbf{n}_B \cdot \mathbf{n}_A}{|\mathbf{n}_B \cdot \mathbf{n}_A|} \\ &= \frac{\mathbf{n}_B \cdot \mathbf{n}_A}{(\mathbf{n}_B \cdot \mathbf{n}_A)^2} \mathbf{x}_{0,s_{n+\frac{1}{2}}} \times \left( (\mathbf{n}_B \times \frac{\partial \mathbf{n}_A}{\partial \mathbf{u}_{s}^{A}}) \mathbf{e}_3 \right). \end{aligned} \quad (79) $$

Using similar procedure to get the results from (70) to (75) by replacing $N$ by $\mathbf{n}_B$, we proved that

$$ \frac{\mathbf{n}_B \cdot \mathbf{n}_A}{(\mathbf{n}_B \cdot \mathbf{n}_A)^2} \mathbf{x}_{0,s_{n+\frac{1}{2}}} \times \left( \left( \mathbf{n}_B \times \frac{\partial \mathbf{n}_A}{\partial \mathbf{u}_{s}^{A}} \right) \mathbf{e}_3 \right) = \mathbf{e}_3. \quad (80) $$

Similarly for $\mathbf{x}_{0,s_{n+\frac{1}{2}}} \times \frac{\partial \omega}{\partial \mathbf{u}_{s}^{B}}$, we get $\mathbf{x}_{0,s_{n+\frac{1}{2}}} \times \frac{\partial \omega}{\partial \mathbf{u}_{s}^{B}} = -\mathbf{e}_3$. Hence, (76) is proven.

Having Eqs. (75) and (76) at hand, Eq. (67) can be simplified to give

$$ \begin{aligned} & \int_L \mathbf{x}_{0,n+\frac{1}{2}} \times \rho A \ddot{\mathbf{u}}_{n+\frac{1}{2}} \, ds + \int_L \mathbf{x}_{0,s_{n+\frac{1}{2}}} \times \rho I \left( \ddot{\mathbf{u}}_{n+\frac{1}{2}} \cdot \frac{\partial \mathbf{n}}{\partial \mathbf{x}_{0,s_{n+\frac{1}{2}}}} \right) \, ds \\ & + \int_L \mathbf{x}_{0,s_{n+\frac{1}{2}}} \times EA \left( \varepsilon_{n+\frac{1}{2}} \frac{\partial \varepsilon}{\partial \mathbf{x}_{0,s_{n+\frac{1}{2}}}} \right) \, ds + \int_L \mathbf{x}_{0,s_{n+\frac{1}{2}}} \times EI \left( \kappa_{n+\frac{1}{2}} \frac{\partial \kappa}{\partial \mathbf{x}_{0,s_{n+\frac{1}{2}}}} \right) \, ds \\ & + \int_L \mathbf{x}_{0,ss_{n+\frac{1}{2}}} \times EI \left( \kappa_{n+\frac{1}{2}} \frac{\partial \kappa}{\partial \mathbf{x}_{0,ss_{n+\frac{1}{2}}}} \right) \, ds \\ &= \int_L \mathbf{x}_{0,n+\frac{1}{2}} \times \mathbf{q}(s) \, ds + \sum_{i=1}^{N} \mathbf{x}_{0,n+\frac{1}{2}} \times \mathbf{P}_i + \sum_{j=1}^{M} M_{3j}. \end{aligned} \quad (81) $$

To further simplify this expression, we first are going to prove that

$$ \mathbf{x}_{0,s} \times \left( \ddot{\mathbf{n}} \cdot \frac{\partial \mathbf{n}}{\partial \mathbf{x}_{0,s}} \right) = \ddot{\mathbf{n}} \times \left( \mathbf{x}_{0,s} \cdot \frac{\partial \mathbf{n}}{\partial \mathbf{x}_{0,s}} \right). \quad (82) $$
Let's consider the second term in (81). After some lengthy algebraic manipulations, Eq. (69) deliver

$$ \frac{\partial n_k}{\partial x_{0,i,s}} = \frac{\epsilon_{3ik}}{|\mathbf{x}_{0,s}|} - \epsilon_{3jk} \frac{x_{0,j,s} x_{0,i,s}}{|\mathbf{x}_{0,s}|^3}, \quad (83) $$

$$ \left( \ddot{\mathbf{n}} \cdot \frac{\partial \mathbf{n}}{\partial \mathbf{x}_{0,s}} \right) \cdot \mathbf{e}_i = \frac{\partial n_k}{\partial x_{0,i,s}} \ddot{\mathbf{n}}_k = \frac{\epsilon_{3ik}}{|\mathbf{x}_{0,s}|} \ddot{\mathbf{n}}_k - \epsilon_{3jk} \frac{x_{0,j,s} x_{0,i,s} \ddot{\mathbf{n}}_k}{|\mathbf{x}_{0,s}|^3}, \quad (84) $$

$$ \left( \mathbf{x}_{0,s} \times \left( \ddot{\mathbf{n}} \cdot \frac{\partial \mathbf{n}}{\partial \mathbf{x}_{0,s}} \right) \right) \cdot \mathbf{e}_a = \epsilon_{abi} \mathbf{x}_{0,b,s} \frac{\partial n_k}{\partial x_{0,i,s}} \ddot{\mathbf{n}}_k = -\delta_{a3} \frac{x_{0,b,s} \ddot{\mathbf{n}}_b}{|\mathbf{x}_{0,s}|}, \quad (85) $$

where use of $x_{0,3,s} = 0$ has been made. Consider further the expression $\ddot{\mathbf{n}} \times (\mathbf{x}_{0,s} \cdot \frac{\partial \mathbf{n}}{\partial \mathbf{x}_{0,s}})$, which is equivalent to

$$ \left( \ddot{\mathbf{n}} \times \left( \mathbf{x}_{0,s} \cdot \frac{\partial \mathbf{n}}{\partial \mathbf{x}_{0,s}} \right) \right) \cdot \mathbf{e}_a = -\delta_{a3} \frac{x_{0,b,s} \ddot{\mathbf{n}}_b}{|\mathbf{x}_{0,s}|}. \quad (86) $$

From (85) and (86), we conclude:

$$ \mathbf{x}_{0,s} \times \left( \ddot{\mathbf{n}} \cdot \frac{\partial \mathbf{n}}{\partial \mathbf{x}_{0,s}} \right) = \ddot{\mathbf{n}} \times \left( \mathbf{x}_{0,s} \cdot \frac{\partial \mathbf{n}}{\partial \mathbf{x}_{0,s}} \right). \quad (87) $$

Hence (82) is proven.

Consider next

$$ \left( \mathbf{x}_{0,s} \cdot \frac{\partial \mathbf{n}}{\partial \mathbf{x}_{0,s}} \right) \cdot \mathbf{e}_i = \frac{\partial n_k}{\partial x_{0,i,s}} \mathbf{x}_{0,k,s} = \frac{\epsilon_{3ik}}{|\mathbf{x}_{0,s}|} \mathbf{x}_{0,k,s} = -n_i. \quad (88) $$

From (87) and (88) we can write:

$$ \mathbf{x}_{0,s_{n+\frac{1}{2}}} \times \left( \ddot{\mathbf{n}}_{n+\frac{1}{2}} \cdot \frac{\partial \mathbf{n}}{\partial \mathbf{x}_{0,s}} \right) = -\ddot{\mathbf{n}}_{n+\frac{1}{2}} \times \mathbf{n}_{n+\frac{1}{2}}. \quad (89) $$

Consider further the third term in (81) and recall the expression for $\epsilon$ defined in (6), one has

$$ \frac{\partial \epsilon}{\partial \mathbf{u}_{,s}} = \mathbf{x}_{0,s}. \quad (90) $$

Hence,

$$ \mathbf{x}_{0,s} \times \mathrm{EA} \left( \epsilon \frac{\partial \epsilon}{\partial \mathbf{u}_{,s}} \right) = 0. \quad (91) $$

Finally, we are going to prove the following statement

$$ \mathbf{x}_{0,s} \times \left( \kappa \frac{\partial \kappa}{\partial \mathbf{x}_{0,s}} \right) + \mathbf{x}_{0,ss} \times \left( \kappa \frac{\partial \kappa}{\partial \mathbf{x}_{0,ss}} \right) = 0. \quad (92) $$

Consider the last two integrals in (81). The expression for $\kappa$, found in (6), can be rewritten as follows:

$$ \kappa = \epsilon_{3ji} \frac{x_{0,i,s} x_{0,j,ss}}{|\mathbf{x}_{0,s}|} - N_{i,s} X_{0,i,s}, \quad (93) $$

which gives

$$ \frac{\partial \kappa}{\partial x_{0,i,s}} = \frac{\epsilon_{3jl} x_{0,j,ss}}{|\mathbf{x}_{0,s}|} - \frac{\epsilon_{3ji} x_{0,i,s} x_{0,j,ss}}{|\mathbf{x}_{0,s}|^3} x_{0,i,s} \quad \text{and} \quad \frac{\partial \kappa}{\partial x_{0,ss,l}} = \frac{\epsilon_{3li} x_{0,i,s}}{|\mathbf{x}_{0,s}|}. \quad (94) $$
The term $\mathbf{x}_{0,s} \times \left(\kappa \frac{\partial \kappa}{\partial \mathbf{x}_{0,s}}\right)$ results in

$$ \left(\mathbf{x}_{0,s} \times \left(\kappa \frac{\partial \kappa}{\partial \mathbf{x}_{0,s}}\right)\right) \cdot \mathbf{e}_a = \mathbf{x}_{0,m,s} \frac{\partial \kappa}{\partial \mathbf{x}_{0,l,s}} \epsilon_{mla} = \delta_{a3} \frac{\mathbf{x}_{0,j,s} \mathbf{x}_{0,j,ss}}{|\mathbf{x}_{0,s}|}. $$

Further, for $\mathbf{x}_{0,ss} \times \left(\kappa \frac{\partial \kappa}{\partial \mathbf{x}_{0,ss}}\right)$, it follows

$$ \left(\mathbf{x}_{0,ss} \times \left(\kappa \frac{\partial \kappa}{\partial \mathbf{x}_{0,ss}}\right)\right) \cdot \mathbf{e}_a = \mathbf{x}_{0,ss,m} \frac{\partial \kappa}{\partial \mathbf{x}_{0,ssl}} \epsilon_{mla} = -\delta_{3a} \frac{\mathbf{x}_{0,ss,i} \mathbf{x}_{0,i,s}}{|\mathbf{x}_{0,s}|}. \quad (95) $$

Hence, Eqs. (95) and (95) mean nothing but

$$ \mathbf{x}_{0,s} \times \left(\kappa \frac{\partial \kappa}{\partial \mathbf{x}_{0,s}}\right) + \mathbf{x}_{0,ss} \times \left(\kappa \frac{\partial \kappa}{\partial \mathbf{x}_{0,ss}}\right) = 0. \quad (96) $$

Using (96), (91), (89) in (81), we get

$$ \begin{aligned} & \rho A \int_L \mathbf{x}_{0n+\frac{1}{2}} \times \dot{\mathbf{u}}_{n+\frac{1}{2}} \, ds + \rho I \int_L \mathbf{n}_{n+\frac{1}{2}} \times \dot{\mathbf{u}}_{n+\frac{1}{2}} \, ds \\ & = \int_L \mathbf{x}_{0n+\frac{1}{2}} \times \mathbf{q}(s) \, ds + \sum_{i=1}^N \mathbf{x}_{0n+\frac{1}{2}} \times \mathbf{P}_i + \sum_{j=1}^M M_{3j}. \end{aligned} \quad (97) $$

which is equivalent to

$$ \begin{aligned} & \rho A \int_L \left( \mathbf{x}_{0n+\frac{1}{2}} \times \ddot{\mathbf{u}}_{n+\frac{1}{2}} + \dot{\mathbf{x}}_{0n+\frac{1}{2}} \times \dot{\mathbf{u}}_{n+\frac{1}{2}} \right) \, ds \\ & \quad + \rho I \int_L \left( \mathbf{n}_{n+\frac{1}{2}} \times \ddot{\mathbf{u}}_{n+\frac{1}{2}} + \dot{\mathbf{n}}_{n+\frac{1}{2}} \times \dot{\mathbf{u}}_{n+\frac{1}{2}} \right) \, ds \\ & = \int_L \mathbf{x}_{0n+\frac{1}{2}} \times \mathbf{q}(s) \, ds + \sum_{i=1}^N \mathbf{x}_{0n+\frac{1}{2}} \times \mathbf{P}_i + \sum_{j=1}^M M_{3j}. \end{aligned} \quad (98) $$

since $\dot{\mathbf{x}}_{0n+\frac{1}{2}} \times \dot{\mathbf{u}}_{n+\frac{1}{2}} = \mathbf{0}$, $\dot{\mathbf{n}}_{n+\frac{1}{2}} \times \dot{\mathbf{n}}_{n+\frac{1}{2}} = \mathbf{0}$, together with the use of the midpoint rule, the expression reads

$$ \begin{aligned} & \rho A \int_L \left( \frac{\mathbf{x}_{0n+1} + \mathbf{x}_{0n}}{2} \times \frac{\dot{\mathbf{u}}_{n+1} - \dot{\mathbf{u}}_n}{\Delta t} + \frac{\mathbf{x}_{0n+1} - \mathbf{x}_{0n}}{\Delta t} \times \frac{\dot{\mathbf{u}}_{n+1} + \dot{\mathbf{u}}_n}{2} \right) \, ds \\ & \quad + \rho I \int_L \left( \frac{\mathbf{n}_{n+1} + \mathbf{n}_n}{2} \times \frac{\dot{\mathbf{u}}_{n+1} - \dot{\mathbf{u}}_n}{\Delta t} \frac{\mathbf{n}_{n+1} - \mathbf{n}_n}{\Delta t} \times \frac{\dot{\mathbf{u}}_{n+1} + \dot{\mathbf{u}}_n}{2} \right) \, ds \\ & = \int_L \mathbf{x}_{0n+\frac{1}{2}} \times \mathbf{q}(s) \, ds + \sum_{i=1}^N \mathbf{x}_{0n+\frac{1}{2}} \times \mathbf{P}_i + \sum_{j=1}^M M_{3j}. \end{aligned} \quad (99) $$

The expression is equivalent to the concise one

$$ \begin{aligned} & \rho A \int_L (\mathbf{x}_{0n+1} \times \dot{\mathbf{u}}_{n+1} - \mathbf{x}_{0n} \times \dot{\mathbf{u}}_n) \, ds + \rho I \int_L (\mathbf{n}_{n+1} \times \dot{\mathbf{u}}_{n+1} - \mathbf{n}_n \times \dot{\mathbf{u}}_n) \, ds \\ & = \int_L \mathbf{x}_{0n+\frac{1}{2}} \times \mathbf{q}(s) \, ds + \sum_{i=1}^N \mathbf{x}_{0n+\frac{1}{2}} \times \mathbf{P}_i + \sum_{j=1}^M M_{3j}. \end{aligned} \quad (100) $$

which indicates that in the case of vanishing moments, we immediately infer

$$ \mathcal{J}_{n+1} - \mathcal{J}_n = \mathbf{0}. \quad (101) $$

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. (21), one can choose the variation to coincide with the velocity vector at time $n + \frac{1}{2}$. Functional (21) takes now the form

$$ 
\begin{aligned}
& \int_L \rho A \ddot{u}_{n+\frac{1}{2}} \cdot \dot{u}_{n+\frac{1}{2}} \, ds + \int_L \rho I \ddot{u}_{n+\frac{1}{2}} \cdot \dot{u}_{n+\frac{1}{2}} \, ds \\
& \qquad + \int_L EA \varepsilon_{n+\frac{1}{2}} \dot{\varepsilon}_{n+\frac{1}{2}} \, ds + \int_L EI \kappa_{n+\frac{1}{2}} \dot{\kappa}_{n+\frac{1}{2}} \, ds + \sum_{k=1}^K C_k \omega_{k n+\frac{1}{2}} \dot{\omega}_{k n+\frac{1}{2}} \\
& = \int_L \boldsymbol{q}(s) \cdot \dot{u}_{n+\frac{1}{2}} \, ds + \sum_{i=1}^N \boldsymbol{P}_i \cdot \dot{u}_{n+\frac{1}{2}} + \sum_{j=1}^M M_{3j} \dot{\theta}_{n+\frac{1}{2}},
\end{aligned} 
\quad (102) $$

which, with the use of the midpoint rule, is equivalent to

$$ 
\begin{aligned}
& \int_L \rho A \frac{(\dot{u}_{n+1} - \dot{u}_n) \cdot (\dot{u}_{n+1} + \dot{u}_n)}{2 \Delta t} \, ds + \int_L \rho I \frac{(\dot{u}_{n+1} - \dot{u}_n) \cdot (\dot{u}_{n+1} + \dot{u}_n)}{2 \Delta t} \, ds \\
& \qquad + \int_L EA \frac{(\varepsilon_{n+1} - \varepsilon_n)(\varepsilon_{n+1} + \varepsilon_n)}{2 \Delta t} \, ds + \int_L EI \frac{(\kappa_{n+1} - \kappa_n)(\kappa_{n+1} + \kappa_n)}{2 \Delta t} \, ds \\
& \qquad + \sum_{k=1}^K C_k \frac{(\omega_{k n+1} - \omega_{k n})(\omega_{k n+1} + \omega_{k n})}{2 \Delta t} \\
& = \int_L \boldsymbol{q}(s) \cdot \frac{\boldsymbol{u}_{n+1} - \boldsymbol{u}_n}{\Delta t} \, ds + \sum_{i=1}^N \boldsymbol{P}_i \cdot \frac{\boldsymbol{u}_{i n+1} - \boldsymbol{u}_{i n}}{\Delta t} + \sum_{j=1}^M M_{3j} \frac{\theta_{j n+1} - \theta_{j n}}{\Delta t}.
\end{aligned} 
\quad (103) $$

The last equation leads to

$$ 
\begin{aligned}
& \frac{1}{2} \int_L \rho A (\dot{u}_{n+1}^2 - \dot{u}_n^2) \, ds + \frac{1}{2} \int_L \rho I (\dot{u}_{n+1}^2 - \dot{u}_n^2) \, ds \\
& \qquad + \frac{1}{2} \int_L EA (\varepsilon_{n+1}^2 - \varepsilon_n^2) \, ds + \frac{1}{2} \int_L EI (\kappa_{n+1}^2 - \kappa_n^2) \, ds + \frac{1}{2} \sum_{k=1}^K C_k (\omega_{k n+1}^2 - \omega_{k n}^2) \\
& = \int_L \boldsymbol{q}(s) \cdot (\boldsymbol{u}_{n+1} - \boldsymbol{u}_n) \, ds + \sum_{i=1}^N \boldsymbol{P}_i \cdot (\boldsymbol{u}_{i n+1} - \boldsymbol{u}_{i n}) + \sum_{j=1}^M M_{3j} (\theta_{j n+1} - \theta_{j n}).
\end{aligned} 
\quad (104) $$

Assuming that the external loading is conservative, from Eq. (104) we finally obtain

$$ [T + \Psi_{int} + \Psi_{ext}]_{n+1} = [T + \Psi_{int} + \Psi_{ext}]_n. \quad (105) $$

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.
<figure><img src="image_8.png" /><figcaption>**Fig. 4** Flexible two-body system</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 5** Loading history</figcaption></figure>

## 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.

### Parameters:

Beam length $L = 5$ m, Height $H = 4$ m

Cross section area $A = 0.3 \text{ m}^2$, Cross section inertia $I = 0.0022 \text{ m}^4$

Young's Modulus $E = 200 \text{ GPa}$, Density $\rho = 7850 \text{ kg/m}^3$

Number of element = 6, Number of steps = $10^6$

Time increment $\Delta t = 10^{-3} \text{ s}$

Figure 6 shows the number of iteration necessary in order for the Newton-Raphson method to converge. For $\Delta t = 10^{-2} \text{ s}$ no more than 5 iterations are needed to achieve convergence whereas for $\Delta t = 10^{-3} \text{ s}$ 3 iterations are sufficient to reach dynamic equilibrium.
<figure><img src="image_10.png" /><figcaption>Fig. 6 Number of iterations until convergence</figcaption></figure>

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: $\Delta t \le 10^{-3}$, $\Delta t = 10^{-2}$, $\Delta t = 0.1$ and $\Delta t = 0.2$. For $\Delta t \le 10^{-3}$ and $\Delta t = 1E-2$, the results are similar and only a slight deviation can be observed for $\Delta t = 0.1$ and $\Delta t = 0.2$. However, in terms of energy conservation, Fig. 8 shows that numerically the energy conservation is achieved with high degree of accuracy for $\Delta t \le 10^{-3}$ and $\Delta t = 10^{-2}$ whereas for larger time steps ($\Delta t \ge 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 $\Delta 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 $\Delta t = 10^{-3}$ s.

### Parameters:

Beam length $L = 5$ m, Height $H = 4$ m

Cross section area $A = 0.1 \text{ m}^2$, Cross section inertia $I = 8.33 \times 10^{-5} \text{ m}^4$
<figure><img src="image_11.png" /><figcaption>Fig. 7 Displacements for different Δt</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 8 Energy history for different dt</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Fig. 9 Vertical displacement of second body's edge for $\Delta t = 10^{-2}$ s</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 10 Energy for $\Delta t = 10^{-2}$ s</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>**Fig. 11** Multi-body figure</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Fig. 12** Loading history</figcaption></figure>

Young's Modulus $E$ = 200 GPa, Density $ρ$ = 7850 kg/m³

Number of elements = 6, Time increment $Δt$ = $10^{-3}$ s

Number of steps = $10^6$, Spring constant $C$ = $10^4$ Nm/rad

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
<figure><img src="image_17.png" /><figcaption>Fig. 13 Horizontal displacement versus time</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Fig. 14 Vertical displacement versus time</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>Fig. 15 Energy time history</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Fig. 16 Spring rotation angle time history</figcaption></figure>

integration method (Eqs. (44) and (48)) where only the increment of $ω$ is of importance. The energy also suffers from the blow-up phenomena as shown Fig. 20.

## Free moving multi-body system

To further assess the performance of the formulation, we consider a free moving multi-body system. More specifically, we investigate the motion of a two-elements mechanism
<figure><img src="image_21.png" /><figcaption>**Fig. 17** Multi-body snap-shots</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>**Fig. 18** Displacement time history—classical midpoint rule for $\Delta t = 10^{-3}$ s</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>Fig. 19 Spring rotation angle time history—classical midpoint rule for $\Delta t = 10^{-3}$ s</figcaption></figure>

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.

P

### Parameters:

Beam length $L = 5$ m, Height $H = 4.8$ m

Cross section area $A = 0.05$ m², Cross section inertia $I = 1.042 \times 10^{-5}$ m⁴

Young's Modulus $E = 200$ GPa, Density $\rho = 7850$ kg/m³

Number of elements = 6, Time increment $\Delta t = 10^{-4}$ s

Number of steps = $10^6$, Spring constant $C = 10^4$ Nm/rad

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
<figure><img src="image_24.png" /><figcaption>Fig. 20 Energy time history—classical midpoint rule for $\Delta t = 10^{-3}$ s</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>Fig. 21 Multi-body figure</figcaption></figure>
<figure><img src="image_26.png" /><figcaption>**Fig. 22 Loading history**</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>**Fig. 23 Energy history**</figcaption></figure>

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 out-of-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.

### Parameters:

Beam length $L_1$ = 3 m, $L_2$ = 10 m, Height $H$ = 2 m

Cross section area $A$ = 0.3 m², Cross section inertia $I$ = 0.0022 m⁴

Young's Modulus $E$ = 200 GPa, Density $ρ$ = 7850 kg/m³

Number of elements = 6, Time increment $Δt$ = 10⁻⁴ s

Number of steps = 10⁶
<figure><img src="image_28.png" /><figcaption>Fig. 24 Momentum history</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>Fig. 25 Angular momentum history</figcaption></figure>
<figure><img src="image_30.png" /><figcaption>**Fig. 26** Spring rotation angle history</figcaption></figure>

<figure><img src="image_31.png" /><figcaption>**Fig. 27** Multi-body snap-shots</figcaption></figure>
<figure><img src="image_32.png" /><figcaption>Fig. 28 Multi-body system</figcaption></figure>

<figure><img src="image_33.png" /><figcaption>Fig. 29 Loading history</figcaption></figure>

<figure><img src="image_34.png" /><figcaption>Fig. 30 Energy history</figcaption></figure>
<figure><img src="image_35.png" /><figcaption>**Fig. 31** Horizontal displacement of the slider</figcaption></figure>

<figure><img src="image_36.png" /><figcaption>**Fig. 32** Multi-body snap-shots</figcaption></figure>

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.
<figure><img src="image_37.png" /><figcaption>Fig. 33 Horizontal displacement</figcaption></figure>

<figure><img src="image_38.png" /><figcaption>Fig. 34 Vertical displacement</figcaption></figure>
## Conclusions

A new approach for geometrically exact Euler-Bernoulli beam in planar multi-body dynamics has been proposed. A new energy-momentum method was designed for long-term dynamics. The advantage of a displacement-only formulation was clearly demonstrated. While the present paper considers only revolute joints, the methodology is general and could be applied to other types of joints as well. Also, the planar nature of the formulations in this paper is not a restriction. A fully three-dimensional formulation is being currently developed. Also a possible modification of the present formulation to incorporate high-frequency damping can be easily implemented, an aspect which was out of the scope in this paper (Figs. 33 and 34).

## Authors' contributions

CS developed the theoretical formulation as well as the time integration scheme. He contributed to the formal proof of the conservation properties of the scheme. TLN and SC implemented the computational procedures and ran several challenging examples. MH has substantially contributed to the formulation. TLN and HM have been involved in drafting the paper. All authors read and approved the final manuscript.

## Availability of data and materials

Data sharing is not applicable to this article as no datasets were generated or analysed during the current study.

## Competing interests

The authors declare that they have no competing interests.

## Funding

The second author received financial support from the Institut des Sciences Appliquées de Rennes through the Dean Ph.D. fellowship.

## Appendix A. Linearisation of functional

In order to solve (25), the discretisation in time and space as presented in “Finite element formulation and time integration scheme” section and the numerical linearisation have been used. To compute the integral, we use the classical Gaussian quadrature. Consider one beam element, Eq. (25) is evaluated at $t_{n+\frac{1}{2}}$

$$ G = \int_{l} \rho A \ddot{u} \cdot \delta \mathbf{u} \, ds + \int_{l} \rho I \ddot{\mathbf{n}} \cdot \frac{\partial \mathbf{n}}{\partial \mathbf{u}_s} \delta \mathbf{u}_s \, ds + \int_{l} EA \varepsilon_{11} \frac{\partial \varepsilon_{11}}{\partial \mathbf{u}_s} \delta \mathbf{u}_s \, ds \\ + \int_{l} EI\kappa \left( \frac{\partial \kappa}{\partial \mathbf{u}_s} \delta \mathbf{u}_s + \frac{\partial \kappa}{\partial \mathbf{u}_{ss}} \delta \mathbf{u}_{ss} \right) \, ds + \sum_{k=1}^{K} C_k \omega_k \frac{\partial \omega}{\partial \mathbf{u}_s} \delta \mathbf{u}_s \\ - \int_{L} \mathbf{q}(s) \cdot \delta \mathbf{u}(s) \, ds - \sum_{i=1}^{N} P_i \cdot \delta \mathbf{u}_i \\ - \sum_{j=1}^{M} M_{3j} \frac{\partial ((N \times \mathbf{n}) \cdot \mathbf{e}_3)}{\partial \mathbf{u}_s} \frac{1}{N \cdot \mathbf{n}} \cdot \delta \mathbf{u}_s \Big|_{n+\frac{1}{2}} = 0. \qquad (A.1) $$
where *l* is the initial length of the element. The latter equation can be rewritten as follows

$$
\begin{equation}
\begin{aligned}
G = & \int_l \rho A \ddot{u} \cdot \delta \mathbf{U} N \, ds + \int_l \rho I \ddot{n} \cdot \frac{\partial \mathbf{n}}{\partial \mathbf{u}_s} \delta \mathbf{U} N_{,s} \, ds + \int_l EA \varepsilon_{11} \frac{\partial \varepsilon_{11}}{\partial \mathbf{u}_s} \delta \mathbf{U} N_{,s} \, ds \\
& + \int_l EI \kappa \left( \frac{\partial \kappa}{\partial \mathbf{u}_s} \delta \mathbf{U} N_{,s} + \frac{\partial \kappa}{\partial \mathbf{u}_{ss}} \delta \mathbf{U} N_{,ss} \right) \, ds + \sum_{k=1}^K C_k \omega_k \frac{\partial \omega}{\partial \mathbf{u}_s} \delta \mathbf{U} N_{,s} \\
& - \int_l \mathbf{q}(s) \cdot \delta \mathbf{U} N \, ds - \sum_{i=1}^N \mathbf{P}_i \cdot \delta \mathbf{U} N \\
& - \sum_{j=1}^M M_{3j} \frac{\partial \left( (N \times \mathbf{n}) \cdot \mathbf{e}_3 \right)}{\partial \mathbf{u}_s} \frac{1}{N \cdot \mathbf{n}} \cdot \delta \mathbf{U} N_{,s} \bigg|_{n+\frac{1}{2}} = 0.
\end{aligned}
\tag{A.2}
\end{equation}
$$

Let's decompose $G$ as follows $G = \rho m + k - f = 0$, where

$$
m = \int_l A N^T \left( 2 \frac{\Delta \mathbf{u}}{\Delta t^2} - 2 \frac{\dot{\mathbf{u}}_n}{\Delta t} \right) ds + \int_l I N_{,s}^T \left( 2 \frac{\Delta \mathbf{n}}{\Delta t^2} - 2 \frac{\dot{\mathbf{n}}_n}{\Delta t} \right) \cdot \frac{\partial \mathbf{n}}{\partial \mathbf{u}_s} ds. \quad (A.3)
$$

$$
\begin{equation}
\begin{aligned}
k = & \int_l EA N_{,s}^T \varepsilon_{11} \frac{\partial \varepsilon_{11}}{\partial \mathbf{u}_s} \, ds + \int_l EI \left( N_{,s}^T \kappa \frac{\partial \kappa}{\partial \mathbf{u}_s} + N_{,ss}^T \kappa \frac{\partial \kappa}{\partial \mathbf{u}_{ss}} \right) \, ds \\
& + \sum_{k=1}^K N_{,s}^T C_k \omega_k \frac{\partial \omega}{\partial \mathbf{u}_s}.
\end{aligned}
\tag{A.4}
\end{equation}
$$

$$
f = \int_l N^T \mathbf{q}(s) \cdot \mathrm{d}s + \sum_{i=1}^{N} N^T \mathbf{P}_i + \sum_{j=1}^{M} M_{3j} N_{,s}^T \frac{\partial \left( (N \times \mathbf{n}) \cdot \mathbf{e}_3 \right)}{\partial \mathbf{u}_s} \frac{1}{N \cdot \mathbf{n}}. \quad (A.5)
$$

In Eq. (A.2), $\Delta \mathbf{U}$ is the unknown as every thing is known at time step $n$. The linerization is made by taking the variation of $G$ with respect to $\Delta \mathbf{U}$ as follows

$$
\Delta \mathbf{m} = \int_l A N^T 2 \frac{1}{dt^2} N \delta \Delta \mathbf{U} \, ds \\
+ \int_l I N_{,s}^T \frac{\partial}{\partial \Delta \mathbf{u}_s} \left( \left( 2 \frac{\Delta \mathbf{n}}{\Delta t^2} - 2 \frac{\dot{\mathbf{n}}_n}{\Delta t} \right) \cdot \frac{\partial \mathbf{n}}{\partial \mathbf{u}_s} \right) N_{,s} \delta \Delta \mathbf{U} \, ds = \mathbf{M} \delta \Delta \mathbf{U}. \quad (A.6)
$$

$$
\begin{align}
\Delta \mathbf{k} = & \int_l EA N_{,s}^T \frac{\partial}{\partial \Delta \mathbf{u}_{,s}} \left( \varepsilon_{11} \frac{\partial \varepsilon_{11}}{\partial \mathbf{u}_{,s}} \right) N_{,s} \delta \Delta \mathbf{U} \, ds \nonumber \\
& + \int_l EI \left( N_{,s}^T \frac{\partial}{\partial \Delta \mathbf{u}_{,s}} \left( \kappa \frac{\partial \kappa}{\partial \mathbf{u}_{,s}} \right) N_{,s} + N_{,s}^T \frac{\partial}{\partial \Delta \mathbf{u}_{,s}} \left( \kappa \frac{\partial \kappa}{\partial \mathbf{u}_{,ss}} \right) N_{,ss} \right. \nonumber \\
& \left. + N_{,ss}^T \frac{\partial}{\partial \Delta \mathbf{u}_{,s}} \left( \kappa \frac{\partial \kappa}{\partial \mathbf{u}_{,ss}} \right) N_{,s} + N_{,ss}^T \frac{\partial}{\partial \Delta \mathbf{u}_{,ss}} \left( \kappa \frac{\partial \kappa}{\partial \mathbf{u}_{,ss}} \right) N_{,ss} \right) \delta \Delta \mathbf{U} \, ds \nonumber \\
& + \sum_{k=1}^K N_{,s}^T \frac{\partial}{\partial \Delta \mathbf{u}_{,ss}} \left( C_k \omega_k \frac{\partial \omega}{\partial \mathbf{u}_{,s}} \right) N_{,s} \delta \Delta \mathbf{U} \nonumber \\
& = \mathbf{K} \delta \Delta \mathbf{U}. \tag{A.7}
\end{align}
$$

$$
\Delta \mathbf{f} = \sum_{j=1}^{M} M_{3j} N_{,s}^T \frac{\partial \left( \frac{\partial \left( (N \times \mathbf{n}) \cdot \mathbf{e}_3 \right)}{\partial \mathbf{u}_s} \frac{1}{N \cdot \mathbf{n}} \right)}{\partial \Delta \mathbf{u}_{,ss}} N_{,s} \delta \Delta \mathbf{U} = \mathbf{F} \delta \Delta \mathbf{U}. \quad (A.8)
$$

Thus, the tangent matrix is $W_{element} = \rho M + K - F$. Finally, one can write $W_{element} \delta \Delta \mathbf{U} = (\rho \mathbf{m} + \mathbf{k} - \mathbf{f})$. Practically, the residual vector and tangent matrix containing the torsional springs could be calculated directly at global level not at element
level as follows

$$
\begin{align}
\sum_{k=1}^{K} C_k \omega_k \delta \omega_k &= \sum_{k=1}^{K} C_k \omega_k \left( \frac{\partial \omega}{\partial \mathbf{u}_{,s}^A} \delta \mathbf{u}_{,s}^A + \frac{\partial \omega}{\partial \mathbf{u}_{,s}^B} \delta \mathbf{u}_{,s}^B \right) \nonumber \\
&= \sum_{k=1}^{K} C_k \omega_k \left( \frac{\partial \omega}{\partial \mathbf{u}_{,s}^A} N_{,s} T_1 \delta \mathbf{\Phi} + \frac{\partial \omega}{\partial \mathbf{u}_{,s}^B} N_{,s} T_2 \delta \mathbf{\Phi} \right). \tag{A.9}
\end{align}
$$

where $\mathbf{\Phi}$ 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

$$
s_p = \sum_{k=1}^{K} \left( T_1^T N_{,s}^T C_k \omega_k \frac{\partial \omega}{\partial \mathbf{u}_{,s}^A} + T_2^T N_{,s}^T C_k \omega_k \frac{\partial \omega}{\partial \mathbf{u}_{,s}^B} \right). \quad (A.10)
$$

Thus, the corresponding global tangent matrix reads

$$
\begin{equation}
\begin{split}
S_p = \sum_{k=1}^{K} & \left( T_1^T N_{,s}^T C_k \frac{\partial}{\partial \Delta \mathbf{u}_{,s}^A} \left( \omega_k \frac{\partial \omega}{\partial \mathbf{u}_{,s}^A} \right) N_{,s} T_1 + T_1^T N_{,s}^T C_k \frac{\partial}{\partial \Delta \mathbf{u}_{,s}^B} \left( \omega_k \frac{\partial \omega}{\partial \mathbf{u}_{,s}^B} \right) N_{,s} T_2 \right. \\
& \left. + T_2^T N_{,s}^T C_k \frac{\partial}{\partial \Delta \mathbf{u}_{,s}^A} \left( \omega_k \frac{\partial \omega}{\partial \mathbf{u}_{,s}^B} \right) N_{,s} T_1 + T_2^T N_{,s}^T C_k \frac{\partial}{\partial \Delta \mathbf{u}_{,s}^B} \left( \omega_k \frac{\partial \omega}{\partial \mathbf{u}_{,s}^B} \right) N_{,s} T_2 \right).
\end{split}
\tag{A.11}
\end{equation}
$$

**Appendix B. Proof of conservation of linear momentum**

To prove the conservation of linear momentum, we start by integrating the field Eq. (55)
which provides us with the statement

$$
\begin{equation}
\begin{split}
& \int_L \rho A \ddot{u}_{n+\frac{1}{2}} \, ds - \int_L \left( \rho I \left( \ddot{u}_{n+\frac{1}{2}} \cdot \frac{\partial \mathbf{n}}{\partial \mathbf{u}_{,s}} \right)_{,s} - EA \left( \varepsilon_{n+\frac{1}{2}} \frac{\partial \varepsilon}{\partial \mathbf{u}_{,s}} \right)_{,s} \right) \, ds \\
& \quad - \int_L \left( EI \left( \kappa_{n+\frac{1}{2}} \frac{\partial \kappa}{\partial \mathbf{u}_{,s}} \right)_{,s} + EI \left( \kappa_{n+\frac{1}{2}} \frac{\partial \kappa}{\partial \mathbf{u}_{,ss}} \right)_{,ss} \right) \, ds = \int_L \mathbf{q}(s) \, ds.
\end{split}
\tag{B.1}
\end{equation}
$$

The subsequent application of Gauss's theorem results in

$$
\begin{equation}
\begin{split}
& \int_L \rho A \ddot{u}_{n+\frac{1}{2}} \, ds - \left( \rho I \ddot{u}_{n+\frac{1}{2}} \cdot \frac{\partial \mathbf{n}}{\partial \mathbf{u}_{,s}} - EA \varepsilon_{n+\frac{1}{2}} \frac{\partial \varepsilon}{\partial \mathbf{u}_{,s}} \right. \\
& \qquad \left. - EI \kappa_{n+\frac{1}{2}} \frac{\partial \kappa}{\partial \mathbf{u}_{,s}} + EI \left( \kappa_{n+\frac{1}{2}} \frac{\partial \kappa}{\partial \mathbf{u}_{,ss}} \right)_{,s} \right) \Big|_0^L = \int_L \mathbf{q}(s) \, ds.
\end{split}
\tag{B.2}
\end{equation}
$$

A view on the boundary condition (56), reveals the general momentum equation

$$
\int_L \rho A \ddot{u}_{n+\frac{1}{2}} \, ds = \int_L \mathbf{q}(s) \, ds + \sum_{i=1}^{N} \mathbf{P}_i = 0, \quad (B.3)
$$

which indicates that in the case of vanishing external forces

$$
\int_L \mathbf{q}(s) \, ds + \sum_i \mathbf{P}_i = 0, \qquad (B.4)
$$
the midpoint rule provides us with the statement

$$ \int_L \rho A \frac{\dot{u}_{n+1} - \dot{u}_n}{\Delta t} ds = 0, \qquad (B.5) $$

leading to

$$ \int_L \rho A (\dot{u}_{n+1} - \dot{u}_n) ds = 0, \qquad (B.6) $$

$$ p_{n+1} - p_n = 0. \qquad (B.7) $$

Hence, the conservation of linear momentum is now proven.

Received: 13 January 2019 Accepted: 30 May 2020

Published online: 07 July 2020

References

1. Shabana AA, Sany JR. A survey of rail vehicle track simulations and flexible multibody dynamics. Nonlinear Dyn. 2001;26:179–212.

2. Shabana AA. Dynamics of multibody systems. Cambridge: Cambridge University Press; 2005.

3. Bauchau OA. Flexible multibody dynamics. Netherlands: Springer; 2011.

4. Escalona JL, Recuero AM. A bicycle model for education in multibody dynamics and real-time interactive simulation. Multibody Syst Dyn. 2012;27:383–402.

5. Xu WF, Meng DS, Chen YQ, Qian HH, Xu YS. Dynamics modeling and analysis of a flexible-base space robot for capturing large flexible spacecraft. Multibody Syst Dyn. 2014;23:357–401.

6. Betsch P. Energy-momentum integrators for elastic cosserat points, rigid bodies, and multibody systems. In: Betsch P, editor. Structure-preserving Integrators in nonlinear structural dynamics and flexible multibody dynamics, vol. 565. 1st ed. Springer: Cham; 2016. p. 31–89.

7. Arnold M, Cardona A, Brüls O. A lie algebra approach to lie group time integration of constrained systems. In: Betsch P, editor. Structure-preserving Integrators in Nonlinear Structural Dynamics and Flexible Multibody Dynamics, vol. 565. 1st ed. Cham: Springer; 2016. p. 91–158.

8. Meng D, She Y, Xu W, Lu W, Liang B. Dynamic modeling and vibration characteristics analysis of flexible-link and flexible-joint space manipulator. Multibody Syst Dyn. 2018;43:321–47.

9. Pasciuto I, Ausejo S, Celigüeta JT, Suescun A, Cazón A. A hybrid dynamic motion prediction method for multibody digital human models based on a motion database and motion knowledge. Multibody Syst Dyn. 2014;32:27–53.

10. Habachi AE, Duprey S, Cheze L, Dumas R. A parallel mechanism of the shoulder-application to multi-body optimisation. Multibody Syst Dyn. 2015;33:439–51.

11. Pettersson R, Nordmark A, Eriksson A. Optimisation of multiple phase human movements. Multibody Syst Dyn. 2013;30:461–84.

12. Schiehlen W. Multibody system dynamics: roots and perspectives. Multibody Syst Dyn. 1997;1:149–88.

13. Wasfy TM, Noor AK. Computational strategies for flexible multibody systems. Appl Mech Rev. 2003;56:553–613.

14. Reissner E. On one-dimensional finite-strain beam theory: the plane problem. J Appl Math Phys. 1972;23:795–804.

15. Simo JC, Vu-Quoc L. On the dynamics of flexible beams under large overall motions-the plane case: part I and II. J Appl Mech. 1986;53:849–63.

16. Shabana AA. Definition of the slopes and the finite element absolute nodal coordinate formulation. Multibody Syst Dyn. 1997;1:339–48.

17. Gerstmayr J, Irschik H. On the correct representation of bending and axial deformation in the absolute nodal coordinate formulation with an elastic line approach. J Sound Vib. 2008;318:461–87.

18. Shabana AA, Berzeri M. Development of simple models for the elastic forces in the absolute nodal coordinate formulation. J Sound Vib. 2000;235:539–65.

19. Ding JY, Wallin M, Wei C, Recuero AM, Shabana AA. Use of independent rotation field in the large displacement analysis of beams. Nonlinear Dyn. 2014;76:1829–43.

20. Omar MA, Shabana AA. A two-dimensional shear deformable beam for large rotation and deformation problems. J Sound Vib. 2001;243:565–76.

21. Gerstmayr J, Matikainen MK, Mikkola AM. A geometrically exact beam element based on the absolute nodal coordinate formulation. Multibody Syst Dyn. 2008;20:359–84.

22. Gerstmayr J, Sugiyama H, Mikkola A. An overview on the developments of the absolute nodal coordinate formulation. Applied mechanics reviews. In: Proceedings of the second joint international conference on multibody system dynamics, Stuttgart, Germany; May 2012.

23. Shabana AA, Christensen AP. Three-dimensional absolute nodal co-ordinate formulation: plate problem. Int J Numer Methods Eng. 1997;40:2775–90.

24. Mikkola AM, Shabana AA. A new plate element based on the absolute nodal coordinate formulation. Paper No DETC20001/VIB-21341, Proc of ASME 2001 DETC, Pittsburgh PA.

25. Gerstmayr J, Schöberl J. A 3D finite element method for flexible multibody systems. Multibody Syst Dyn. 2006;15:309–24.

26. Lang H, Linn J, Arnold M. Multi-body dynamics simulation of geometrically exact Cosserat rods. Multibody Syst Dyn. 2011;25:285–312.

27. Boer SE, Aarts RGKM, Hakvoort WBJ. Model reduction for efficient time-integration of spatial flexible multibody models. Multibody Syst Dyn. 2014;31:69–91.

28. Olshevskiy A, Dmitrochenko O, Kim C. Three- and four-noded planar elements using absolute nodal coordinate formulation. Multibody Syst Dyn. 2013;29:255–69.
29. Nachbagauer K, Pechstein AS, Irschik H, Gerstmayr J. A new locking-free formulation for planar, shear deformable, linear and quadratic beam finite elements based on the absolute nodal coordinate formulation. Multibody Syst Dyn. 2011;26:245-63.

30. Dmitrochenko O, Mikkola A. A formal procedure and invariants of a transition from conventional finite elements to the absolute nodal coordinate formulation. Multibody Syst Dyn. 2009;22:323-39.

31. Borri M, Bottasso C. An intrinsic beam model based on helicoidal approximation—part 1: formulation. Int J Numer Methods Eng. 1994;37:2267-89.

32. Borri M, Bottasso C. An intrinsic beam model based on helicoidal approximation—part 2: linearization and finite element implementation. Int J Numer Methods Eng. 1994;37:2291-309.

33. Bottasso C, Borri M. Energy preserving/decaying schemes for non-linear beam dynamics using the helicoidal approximation. Comput Methods Appl Mech Eng. 1997;143:393-415.

34. Merlini T, Morandini M. The helicoidal modeling in computational finite elasticity. Part I: variational formulation. Int J Solids Struct. 2004;41:5351-81.

35. Merlini T, Morandini M. The helicoidal modeling in computational finite elasticity. Part II: multiplicative interpolation. Int J Solids Struct. 2004;41:5383-409.

36. Betsch P, Uhlar S. Energy-momentum conserving integration of multibody dynamics. Multibody Syst Dyn. 2007;17:243-89.

37. Bauchau OA, Bottasso CL. On the design of energy preserving and decaying schemes for flexible, nonlinear multi-body systems. Comput Methods Appl Mech Eng. 1999;169:61-79.

38. Simo JC, Tarnow N. The discrete energy-momentum method. Conserving algorithms for nonlinear elastodynamics. J Appl Math Phys. 1992;43:757-92.

39. Goicolea JM, Garcia Orden JC. Dynamic analysis of rigid and deformable multibody systems with penalty methods and energy-momentum schemes. Comput Methods Appl Mech Eng. 2000;188:789-804.

40. Leyendecker S, Betsch P, Steinmann P. The discrete null space method for the energy consistent integration of constrained mechanical systems. Part III: flexible multibody dynamics. Multibody Syst Dyn. 2008;19:45-72.

41. Bottasso CL, Croce A. Optimal control of multibody systems using an energy preserving direct transcription method. Multibody Syst Dyn. 2004;12:17-45.

42. Géradin M, Cardona A. Flexible multibody dynamics: a finite element approach. New York: Wiley; 2001.

43. Sansour C, Wriggers P, Sansour J. Nonlinear dynamics of shells: theory, finite element formulation, and integration schemes. Nonlinear Dyn. 1997;13:279-305.

44. Sansour C, Wagner W, Wriggers P, Sansour J. An energy-momentum integration scheme and enhanced strain finite elements for the non-linear dynamics of shells. Int J Non-Linear Mech. 2002;37:951-66.

45. Sansour C, Nguyen TL, Hjiaj M. An Energy-momentum method for in-plane geometrically exact Euler-Bernoulli beam dynamics. Int J Numer Methods Eng. 2015;102:99-134.

**Publisher's Note**

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Submit your manuscript to a **SpringerOpen**® journal and benefit from:

► Convenient online submission

► Rigorous peer review

► Open access: articles freely available online

► High visibility within the field

► Retaining the copyright to your article

Submit your next manuscript at ► springeropen.com