A novel director-based Bernoulli–Euler beam finite element in absolute nodal coordinate formulation free of geometric singularities

A three-dimensional nonlinear finite element for thin beams is proposed within the absolute nodal coordinate formulation (ANCF). The deformation of the element is described by means of displacement vector, axial slope and axial rotation parameter per node. The element is based on the Bernoulli–Euler theory and can undergo coupled axial extension, bending and torsion in the large deformation case. Singularities – which are typically caused by such parameterizations – are overcome by a director per element node. Once the directors are properly defined, a cross sectional frame is defined at any point of the beam axis. Since the director is updated during computation, no singularities occur. The proposed element is a three-dimensional ANCF Bernoulli–Euler beam element free of singularities and without transverse slope vectors. Detailed convergence analysis by means of various numerical static and dynamic examples and comparison to analytical solutions shows the performance and accuracy of the element.


Introduction
In the present paper, a Bernoulli-Euler beam finite element based on the absolute nodal coordinate formulation (ANCF) is introduced.The ANCF has been developed by Shabana (1997Shabana ( , 2005) ) as an alternative to the classical large rotation vector formulation for the modeling of large deformation structural problems in two and three dimensions.Originally, thin (Bernoulli-Euler) ANCF elements led to a different solution as the classical Euler Elastica, however in the two-dimensional case, Gerstmayr and Irschik (2008) showed how to write the work of elastic forces, such that classical solutions can be retrieved by ANCF elements.In the threedimensional case, however, there are no singularity-free Bernoulli-Euler ANCF elements so far (in contrast to fully parametrized shear deformable ANCF elements, see Schwab and Meijaard (2010), who present and validate a locking free 3-D beam element with the usage of an elastic lien approach and the Whu-Washizu variational principle).Among earlier works on Bernoulli-Euler ANCF elements we mention von Dombrowski (2003); Dmitrochenko (2005); Gerstmayr and Shabana (2006).The formulation of von Dombrowski (2003) contains a time-dependent mass matrix and due to the parameterization of rotations it suffers from singularities.A similar approach as von Dombrowski ( 2003) is chosen by Dmitrochenko (2005); Dmitrochenko and Pogorelov (2003).In the latter references, the usage of a Frenet basis leads to singularities in the inflection points, because the torsion angle is not uniquely defined.Gerstmayr and Shabana (2006) presented an efficient approach for thin structural problems, however, torsion has not been included in their formulation.Note that in von Dombrowski (2003); Dmitrochenko (2005); Gerstmayr and Shabana (2006), the material measure of curvature has not been utilized, therefore the axial extension and bending are not decoupled, which may lead to erroneous results under large axial forces, see Gerstmayr and Irschik (2008).The original three-dimensional approach by Yakoub and Shabana (2001) describes a shear-deformable beam element, which involves a constant mass matrix and does not suffer from singularities.
Published by Copernicus Publications.The present ANCF formulation features strain measures for axial extension, bending and torsion.For the consistent derivation of such strain measures see Eliseev (2003).The principle of virtual work includes the kinematic relations for the strain measures of the rod, evaluated at a material line.For similar formulations of the direct approach to a material line see e.g. the works of Reissner (1973) or Antman (1972).
The parameterization of the kinematics of a spatial curved rod without shear was presented by Krommer and Vetyukov (2009), however there, the calculation was performed using a global Ritz approach.The kinematic description and the virtual work for the proposed element is chosen in a way similar to Vetyukov and Eliseev (2010) or rather Vetyukov (2008).In the latter paper, only static problems are regarded, while in the present paper, also dynamic behavior is analyzed.In contrast to standard ANCF elements, in which the mass matrix is constant, here the exact mass matrix is time-dependent, since a rigid cross section is assumed.In addition, a constant mass matrix, similar to Dmitrochenko (2005), can be applied for linear problems or small deformation problems, e.g.small oscillations, while for nonlinear problems a nonconstant mass matrix is utilized.The geometry of the beam is described by a curve, representing the beam's axis, as well as by an orientation of the cross section at each point of this axis, see Fig. 1.The orientation of the cross section is defined by an angle of axial rotation relative to a certain reference direction, which is introduced in order to perform the rotation in a correct manner and to avoid singularities originating from the usage of the rotation parameter.There exist different interpolation procedures in literature, in which displacements and rotations (see Simo and Vu-Quoc, 1986), displacements and slopes (see Shabana, 1997) or strains (see Gams et al., 2007;Zupan and Saje, 2003) are used as basic interpolated variables.In contrast to the above-mentioned formulations, the present approach is based on the interpolation of displacements, slopes and a rotation around the beam axis.For the interpolation of displacements and slopes, cubic shape functions are chosen, while the angle of axial rotation is interpolated linearly with respect to the beam's axis.
In the following, the geometric description of the finite element, the choice of degrees of freedom, as well as the definition of the strain energy for the ANCF beam element are  presented.In contrast to previous Bernoulli Euler ANCF elements, the proposed element is investigated using a large number of numerical examples, many of them suggested by Schwab and Meijaard (2009).In this work, the examples are restricted to small and large deformation static problems, buckling and linearized dynamic (eigenvalue) problems.In addition to the content of Nachbagauer et al. (2011), the exact dynamic terms of the beam element and a nonlinear dynamic example are provided in this work.

