# A Hybrid Interpolation Method for Geometric Nonlinear Spatial Beam Elements with Explicit Nodal Force.

1. Introduction

In engineering practice, a majority of three-dimensional beams can be considered slender beams and can thus be regarded as Euler-Bernoulli beams. Much research in recent years has been focused on modeling geometrical nonlinearity for 3D Euler-Bernoulli beams subjected to large deformations. Because spatial large rotations are physically nonadditive, an improper discretization of spatial rotations may lead to nonobjectivity of strain measures [1]. Different from the linear beam theory, large rotations of beam cross sections cannot be totally determined by the beam centerline. However, displacement and rotation of the cross section of a slender beam must satisfy the Kirchhoff constraints. Therefore, there is a critical need to develop an efficient interpolation method to solve these problems.

Geometrically exact beam theory, developed by Reissner [2,3] and Simo [4,5], has been widely adopted for modeling of geometrical nonlinear beams [6-11]. For a geometrically exact beam, the configuration of a spatial beam can be described by the position vector of the beam centerline and an orthogonal matrix, which specifies rotation of a cross section. A curvature vector is used to describe the rotational strain measure, which is anisotropic and proportional to the stress resultants. Relationship between rotation and rotational strain was found to be consistent with the virtual power principle and is valid for configurations with arbitrary large displacements and rotations. Due to independent interpolations of displacement and rotation, use of geometrically exact beam theory may lead to shear-locking problem for slender beams and nonobjective interpolation of strain measure [1]. In order to avoid these problems, many researchers presented different formulations of geometrically exact beam elements. The core technology in developing these formulations has been on the discretization of rotational degrees of freedom or the incremental rotation [6,11-13]. In spite of the efforts, a problem still exists with accurate determination of an analytical expression of elastic forces using the displacement-based geometrically exact beam theory.

To avoid the interpolation of rotation parameters, an absolute nodal coordinate formulation (ANCF) was proposed by Shabana [14-18]. Three polynomials were used as the assumed displacement field in the ANCF-based beam element. In ANCF, 12 variables were used for each node of a 3D beam element, which includes the position vector and 9 slopes. Therefore, dimension of governing equations is significantly larger in comparison with traditional beam elements. Based on the interpolation of beam curvatures, a new beam element was proposed by Zupan [7,19], in which displacement and rotational vector were not interpolated at all. In this formulation, strain measure-interpolation based element was demonstrated to preserve the objectivity, namely, invariant with rigid body motions, and path independence of strain measures. However, use of the beam element based on the curvature interpolation was not convenient when used in calculating nodal displacement and rotation, since they cannot be analytically integrated by strain measures.

In this paper, a hybrid interpolation scheme for geometrically exact Euler-Bernoulli beams is proposed. Global nodal displacement and the rotational vector are the basic unknown variables for the formulation. Interpolation of the beam centerline was used to determine nodal value and its derivatives of the curvatures. These new parameters were also interpolated. Using combined interpolation of the nodal displacement and strain measures, an analytical expression for nodal forces was formulated and objectivity and path independence of strain measures were satisfied.

The outline of this work is listed as follows. In Section 2, basic theory of geometrically exact Euler-Bernoulli beams with large deformation is briefly reviewed. In Section 3, based on the principle of virtual power, weak form of the dynamic equilibrium equations for 3D Euler-Bernoulli beams was deduced using the novel method proposed here, which is concise and convenient for large rotations. Next, hybrid interpolations of nodal displacements and curvatures are introduced in Section 4. In Section 5, governing equations for the discrete finite beam element model are given. Property tests and applications of the beam element are presented in Sections 6 and 7, respectively.

2. General Theory of Geometrically Exact Euler-Bernoulli Beams

In this section, geometrically exact Euler-Bernoulli beam theory is reviewed, including kinematics of large deformation beams, strain measures, stress resultants, and constitutive relations.

2.1. Kinematics and Strain Measures. The initial and deformed configurations of a spatial Euler-Bernoulli beam in the inertial reference frame (X-Y-Z) are shown in Figure 1. The main assumption of 3D Euler-Bernoulli beam theory is that all cross sections keep rigid and perpendicular to the tangent vector of the centerline during deformation. The configuration can be fully described by the following.

(1) The position vector r([s.sub.0]) of beam centerline, where [s.sub.0] [member of] [0, [L.sub.0]] is the arc-length coordinate of the beam centerline at the initial configuration and [L.sub.0] is the initial length of the beam.

(2) The cross-sectional reference frame {[e.sub.s], [e.sub.t], [e.sub.b]}: the orientation of a cross section is denoted by the right-handed orthonormal triads attached to the cross section centroid. With the Kirchhoff constraint, the normal vector of the cross section [e.sub.s] is parallel to the tangent vector of the beam centerline r', which yields

[e.sub.s] = r'/[parallel]r'[parallel], [parallel]r'[parallel] = 1 + [[epsilon].sub.s], (1)

where a symbol with a superscript " '" denotes its derivative with respect to the arc-length coordinate [s.sub.0]. [[epsilon].sub.s] is the axial strain. The second and third axial vectors [e.sub.t] and [e.sub.b] are in the directions of two principal axes of inertia of the cross section. Since [e.sub.s] = [e.sub.t] x [e.sub.b], it is obvious that the tangent vector of the beam centerline r' does not have any components in the directions [e.sub.t] and [e.sub.b]; that is,

[e.sup.T.sub.t] r' = [e.sup.T.sub.b] r' = 0. (2)

The equation means that the beam is shear-free. In consequence, orientation of a cross section is in some sense determined by the beam centerline and the position vector derivatives seeing (1) and (2). That is to say, the so-called Kirchhoff constraint is carried out accurately in a strong manner.

Under the definition of the local coordinate system of the beam cross section, the rotational matrix of the cross section with respect to the inertial one is given by

R = [[e.sub.s] [e.sub.t] [e.sub.b]]. (3)

The cross-sectional rotation is parameterized by rotational vector [[theta].sub.a]. Hence, the rotational matrix can be expressed as

[mathematical expression not reproducible] (4)

here with [[rho].sub.a] = [parallel][[theta].sub.a][parallel]. In this paper, a skew-symmetric tensor with symbol "~" denotes the cross product of two vectors; that is, a x b = [??]b. The skew-symmetric tensor is given by