Geometric description
The geometry of the beam is described by a curve, representing the beam's axis, and a cross section at each point of this axis, see Fig. 1.In the present paper, a Bernoulli-Euler beam is considered, which means, that the cross section of the beam is thin, undeformed, and orthogonal to the beam axis (since shear deformation can be neglected) at any time and for any point of the beam axis.Moreover we assume, that the beam's axis intersects each cross section exactly at the cross section's centroid.The position p of an arbitrary point (or particle) in the cross section may thus be computed by . (1) Here, r(ξ) denotes the axial position (i.e., where the cross section belonging to p meets the beam's axis).The rotation tensor A(ξ) is defined by in which, e i (ξ) denotes the i-th base vector of the local axis frame at ξ (see Fig. 2).These vectors can be defined in terms of the slope vector r = ∂r ∂x , the reference direction (director) d and an angle θ as follows (for easier reading, the dependencies on ξ are omitted): in which e 30 denotes the normalized projection of the director d into the normal plane of the slope r , i. e., and e 20 is defined by the cross product Note that a similar parameterization of the cross section orientation has already been used by von Dombrowski (2003), however, not regarding a director.Summarizing, the geometry of a Bernoulli-Euler beam can be fully described by the axial position r, slope r , and cross section orientation angle θ, for a fixed choice of a director d.Equation ( 6) says, that geometric singularities occur, if (and only if) the slope r is collinear with the director d.A safe strategy for preventing singularities throughout the deformation process using an updating procedure is addressed in Sect. 4. Let us first turn to the finite element discretization.

Finite element discretization
The proposed beam finite element is defined by a pair of axial position vectors, axial slope vectors, and axial rotation angles, i.e., All together those are 14 degrees of freedom, which we collect in the vector of generalized coordinates see Fig. 3 for a sketch of the element.The element coordinates are chosen similar to von Dombrowski (2003) and Dmitrochenko (2005), however due to the use of the director, the meaning of the axial rotation angle is different in the present approach.For ξ ∈ [−L/2, L/2] the axial position vector r(ξ) is interpolated by cubic polynomials, with denoting the positional shape functions, and the rotation angle θ(ξ) is interpolated by linear polynomials, with denoting the rotational shape functions.The position and orientation of an arbitrary point is computed by with the shape function matrix The sub matrices S r and S θ are defined by with I denoting the identity in R 3 .The director d is also defined by the pairs,  In large deformation problems it is crucial to prevent the director d from becoming collinear with the beam axis direction e 1 .Therefore a director update in form of a projection into the normal plane of e 1 is suggested at every time or load step of the simulation.Note, that an update of this kind does not change the configuration of the local frame, see also Eq. ( 6) and Fig. 2.

Director update
Due to Eq. ( 6) the local frame is well defined if and only if the director d is not (numerically) collinear with the axial slope r at any point on the beam axis.A safe strategy to avoid this situation is to successively update the director at the FEnodes at every load/time step.As experienced in the scope of this work, a simple projection of the directors at the FE-nodes into the cross section of the beam seems to be sufficient, see Fig. 4. In case the axial direction e 1 becomes collinear during a load or time step of the simulation, it helps to repeat this same iteration step with a reduced load or time increment, and to update the director as explained at the intermediate step.Moreover, a subsequent rotation of the projected director around the beam axis by the torsional angle θ, and subsequently resetting θ to zero, provides an easy way to decrease the complexity of the terms appearing in the equations of motion and thereby speed up the static/dynamic calculations.As it seems, a simple linear interpolation of the director between the FE-nodes is sufficient for even complex scenarios and large deformations, given that a reliable initial configuration was chosen.For instance, two neighbouring nodal directors in opposite directions would lead to a vanishing director at the center of the finite element.Let be finally mentioned, that the proposed Bernoulli-Euler beam finite element provides C 1 -continuity along element borders only for the beam axis, whereas the torsion of the cross section, i.e. angle θ, is just C 0 -continuous.Hence, a 4th-order convergence will be prohibited particularly in problems with significant torsional effects.A fully C 1 continuous setting, requiring the rate of the torsional angle θ at the FE-nodes to serve as a generalized coordinate, together with a conforming interpolation of the director along the beam axis is left for further investigation.

Equations of motion
The equations of motion are written in the form of the Lagrange equations of second kind: where T denotes the kinetic energy, Π the strain energy, Q i are generalized external forces, q i are the generalized coordinates, as defined in Eq. ( 11), and qi denotes the partial derivative of q i with respect to the time variable t.In the case of static problems, the kinetic energy T vanishes, such that the equations of the static equilibrium are

Strain energy
With the principle of virtual work applied to the onedimensional Cosserat continuum (see, e.g., Antman, 1972;Eliseev, 2003;Reissner, 1973) one can prove that the strain energy per unit length is a function of two strain measures, one responsible for shear and axial extension, and the second one responsible for bending and torsion.In the considered model, the shear is kinematically constrained: the rotation of particles is in correspondence with the variation of the tangential direction r .Conforming to Simo (1985, Eq. 4.10), the axial extension is described by the axial strain in which ds = |r (ξ)| dξ and dS = |r 0 (ξ)| dξ denote the arc length in the deformed and undeformed state, respectively (see Fig. 1).Utilizing Einstein's summation convention, the vector of bending and torsional strain (see Eliseev, 2003) reads: in which the torsional strain κ 1 and bending strains κ 2 , κ 3 are expressed as the difference of the components of the vector of twist and curvature in deformed and undeformed state, As an important geometrical property the derivative of the basis with respect to the beam's axis may be determined by e i = k×e i and e 0 i = k 0 ×e 0 i , respectively.Notice, the components κ i are considered in the local basis e i , i.e., there holds Interpolated values for the axial strain are computed by r (ξ) = S r (ξ)q r , whereas curvature components are functions Mech.Sci., 4, 279-289, 2013 www.mech-sci.net/4/279/2013/ of the local frame (e i ) 3 i=1 , and thus of S(ξ)q, see Eqs. ( 3)-( 7) and Eq. ( 16).
If the local basis vectors e 2 , e 3 are chosen in the directions of the principal axis of the cross section, then the following quadratic approximation of the strain energy can be adopted: with EA, GJ, EI y and EI z being correspondingly the stiffnesses for axial extension, torsion and bending in the two principal directions.The form of the strain energy in Eq. ( 27) corresponds to the numerical formulation of Simo and Vu-Quoc (1986) if neglecting (cross sectional) transverse shear.
The strain energy in Eq. ( 27) is also similar to that suggested by von Dombrowski (2003) and Dmitrochenko (2005), however the curvature components are different in the present approach.Additional theoretical argumentation for the strain energy in the form in Eq. ( 27) and the kinematical definitions in Eqs.(23,24) are provided by the asymptotic analysis of the three-dimensional problem of the theory of elasticity for a naturally curved and twisted rod, presented by Yeliseyev and Orlov (1999).The variation the strain energy in Eq. ( 27) follows as The variation of ε in Eq. ( 23) is easily computed as The variation of κ i in Eq. ( 26) has to be computed columnwise by Here, ∂v ∂q can be easily computed by help of the shape functions in Eqs.(13, 15).

Kinetic energy in dynamic problems
In case of dynamic problems, the full Lagrange Equations in Eq. ( 21) have to be solved.The kinetic energy T is defined as in which ρ denotes the material density, and p, as defined in Eq. ( 1), points to an arbitrary point of the volume V of the element in undeformed state, see Fig. 1.Since ṗ depends linearly on q, i. e., ṗ(q, q) = ∂r ∂q (q) + ∂e 2 ∂q (q)η + ∂e 3 ∂q (q)ζ q, (32) Equation ( 31) turns into where the mass matrix M is defined by the integral and the matrices D and L 0 are defined in which the moments of inertia read I z = h 3 w 12 and I y = w 3 h 12 .The kinetic energy T appears in the following two terms in the Lagrange's equations Eq. ( 21): which means, that the following two terms have to be implemented regarding the kinetic energy:

Numerical examples
The proposed ANCF element has been implemented in the framework of the multibody and finite element research code HOTINT1 .In order to show the performance and accuracy of the proposed ANCF element, several numerical examples are considered.Many of the following examples are based on the beam benchmark problems proposed by Schwab and Meijaard (2009).As a first example, a cantilever beam under a tip load combined with a bending and torsional moment leading to small deformations is investigated.Further, we discuss large bending and torsion problems, Euler and lateral buckling, as well as the linear dynamic problem of an eigenfrequency analysis in the case of a simply supported beam with pre-stress as well as a nonlinear dynamic pendulum.If not mentioned differently in the description of the examples, we assume a cantilever beam, see Fig. 5 with length L in x direction, width w in y direction, and height h in z direction.By default all components of the tip resultant force f = F x , F y , F z T and the tip resultant moment M = M x , M y , M z T are zero, Young's modulus is assumed to be E = 2.1 × 10 11 N m −2 , the Poisson ratio is ν = 0.3, and the director points into z direction, d = [0, 0, 1] T , on the whole axis of the cantilever.

Small bending and torsion
As a first static example, a cantilever beam loaded at the tip by a vertical force, by a bending and torsional moment is investigated according to Schwab and Meijaard (2009, Sect.2.1).The chosen load case leads to small displacements and small rotations.The cantilever beam has length L = 1 m, width w = 0.01 m and height h = 0.02 m.The vertical tip load is chosen as F z = 1 × 10 −4 N, and the bending moment M y and the torsional moment M x are of the same magnitude M x = M y = 1×10 −4 Nm.The torsional stiffness of the rectangular cross section is set to GJ = Ghw 3 /3, where G denotes the standard shear modulus G = E/(2 + 2ν).The theoretical value of the tip displacement in z direction u z , as well as the axial rotation θ and the rotation around the y axis, ϕ y , can be computed by the following formulas: In our tests, the difference between theoretical and numerical solution was less than 10 −6 m for one element, and less than 10 −12 for two and more elements.

Large bending
A cantilever beam with length L = 2 m, width and height w = h = 0.1 m under a transversal tip load F y = 12EI y /L 2 is investigated.The obtained tip displacement is compared to a solution computed with arbitrary precision given in Gerstmayr and Irschik (2008), see Table 1.The convergence rate is of order 4, as can be seen in Fig. 6 Table 1.Convergence analysis (order 4) of the axial and transverse displacement at the tip of the cantilever beam in Sect.6.2 compared to the exact solution stated in Gerstmayr and Irschik (2008).Fig. 5. Cantilever beam with tip force and moment resultants.

Small bending and torsion
230 As a first static example, a cantilever beam loaded at the tip by a vertical force, by a bending and torsional moment is investigated according to (Schwab and Meijaard, 2009, Sec.2.1).The chosen load case leads to small displacements and small rotations.The cantilever beam has length L = 1 m, width w = 0.01 m and height h = 0.02 m.The vertical tip load is chosen as F z = 1 • 10 −4 N, and the bending moment M y and the torsional moment M x are of the same magnitude M x = M y = 1 • 10 −4 Nm.The torsional stiffness of the rectangular cross section is set to GJ = Ghw 3 /3, where G denotes the standard shear modulus G = E/(2 + 2ν).The theoretical value of the tip displacement in z-direction u z , as well as the axial rotation θ and the rotation around the y-axis, ϕ y , can be computed by the following formulas: In our tests, the difference between theoretical and numerical solution was less than 10 −6 m for one element, and less than 10 −12 for two and more elements.and Meijaard, 2009, Sec.2.3).Since it was not possible for the software package ANSYS to calculate a numerical solution for the full amount of loading, only 50 percent were pre-270 scribed.Thereafter, the solution of ANSYS was compared to the solution by using the proposed beam finite element.Note, that the solution of the proposed element coincides with the solution of ANSYS up to 5 digits, which significantly behaves better than an approach by (von Dombrowski, 2003, 275 Fig. 3).

Full circle bending problem
A cantilever beam with length L = 2 m, width and height w = h = 0.1 m under a tip bending moment M z = 2πE I z /L is analyzed.Table 3 reports on the convergence behavior.
280 Figure 6.The error |u k −u| between computed and exact solution is plotted versus the degrees of freedom for the large bending example Sect.6.2.Actually, fourth order convergence may be observed only in plane problems.