[mathematical expression not reproducible]. (5)

The current cross-sectional orientation matrix in (4) could be decomposed into the original orientation of the cross section and the rotation of the local coordinate system; that is,

[[e.sub.s] [e.sub.t] [e.sub.b]] = [bar.R] [[e.sub.0s] [e.sub.0t] [e.sub.0b]], [bar.R] = R([theta]), (6)

where [theta] is the rotation of the cross section. The time derivative of base vectors of the cross-sectional reference frame can be obtained by the cross product of its angular velocity and the base vectors:

[[??].sub.s] = [omega] x [e.sub.s], [[??].sub.t] = [omega] x [e.sub.t], [[??].sub.b] = [omega] x [e.sub.b], (7)

where the superscript dot denotes the derivative with respect to time and the angular velocity of the cross section can be expressed by

[omega] = T([[theta].sub.a]) x [[??].sub.a] (8)

here with

[mathematical expression not reproducible]. (9)

The skew-symmetric tensor of the angular velocity can be expressed in terms of the time derivative of the rotational matrix referring to (7):

[mathematical expression not reproducible]. (10)

Similarly, the derivative of base vector with respect to the arc-length coordinate could be given by the cross product of the curvature and base vectors:

[e'.sub.s] = [kappa] x [e.sub.s], [e'.sub.t] = [kappa] x [e.sub.t], [e'.sub.b] = [kappa] x [e.sub.b]. (11)

Curvature is used as the strain measure in Reissner's work [2,3]. Simo [4,5] has also discussed it in detail in geometrically exact beam theory. The skew-symmetric tensor of the curvature vector can be formulated by the derivative of rotational matrix with respect to the arc-length parameter referring to (11):

[??] = R'[R.sup.T]. (12)

Furthermore, the curvature vector is also represented by the derivative of rotational vector with respect to the arc-length coordinate:

[kappa] = T([[theta].sub.a]) x [[theta]'.sub.a]. (13)

Angular velocity and curvature vector can be expressed in the cross-sectional coordinate system as follows:

[omega] = [[omega].sub.s][e.sub.s] + [[omega].sub.t][e.sub.t] + [[omega].sub.b][e.sub.b], [kappa] = [[kappa].sub.s][e.sub.s] + [[kappa].sub.t][e.sub.t] + [[kappa].sub.s][e.sub.b], (14)

where [[kappa].sub.s] denotes the cross-sectional torsional curvature and [[kappa].sub.t] and [[kappa].sub.b] are the cross-sectional bending curvatures in the two perpendicular directions.

The axial strain [[epsilon].sub.s] can be obtained by the elongation referring to (1):

[[epsilon].sub.s] = [e.sup.T.sub.s] r' - 1. (15)

2.2. Stress Resultants and Constitutive Relations. The stress-resultant force and moment of a cross section can be expressed in the cross-sectional coordinate system as follows:

f = [f.sub.s][e.sub.s] + [f.sub.t][e.sub.t] + [f.sub.b][e.sub.b], m = [m.sub.s][e.sub.s] + [m.sub.t][e.sub.t] + [m.sub.b][e.sub.b]. (16)

The constitutive relation between the strain measure and the resultant force and moment of stress for the beam of linear elastic material is given by [3,4]

[f.sub.s] = EA[[epsilon].sub.s], [m.sub.s] = GJ[[kappa].sub.s], [m.sub.t] = [EI.sub.t][[kappa].sub.t], [m.sub.b] = [EI.sub.b][[kappa].sub.b], (17)

where EA, GJ, [EI.sub.t], [EI.sub.b] are the axial tensile, torsional, and principal bending stiffness, respectively. In some sense, when beams are regarded as one-dimensional structures, their constitutive relations are anisotropy. This feature must be taken into account in the modeling of beam elements.

3. Formulation of 3D Geometrically Exact Euler-Bernoulli Beam

In most formulations of linear finite elements, the principle of virtual work is adopted. However, if rotations are large, it is difficult to obtain expression of the rotational virtual work. Due to the nonadditive characteristic of 3D rotations, product of the virtual rotational vector and moment do not constitute virtual work. There are kinds of parameters describing large rotation. Different parameters have different meanings. The principle of virtual power has a natural advantage in finite element formulation of large rotations since rational virtual power is the product of the virtual angular velocity and moment. Virtual angular velocity is definitely conjugate to internal moment, and angular velocity has clear physical meanings. In this work, principle of virtual power is adopted to deduce the finite element formulation of 3D Euler-Bernoulli beams with large rotations.

A beam can be regarded as a series of rigid bodies under the assumption of rigid cross section. In such case, current arc-length coordinate s of the beam centerline is a function of time t and the initial arc-length coordinate [s.sub.0]. Therefore, all strains, stresses, and internal forces can be regarded as functions of these two variables. Thus, axial strain is given by

[[epsilon].sub.s] = [partial derivative]s/[partial derivative][s.sub.0] - 1. (18)

An infinitesimal beam ds = (1 + [[epsilon].sub.s])[ds.sub.0] is shown in Figure 2, where f and m are the stress-resultant force and moment, respectively. In order to simplify the derivation, a reasonable assumption could be made that there is no external force on the infinitesimal beam. The strong form of the dynamic equilibrium equations of the infinitesimal beam can be expressed as

[mathematical expression not reproducible], (19)

where J is the cross-sectional moment of inertia and [rho] is the mass density per unit length before deformation.

The virtual power produced by applied forces of the infinitesimal beam can be formulated as

[mathematical expression not reproducible]. (20)

The virtual power of inertia forces of the infinitesimal beam is given by

[mathematical expression not reproducible]. (21)

According to the principle of virtual power, the virtual power of internal forces can be given by the difference of the virtual powers of applied forces and inertia forces, which yields

[mathematical expression not reproducible]. (22)

From another point of view, the virtual power of inertia forces also can be obtained through the integration of (22) along the beam length:

[mathematical expression not reproducible]. (23)

The components in (14) could be given by

[mathematical expression not reproducible], (24)

[mathematical expression not reproducible]. (25)

The first-order derivatives of (14) with respect to arc-length coordinate [s.sub.0] and time t, respectively, are

[omega]' = ([[omega]'.sub.s][e.sub.s] + [[omega]'.sub.t][e.sub.t] + [[omega]'.sub.b][e.sub.b]) + [kappa] x [omega], (26)

[??] = [[??].sub.s][e.sub.s] + [[??].sub.t][e.sub.t] + [[??].sub.b][e.sub.b] + [omega] x [kappa]. (27)

As arc-length [s.sub.0] and time are independent, the derivation orders of base vectors can be exchanged. Therefore, the following relationships can be obtained referring to (24):

[mathematical expression not reproducible]. (28)

Considering (28) and (27), (26) could be rewritten as

[mathematical expression not reproducible]. (29)

Taking the time derivative of (1), one yields

[??]' - [omega] x r' = [[??].sub.s][e.sub.s]. (30)

With (29)-(30) and (23), the virtual power of internal forces can be expressed as

[mathematical expression not reproducible]. (31)

All parameters in the virtual power equations are defined in the cross-sectional coordinate system, which do not contain rigid body motions, and the beam element formulated using the present method cannot produce self-strains, even for large deformations.

4. Discretization

It can be concluded from Section 3 that the internal virtual power of 3D Euler-Bernoulli beams is directly dependent on the strain measure rather than the displacement and rotation. Unlike the rotational vector, the strain measure is additive. It cannot cause nonobjective strain measure [7] by the strain interpolation instead of the displacement or rotational interpolation. Nodal coordinates which are global nodal displacement and the rotational vector are the basic unknown variables. Based on a Hermitian interpolation of the beam centerline, nodal value and its derivatives of the bending curvature are consistently determined which are afterwards reinterpolated again by Hermitian polynomials. A constant axial strain field as well as a constant torsion field is constructed from the nodal degrees of freedom based on the geometrical approximations.

4.1. Nodal Coordinates. Basis vectors of cross sections at the two ends are {[e.sub.1s], [e.sub.1t], [e.sub.1b]} and {[e.sub.2s], [e.sub.2t], [e.sub.2b]}, respectively. The rotational vectors on the two ends are [[theta].sub.a1] and [[theta].sub.a2], respectively. A beam is divided into a series of beam elements. The left and right cross-sectional orientations of a node are not the same when two elements are joined at the node in an intersection angle different from 180[degrees]. The left and right original orientations of the cross sections of a node may be different in a nonstraight beam. However, the rotations of these two cross sections [theta] are the same. Nodal displacements [r.sub.1], [r.sub.2] and rotational vectors [[theta].sub.1], [[theta].sub.2] are adopted as the basic unknown functions of a two-node beam element, namely, nodal coordinates:

q = [[[r.sup.T.sub.1] [[theta].sup.T.sub.1] [r.sup.T.sub.2] [[theta].sup.T.sub.2]].sup.T]. (32)

The corresponding generalized velocities [??] include the nodal velocities [[??].sub.1], [[??].sub.2] and the angular velocities [[omega].sub.1], [[omega].sub.2]:

[mathematical expression not reproducible]. (33)

4.2. Interpolation of the Beam Centerline. The relative displacement between two nodes of an element after deformation can be expressed as

[r.sub.a] = [r.sub.2] - [r.sub.1] = [chi][e.sub.a], [chi] = [parallel][r.sub.2] - [r.sub.1][parallel]. (34)

If the beam element is small enough, the centerline after deformation could be regarded as a circular arc, as shown in Figure 3. The angle between the two base vectors [e.sub.1s] and [e.sub.2s] is denoted by 2[theta]; that is, cos(2[theta]) = [e.sup.T.sub.1s] [e.sub.2s]. The current element length L is related to the chord length x, which yields

L = [eta][chi], (35)

where [eta] = [theta]/ sin([theta]). In an element, the angular 2[theta] cannot be greater than 2[pi]. Thus, [eta] is definitely positive and approaches 1 when 2[theta] tends to 0. This conclusion conforms to the reality.

To fulfill C1 continuity at the element boundaries, the beam centerline of an element is assumed by the Hermite interpolation. Due to the situation that beams undergo large displacements and rotations, three components of the position vector on the beam centroid line need to be treated in the same way:

r = [N.sub.1][r.sub.1] + [N.sub.2][r.sub.2] + L ([N.sub.3][e.sub.1s] + [N.sub.4][e.sub.2s]), (36)

where [r.sub.1], [r.sub.2] are nodal position vectors and [e.sub.1s], [e.sub.2s] are unit tangent vectors of beam centerline at the two nodes. The influence of axial deformation has been considered in the expression of the current element length L, as shown in (38). The Hermite shape functions are

[mathematical expression not reproducible]. (37)

Equation (36) for the interpolation of the beam centerline is of great help for calculating the nodal bending curvatures.

4.3. Resolving Nodal Strain Measures. As a kind of slender structure, the beam axial tensile and torsional strains are small. The axial strain could be assumed as a constant along the element centerline; that is,

[epsilon] = L/[L.sub.0] - 1. (38)

Based on Euler-Bernoulli beam theory, the cross-sectional normal vector [e.sub.s] can be derived by taking the derivative of (36) with respect to [s.sub.0]. The first- and second-order derivatives of [e.sub.s] with respect to arc-length [s.sub.0] are

[mathematical expression not reproducible]. (39)

The bending curvatures are expressed by cross-sectional base vectors and their derivatives with respect to arc-length [s.sub.0] referring to (11):

[[kappa].sub.t] = -[e.sup.T.sub.b][e'.sub.s], [[kappa].sub.b] = [e.sup.T.sub.t][e'.sub.s]. (40)

Referring to (11) and (25), the derivatives of bending curvatures in (40) with respect to arc-length [s.sub.0] are given as follows:

[mathematical expression not reproducible]. (41)

Taking the derivatives of shape functions, the derivatives of base vector on the beam ends can be obtained referring to (39):

[mathematical expression not reproducible]. (42)

Nodal curvature parameters in two perpendicular bending directions are resolved by substituting (42) with (40)-(41):

[mathematical expression not reproducible]. (43)

Nodal curvature parameters of a beam element are rewritten as

[mathematical expression not reproducible]. (44)

With the given process, nodal curvature parameters and their derivatives with respect to arc-length parameter [s.sub.0] are expressed by nodal coordinates.