Large bending and torsion
A cantilever beam with length L = 1 m, width w = 0.005 m, and height h = 0.02 m is investigated.At the tip, a bending moment M y = 50 Nm, and a torsional moment M x = 12.5 Nm is applied.The torsional stiffness of the cross section is supposed to be GJ = G h w 3 /3 with G again denoting the shear modulus.The obtained tip displacement is compared in Table 2 to the numerical solution based on an ANSYS implementation with 40 beam elements.Notice, this example is almost identical with a benchmark problem in Schwab and Meijaard (2009, Sect. 2.3).Since it was not possible for the software package ANSYS to calculate a numerical solution for the full amount of loading, only 50 percent were prescribed.Thereafter, the solution of ANSYS was compared to the solution by using the proposed beam finite element.Note, that the solution of the proposed element coincides with the solution of ANSYS up to 5 digits, which significantly

Euler buckling
As a first buckling problem, a Euler buckling of a cantilever beam under a compressive normal force F x is analyzed.The dimensions of the beam are chosen as length L = 1 m, width w = 0.01 m and height h = 0.02 m, according to the benchmark problem in (Schwab and Meijaard, 2009, Sec. 2.5).
The theoretical buckling load can be calculated by F th = (π/4)(E I min /L 2 ), in which I min denotes the minimum of the bending stiffnesses I min = min(I y ,I z ).Notice, that differently to Schwab and Meijaard (2009) the eigenvalues are 290 numerically calculated in the deformed state which corre-  4. The ratio of the numerical to the theoretical buckling load Fnum/Fth in case of Euler buckling, Sect.6.5, compared to the results of Schwab and Meijaard (2009).Note, that the ratio Fnum/Fth goes up after using 8 elements when using the original value for EA, which is due to the fact, that the numerical buckling loads could only be determined by solving an eigenvalue problem for the deformed system.A second test for a stiffer value of EA showed order 4 convergence.sponds to the numerical compressive normal force F num .In other words, the buckling loads were determined by observing the lowest bending mode of the system, which corresponds to the deformed state of the beam (since the simu-295 lation sofware we used does not provide linearized equations directly for the eigenmode analysis, but the system matrix is assembled automatically for the deformed state of the beam).Thus, in order to allow a comparison to the benchmark results, another test was done, in which the axial stiffness of the 300 beam is penalized (multiplied by the factor 100), such that buckling occurs at a less deformed state of the beam.Both cases (original and penalized axial stiffnesses) are compared to a benchmark example by Schwab and Meijaard (2009) in Tab. 4. 305

Lateral buckling
As a further buckling problem, the lateral buckling of a cantilever beam under a lateral force F is investigated.The dimensions of the beam are chosen as length L = 1 m, width w = 0.002 m and height h = 0.02 m, as in 310 behaves better than an approach by von Dombrowski (2003, Fig. 3).

Full circle bending problem
A cantilever beam with length L = 2 m, width and height w = h = 0.1 m under a tip bending moment M z = 2 π E I z /L is analyzed.Table 3 reports on the convergence behavior.

Euler buckling
As a first buckling problem, a Euler buckling of a cantilever beam under a compressive normal force F x is analyzed.The dimensions of the beam are chosen as length  4. The ratio of the numerical to the theoretical buckling load F num /F th in case of Euler buckling, Sect.6.5, compared to the results of Schwab and Meijaard (2009).Note, that the ratio F num /F th goes up after using 8 elements when using the original value for EA, which is due to the fact, that the numerical buckling loads could only be determined by solving an eigenvalue problem for the deformed system.A second test for a stiffer value of EA showed order 4 convergence.L = 1 m, width w = 0.01 m and height h = 0.02 m, according to the benchmark problem in Schwab and Meijaard (2009, Sect. 2.5).The theoretical buckling load can be calculated by F th = (π/4)(E I min /L 2 ), in which I min denotes the minimum of the bending stiffnesses I min = min(I y , I z ).Notice, that differently to Schwab and Meijaard (2009) the eigenvalues are numerically calculated in the deformed state which corresponds to the numerical compressive normal force F num .In other words, the buckling loads were determined by observing the lowest bending mode of the system, which corresponds to the deformed state of the beam (since the simulation sofware we used does not provide linearized equations directly for the eigenmode analysis, but the system matrix is assembled automatically for the deformed state of the beam).Thus, in order to allow a comparison to the benchmark results, another test was done, in which the axial stiffness of the beam is penalized (multiplied by the factor 100), such that buckling occurs at a less deformed state of the beam.Both cases (original and penalized axial stiffnesses) are compared to a benchmark example by Schwab and Meijaard (2009)  Table 5.The ratio of the numerical to theoretical buckling load F num /F th in case of lateral buckling, Sect.6.6, compared to the results of Schwab and Meijaard (2009).Alike in Sect.6.5, the numerical buckling loads could only be determined by solving an eigenvalue problem for the deformed system.In order to compare the results, a thinner beam has been utilized, such that buckling occured already at a less deformed state.

Lateral buckling
As a further buckling problem, the lateral buckling of a cantilever beam under a lateral force F z is investigated.The dimensions of the beam are chosen as length L = 1 m, width w = 0.002 m and height h = 0.02 m, as in the benchmark example by Schwab and Meijaard (2009, Sect. 2.6).
The theoretical buckling load can be calculated by F th = 4.012599344 √ E I min GJ/L 2 , in which E I min is the smaller flexural stiffness with I min = min(I y , I z ) and GJ = G hw 3 /3 the torsional stiffness of the rectangular cross section, in which G denotes the shear modulus.The numerical solution is compared to the convergence rate of the nonlinear beam finite element proposed by Schwab and Meijaard (2009), see Table 5. Differently to Schwab and Meijaard (2009), the eigenvalues are calculated in the deformed state, which corresponds to the applied numerical buckling load F num .Since buckling occurs already at a less deformed state as the ratio min{w, h}/ max{w, h} becomes smaller, the width w of the beam was chosen to be w = 0.0002 m (which is 100 times smaller than in the benchmark example by Schwab and Meijaard, 2009).The ratio F num /F th for both the cases w = 0.0002 m and w = 0.002 m are outlined in Table 5.

Eigenfrequencies of a simply supported beam
According to Schwab and Meijaard (2009, Sect. 2.9), the eigenfrequencies of a simply supported beam with dimensions length L = 1 m, width w = 0.02 m and height h = 0.02 m with pre-stress F x are computed.For a sketch of the problem setup, see Fig. 8.The exact values for the zero-load frequency w 0 , the non-dimensional pre-stress α and the first non-dimensionalized frequency w B can be calculated by the formulas The convergence analysis in dependency on the nondimensional pre-stress α is outlined in Table 6.Here, the numerical eigenvalue w num is divided by the theoretical eigenvalue w B .Differently to Schwab and Meijaard (2009), the numerical eigenvalues are calculated in deformed state (by using the stiffness matrix of the deformed geometry), which corresponds to the load F = α π 2 E I y L −2 .Thus, the numerical solution is different from the theoretical solution w B , but converging as α tends to zero.

Dynamic example of a double pendulum
In this example the nonlinear dynamic behavior of the proposed elements is tested by a rigid flexible pendulum as in Sugiyama et al. (2003).Note, that the original example is used for verifying thick beam finite elements including shear terms in the variational formulation according to Timoshenko's beam theory.In contrast to that, in this work the example serves for the experimental verification of the dynamic solution of the proposed Bernoulli Euler beam finite element, by comparing it to another solution obtained by already verified thick ANCF-beam finite elements based on Timoshenko's theory, see Nachbagauer andGerstmayr (2012, 2013).The theory says, that both solutions (according to Bernoulli-Euler's and Timoshenko's theory, respectively) converge to each other, as the spatial (thickness to length) ratio of the beams goes to zero.Two tests are done, one for the original experiment setup, where the beams' spatial ratio is 1 : 5, and another test for the ratio of 1 : 50 and an appropriate adjustment of the material parameters, which is necessary to obtain a comparable solution.For both tests a sufficiently refined spatial and temporal discretization is chosen, such that the discretization error has no significant influence.