4.4. Discretization of Curvature Components in Space. In the original theory of geometrically exact beam, the analytical expression of centerline curvatures is almost impossible to be calculated by the nodal rotational parameters. Thus, a numerical integration technique is implemented. Unfortunately, the numerical integration needs to be taken for all elements at every step. In this work, the internal centerline curvatures are given directly by Hermitian interpolation, which makes the curvature expression explicit and analytically integrable. This property can benefit the calculation of explicit nodal forces.

The torsional curvature could be assumed as a constant in the beam element. Basis vector [e.sub.2t] has three components in the left end cross-sectional frame after deformation, as shown in Figure 4:

[e.sub.2t] = [c.sub.1][e.sub.1s] + [c.sub.2][e.sub.1t] + [c.sub.3][e.sub.1b], (45)

where [c.sub.2] = [([e.sub.2t]).sup.T][e.sub.1t] and [c.sub.3] = [([e.sub.2t]).sup.T][e.sub.1b].

The relative torsional angle between the left and right end sections of a beam element can be calculated by the arcsine function according to the geometrical relationship:

[alpha] = arcsin [c.sub.3]/[square root of ([c.sup.2.sub.2 + [c.sup.2.sub.3])]. (46)

Torsional curvature in the element can be approximated as the average:

[[kappa].sub.s] = [alpha]/[L.sub.0]. (47)

Nodal bending curvatures and their arc-length derivatives could be consistently derived from the displacement field in Sections 4.2 and 4.3. It affords a possibility to employ cubic Hermite polynomials for bending curvatures [[kappa].sub.t] and [[kappa].sub.b]:

[mathematical expression not reproducible]. (48)

Cubic Hermite interpolation for bending curvatures has third-order accuracy, which can make the interpolation curve comparatively smooth. With the above description, four independent strain measures in the cross-sectional reference frame are discretized in space.

4.5. Time Derivatives of Strain Measures. Rates of the strains can also be formulated by nodal parameters. The time derivative of the current beam element length in (35) is expressed as

[mathematical expression not reproducible], (49)

where

[mathematical expression not reproducible]. (50)

The time derivative of the axial strain can be given referring to (49):

[mathematical expression not reproducible]. (51)

The time derivative of the torsional curvature equation (47) is given in the following equation:

[[??].sub.s] = 1/[L.sub.0]([[tau].sup.T]([[omega].sub.2] - [[omega].sub.1])), (52)

where [tau] = (1/([c.sup.2.sub.2] + [c.sup.2.sub.3]))[[??].sub.2t]([c.sub.2][e.sub.1b] - [c.sub.3][e.sub.1t]). Taking the time derivative of nodal curvatures (43), the rates of the nodal curvature parameters can be given as follows:

[mathematical expression not reproducible]. (53)

5. Discrete Governing Equations

In this section, the discrete governing equations are formulated based on the dynamic equilibrium equations introduced in Section 3. Virtual power of a beam element is determined by the sum of the virtual power of internal, inertia, and external forces. These forces are introduced separately.

5.1. Elastic Nodal Forces. Using the interpolations of the strain measures (38) and (47)-(48), the internal virtual power of an element is calculated by integrating (31) along the beam domain as follows:

[mathematical expression not reproducible], (54)

where the coefficient matrix is

[mathematical expression not reproducible]. (55)

By substituting the virtual variations of the strain rates in (51)-(53) into (54), the virtual power produced by internal forces of the element can be rewritten in the following form:

[mathematical expression not reproducible], (56)

where [f.sub.1], [f.sub.2], [m.sub.1], and [m.sub.2] are corresponding nodal forces and moments, which are power-conjugated to [[??].sub.1], [[??].sub.2], [[omega].sub.1], and [[omega].sub.2], respectively. Expression of the left end nodal force is

[mathematical expression not reproducible]. (57)

The right end nodal force is

[mathematical expression not reproducible]. (58)

The left end nodal moment of the element is

[mathematical expression not reproducible]. (59)

The right end nodal moment of the element is

[mathematical expression not reproducible], (60)

where the nodal parameters are defined in matrix form as

[mathematical expression not reproducible]. (61)

Examining these expressions of nodal forces, it is easy to verify that

[f.sub.1] + [f.sub.2] = 0, [m.sub.1] + [m.sub.2] + ([[??].sub.2] - [[??].sub.1]) [f.sub.2] = 0. (62)

It can be proved that nodal forces described in (62) constitute an equilibrium force system. Hence, the presented beam element in this work is applicable to dynamic analysis of flexible multibody systems. Furthermore, the stiffness matrix of the proposed beam element is symmetric. Nodal forces were derived in the current configuration taking into consideration the coupled effects of the axial tensile, torsional, and bending deformations. Finally, the proposed beam element allows the explicit expression of nodal forces to be determined without using the numerical integration. This has a distinct advantage in the numerical solution of governing equations.

5.2. Mass Matrix. The displacement and rotation of the beam element are not directly discretized in space. Hence, the mass matrix and inertia force of a beam element are given by independent interpolations in this section. The position vector of an arbitrary point [r.sub.a] in a rigid cross section of Euler-Bernoulli beam could be expressed as

[r.sub.a] = r + [r.sub.p], (63)

where r is the position vector of the centroid and [r.sub.p] is the relative position vector with respect to the centroid. The virtual power of the inertia forces of the beam element can be written as

[mathematical expression not reproducible]. (64)

Using global nodal coordinates, the equation of flexible multibody system dynamics contains no centripetal force or Coriolis force.

If the acceleration caused by the axial extension is ignored, the velocity and acceleration of the beam centerline are given referring to (36):

[mathematical expression not reproducible]. (65)

The change of the cross-sectional angular velocity in a beam element is usually small. Hence, a linear interpolation is employed; that is, the angular velocity of the beam cross section is expressed in the cross-sectional frame as follows:

[??] = (1 - [xi]) [[??].sub.1] + [xi][[??].sub.2], (66)

where [[??].sub.1] and [[??].sub.2] are angular velocities of the nodes. Angular velocities in inertia coordinate system are

[mathematical expression not reproducible]. (67)

The inertia tensor of an arbitrary cross section in its cross-sectional frame can be expressed as [??] = diag([J.sub.s], [I.sub.t], [I.sub.b]). The virtual power of inertia forces in (64) can be rewritten as

[mathematical expression not reproducible], (68)

where [mathematical expression not reproducible], M is the mass matrix of the beam element, and [F.sub.i] is the inertia force:

[mathematical expression not reproducible]. (69)

5.3. The Governing Equations. Based on the principle of virtual power, the virtual power produced by external forces is equal to the sums of the virtual powers of the internal and inertia forces:

[mathematical expression not reproducible], (70)

where [F.sub.a] is the generalized external force. The governing equation of the beam element can be written as

M[??] = [F.sub.a] - [F.sub.i] - [F.sub.e]. (71)

The state variables of the governing equation are the nodal coordinates q in (32) and the general velocity [??] in (33):

[mathematical expression not reproducible]. (72)

The governing equation can be rewritten as a set of first-order differential equations, which can be solved by the well-developed numerical methods for ordinary differential equations:

[mathematical expression not reproducible]. (73)

6. Tests of the Proposed Beam Element

As mentioned in Section 5.1, an explicit expression of nodal forces without using the numerical integration can be given by the proposed large deformation beam element. Therefore, a linear beam element can be evolved. Because of the additivity of strain measures, the path independence of numerical solutions has been proved by Zupan and Saje [7]. In addition, the objectivity of strain measures and the convergence of the beam element are tested.

6.1. Degenerate Linear Beam Element. With the assumption of small displacements, the equilibrium equation is formulated in the initial configuration. Simplifications for the basis vectors of cross-sectional frames on the beam ends are given by

[mathematical expression not reproducible], (74)

where the frame x-y-z is the initial cross-sectional reference frame of a straight uniform cross-sectional element. Hence, [eta] [approximately equal to] 1. Linearizing the relative displacement in (34), one obtains

[r.sub.2] - [r.sub.1] [approximately equal to] [L.sub.0][e.sub.x], [chi] [approximately equal to] [L.sub.0]. (75)

As [c.sub.2] = [([e.sub.2t]).sup.T] [e.sub.1t] [approximately equal to] 1 and [c.sub.3] = [([e.sub.2t]).sup.T] [e.sub.1b] [approximately equal to] 0, the following approximations could be made for the intermediate parameters:

[l.sub.w] = 0, [l.sub.s] = [e.sub.x], [tau] = [e.sub.x]. (76)

Neglecting the [(o).sup.2] order of strains, the linearized nodal forces of the beam element are expressed as follows:

[mathematical expression not reproducible]. (77)

Under small deformation, the linearized nodal forces still satisfy the equilibrium equations (62). The parameters [[psi].sub.ij], [[PHI].sub.ij] (i = 1,2, j = t, b), [[kappa].sub.s], and [epsilon] are corresponding bending, torsional, and axial strains, respectively. The proposed element can degenerate to classical linear beam element.

6.2. Objectivity of Strain Measures. In order to test the objectivity of strain measures, the rigid body motions, including the translation [r.sub.r] and the rotation [R.sub.r], are superimposed onto a deformed beam element. Since [R.sub.r][R.sub.r.sup.T] = E is a unit matrix, the increment of intermediate variables and strain measures are zeros referring to (47), (38), and (43):

[mathematical expression not reproducible]. (78)

The discrete bending curvatures given by (48) are also unaffected by the rigid body motions:

[mathematical expression not reproducible]. (79)

It is proved that the strain measures in material reference frame are invariant under the rigid body motions.

6.3. Convergence and Patch Test of the Beam Element. In this section, the convergence of the proposed beam element is tested. The proposed beam formulation is displacement-based, considering its primary variables. Since the interpolation curve of the beam centerline has C1 continuity, the finite element solution converges to the exact solution when the length of a beam element tends to zero.

Instead of displacement and rotation parameters, strain measures are the discretized parameters in the proposed beam element. Hence, patch test is needed to verify the convergence of the element. In this section, patch tests of the three kinds of constant strains are given below.

Without loss of generality, an arbitrary finite element mesh for a beam is given in Figure 5. The node i connects the elements a and b. [L.sub.a] and [L.sub.b] are lengths of the elements. The material properties are uniform. Loads and displacements of these nodes are chosen to make sure that strains of the elements are constant. No body force is applied on the elements and no external force is acted on the node i. The element will be convergent if the elastic force on the node i satisfies the equilibrium equation.

(1) The axial strain in elements [epsilon] = [c.sub.1] is a constant, and other strain measures are all zeros. A pair of equal and opposite forces along X-axis is applied on nodes j and k.

Displacements and rotations of the nodes should be chosen from the following modes:

[mathematical expression not reproducible], (80)

where r with the corresponding subscript are displacements of the nodes, the first subscript of basis vectors [e.sub.is] is the number of node, and the second subscript is the number of base vector. Parameters [c.sub.i] (i = 1,2,3,...) are constants. The curvature parameters are zeros in the aforementioned displacement modes. By substituting (80) into the expression of nodal forces, the node forces on node i have only one nonzero component which is induced by the axial strains:

[mathematical expression not reproducible]. (81)

Thus, the following equilibrium equation of node i is deduced:

[F.sub.i] = [f.sub.ai] + [f.sub.bi] = 0, [M.sub.i] = [m.sub.ai] + [m.sub.bi] = 0. (82)

(2) The torsional curvature along the elements is constant; that is, [[kappa].sub.b] = [c.sub.2], and other strain measures are all zeros. A pair of equal and opposite torques along X-axis is applied on nodes j and fc. Displacements and rotations of the nodes are chosen from the following modes:

[mathematical expression not reproducible], (83)

where [angle]([e.sub.it], [e.sub.jt]) is the angle between two base vectors [e.sub.it] and [e.sub.jt]. The base vectors [e.sub.is], [e.sub.js], and [e.sub.ks] are along one straight line, which yields [tau] = [e.sub.is]. Therefore,

[m.sub.ai] = GJ[[kappa].sub.s][tau], [m.sub.bi] = -GJ[[kappa].sub.s][tau], [f.sub.ai] = [f.sub.bi] = 0. (84)

The equilibrium equations of node i in (82) are satisfied.

(3) The bending curvature along the elements is constant; that is, [[kappa].sub.b] = [c.sub.3], and other strain measures are all zero. A pair of equal and opposite moments [M.sub.ex] along Z-axis is acted on nodes j and k. Displacements and rotations of the nodes are chosen from the following modes:

[mathematical expression not reproducible]. (85)

The derivatives of the curvature components with respect to arc-length coordinate tend to zero when the angles [[theta].sub.1], [[theta].sub.2] are infinitesimal. The nonzero parameters for the pure bending have the following relationships:

[mathematical expression not reproducible]. (86)

By substituting (86) into node forces formulas, it obtains

[mathematical expression not reproducible]. (87)

The equilibrium equations of node i in (82) are fulfilled.

The equilibrium equations of node i verify that the proposed element passes the patch test when the element size becomes infinitesimal. Hence, the proposed element is convergent.

The Hermitian interpolation of the beam centerline is suitable for linear beams. However, it gives lower precision for large curvatures. The projection of the vector [e'.sub.s] on the base vector [e.sub.s] for Euler-Bernoulli beams is theoretically zero. However, the projection at the beam end node [e.sup.T.sub.1s][e'.sub.1s] = (6/[L.sub.0]L)[e.sup.T.sub.1s][r.sub.a] + (1/[L.sub.0])[e.sup.T.sub.1s](-4[e.sub.1s] - 2[e.sub.2s]) equals zero only when their included angles are infinitesimal. More elements need to be employed for beams with large deformation to derive the real curvatures.

7. Numerical Examples

Four numerical examples are examined to establish the validity, accuracy, and the convergence performance of the proposed beam element. First three examples deal with static analysis of a cantilever beam, an initially deformed beam, and a spatial beam structure. The last example is a dynamic problem of a pendulum with the stiffening effect taken into consideration.

The procedure is written in Matlab computing environment. The fsolve function in Optimization Toolbox of Matlab was employed to solve static problems. TolFun was set to [10.sup.-18] as the size of the latest change in sum of squares of function values and gradient of the sum of squares. TolX was set to [10.sup.-18], which is the lower bound of change in a step size. Iterations were terminated when the last step was smaller than TolFun or TolX. The ode15s function was employed to solve dynamic equations. The initial step was set as [10.sup.-4]. Error control properties were chosen as relative error tolerance RelTol = [10.sup.-4] and absolute error tolerance AbsTol = [10.sup.-4].

Example 1 (a cantilever beam subjected to an end moment). An initially straight cantilever beam is subjected to a concentrated moment M at the free end, as shown in Figure 6. The cantilever beam undergoing large deformation, which can be described as of pure bending, was investigated by many researchers [20-22], because an analytical solution exists in this example. Here, curvature radius 1/R is related to the applied moment M and the bending stiffness 1/R = M/EI.

When M satisfies the equation [M.sub.i] = ML/2[pi]EI = 1, the beam is curved into a complete circle and rotational angle of the end section is 2[pi]. The cantilever beam is divided into N elements. The rotational angle of the end section is given in Table 1. Furthermore, the errors with respect to the analytical solution are also given in Table 1, which shows that errors reduce when number of elements increases. Figure 7 shows deformed shapes of the cantilever beam when subjected to different end moments [M.sub.i]. In this particular case, the cantilever beam is divided into 10 elements. Figure 8 shows results for dimensionless vertical u/L and horizontal v/L tip displacements as a function of dimensional load ML/2[pi]EI, which are also compared with the analytical results. It can be seen from the comparison that numerical and analytical results are in good agreement.

Example 2 (cantilever beam initially curved to [pi]/4). A cantilever beam initially 45-degree bending subjected to a concentrated end load is shown in Figure 9. This example was first proposed by Bathe and Bolourchi [23] and has been employed by numerous researchers [5,10,11,21,24].

The bending beam lies in the X-Y plane and has an average radius of 100 m, as shown in Figure 9. Initial tip position is (70.71 m, 29.29 m, 0). A vertical concentrated load P = 600 N was applied at the tip along the Z direction. The beam cross section was a square with an area of 1 [m.sup.2]. The material was linear elastic. Elastic modulus and Poisson's ratio were E = [10.sup.7] Pa and v = 0, respectively.

Initially, the structure can be regarded as a curved beam with constant curvature, which maintains equilibrium under load and displacement boundary conditions. Eight equal straight beam elements were used to approximate the curved beam. An increasing magnitude of concentrated load on the beam end tip was applied in 40 load step increments. Figure 10 shows the evolution of tip displacement as a function of the applied load modulus. This result was compared to results reported by Bathe and Bolourchi [23], in which beam shear deformations were neglected. In Bathe's results, the Newton-Cotes formula of order 3 * 3 * 3 and Gauss integration of order 2 * 2 * 2 were employed for the isoparametric elements. Our results coincide with results reported by Bathe, as shown in Figure 8. Relative error of the numerical solution was found to increase with increase in displacement amplitude. Displacement magnitude of the beam end tip under 600 N force was determined as 73.5308 m in this work whereas the reference value is 73.022 m in Bathe's work, as shown in Table 2. Relative error was 0.7%. In the proposed beam element, assumptions about displacement model have been given up. The proposed beam element avoids overstiff conditions and tends to be flexible.

Example 3 (spatial beam structure under coupling deformation). A spatial beam structure was anchored on the support, as shown in Figure 11. Three straight beams of 1 m were connected perpendicular to each other. Cross section was a square with an area of 0.01 [m.sup.2]. Two forces were loaded on the tip along the X and Z directions. Elastic modulus and Poisson's ratio were E = [10.sup.6] Pa and v =0, respectively.

Every straight beam was divided into 4 elements. 80 load steps were used to observe load-deflection curves. As shown in Figure 12, load-deflection curves of the beam tip were compared to results from absolute nodal coordinate formulation (ANCF) [25]. Force-displacement curves were found to agree with results of ANCF. Displacement error compared to results in [25] was found to reduce with more elements in every straight beam, as shown in Table 3.

Example 4 (flexible pendulum subjected to spin-up motion). This example considers a clockwise spinning motion of a flexible pendulum, as shown in Figure 13. The pendulum was attached to a shaft orthogonal to the beam axis and driven by a prescribed rotation. Left end of the beam was on the origin of the inertia coordinate system. Angular velocity is given by

[mathematical expression not reproducible]. (88)

This is a popular example to test geometric stiffening effect and has been studied by many researchers [20,26, 27]. This spin-up model is usually motivated by a spacecraft antenna and helicopter blade. The beam has a length of L = 8 m. Cross-sectional area is A = 7.3 x [10.sup.-5] [m.sup.2] and the second-area moment of inertia is [I.sub.t] = 8.218 x [10.sup.-9] [m.sup.4]. Elastic modulus used was E = 6.895 x [10.sup.10] Pa. Mass density was [rho] = 2766.67 kg/[m.sup.3]. First three natural frequencies were determined as 2.9 Hz, 18.2 Hz, and 51.1 Hz. The pendulum reached a steady-state angular velocity after the spin-up time T = 15 s. For this simulation, gravity was not considered.

The pendulum was divided into 16-beam element. Steady-state angular velocity of [OMEGA] = 4 rad/s was used for simulation. Transversal deformation of the tip is given in Figure 14(a). A comparison of the deformation curve agrees well with the results of [27]. In the first 15 seconds, angular acceleration due to the inertia forces caused bending deformation of the flexible pendulum. When angular velocity was in steady state, centrifugal force straightened the flexible beam. Steady-state angular velocity [OMEGA] = 6 rad/s was also tested. Transversal deformation of the tip is given in Figure 14(b). Transversal deformation decreases under higher rotational speed, which can be observed from the results shown in Figures 14(a) and 14(b). Geometric stiffening effects could be reflected in this dynamic simulation.

8. Conclusions

In this paper, a new 3D geometrically exact Euler-Bernoulli beam element was proposed based on a hybrid interpolation of the beam centerlines and curvatures. The interpolation of the beam centerlines determined nodal curvatures, where the C1 continuity was satisfied, and nodal strain measures were consistently formulated. Analytical expression of nodal forces can be obtained using this method, which has an advantage in numerical simulation. While deformation is small, the proposed beam element can degenerate into classical linear beam element. Furthermore, objectivity and path independence of strain measures were satisfied by the cubic Hermitian interpolation of the beam curvatures. Through four numerical examples, validity of the proposed beam element was proved by the static analysis of 3D Euler-Bernoulli elements with large deformation and dynamic analysis of flexible multibody systems.

http://dx.doi.org/10.1155/2016/8980676

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

The authors are grateful for the National Natural Science Foundation of China (11372057). Their sincere thanks are due to Gang Wang for assistance with the revision of this paper.

References

[1] M. A. Crisfield and G. Jelenic, "Objectivity of strain measures in the geometrically exact three-dimensional beam theory and its finite-element implementation," The Royal Society of London: Proceedings Series A: Mathematical, Physical and Engineering Sciences, vol. 455, no. 1983, pp. 1125-1147, 1999.

[2] E. Reissner, "On one-dimensional large-displacement finite-strain beam theory," Studies in Applied Mathematics, vol. 52, no. 2, pp. 87-95, 1973.

[3] E. Reissner, "On one-dimensional finite-strain beam theory: the plane problem," Zeitschrift fur angewandte Mathematik und Physik ZAMP, vol. 23, no. 5, pp. 795-804, 1972.

[4] J. C. Simo, "A finite strain beam formulation. The three-dimensional dynamic problem. Part I," Computer Methods in Applied Mechanics and Engineering, vol. 49, no. 1, pp. 55-70, 1985.

[5] J. C. Simo and L. Vu-Quoc, "A three-dimensional finite-strain rod model. Part II: computational aspects," Computer Methods in Applied Mechanics and Engineering, vol. 58, no. 1, pp. 79-116, 1986.

[6] I. Romero, "The interpolation of rotations and its application to finite element models of geometrically exact rods," Computational Mechanics, vol. 34, no. 2, pp. 121-133, 2004.

[7] D. Zupan and M. Saje, "Finite-element formulation of geometrically exact three-dimensional beam theories based on interpolation of strain measures," Computer Methods in Applied Mechanics and Engineering, vol. 192, no. 49-50, pp. 5209-5248, 2003.

[8] F. Armero and J. Valverde, "Invariant Hermitian finite elements for thin Kirchhoff rods. I: the linear plane case," Computer Methods in Applied Mechanics and Engineering, vol. 213-216, pp. 427-457, 2012.

[9] C. Meier, A. Popp, and W. A. Wall, "An objective 3D large deformation finite element formulation for geometrically exact curved Kirchhoff rods," Computer Methods in Applied Mechanics and Engineering, vol. 278, pp. 445-478, 2014.

[10] J. B. Jonker and J. P. Meijaard, "A geometrically non-linear formulation of a three-dimensional beam element for solving large deflection multibody system problems," International Journal of Non-Linear Mechanics, vol. 53, pp. 63-74, 2013.

[11] G. Jelenic and M. A. Crisfield, "Geometrically exact 3D beam theory: implementation of a strain-invariant finite element for statics and dynamics," Computer Methods in Applied Mechanics and Engineering, vol. 171, no. 1-2, pp. 141-171, 1999.

[12] C. Meier, A. Popp, and W. A. Wall, "A locking-free finite element formulation and reduced models for geometrically exact Kirchhoff rods," Computer Methods in Applied Mechanics and Engineering, vol. 290, pp. 314-341, 2015.

[13] S. Ghosh and D. Roy, "A frame-invariant scheme for the geometrically exact beam using rotation vector parametrization," Computational Mechanics, vol. 44, no. 1, pp. 103-118, 2009.

[14] R. Y. Yakoub and A. A. Shabana, "Three dimensional absolute nodal coordinate formulation for beam elements: implementation and applications," Journal of Mechanical Design, vol. 123, no. 4, pp. 614-621, 2001.

[15] J. T. Sopanen and A. M. Mikkola, "Description of elastic forces in absolute nodal coordinate formulation," Nonlinear Dynamics, vol. 34, no. 1, pp. 53-74, 2003.

[16] A. A. Shabana and R. Y. Yakoub, "Three dimensional absolute nodal coordinate formulation for beam elements: theory," Journal of Mechanical Design, vol. 123, no. 4, pp. 606-613, 2001.

[17] A. A. Shabana, "Computer implementation of the absolute nodal coordinate formulation for flexible multibody dynamics," Nonlinear Dynamics, vol. 16, no. 3, pp. 293-306, 1998.

[18] O. Dmitrochenko and A. Mikkola, "A formal procedure and invariants of a transition from conventional finite elements to the absolute nodal coordinate formulation," Multibody System Dynamics, vol. 22, no. 4, pp. 323-339, 2009.

[19] D. Zupan and M. Saje, "A new finite element formulation of three-dimensional beam theory based on interpolation of curvature," CMES--Computer Modeling in Engineering, vol. 4, no. 2, pp. 301-318, 2003.

[20] J. M. Mayo, D. Garcia-Vallejo, and J. Dominguez, "Study of the geometric stiffening effect: comparison of different formulations," Multibody System Dynamics, vol. 11, no. 4, pp. 321-341, 2004.

[21] M. A. Crisfield, "A consistent co-rotational formulation for nonlinear, three-dimensional, beam-elements," Computer Methods in Applied Mechanics and Engineering, vol. 81, no. 2, pp. 131-150, 1990.

[22] P. Nanakorn and L. N. Vu, "A 2D field-consistent beam element for large displacement analysis using the total Lagrangian formulation," Finite Elements in Analysis and Design, vol. 42, no. 14-15, pp. 1240-1247, 2006.

[23] K. Bathe and S. Bolourchi, "Large displacement analysis of three-dimensional beam structures," International Journal for Numerical Methods in Engineering, vol. 14, no. 7, pp. 961-986, 1979.

[24] H.-G. Kwak, D.-Y. Kim, and H.-W. Lee, "Effect of warping in geometric nonlinear analysis of spatial beams," Journal of Constructional Steel Research, vol. 57, no. 7, pp. 729-751, 2001.

[25] I. Romero, "A comparison of finite elements for nonlinear beams: the absolute nodal coordinate and geometrically exact formulations," Multibody System Dynamics, vol. 20, no. 1, pp. 51-68, 2008.

[26] A. K. Banerjee and J. M. Dickens, "Dynamics of an arbitrary flexible body in large rotation and translation," Journal of Guidance, Control, and Dynamics, vol. 13, no. 2, pp. 221-227, 1990.

[27] J. Ryu, S.-S. Kim, and S.-S. Kim, "A criterion on inclusion of stress stiffening effects in flexible multibody dynamic system simulation," Computers & Structures, vol. 62, no. 6, pp. 1035-1048, 1997.

[28] A. Cardona and M. Geradin, "A beam finite element nonlinear theory with finite rotations," International Journal for Numerical Methods in Engineering, vol. 26, no. 11, pp. 2403-2438, 1988.

Huiqing Fang and Zhaohui Qi

State Key Laboratory of Structural Analysis for Industrial Equipment, Faculty of Vehicle Engineering and Mechanics, Dalian University of Technology, Dalian 116024, China

Correspondence should be addressed to Zhaohui Qi; zhaohuiq@dlut.edu.cn

Received 1 November 2015; Revised 20 January 2016; Accepted 4 February 2016

Caption: Figure 1: A spatial Euler-Bernoulli beam in the inertial system.

Caption: Figure 2: An infinitesimal beam.

Caption: Figure 3: The arc and chord lengths of a beam element.

Caption: Figure 4: Torsional angle of the right end section relative to the left one.

Caption: Figure 5: The mesh for patch test.

Caption: Figure 6: A cantilever beam subjected to concentrated end moment.

Caption: Figure 7: Deformed geometric shapes of the cantilever beam subjected to different end moments.

Caption: Figure 8: Dimensionless deflection of tip-loaded cantilever beam.

Caption: Figure 9: Cantilever beam initially curved to 45[degrees] subjected to a concentrated tip force.

Caption: Figure 10: Load-displacement curves of the tip.

Caption: Figure 12: Tip force-displacement curves.

Caption: Figure 13: Flexible pendulum subjected to spin-up motion.

Caption: Figure 14: Transversal deformation of the tip driven by different rotational velocities.
```Table 1: Rotational angle of the end section for a cantilever beam
with different numbers of element divisions.

Number of             Rotation of
elements                the tip       Absolute error

N = 6                    6.384            -0.1011
N = 8                    6.313            -0.0305
N =10                    6.295            -0.0123
N =12                    6.289            -0.0059
N = 20                   6.283      -7.6 * [10.sup.-4]
N = 30                   6.283      -1.5 * [10.sup.-4]
N = 40                   6.283      -4.7 * [10.sup.-5]
Analytical solution      6.283

Number of
elements                Relative error

N = 6                       -0.0161
N = 8                       -0.0048
N =10                       -0.002
N =12                 -9.4 * [10.sup.-4]
N = 20                -1.2 * [10.sup.-4]
N = 30                -2.4 * [10.sup.-5]
N = 40                -7.6 * [10.sup.-6]
Analytical solution

Table 2: Initial curved cantilever. Position components of the tip for
end force of 600 N with 8 elements. The initial position was
(70.71, 29.29, 0) (m).

Formulation                                           rx (m)   ry (m)

Bathe and Bolourchi [23]                               47.2    -15.9
Simo and Vu-Quoc [5]       Geometrically exact beam   47.23    -15.79
Cardona and Geradin [28]   Geometrically exact beam   47.04    -15.56
Romero [25]                Geometrically exact beam   46.98    -15.69
Romero [25]                Absolute nodal method      50.16    -17.16
Present                                               47.87    -14.07

Formulation                                           rz (m)

Bathe and Bolourchi [23]                               53.4
Simo and Vu-Quoc [5]       Geometrically exact beam   53.37
Cardona and Geradin [28]   Geometrically exact beam   53.50
Romero [25]                Geometrically exact beam   53.50
Romero [25]                Absolute nodal method      50.26
Present                                               54.01

Table 3: Spatial beam structure. Displacement components of the
tip under force = 5 N with different numbers of element divisions.

Ux (m)    Uy (m)   Uz (m)

Result in [25]   -1.3568   0.197     -1.8
N = 4            -1.3458   0.1983   -1.8563
N = 5            -1.3478   0.1979   -1.8453
N = 6            -1.3490   0.1975   -1.8380
```