Test A
The considered multibody system consists of a rigid body connected to a very flexible right angle frame, see Fig. 9.
The involved components are a rigid and a flexible body, and two revolute joints, which connect the rigid body both to the ground and to the right angle frame.The rigid body has a length of L r = 0.2 m, a cross section area of A r = 2.5 × 10 −3 m 2 , and a material density of ρ r = 7200 kg m −3 .The flexible right angle frame is modeled by two beams connected by a rigid joint.Each beam in the original example (see Sugiyama et al., 2003;Nachbagauer and Gerstmayr, 2013) has a length of L f = 0.5 m, a cross section area of A f = 1×10 −2 m 2 , a material density of ρ f = 7200 kg m −3 , a Young's modulus of E f = 2.0×10 6 Pa, a Poisson ratio of ν f = 0.3, and a second moment of area I f A = 8.33 × 10 −6 m 4 .Apart from gravity, there is no external force acting on the multibody system.The rigid body is rotated 30 degrees around the y-axis.Since the revolute joint which connects the rigid body to the ground allows rotation around the global y-axis, the motion of the rigid body is executed in the global (x, z)-plane.The flexible beam is rotated 30 degrees around the z-axis, and coherently the axis of the second revolute joint is chosen 30 degrees off the global y-axis.Therefore, the flexible body performs an in-and out-of-plane motion, which causes torsion.Two Finite Element simulations are performed for this Test: one simulation with an FE-approximation by 32 of the proposed Bernoulli-Euler beam finite elements, and another simulation with an FE-approximation by 32 already verified Timoshenko beam finite elements (see Nachbagauer

Test B
The thickness of the pendulum is divided by the factor 10 compared to Test A, such that the cross section area of the rigid part measures A r = 2.5 × 10 −5 m 2 and the cross section area of the flexible part A f = 1 × 10 −4 m 2 .In order to obtain a similar dynamic solution as in Test A, Young's modulus is chosen E f = 2.0 × 10 10 Pa, whereas the density of both the rigid and the flexible part is multiplied by the factor 10 2 , such that ρ r = ρ f = 720 000 kg m −3 .Again, as in Test A, two Finite Element simulations are performed.The dynamic behavior of the double pendulum is outlined in the Plot (b) of Fig. 10.In the first simulation (blue curves) Bernoulli-Euler beam finite elements are used, whereas the second simulation utilizes Timoshenko beam finite elements (black curves).

Large bending and torsion: director update
In this final example the impact of the director update, as suggested in Sect.4, is studied.The geometry of the beam and its material properties are chosen as in Sect.6.3 (large bending and torsion).In difference to Sect.6.3 a tip load is applied with components F y = 25 N and F z = 250 N instead of moments.In the simulation exactly two load steps are performed, the first step with the half, and the second step with the full loading.For the convergence analysis at each FEnode the director is chosen to point in z direction.Two tests are performed: In Test A the directors are kept constant, in Test B the directors are updated at each node (see Sect.

A remark on the use of constraints
The considered example includes a right angle joint of two beam segments.The demand on identity of rotations on the end points of both segments is trivial in the framework of Timoshenko theory, in which rotations appear as indepen-440 dent degrees of freedom.An efficient implementation of this condition for Bernoulli-Euler elements under consideration has been realized via Lagrange multipliers for six scalar conditions of equivalence of translations and rotations in the adjacent cross-sections of the two rod segments.

7 Conclusions
In the present paper, a spatial thin beam element with axial, bending, and torsional deformation is presented.The proposed element underlies the absolute nodal coordinate formulation, which is designed for large displacements, large 450 rotations and large deformations and multibody dynamics problems.The element kinematics and deformation energy are chosen according to Vetyukov and Eliseev (2010), including an additional rotation director, which is updated during computation process.The investigation of numerical exam-455 ples, both in the static and dynamic case, shows the expected accuracy and a fast performance of the proposed element.Moreover, in contrast to already existing 3D ANCF elements, it does not suffer from geometrical singularities.their projection into the cross section.Table 7 shows clearly, that the update of the directors does not add any additional inaccuracy to the discretized solutions.
A remark on the use of constraints The considered example includes a right angle joint of two beam segments.The demand on identity of rotations on the end points of both segments is trivial in the framework of Timoshenko theory, in which rotations appear as independent degrees of freedom.An efficient implementation of this condition for Bernoulli-Euler elements under consideration has been realized via Lagrange multipliers for six scalar con-ditions of equivalence of translations and rotations in the adjacent cross-sections of the two rod segments.

Conclusions
In the present paper, a spatial thin beam element with axial, bending, and torsional deformation is presented.The proposed element underlies the absolute nodal coordinate formulation, which is designed for large displacements, large rotations and large deformations and multibody dynamics problems.The element kinematics and deformation energy are chosen according to Vetyukov and Eliseev (2010), including an additional rotation director, which is updated during computation process.The investigation of numerical examples, both in the static and dynamic case, shows the expected accuracy and a fast performance of the proposed element.Moreover, in contrast to already existing 3-D ANCF elements, it does not suffer from geometrical singularities.

Figure 2 .
Figure 2. In order to have a uniquely defined orientation of the cross section about the beam axis direction e 1 = r |r | a non-collinear director d is utilized.The local frame (e 1 , e 2 , e 3 ) is defined by the normalized projection e 30 of d into the normal plane of e 1 , and a subsequent rotation around a torsional angle θ.

e 2 =Figure 3 .
Figure 3. Degrees of freedom for one beam finite element in the nodes α and β: axial position vector r, slope vector r , and cross section orientation angle θ.

Figure 4 .
Figure 4.In large deformation problems it is crucial to prevent the director d from becoming collinear with the beam axis direction e 1 .Therefore a director update in form of a projection into the normal plane of e 1 is suggested at every time or load step of the simulation.Note, that an update of this kind does not change the configuration of the local frame, see also Eq. (6) and Fig.2.

Figure 5 .
Figure 5. Cantilever beam with tip force and moment resultants.
with length L = 2 m, width and height w = h = 0.1 m under a transversal tip load F y = 12EI y /L 2 is investigated.The obtained tip displacement is compared to a solution computed with arbitrary precision given in Ger-255 stmayr and Irschik (2008), see Tab. 1.The convergence rate is of order 4, as can be seen in Fig. 66.3 Large bending and torsionA cantilever beam with length L = 1 m, width w = 0.005 m, and height h = 0.02 m is investigated.At the tip, a bend-260 ing moment M y = 50 Nm, and a torsional moment M x = 12.5 Nm is applied.The torsional stiffness of the cross section is supposed to be GJ = Ghw 3 /3 with G again denoting the shear modulus.The obtained tip displacement is compared in Tab. 2 to the numerical solution based on an ANSYS implementation with 40 beam elements.Notice, this example is almost identical with a benchmark problem in(Schwab

Fig. 6 .
Fig.6.The error|u k − u| between computed and exact solution is plotted versus the degrees of freedom for the large bending example in Sect.6.2.Actually, fourth order convergence may be observed only in plane problems.

Fig. 7 .
Fig. 7. Deformation of the beam in Sect.6.3 with color-plot of the torsional moment.

Figure 7 .
Figure 7. Deformation of the beam in Sect.6.3 with color-plot of the torsional moment.

Figure 9 .
Figure9.Geometry of a double pendulum consisting of a rigid and a flexible body, as proposed inSugiyama et al. (2003).

Fig. 10 .
Fig. 10.Trajectories of the tip of the double pendulum in Test A (a) and in Test B (b) in Sect.6.8.For each of the two tests either two simulations are performed: One by using an FE-discretization with 32 Bernoulli-Euler Finite Elements (3 blue curves for x, y, and z position components), and another by using 32 Timoshenko beam finite elements (3 black curves).The three black curves cannot be seen in Plot (b), since the blue curves are on top of them.In other words, the FE-solutions converge as the ratio of thickness to length goes to zero.

Figure 10 .
Figure 10.Trajectories of the tip of the double pendulum in Test A (a) and in Test B (b) in Sect.6.8.For each of the two tests either two simulations are performed: one by using an FEdiscretization with 32 Bernoulli-Euler Finite Elements (3 blue curves for x, y, and z position components), and another by using 32 Timoshenko beam finite elements (3 black curves).The three black curves cannot be seen in Plot (b), since the blue curves are on top of them.In other words, the FE-solutions converge as the ratio of thickness to length goes to zero.

Table 1 .
Convergence analysis (order 4) of the axial and transverse displacement at the tip of the cantilever beam in Sect.6.2 compared to the exact solution stated inGerstmayr and Irschik (2008).

Table 2 .
Displacement of the tip in case of large bending and torsion in Sect.6.3 compared to a numerical solution with 40 beam elements in the software package ANSYS.

Table 3 .
Displacements ux and uy for the full circle bending in Sect.6.4 converge at order 4.

Table 3 .
Displacements u x and u y for the full circle bending in Sect.6.4 converge at order 4.

Table 6 .
First non-dimensionalized eigenfrequency analysis for a simply supported beam with pre-stress in Sect.6.7 versus number (#) of finite elements.

Table 7 .
Convergence Table Sect.6.9.In the first column the number of Finite Elements is displayed.The next two columns show the error of the Finite Element iterates compared the converged solution |u FE A/B − u * A/B | in Test A and Test B, whereas the third column reports on the difference of the iterates |u FE A −u FE B |.As converged solution, the FE-solution for 512 Elements was used.The converged solutions in Test A and B differ by |u * A − u * B | = 2.4050 × 10 −7 .