# A Spatial Euler-Bernoulli Beam Element for Rigid-Flexible Coupling Dynamic Analysis of Flexible Structures

Authors: Zhigang Zhang; Zhaohui Qi; Zhigang Wu; Huiqing Fang
Journal: Shock and Vibration
 (2015)

---

## Abstract

A two-node spatial beam element with the Euler-Bernoulli assumption is developed for the nonlinear dynamic analysis of slender beams undergoing arbitrary rigid motions and large deformations. During the analysis, the global displacement and rotation vectors with six degrees of freedom are selected as the nodal coordinates. In addition, the “shear locking” problem is avoided successfully since the beam cross-sections are always perpendicular to the current neutral axes by employing a special coupled interpolation of the centroid position and the cross-section orientation. Then a scheme is presented where the original transient strains representing the nodal forces are replaced by proposed average strains over a small time interval. Thus all the high frequencies can be filtered out and a corresponding equivalent internal damping will be produced in this new formulation, which can improve the computation performance of the proposed element for solving the stiff problem and evaluate the governing equations even by using the nonstiff ordinary differential equation solver. Finally, several numerical examples are carried out to verify the validation and efficiency of this proposed formulation by comparison with the analytical solutions and other research works.


---

## Body

## 1. Introduction

Lots of slender structures in engineering are composed by beams, such as framed structures, robot arms, large deployable space structures, and turbine propellers. Motivated by the higher standards and increasing demands in the field of engineering design, the rigid-flexible coupling dynamic analysis of the spatial beams undergoing large deformation has drawn great attentions in the past few decades [1–10].According to the rigid cross-section assumption, the dimension of beam model is reduced. Meanwhile, it needs paying careful attention to the parameterization and discretization of the finite rotation. As well known, the geometrically exact beam theory [4, 5, 11] allows formulating problems involving arbitrarily large displacements, rotations, and strains, which has provided the basis of many recent finite element formulations. In the original articles [5, 6], Simo and Vu-Quoc employed the global position and rotation vector as the nodal coordinates to independently interpolate the incremental displacement and rotation fields. Since then, many innovative numerical formulations of spatial beams with different discretization of the rotation fields have been proposed. Cardona and Geradin [1] advocated an interpolation of the material incremental rotations. Ibrahimbegović et al. [12] used a procedure to interpolate the total rotation vector directly. Jelenić and Crisfield [13] proposed a new formulation based on the interpolation of incremental local rotations. An alternative approach was presented by Betsch and Steinmann [10], in which the base vectors of the cross-section were interpolated. Such beam elements, in which the fields of displacement and rotation are interpolated independently, are generally referred to as the “Timoshenko beam” elements, however, accompanied by the “shear locking” phenomenon. Moreover, the rotation vectors are physically nonadditive quantities, and the direct discretization of rotations without spoiling the objectivity of rotational strain measures is also a challenging task [14, 15].With one of the purposes of sidestepping the problems caused by the finite rotation, the so-called absolute nodal coordinate formulation of beam was originally developed by Shabana and Yakoub [16, 17]. By selecting absolute coordinates and their slops, totally 12 degrees of freedom per node, the particularly designed shape function can represent arbitrary large rigid body motions exactly. This formulation produces a constant mass matrix but suffers Poisson locking problem [18].The Euler-Bernoulli beam model inherently requires that the cross-section should keep perpendicular to the tangent of the centerline during the deformation processes, which means that the translation displacement and rotation fields are coupled. It is easy to construct a set of field-consistent interpolation functions under the assumption of small displacements. Thus the Euler-Bernoulli beam elements are commonly used in the corotational technique [19–22]. However, few Euler-Bernoulli beam formulations are proposed for the structures with arbitrary rigid motions and large deformations. Nanakorn and Vu [23] developed a planar Euler-Bernoulli beam element for large displacement analysis using the total Lagrangian formulation. More recently, Zhao and Ren [7] proposed a singularity-free Euler-Bernoulli beam element, where each node includes totally eight degrees of freedom (three global positions, four Euler parameters, and one normal strain), to coupled interpolate the position and orientation. However, the resulting governing equations are a set of differential algebraic equations [7], which may lead to additional challenges during calculating the numerical solutions.Although Euler-Bernoulli hypothesis makes the construction of the interpolation function of the beam element more difficult, it provides an idea to interpolate the finite rotation field without direct discretization, which will benefit from shear-locking-free and objectivity of rotational strains. Based on this idea, a spatial Euler-Bernoulli element is developed in this paper for the rigid-flexible coupling dynamic analysis. In this element, each node has six generalized nodal coordinates, including the global displacement vector and the rotation vector. Considering the normal strain being a small quantity, we ignore the influence of the normal strain at two end points on the shape of the centerline and get the centroid position field through the nodal position vectors and the normal vectors of the end sections with Hermite interpolation. So the perpendicularity assumption at two end points is guaranteed and the tangent field of centerline is determined. Then the orientation of the cross-section is achieved by the sequential rotations to keep the cross-section perpendicular to the tangent vector. In this way, the Euler-Bernoulli beam assumption is satisfied and the “shear locking” problem will be avoided.The governing dynamic equations of the system with a large rigid body motion and small elastic vibration are nonlinear stiff differential equations. Because of the presence of high frequencies, whose response is an artifact of the spatial discretization that contains no information about the physical behavior of the system, the time integration of stiff equations becomes very difficult. Many researches therefore try to introduce a numerical dissipation in the high-frequency into the time integration algorithms, such as Newmark method [24], Wilson θ-method [25], and generalized α-method [26]. In this paper, an effective scheme is proposed for reducing the difficulty of solving the stiff governing equation from the finite element modeling. Based on this spatial beam element, the average strains over a small time interval are introduced to replace the original transient strains in the expression of the nodal forces. Then an innovative finite element model is derived, which can filter out all the high frequencies and introduce the equivalent internal damping. The numerical examples show that the governing equations can be solved even by using the nonstiff ordinary differential equation solver with an appropriate time interval parameter.This paper is organized as follows. In Section2 a brief review of the parameterization of finite rotations is introduced and the internal virtual power of a spatial Euler-Bernoulli beam is derived. Then the finite element implementation for this beam is provided in Section 3. In addition, several illustrative numerical examples are carried out in Section 4 to verify the validation of the developed beam element. Finally, some conclusions are presented in Section 5.


## 2. Geometry and Internal Virtual Power of Euler-Bernoulli Beam

### 2.1. Parameterization of the Finite Rotation, Angular Velocity, and Curvature Vector

According to the Bernoulli hypothesis, which states that “plane sections remain plane,” the geometry of the spatial beam is described by the centroid line and a family of the corresponding cross-sections as shown in Figure1, where a fixed global frame with orthonormal base vectors (
g
1
,
g
2
,
g
3
) is defined and the centroid line is described by the position vector r
(
s
) with s representing the arc-length coordinate of the reference configuration.Figure 1
Geometric configuration of the spatial beam.The cross-section, whose geometric shape is assumed to be arbitrary and unchanged along the beam, executes rotational motion during deformation. A family of orthonormal basis(
e
1
,
e
2
,
e
3
), called the cross-section basis or the material basis, is employed to describe the orientation of the cross-sections. The base vectors e
2 and e
3 are directed along the principal axes of inertia of the cross-section, and e
1 is the normal vector of the cross-section; that is, e
1
=
e
2
×
e
3. The relationship between the global base vector g
i and cross-section base vector e
i can be expressed as(1)
e
i
=
R
g
i
,
i
=
1,2
,
3
,where R is the rotation matrix with the following properties; that is,(2)
R
R
T
=
I
,
det
⁡
R
=
1
,where I is the identity matrix.The time derivative of (1) leads to(3)
e
˙
i
=
R
˙
g
i
=
R
˙
R
T
R
g
i
=
R
˙
R
T
e
i
.In this paper, the symbol (
˙
) denotes the time derivative. For (2), we have(4)
R
˙
R
T
=
-
R
R
˙
T
=
-
R
˙
R
T
T
.This means that R
˙
R
T is a skew-symmetric matrix and can be denoted as(5)
ω
~
=
R
˙
R
T
.Then substituting (5) into (3) results in(6)
e
˙
i
=
ω
~
e
i
=
ω
×
e
i
,where ω is the axial vector of skew-symmetric matrix ω
~, and it denotes the angular velocity vector.Similarly, taking derivative of (1) with respect to the arc-length coordinate s leads to(7)
e
i
′
=
κ
~
e
i
=
κ
×
e
i
,where the symbol ′ denotes the derivative with respect to the arc-length coordinate s, and the skew-symmetric matrix κ
~ is given below(8)
κ
~
=
R
′
R
T
.The axial vector κ corresponding to κ
~ is defined as the curvature vector.The curvature vectorκ can be expressed with respect to the material basis(9)
κ
=
κ
e
1
e
1
+
κ
e
2
e
2
+
κ
e
3
e
3
=
κ
e
i
e
i
,where Einstein notation is introduced. By substituting (9) into (7), we have(10)
e
j
′
=
κ
e
i
e
i
×
e
j
=
e
i
j
k
κ
e
i
e
k
,where e
i
j
k is the permutation symbol in 3D. For (10), we have the following relationship(11)
e
i
j
k
κ
e
i
=
e
j
′
·
e
k
.The time derivative of (11) could be expressed as(12)
e
i
j
k
κ
˙
e
i
=
e
˙
j
′
·
e
k
+
e
j
′
·
e
˙
k
=
ω
′
×
e
j
+
ω
×
e
j
′
·
e
k
+
e
j
′
·
ω
×
e
k
=
ω
′
·
e
j
×
e
k
=
e
i
j
k
ω
′
·
e
i
.Therefore, we can obtain the following equation:(13)
ω
′
=
κ
˙
e
i
e
i
.Similarly the angular velocity ω can be given as(14)
ω
=
ω
e
i
e
i
.By substituting (14) into (6), finally, we get(15)
e
i
j
k
ω
e
i
=
e
˙
j
·
e
k
,
κ
˙
=
ω
e
i
′
e
i
.The rotation matrix R could be written as a function of rotational vector θ; that is,(16)
R
θ
=
I
+
sin
⁡
ϑ
ϑ
θ
~
+
1
-
cos
⁡
ϑ
ϑ
2
θ
~
2
=
exp
⁡
θ
~
,where θ
~ denotes the skew-symmetric matrix corresponding to θ, and ϑ
≔
θ is the norm of the rotational vector. Formula (16) is also called the Rodrigues formula.The relationship between the angular velocity vectorω and the time derivative of the rotation vector θ
˙ can be expressed by(17)
ω
=
T
θ
θ
˙
.Similarly the curvature vector κ can be also expressed by the arc-length derivative of the rotational vector θ
′; that is,(18)
κ
=
T
θ
θ
′
,where T
θ denotes a matrix given by(19)
T
θ
=
I
+
1
-
cos
⁡
ϑ
ϑ
2
θ
~
+
ϑ
-
sin
⁡
ϑ
ϑ
3
θ
~
2
.According to (17), the time derivative of the rotational vector can be obtained by(20)
θ
˙
=
T
θ
-
1
ω
,where(21)
T
θ
-
1
=
I
-
1
2
θ
~
+
1
ϑ
2
1
-
ϑ
2
tan
⁡
ϑ
/
2
θ
~
2is the inverse matrix of T
θ. According to (20), the singularity problem will occur at ϑ
=
2
k
π
,
k
∈
Z
∖
0, because 2
k
π
/
2
tan
⁡
k
π tends to infinity. That is the major shortcoming while using the rotation vector to represent the finite rotation.Let’s consider a rotation vectorφ with a rotation angle which is greater than π and less than 2
π; that is, π
<
φ
<
2
π; then its complement rotation vector is defined by (22)
φ
∗
=
φ
-
2
π
φ
φ
.It can be seen that the magnitude of the new rotation vector is φ
∗
=
2
π
-
φ, which means 0
<
φ
∗
<
π. Note that the rotation vector φ and its complement φ
∗ represent the same rotation; that is, R
(
φ
∗
)
=
R
(
φ
). So the rotation could be represented by these two rotation vectors. When a rotation angle is greater than π, the change of parameterization according to (22) will be achieved. Thus, the singularity problems at the rotation angle 2
π will never be encountered. The similar procedure of passing the singularity of the rotation vector can be found in [1, 3, 12].


### 2.2. Internal Virtual Power of Euler-Bernoulli Beam

An infinitesimal element of a spatial beam is shown in Figure2, where the beam is subjected to the external distributed forces and moment vectors n
(
s
) and m
(
s
) per unit length of the reference line of centroids. And, in that figure, N
(
s
) and M
(
s
) are the stress-resultant force and resulting moment vectors over the cross-section, respectively.Figure 2
Equilibrium of an infinitesimal element of a beam.The equilibrium equations of the infinitesimal element can be written as(23)
A
ρ
r
¨
=
n
+
N
′
,
J
ρ
ω
˙
+
ω
×
J
ρ
ω
=
m
+
M
′
+
r
′
×
N
,with A
ρ and J
ρ representing the linear mass density and moment of inertia, respectively. The virtual kinetic power of the infinitesimal element is given by(24)
δ
p
i
n
e
=
δ
r
˙
·
r
¨
A
ρ
d
s
+
δ
ω
·
J
ρ
ω
˙
+
ω
×
J
ρ
ω
d
s
.By substituting (23) into (24), we can obtain(25)
δ
p
i
n
e
=
δ
r
˙
·
n
+
N
′
d
s
+
δ
ω
·
m
+
M
′
+
r
′
×
N
d
s
.The external virtual power of the infinitesimal element can be expressed as(26)
δ
p
e
x
t
=
δ
r
˙
·
n
d
s
+
δ
r
˙
·
N
′
d
s
+
δ
r
˙
′
·
N
d
s
+
δ
ω
·
m
d
s
+
δ
ω
·
M
′
d
s
+
δ
ω
′
·
M
d
s
.By using the principle of virtual power, which states that at any time t, the total virtual power of the external, internal, and inertia forces is zero in any admissible virtual state of motion, the internal virtual power of the infinitesimal element is given by(27)
δ
p
i
n
t
=
δ
p
e
x
t
-
δ
p
i
n
e
=
δ
r
˙
′
·
N
d
s
-
δ
ω
·
r
′
×
N
d
s
+
δ
ω
′
·
M
d
s
=
δ
r
˙
′
-
ω
×
r
′
·
N
d
s
+
δ
ω
′
·
M
d
s
.For the Euler-Bernoulli beam, the cross-section remains perpendicular to the tangent of the centerline during deformation; that is,(28)
r
′
=
1
+
ε
e
1
,where ε is the normal strain, and we also have(29)
r
˙
′
-
ω
×
r
′
=
ε
˙
e
1
.Substituting (13) and (29) into the internal virtual power equation (27) and calculating the integral over the range [
0
,
L
] leads to(30)
δ
p
i
n
t
=
∫
0
L
δ
ε
˙
N
e
1
+
δ
κ
˙
e
1
M
e
1
+
δ
κ
˙
e
2
M
e
2
+
δ
κ
˙
e
3
M
e
3
d
s
,where N
e
i and M
e
i are components of N and M given relative to the material basis; that is,(31)
N
=
N
e
1
e
1
+
N
e
2
e
2
+
N
e
3
e
3
,
M
=
M
e
1
e
1
+
M
e
2
e
2
+
M
e
3
e
3
.The constitutive equation of the Euler-Bernoulli beam can be described as(32)
N
e
1
=
E
A
ε
,
M
e
1
=
G
J
κ
e
1
,
M
e
2
=
E
I
e
2
κ
e
2
,
M
e
3
=
E
I
e
3
κ
e
3
.By substituting (32) into (30), the internal virtual power can be expressed as(33)
δ
p
i
n
t
=
∫
0
L
δ
ε
˙
E
A
ε
+
δ
κ
˙
e
1
G
J
κ
e
1
+
δ
κ
˙
e
2
E
I
e
2
κ
e
2
+
δ
κ
˙
e
3
E
I
e
3
κ
e
3
d
s
,where E
A and G
J are the axial stiffness and the torsion stiffness, respectively. E
I
e
2 and E
I
e
3 are the principal bending stiffness.This above introduced principle is usually regarded as the geometrically exact beam theory, which is firstly proposed by Reissner [11] and then has been mainly developed by Simo and Vu-Quoc [4, 5].


## 2.1. Parameterization of the Finite Rotation, Angular Velocity, and Curvature Vector

According to the Bernoulli hypothesis, which states that “plane sections remain plane,” the geometry of the spatial beam is described by the centroid line and a family of the corresponding cross-sections as shown in Figure1, where a fixed global frame with orthonormal base vectors (
g
1
,
g
2
,
g
3
) is defined and the centroid line is described by the position vector r
(
s
) with s representing the arc-length coordinate of the reference configuration.Figure 1
Geometric configuration of the spatial beam.The cross-section, whose geometric shape is assumed to be arbitrary and unchanged along the beam, executes rotational motion during deformation. A family of orthonormal basis(
e
1
,
e
2
,
e
3
), called the cross-section basis or the material basis, is employed to describe the orientation of the cross-sections. The base vectors e
2 and e
3 are directed along the principal axes of inertia of the cross-section, and e
1 is the normal vector of the cross-section; that is, e
1
=
e
2
×
e
3. The relationship between the global base vector g
i and cross-section base vector e
i can be expressed as(1)
e
i
=
R
g
i
,
i
=
1,2
,
3
,where R is the rotation matrix with the following properties; that is,(2)
R
R
T
=
I
,
det
⁡
R
=
1
,where I is the identity matrix.The time derivative of (1) leads to(3)
e
˙
i
=
R
˙
g
i
=
R
˙
R
T
R
g
i
=
R
˙
R
T
e
i
.In this paper, the symbol (
˙
) denotes the time derivative. For (2), we have(4)
R
˙
R
T
=
-
R
R
˙
T
=
-
R
˙
R
T
T
.This means that R
˙
R
T is a skew-symmetric matrix and can be denoted as(5)
ω
~
=
R
˙
R
T
.Then substituting (5) into (3) results in(6)
e
˙
i
=
ω
~
e
i
=
ω
×
e
i
,where ω is the axial vector of skew-symmetric matrix ω
~, and it denotes the angular velocity vector.Similarly, taking derivative of (1) with respect to the arc-length coordinate s leads to(7)
e
i
′
=
κ
~
e
i
=
κ
×
e
i
,where the symbol ′ denotes the derivative with respect to the arc-length coordinate s, and the skew-symmetric matrix κ
~ is given below(8)
κ
~
=
R
′
R
T
.The axial vector κ corresponding to κ
~ is defined as the curvature vector.The curvature vectorκ can be expressed with respect to the material basis(9)
κ
=
κ
e
1
e
1
+
κ
e
2
e
2
+
κ
e
3
e
3
=
κ
e
i
e
i
,where Einstein notation is introduced. By substituting (9) into (7), we have(10)
e
j
′
=
κ
e
i
e
i
×
e
j
=
e
i
j
k
κ
e
i
e
k
,where e
i
j
k is the permutation symbol in 3D. For (10), we have the following relationship(11)
e
i
j
k
κ
e
i
=
e
j
′
·
e
k
.The time derivative of (11) could be expressed as(12)
e
i
j
k
κ
˙
e
i
=
e
˙
j
′
·
e
k
+
e
j
′
·
e
˙
k
=
ω
′
×
e
j
+
ω
×
e
j
′
·
e
k
+
e
j
′
·
ω
×
e
k
=
ω
′
·
e
j
×
e
k
=
e
i
j
k
ω
′
·
e
i
.Therefore, we can obtain the following equation:(13)
ω
′
=
κ
˙
e
i
e
i
.Similarly the angular velocity ω can be given as(14)
ω
=
ω
e
i
e
i
.By substituting (14) into (6), finally, we get(15)
e
i
j
k
ω
e
i
=
e
˙
j
·
e
k
,
κ
˙
=
ω
e
i
′
e
i
.The rotation matrix R could be written as a function of rotational vector θ; that is,(16)
R
θ
=
I
+
sin
⁡
ϑ
ϑ
θ
~
+
1
-
cos
⁡
ϑ
ϑ
2
θ
~
2
=
exp
⁡
θ
~
,where θ
~ denotes the skew-symmetric matrix corresponding to θ, and ϑ
≔
θ is the norm of the rotational vector. Formula (16) is also called the Rodrigues formula.The relationship between the angular velocity vectorω and the time derivative of the rotation vector θ
˙ can be expressed by(17)
ω
=
T
θ
θ
˙
.Similarly the curvature vector κ can be also expressed by the arc-length derivative of the rotational vector θ
′; that is,(18)
κ
=
T
θ
θ
′
,where T
θ denotes a matrix given by(19)
T
θ
=
I
+
1
-
cos
⁡
ϑ
ϑ
2
θ
~
+
ϑ
-
sin
⁡
ϑ
ϑ
3
θ
~
2
.According to (17), the time derivative of the rotational vector can be obtained by(20)
θ
˙
=
T
θ
-
1
ω
,where(21)
T
θ
-
1
=
I
-
1
2
θ
~
+
1
ϑ
2
1
-
ϑ
2
tan
⁡
ϑ
/
2
θ
~
2is the inverse matrix of T
θ. According to (20), the singularity problem will occur at ϑ
=
2
k
π
,
k
∈
Z
∖
0, because 2
k
π
/
2
tan
⁡
k
π tends to infinity. That is the major shortcoming while using the rotation vector to represent the finite rotation.Let’s consider a rotation vectorφ with a rotation angle which is greater than π and less than 2
π; that is, π
<
φ
<
2
π; then its complement rotation vector is defined by (22)
φ
∗
=
φ
-
2
π
φ
φ
.It can be seen that the magnitude of the new rotation vector is φ
∗
=
2
π
-
φ, which means 0
<
φ
∗
<
π. Note that the rotation vector φ and its complement φ
∗ represent the same rotation; that is, R
(
φ
∗
)
=
R
(
φ
). So the rotation could be represented by these two rotation vectors. When a rotation angle is greater than π, the change of parameterization according to (22) will be achieved. Thus, the singularity problems at the rotation angle 2
π will never be encountered. The similar procedure of passing the singularity of the rotation vector can be found in [1, 3, 12].


## 2.2. Internal Virtual Power of Euler-Bernoulli Beam

An infinitesimal element of a spatial beam is shown in Figure2, where the beam is subjected to the external distributed forces and moment vectors n
(
s
) and m
(
s
) per unit length of the reference line of centroids. And, in that figure, N
(
s
) and M
(
s
) are the stress-resultant force and resulting moment vectors over the cross-section, respectively.Figure 2
Equilibrium of an infinitesimal element of a beam.The equilibrium equations of the infinitesimal element can be written as(23)
A
ρ
r
¨
=
n
+
N
′
,
J
ρ
ω
˙
+
ω
×
J
ρ
ω
=
m
+
M
′
+
r
′
×
N
,with A
ρ and J
ρ representing the linear mass density and moment of inertia, respectively. The virtual kinetic power of the infinitesimal element is given by(24)
δ
p
i
n
e
=
δ
r
˙
·
r
¨
A
ρ
d
s
+
δ
ω
·
J
ρ
ω
˙
+
ω
×
J
ρ
ω
d
s
.By substituting (23) into (24), we can obtain(25)
δ
p
i
n
e
=
δ
r
˙
·
n
+
N
′
d
s
+
δ
ω
·
m
+
M
′
+
r
′
×
N
d
s
.The external virtual power of the infinitesimal element can be expressed as(26)
δ
p
e
x
t
=
δ
r
˙
·
n
d
s
+
δ
r
˙
·
N
′
d
s
+
δ
r
˙
′
·
N
d
s
+
δ
ω
·
m
d
s
+
δ
ω
·
M
′
d
s
+
δ
ω
′
·
M
d
s
.By using the principle of virtual power, which states that at any time t, the total virtual power of the external, internal, and inertia forces is zero in any admissible virtual state of motion, the internal virtual power of the infinitesimal element is given by(27)
δ
p
i
n
t
=
δ
p
e
x
t
-
δ
p
i
n
e
=
δ
r
˙
′
·
N
d
s
-
δ
ω
·
r
′
×
N
d
s
+
δ
ω
′
·
M
d
s
=
δ
r
˙
′
-
ω
×
r
′
·
N
d
s
+
δ
ω
′
·
M
d
s
.For the Euler-Bernoulli beam, the cross-section remains perpendicular to the tangent of the centerline during deformation; that is,(28)
r
′
=
1
+
ε
e
1
,where ε is the normal strain, and we also have(29)
r
˙
′
-
ω
×
r
′
=
ε
˙
e
1
.Substituting (13) and (29) into the internal virtual power equation (27) and calculating the integral over the range [
0
,
L
] leads to(30)
δ
p
i
n
t
=
∫
0
L
δ
ε
˙
N
e
1
+
δ
κ
˙
e
1
M
e
1
+
δ
κ
˙
e
2
M
e
2
+
δ
κ
˙
e
3
M
e
3
d
s
,where N
e
i and M
e
i are components of N and M given relative to the material basis; that is,(31)
N
=
N
e
1
e
1
+
N
e
2
e
2
+
N
e
3
e
3
,
M
=
M
e
1
e
1
+
M
e
2
e
2
+
M
e
3
e
3
.The constitutive equation of the Euler-Bernoulli beam can be described as(32)
N
e
1
=
E
A
ε
,
M
e
1
=
G
J
κ
e
1
,
M
e
2
=
E
I
e
2
κ
e
2
,
M
e
3
=
E
I
e
3
κ
e
3
.By substituting (32) into (30), the internal virtual power can be expressed as(33)
δ
p
i
n
t
=
∫
0
L
δ
ε
˙
E
A
ε
+
δ
κ
˙
e
1
G
J
κ
e
1
+
δ
κ
˙
e
2
E
I
e
2
κ
e
2
+
δ
κ
˙
e
3
E
I
e
3
κ
e
3
d
s
,where E
A and G
J are the axial stiffness and the torsion stiffness, respectively. E
I
e
2 and E
I
e
3 are the principal bending stiffness.This above introduced principle is usually regarded as the geometrically exact beam theory, which is firstly proposed by Reissner [11] and then has been mainly developed by Simo and Vu-Quoc [4, 5].


## 3. Finite Element Implementation

### 3.1. Coupled Interpolation for Euler-Bernoulli Beam Element

A two-node straight beam element with initial lengthL and unchanged cross-section is considered as depicted in Figure 3.Figure 3
The kinematic description of the Euler-Bernoulli beam element.In the initial configuration,r
i
0 and R
θ
0 represent the nodal position vector and the orientation matrix of the cross-section, respectively. For simplicity, the cross-section basis is denoted by (
n
,
t
,
b
), which means that n
=
e
1, t
=
e
2, b
=
e
3. The global displacement and rotation vectors are chosen as the nodal coordinates:(34)
p
1
=
u
1
φ
1
,
p
2
=
u
2
φ
2
.And the element generalized coordinates p
e can be defined by(35)
p
e
=
p
1
p
2
.The nodal position vector and the orientation matrix of the cross-section in the current configuration can be expressed as(36)
r
i
=
r
i
0
+
u
i
,
R
i
=
n
i
t
i
b
i
=
R
φ
i
R
θ
0
.For the Euler-Bernoulli beam, the arc-length derivative of position vector remains parallel to the normal direction of the cross-section, that is, (28). At the both ends of the element, we have the relationship r
i
′
=
(
1
+
ε
i
)
n
i. It should be pointed out that the beam structure is mainly subjected to a transverse load and produces a bending deformation. The normal strain is a small quantity; that is, ε
≪
1. Therefore, the influence of the normal strain at two ends on the shape of the centerline can be ignored, and the global position vector of centerline r is obtained by employing the Hermite interpolation; that is,(37)
r
ξ
=
N
1
0
r
1
+
N
1
1
n
1
+
N
2
0
r
2
+
N
2
1
n
2
,where(38)
N
1
0
=
1
-
3
ξ
2
+
2
ξ
3
,
N
1
1
=
L
ξ
3
-
ξ
2
,
N
2
0
=
3
ξ
2
-
2
ξ
3
,
N
2
1
=
L
ξ
-
2
ξ
2
+
ξ
3are the Hermite shape functions with ξ
=
s
/
L
,
ξ
∈
[
0,1
].The unit normal vector of the cross-section is given by(39)
n
=
r
′
r
′
.The other two orthogonal unit vectors of the material basis t and b can be determined by the successive application of the following two rotation vectors to the left end section basis (
n
1
,
t
1
,
b
1
).Firstly, an intermediate reference frame(
n
,
t
_
,
b
_
) arrived when rotating (
n
1
,
t
1
,
b
1
) an angle β with respect to the direction n
1
×
n, which is illustrated in Figure 4. The corresponding rotation vector is expressed as(40)
β
=
λ
β
n
1
×
n
,where(41)
λ
β
=
β
sin
⁡
β
,
β
=
arccos
⁡
n
1
T
n
.Then, the final material basis (
n
,
t
,
b
) is achieved by rotating the (
n
,
t
_
,
b
_
) by an angle α with respect to the vector n. The corresponding rotation vector is(42)
α
=
α
n
.Thus the orientation matrix of the cross-section can be expressed by(43)
R
θ
=
n
t
b
=
R
α
R
β
R
1
.Figure 4
The decomposition of beam section rotation.The above analysis shows that the first rotation vectorβ can be determined by the interpolation of the position vector r
(
s
). Hence, we just need to interpolate the second rotation vector α within the element. Moreover, the direction of rotation vector α is n, and the angle α can be interpolated linearly by(44)
α
=
1
-
ξ
α
1
+
ξ
α
2
=
ξ
α
2
,were α
1
=
0 and the relative torsion angle of the right end section α
2 can be obtained according to the following steps:(a)
The rotation vectorα
2 can be calculated according to R
α
2
=
R
2
R
1
T
R
β
2
T. The extraction of α
2 from R
α
2 can be performed by using the Spurrier [27] algorithm as given by Simo and Vu-Quoc [5].
(b)
Then, the relative torsion angle is given asα
2
=
n
2
T
α
2.It should be mentioned that (40) cannot provide a unique solution if the first rotation angle β is equal to π, which means n
=
-
n
1. In this case, the relative bending angle within one element will be greater than 180°, which is prohibited in this paper while constructing the element.


### 3.2. Average Strains over a Small Time Interval and Nodal Forces

According to (28), the normal strain within the element is given by(45)
ε
ξ
=
n
T
r
′
-
1
.It is necessary and interesting to investigate the influence of the normal strain at two ends on the element deformation. Figure 5 shows the normal strain distribution in which the element is subjected to pure axial tensile and has the extension Δ
L.Figure 5
The normal strain distribution.In this example, the real normal strain distribution is uniform. But the normal strain directly derived from the position interpolation (37) is nonuniform, as shown in Figure 5. The magnitude of the normal strain reaches its maximal value in the middle and tends to zero at two ends. The strain energy is 20% more than that of the uniform case, which means that the beam becomes stiffer. In order to improve this solution, we assume that the normal strain is uniform within the element and equal to the average normal strain, which can be written by(46)
ε
e
=
∫
0
1
ε
ξ
d
ξ
.For the axial tension case, the exact solution can be obtained by using only two Gauss integration points. And the accuracy of this assumption will be illustrated in Section 4.In order to facilitate the following derivation, the time derivative ofε
(
ξ
) is defined by (47)
ε
˙
ξ
=
n
T
r
˙
′
,in which the identity n
˙
T
n
=
0 is used in the simplification. According to the interpolation of the global position vector (37), the velocity of the centroid can be expressed as (48)
r
˙
=
N
1
0
u
˙
1
-
N
1
1
n
~
1
ω
1
+
N
2
0
u
˙
2
-
N
2
1
n
~
2
ω
2
=
N
r
q
˙
e
,where r
˙
i
=
u
˙
i denotes the velocity of the element node and ω
i is the angular velocity. The element generalized velocity is given by(49)
q
˙
e
=
u
˙
1
T
ω
1
T
u
˙
2
T
ω
2
T
T
,
N
r
=
N
ξ
D
e
,
N
ξ
=
N
1
0
I
N
1
1
I
N
2
0
I
N
2
1
I
,
D
e
=
diag
⁡
I
,
-
n
~
1
,
I
,
-
n
~
2
.Substituting (48) into (47) yields(50)
ε
˙
=
Ψ
ε
q
˙
e
,
Ψ
ε
=
n
T
N
r
′
.Therefore, the time derivate of the average normal strain and its virtual variation can be expressed by(51)
ε
˙
e
=
Ψ
e
q
˙
e
,
δ
ε
˙
e
=
Ψ
e
δ
q
˙
e
,where(52)
Ψ
e
=
∫
0
1
Ψ
ε
d
ξ
.Based on the definition of the curvature vector (18) and the decomposition of beam cross-section rotation (43), the curvature vector of the beam is given by(53)
κ
=
T
α
α
′
+
R
α
T
β
β
′
.Furthermore, the components of curvature κ with respect to the material basis are given by (54)
κ
=
κ
n
n
+
κ
t
t
+
κ
b
b
.The two bending curvature components κ
t and κ
b, in terms of (11), are written as(55)
κ
t
=
-
b
T
n
′
,
κ
b
=
t
T
n
′
.The torsion component κ
n is(56)
κ
n
=
n
T
κ
=
n
T
T
α
α
′
+
n
T
R
α
T
β
β
′
.By substituting (40) and (42) into the two items on the right-hand side of (56), respectively, we obtain(57)
n
T
T
α
α
′
=
α
′
=
α
2
L
,
n
T
R
α
T
β
β
′
=
n
T
T
β
β
′
=
c
n
T
n
′
,where(58)
c
n
=
γ
n
n
×
n
1
,
γ
n
=
1
1
+
n
1
T
n
,
(59)
n
′
=
Γ
n
r
′
′
,
Γ
n
=
I
-
n
n
T
r
′
.After substituting (57) into (56), we have(60)
κ
n
=
α
2
L
+
c
n
T
n
′
.The time derivatives of (55) and (60) can be expressed by(61)
κ
˙
b
=
t
T
n
˙
′
+
n
′
T
t
˙
,
κ
˙
t
=
-
b
T
n
˙
′
-
n
′
T
b
˙
,
κ
˙
n
=
α
˙
2
L
+
c
n
T
n
˙
′
+
n
′
T
c
˙
n
,in which the time derivative of the unit normal vector n is given by(62)
n
˙
=
Γ
n
r
˙
′
=
S
n
q
˙
e
,
S
n
=
Γ
n
N
r
′
.Therefore, the derivative of (62) with respect to arc-length coordinate s is(63)
n
˙
′
=
S
n
′
q
˙
e
.According to the decomposition of the cross-section rotation (43), the relationship of the angular velocities at two ends of the element can be expressed as(64)
ω
2
=
T
α
2
α
˙
2
+
R
α
2
T
β
2
β
˙
2
+
R
α
2
R
β
2
ω
1
.So the time derivative of α
2 is(65)
α
˙
2
=
n
2
T
T
α
2
α
˙
2
=
n
2
T
+
c
n
2
T
n
~
2
ω
2
-
n
1
T
ω
1
=
S
α
2
q
˙
e
,where(66)
S
α
2
=
n
2
T
+
c
n
2
T
n
~
2
S
ω
2
-
n
1
T
S
ω
1
,and S
ω
i is constant matrix with which the angular velocity can be denoted by ω
i
=
S
ω
i
q
˙
e.The time derivative of (58) is(67)
c
˙
n
=
γ
˙
n
n
~
n
1
+
γ
n
n
~
n
˙
1
-
γ
n
n
~
1
n
˙
=
C
n
q
˙
e
,where(68)
C
n
=
n
~
Λ
n
n
~
1
S
ω
1
+
n
~
1
Λ
n
T
S
n
,
Λ
n
=
γ
n
2
n
1
n
T
-
γ
n
I
.The angular velocity ω can be expressed with respect to the material basis(69)
ω
=
ω
n
n
+
ω
t
t
+
ω
b
b
.Based on the relationship between the derivatives with respect to time t and arch-length coordinate s, we have(70)
ω
b
=
t
T
n
˙
=
W
b
q
˙
e
,
ω
t
=
-
b
T
n
˙
=
W
t
q
˙
e
,
ω
n
=
ξ
α
˙
2
+
c
n
T
n
˙
+
n
1
T
ω
1
=
W
n
q
˙
e
,where(71)
W
b
=
t
T
S
n
,
W
t
=
-
b
T
S
n
,
W
n
=
ξ
S
α
2
+
c
n
T
S
n
+
n
1
T
S
ω
1
.Then the time derivative of the other base vectors t and b can be given by(72)
t
˙
=
ω
×
t
=
ω
n
b
-
ω
b
n
=
S
t
q
˙
e
,
S
t
=
b
W
n
-
n
W
b
,
b
˙
=
ω
×
b
=
ω
t
n
-
ω
n
t
=
S
b
q
˙
e
,
S
b
=
n
W
t
-
t
W
n
.By substituting (63), (65), (67), and (72) into (61), the time derivatives of the curvature and their virtual variations can be expressed as(73)
κ
˙
b
=
Ψ
b
q
˙
e
,
δ
κ
˙
b
=
Ψ
b
δ
q
˙
e
,
κ
˙
t
=
Ψ
t
q
˙
e
,
δ
κ
˙
t
=
Ψ
t
δ
q
˙
e
,
κ
˙
n
=
Ψ
n
q
˙
e
,
δ
κ
˙
b
=
Ψ
n
δ
q
˙
e
,where(74)
Ψ
b
=
t
T
S
n
′
+
n
′
T
S
t
,
Ψ
t
=
-
b
T
S
n
′
+
n
′
T
S
b
,
Ψ
n
=
1
L
S
α
2
+
c
n
T
S
n
′
+
n
′
T
C
n
.By substituting the normal strain ε
e, torsion κ
n, two bending curvatures κ
t, κ
b, and their virtual time derivatives (51) and (73) into the internal virtual power equation (33), we obtain(75)
δ
p
i
n
t
=
δ
ε
˙
e
E
A
L
ε
e
+
L
∫
0
1
δ
κ
˙
n
G
J
κ
n
+
δ
κ
˙
t
E
I
t
κ
t
+
δ
κ
˙
b
E
I
b
κ
b
d
ξ
=
δ
q
˙
e
T
F
i
n
t
,where F
i
n
t is called the generalized nodal force:(76)
F
i
n
t
=
Ψ
e
T
E
A
L
ε
e
+
L
∫
0
1
Ψ
n
T
G
J
κ
n
+
Ψ
t
T
E
I
t
κ
t
+
Ψ
b
T
E
I
b
κ
b
d
ξ
.The time derivative of (76) can be written as(77)
F
˙
i
n
t
=
Ψ
e
T
E
A
L
ε
˙
e
+
Ψ
˙
e
T
E
A
L
ε
e
+
L
∫
0
1
Ψ
n
G
J
κ
˙
n
+
Ψ
t
E
I
t
κ
˙
t
+
Ψ
b
E
I
b
κ
˙
b
d
ξ
+
L
∫
0
1
Ψ
˙
n
G
J
κ
n
+
Ψ
˙
t
E
I
t
κ
t
+
Ψ
˙
b
E
I
b
κ
b
d
ξ
=
K
t
1
+
K
t
σ
q
˙
e
,where(78)
Ψ
˙
e
=
Ω
e
T
q
˙
e
,
Ψ
˙
n
=
Ω
n
T
q
˙
e
,
Ψ
˙
t
=
Ω
t
T
q
˙
e
,
Ψ
˙
b
=
Ω
b
T
q
˙
e
,
K
t
1
=
E
A
L
Ψ
e
T
Ψ
e
+
L
∫
0
1
G
J
Ψ
n
T
Ψ
n
+
E
I
t
Ψ
t
T
Ψ
t
+
E
I
b
Ψ
b
T
Ψ
b
d
ξ
,
K
t
σ
=
Ω
e
T
E
A
L
ε
e
+
L
∫
0
1
Ω
n
T
G
J
κ
n
+
Ω
t
T
E
I
t
κ
t
+
Ω
b
T
E
I
b
κ
b
d
ξ
.The first part K
t
1 is generally called the material stiffness matrix, and the second part K
t
σ is the “geometric” stiffness matrix.The beams studied in the rigid-flexible coupled dynamic systems undergo a large rigid body motion and small elastic vibration. Due to the influence of high frequencies, the time integration of the system equations becomes very difficult. In order to eliminate their effects, the average strains over a small time interval are proposed to replace the traditional transient strains in the expression of the nodal forces (76). Firstly, a second-order Taylor series expansion is employed to represent the transient strain. Taking the normal strain as an example, it can be given by(79)
ε
e
t
+
h
≈
ε
e
t
+
ε
˙
e
t
h
+
1
2
ε
¨
e
t
h
2
.Then the average normal strain over a small time interval h is defined as(80)
ε
⌢
e
=
1
h
∫
0
h
ε
e
t
+
h
d
h
=
ε
e
t
+
a
1
ε
˙
e
t
+
a
2
ε
¨
e
t
,where the coefficients are expressed as a
1
=
h
/
2 and a
2
=
h
2
/
6. Thus the average rotational strains are given by(81)
κ
⌢
n
=
κ
n
+
a
1
κ
˙
n
+
a
2
κ
¨
n
,
κ
⌢
b
=
κ
b
+
a
1
κ
˙
b
+
a
2
κ
¨
b
,
κ
⌢
t
=
κ
t
+
a
1
κ
˙
t
+
a
2
κ
¨
t
.By replacing the transient strains in (76) with average strains (80) and (81), the modified generalized nodal force can be written as(82)
F
⌢
i
n
t
=
Ψ
e
T
E
A
L
ε
⌢
e
+
L
∫
0
1
Ψ
n
T
G
J
κ
⌢
n
+
Ψ
t
T
E
I
t
κ
⌢
t
+
Ψ
b
T
E
I
b
κ
⌢
b
d
ξ
=
F
i
n
t
+
a
1
K
t
1
+
a
2
K
˙
C
q
˙
e
+
a
2
K
t
1
q
¨
e
,where(83)
K
˙
C
=
E
A
L
Ψ
e
T
Ψ
˙
e
+
L
∫
0
1
G
J
Ψ
n
T
Ψ
˙
n
+
E
I
t
Ψ
t
T
Ψ
˙
t
+
E
I
b
Ψ
b
T
Ψ
˙
b
d
ξ
.


### 3.3. Virtual Power of Inertial Forces and Mass Matrix

In terms of the centroid velocity (48), the virtual velocity of the centroid can be expressed as(84)
δ
r
˙
s
=
N
r
δ
q
˙
e
.The accelerated velocity of the centroid can be obtained by the time derivation of (48):(85)
r
¨
s
=
N
r
q
¨
e
+
a
r
,
a
r
=
N
˙
r
q
˙
e
.Therefore, the virtual kinetic power caused by the translation can be expressed as(86)
δ
p
t
r
a
=
m
∫
0
1
δ
r
˙
T
r
¨
d
ξ
=
δ
q
˙
e
T
M
t
r
a
q
¨
e
+
F
t
r
a
,where(87)
M
t
r
a
=
m
∫
0
1
N
r
T
N
r
d
ξ
=
m
D
e
T
H
N
D
e
,
F
t
r
a
=
m
∫
0
1
N
r
T
a
r
d
ξ
=
m
D
e
T
H
N
D
˙
e
q
˙
e
,
H
N
=
∫
0
1
N
ξ
T
N
ξ
d
ξ
=
1
420
156
I
22
L
I
54
I
-
13
L
I
4
L
2
I
13
L
I
-
3
L
2
I
sym.
156
I
-
22
L
I
4
L
2
I
.The angular velocity of the cross-section is(88)
ω
=
R
θ
ω
-
,
ω
~
=
R
θ
ω
-
~
R
θ
T
.According to (70), the angular velocity in the material form and its virtual variation can be expressed as(89)
ω
-
=
N
ω
q
˙
e
,
δ
ω
-
=
N
ω
δ
q
˙
e
,where(90)
N
ω
=
W
n
T
W
t
T
W
b
T
T
.The angular acceleration is(91)
ω
˙
=
ω
˙
n
n
+
ω
˙
t
t
+
ω
˙
b
b
=
R
θ
ω
-
˙
,where(92)
ω
-
˙
=
N
ω
q
¨
e
+
a
ω
,
a
ω
=
N
˙
ω
q
˙
e
.The virtual kinetic power caused by the rotation may be expressed by(93)
δ
p
r
o
t
=
L
∫
0
1
δ
ω
T
J
ρ
ω
˙
+
ω
~
J
ρ
ω
d
ξ
=
L
∫
0
1
δ
ω
-
T
J
-
ρ
ω
-
˙
+
ω
-
~
J
-
ρ
ω
-
d
ξ
,where J
-
ρ
=
R
θ
T
J
ρ
R
θ is the centroidal inertia matrix of the cross-section with respect to the material basis. Substituting (89) and (92) into (93) yields(94)
δ
p
r
o
t
=
δ
q
˙
e
T
M
r
o
t
q
¨
e
+
F
r
o
t
,where(95)
M
r
o
t
=
L
∫
0
1
N
ω
T
J
-
ρ
N
ω
d
ξ
,
F
r
o
t
=
L
∫
0
1
N
ω
T
J
-
ρ
a
ω
+
ω
-
~
J
-
ρ
ω
-
d
ξ
.In terms of (86) and (94), the virtual power of inertial forces δ
p
i
n
e can be written as(96)
δ
p
i
n
e
=
δ
p
t
r
a
+
δ
p
r
o
t
=
δ
q
˙
e
T
M
i
n
e
q
¨
e
+
F
i
n
e
,where(97)
M
ine
=
M
tra
+
M
rot
⁡
,
F
ine
=
F
tra
+
F
rot
⁡
.


### 3.4. Governing Equation

The virtual power equation of the beam element can be written as(98)
δ
p
i
n
e
+
δ
p
i
n
t
=
δ
p
e
x
t
,in which δ
p
e
x
t
=
δ
q
˙
e
T
F
e
x
t is the external virtual power and F
e
x
t is the generalized external force. Substituting the internal virtual power (75) and the inertial virtual power (96) into the virtual power equation (98) yields(99)
δ
q
˙
e
T
M
ine
q
¨
e
+
F
ine
+
F
⌢
i
n
t
-
F
e
x
t
=
0
.By substituting the modified generalized nodal force (82) into (99), the governing equations of beam element can be expressed as(100)
M
ine
+
a
2
K
t
1
q
¨
e
+
a
1
K
t
1
+
a
2
K
˙
C
q
˙
e
+
F
ine
+
F
i
n
t
-
F
e
x
t
=
0
.Thus, the generalized acceleration is(101)
q
¨
e
=
M
⌢
-
1
F
⌢
,where the mass matrix and generalized force can be written by(102)
M
⌢
=
M
ine
+
a
2
K
t
1
,
F
⌢
=
F
e
x
t
-
F
ine
-
F
i
n
t
-
a
1
K
t
1
+
a
2
K
˙
C
q
˙
e
.The underlined items in (102) are the additional mass matrix and the additional generalized forces, which are introduced by the average strains in the expression of the nodal forces, and will vanish when the time parameter h is equal to zero. Generally, the system frequency will reduce when increasing the mass. Referring to the damping model in [28], it is found that the item a
1
K
t
1
q
˙
e is equivalent to the internal damping force. If q
˙
e represents the rigid motion only, this equivalent damping force will disappear since K
t
1
q
˙
e is equal to zero for the allowed rigid motion.The element generalized velocityq
˙
e and the generalized coordinates p
e are set as the state variables; that is,(103)
X
=
q
˙
e
p
e
,in which the relationship between the time derivative of p
e and q
˙
e becomes(104)
p
˙
e
=
Λ
q
˙
e
,
Λ
=
diag
⁡
I
,
T
φ
1
-
1
,
I
,
T
φ
2
-
1
.The governing equations, which are a set of the second order ordinary differential equations, can be written as(105)
X
˙
=
M
⌢
-
1
F
⌢
Λ
q
˙
e
.Equations (105) are the state equations, which are a set of the first order ordinary differential equations, and can be directly solved by the ordinary differential equation solvers.As mentioned in Section2, the inherent singularity problem of the rotation vector can be avoided according to the procedures as shown in Figure 6. After each successful integration step, the rotation vector of every node will be checked. Once the angle large is larger than π, the parameterization will be changed into its complement rotation vector (22). Then the next integration step can proceed.Figure 6
Flowchart for avoiding the singularity of the rotation vector.


## 3.1. Coupled Interpolation for Euler-Bernoulli Beam Element

A two-node straight beam element with initial lengthL and unchanged cross-section is considered as depicted in Figure 3.Figure 3
The kinematic description of the Euler-Bernoulli beam element.In the initial configuration,r
i
0 and R
θ
0 represent the nodal position vector and the orientation matrix of the cross-section, respectively. For simplicity, the cross-section basis is denoted by (
n
,
t
,
b
), which means that n
=
e
1, t
=
e
2, b
=
e
3. The global displacement and rotation vectors are chosen as the nodal coordinates:(34)
p
1
=
u
1
φ
1
,
p
2
=
u
2
φ
2
.And the element generalized coordinates p
e can be defined by(35)
p
e
=
p
1
p
2
.The nodal position vector and the orientation matrix of the cross-section in the current configuration can be expressed as(36)
r
i
=
r
i
0
+
u
i
,
R
i
=
n
i
t
i
b
i
=
R
φ
i
R
θ
0
.For the Euler-Bernoulli beam, the arc-length derivative of position vector remains parallel to the normal direction of the cross-section, that is, (28). At the both ends of the element, we have the relationship r
i
′
=
(
1
+
ε
i
)
n
i. It should be pointed out that the beam structure is mainly subjected to a transverse load and produces a bending deformation. The normal strain is a small quantity; that is, ε
≪
1. Therefore, the influence of the normal strain at two ends on the shape of the centerline can be ignored, and the global position vector of centerline r is obtained by employing the Hermite interpolation; that is,(37)
r
ξ
=
N
1
0
r
1
+
N
1
1
n
1
+
N
2
0
r
2
+
N
2
1
n
2
,where(38)
N
1
0
=
1
-
3
ξ
2
+
2
ξ
3
,
N
1
1
=
L
ξ
3
-
ξ
2
,
N
2
0
=
3
ξ
2
-
2
ξ
3
,
N
2
1
=
L
ξ
-
2
ξ
2
+
ξ
3are the Hermite shape functions with ξ
=
s
/
L
,
ξ
∈
[
0,1
].The unit normal vector of the cross-section is given by(39)
n
=
r
′
r
′
.The other two orthogonal unit vectors of the material basis t and b can be determined by the successive application of the following two rotation vectors to the left end section basis (
n
1
,
t
1
,
b
1
).Firstly, an intermediate reference frame(
n
,
t
_
,
b
_
) arrived when rotating (
n
1
,
t
1
,
b
1
) an angle β with respect to the direction n
1
×
n, which is illustrated in Figure 4. The corresponding rotation vector is expressed as(40)
β
=
λ
β
n
1
×
n
,where(41)
λ
β
=
β
sin
⁡
β
,
β
=
arccos
⁡
n
1
T
n
.Then, the final material basis (
n
,
t
,
b
) is achieved by rotating the (
n
,
t
_
,
b
_
) by an angle α with respect to the vector n. The corresponding rotation vector is(42)
α
=
α
n
.Thus the orientation matrix of the cross-section can be expressed by(43)
R
θ
=
n
t
b
=
R
α
R
β
R
1
.Figure 4
The decomposition of beam section rotation.The above analysis shows that the first rotation vectorβ can be determined by the interpolation of the position vector r
(
s
). Hence, we just need to interpolate the second rotation vector α within the element. Moreover, the direction of rotation vector α is n, and the angle α can be interpolated linearly by(44)
α
=
1
-
ξ
α
1
+
ξ
α
2
=
ξ
α
2
,were α
1
=
0 and the relative torsion angle of the right end section α
2 can be obtained according to the following steps:(a)
The rotation vectorα
2 can be calculated according to R
α
2
=
R
2
R
1
T
R
β
2
T. The extraction of α
2 from R
α
2 can be performed by using the Spurrier [27] algorithm as given by Simo and Vu-Quoc [5].
(b)
Then, the relative torsion angle is given asα
2
=
n
2
T
α
2.It should be mentioned that (40) cannot provide a unique solution if the first rotation angle β is equal to π, which means n
=
-
n
1. In this case, the relative bending angle within one element will be greater than 180°, which is prohibited in this paper while constructing the element.


## 3.2. Average Strains over a Small Time Interval and Nodal Forces

According to (28), the normal strain within the element is given by(45)
ε
ξ
=
n
T
r
′
-
1
.It is necessary and interesting to investigate the influence of the normal strain at two ends on the element deformation. Figure 5 shows the normal strain distribution in which the element is subjected to pure axial tensile and has the extension Δ
L.Figure 5
The normal strain distribution.In this example, the real normal strain distribution is uniform. But the normal strain directly derived from the position interpolation (37) is nonuniform, as shown in Figure 5. The magnitude of the normal strain reaches its maximal value in the middle and tends to zero at two ends. The strain energy is 20% more than that of the uniform case, which means that the beam becomes stiffer. In order to improve this solution, we assume that the normal strain is uniform within the element and equal to the average normal strain, which can be written by(46)
ε
e
=
∫
0
1
ε
ξ
d
ξ
.For the axial tension case, the exact solution can be obtained by using only two Gauss integration points. And the accuracy of this assumption will be illustrated in Section 4.In order to facilitate the following derivation, the time derivative ofε
(
ξ
) is defined by (47)
ε
˙
ξ
=
n
T
r
˙
′
,in which the identity n
˙
T
n
=
0 is used in the simplification. According to the interpolation of the global position vector (37), the velocity of the centroid can be expressed as (48)
r
˙
=
N
1
0
u
˙
1
-
N
1
1
n
~
1
ω
1
+
N
2
0
u
˙
2
-
N
2
1
n
~
2
ω
2
=
N
r
q
˙
e
,where r
˙
i
=
u
˙
i denotes the velocity of the element node and ω
i is the angular velocity. The element generalized velocity is given by(49)
q
˙
e
=
u
˙
1
T
ω
1
T
u
˙
2
T
ω
2
T
T
,
N
r
=
N
ξ
D
e
,
N
ξ
=
N
1
0
I
N
1
1
I
N
2
0
I
N
2
1
I
,
D
e
=
diag
⁡
I
,
-
n
~
1
,
I
,
-
n
~
2
.Substituting (48) into (47) yields(50)
ε
˙
=
Ψ
ε
q
˙
e
,
Ψ
ε
=
n
T
N
r
′
.Therefore, the time derivate of the average normal strain and its virtual variation can be expressed by(51)
ε
˙
e
=
Ψ
e
q
˙
e
,
δ
ε
˙
e
=
Ψ
e
δ
q
˙
e
,where(52)
Ψ
e
=
∫
0
1
Ψ
ε
d
ξ
.Based on the definition of the curvature vector (18) and the decomposition of beam cross-section rotation (43), the curvature vector of the beam is given by(53)
κ
=
T
α
α
′
+
R
α
T
β
β
′
.Furthermore, the components of curvature κ with respect to the material basis are given by (54)
κ
=
κ
n
n
+
κ
t
t
+
κ
b
b
.The two bending curvature components κ
t and κ
b, in terms of (11), are written as(55)
κ
t
=
-
b
T
n
′
,
κ
b
=
t
T
n
′
.The torsion component κ
n is(56)
κ
n
=
n
T
κ
=
n
T
T
α
α
′
+
n
T
R
α
T
β
β
′
.By substituting (40) and (42) into the two items on the right-hand side of (56), respectively, we obtain(57)
n
T
T
α
α
′
=
α
′
=
α
2
L
,
n
T
R
α
T
β
β
′
=
n
T
T
β
β
′
=
c
n
T
n
′
,where(58)
c
n
=
γ
n
n
×
n
1
,
γ
n
=
1
1
+
n
1
T
n
,
(59)
n
′
=
Γ
n
r
′
′
,
Γ
n
=
I
-
n
n
T
r
′
.After substituting (57) into (56), we have(60)
κ
n
=
α
2
L
+
c
n
T
n
′
.The time derivatives of (55) and (60) can be expressed by(61)
κ
˙
b
=
t
T
n
˙
′
+
n
′
T
t
˙
,
κ
˙
t
=
-
b
T
n
˙
′
-
n
′
T
b
˙
,
κ
˙
n
=
α
˙
2
L
+
c
n
T
n
˙
′
+
n
′
T
c
˙
n
,in which the time derivative of the unit normal vector n is given by(62)
n
˙
=
Γ
n
r
˙
′
=
S
n
q
˙
e
,
S
n
=
Γ
n
N
r
′
.Therefore, the derivative of (62) with respect to arc-length coordinate s is(63)
n
˙
′
=
S
n
′
q
˙
e
.According to the decomposition of the cross-section rotation (43), the relationship of the angular velocities at two ends of the element can be expressed as(64)
ω
2
=
T
α
2
α
˙
2
+
R
α
2
T
β
2
β
˙
2
+
R
α
2
R
β
2
ω
1
.So the time derivative of α
2 is(65)
α
˙
2
=
n
2
T
T
α
2
α
˙
2
=
n
2
T
+
c
n
2
T
n
~
2
ω
2
-
n
1
T
ω
1
=
S
α
2
q
˙
e
,where(66)
S
α
2
=
n
2
T
+
c
n
2
T
n
~
2
S
ω
2
-
n
1
T
S
ω
1
,and S
ω
i is constant matrix with which the angular velocity can be denoted by ω
i
=
S
ω
i
q
˙
e.The time derivative of (58) is(67)
c
˙
n
=
γ
˙
n
n
~
n
1
+
γ
n
n
~
n
˙
1
-
γ
n
n
~
1
n
˙
=
C
n
q
˙
e
,where(68)
C
n
=
n
~
Λ
n
n
~
1
S
ω
1
+
n
~
1
Λ
n
T
S
n
,
Λ
n
=
γ
n
2
n
1
n
T
-
γ
n
I
.The angular velocity ω can be expressed with respect to the material basis(69)
ω
=
ω
n
n
+
ω
t
t
+
ω
b
b
.Based on the relationship between the derivatives with respect to time t and arch-length coordinate s, we have(70)
ω
b
=
t
T
n
˙
=
W
b
q
˙
e
,
ω
t
=
-
b
T
n
˙
=
W
t
q
˙
e
,
ω
n
=
ξ
α
˙
2
+
c
n
T
n
˙
+
n
1
T
ω
1
=
W
n
q
˙
e
,where(71)
W
b
=
t
T
S
n
,
W
t
=
-
b
T
S
n
,
W
n
=
ξ
S
α
2
+
c
n
T
S
n
+
n
1
T
S
ω
1
.Then the time derivative of the other base vectors t and b can be given by(72)
t
˙
=
ω
×
t
=
ω
n
b
-
ω
b
n
=
S
t
q
˙
e
,
S
t
=
b
W
n
-
n
W
b
,
b
˙
=
ω
×
b
=
ω
t
n
-
ω
n
t
=
S
b
q
˙
e
,
S
b
=
n
W
t
-
t
W
n
.By substituting (63), (65), (67), and (72) into (61), the time derivatives of the curvature and their virtual variations can be expressed as(73)
κ
˙
b
=
Ψ
b
q
˙
e
,
δ
κ
˙
b
=
Ψ
b
δ
q
˙
e
,
κ
˙
t
=
Ψ
t
q
˙
e
,
δ
κ
˙
t
=
Ψ
t
δ
q
˙
e
,
κ
˙
n
=
Ψ
n
q
˙
e
,
δ
κ
˙
b
=
Ψ
n
δ
q
˙
e
,where(74)
Ψ
b
=
t
T
S
n
′
+
n
′
T
S
t
,
Ψ
t
=
-
b
T
S
n
′
+
n
′
T
S
b
,
Ψ
n
=
1
L
S
α
2
+
c
n
T
S
n
′
+
n
′
T
C
n
.By substituting the normal strain ε
e, torsion κ
n, two bending curvatures κ
t, κ
b, and their virtual time derivatives (51) and (73) into the internal virtual power equation (33), we obtain(75)
δ
p
i
n
t
=
δ
ε
˙
e
E
A
L
ε
e
+
L
∫
0
1
δ
κ
˙
n
G
J
κ
n
+
δ
κ
˙
t
E
I
t
κ
t
+
δ
κ
˙
b
E
I
b
κ
b
d
ξ
=
δ
q
˙
e
T
F
i
n
t
,where F
i
n
t is called the generalized nodal force:(76)
F
i
n
t
=
Ψ
e
T
E
A
L
ε
e
+
L
∫
0
1
Ψ
n
T
G
J
κ
n
+
Ψ
t
T
E
I
t
κ
t
+
Ψ
b
T
E
I
b
κ
b
d
ξ
.The time derivative of (76) can be written as(77)
F
˙
i
n
t
=
Ψ
e
T
E
A
L
ε
˙
e
+
Ψ
˙
e
T
E
A
L
ε
e
+
L
∫
0
1
Ψ
n
G
J
κ
˙
n
+
Ψ
t
E
I
t
κ
˙
t
+
Ψ
b
E
I
b
κ
˙
b
d
ξ
+
L
∫
0
1
Ψ
˙
n
G
J
κ
n
+
Ψ
˙
t
E
I
t
κ
t
+
Ψ
˙
b
E
I
b
κ
b
d
ξ
=
K
t
1
+
K
t
σ
q
˙
e
,where(78)
Ψ
˙
e
=
Ω
e
T
q
˙
e
,
Ψ
˙
n
=
Ω
n
T
q
˙
e
,
Ψ
˙
t
=
Ω
t
T
q
˙
e
,
Ψ
˙
b
=
Ω
b
T
q
˙
e
,
K
t
1
=
E
A
L
Ψ
e
T
Ψ
e
+
L
∫
0
1
G
J
Ψ
n
T
Ψ
n
+
E
I
t
Ψ
t
T
Ψ
t
+
E
I
b
Ψ
b
T
Ψ
b
d
ξ
,
K
t
σ
=
Ω
e
T
E
A
L
ε
e
+
L
∫
0
1
Ω
n
T
G
J
κ
n
+
Ω
t
T
E
I
t
κ
t
+
Ω
b
T
E
I
b
κ
b
d
ξ
.The first part K
t
1 is generally called the material stiffness matrix, and the second part K
t
σ is the “geometric” stiffness matrix.The beams studied in the rigid-flexible coupled dynamic systems undergo a large rigid body motion and small elastic vibration. Due to the influence of high frequencies, the time integration of the system equations becomes very difficult. In order to eliminate their effects, the average strains over a small time interval are proposed to replace the traditional transient strains in the expression of the nodal forces (76). Firstly, a second-order Taylor series expansion is employed to represent the transient strain. Taking the normal strain as an example, it can be given by(79)
ε
e
t
+
h
≈
ε
e
t
+
ε
˙
e
t
h
+
1
2
ε
¨
e
t
h
2
.Then the average normal strain over a small time interval h is defined as(80)
ε
⌢
e
=
1
h
∫
0
h
ε
e
t
+
h
d
h
=
ε
e
t
+
a
1
ε
˙
e
t
+
a
2
ε
¨
e
t
,where the coefficients are expressed as a
1
=
h
/
2 and a
2
=
h
2
/
6. Thus the average rotational strains are given by(81)
κ
⌢
n
=
κ
n
+
a
1
κ
˙
n
+
a
2
κ
¨
n
,
κ
⌢
b
=
κ
b
+
a
1
κ
˙
b
+
a
2
κ
¨
b
,
κ
⌢
t
=
κ
t
+
a
1
κ
˙
t
+
a
2
κ
¨
t
.By replacing the transient strains in (76) with average strains (80) and (81), the modified generalized nodal force can be written as(82)
F
⌢
i
n
t
=
Ψ
e
T
E
A
L
ε
⌢
e
+
L
∫
0
1
Ψ
n
T
G
J
κ
⌢
n
+
Ψ
t
T
E
I
t
κ
⌢
t
+
Ψ
b
T
E
I
b
κ
⌢
b
d
ξ
=
F
i
n
t
+
a
1
K
t
1
+
a
2
K
˙
C
q
˙
e
+
a
2
K
t
1
q
¨
e
,where(83)
K
˙
C
=
E
A
L
Ψ
e
T
Ψ
˙
e
+
L
∫
0
1
G
J
Ψ
n
T
Ψ
˙
n
+
E
I
t
Ψ
t
T
Ψ
˙
t
+
E
I
b
Ψ
b
T
Ψ
˙
b
d
ξ
.


## 3.3. Virtual Power of Inertial Forces and Mass Matrix

In terms of the centroid velocity (48), the virtual velocity of the centroid can be expressed as(84)
δ
r
˙
s
=
N
r
δ
q
˙
e
.The accelerated velocity of the centroid can be obtained by the time derivation of (48):(85)
r
¨
s
=
N
r
q
¨
e
+
a
r
,
a
r
=
N
˙
r
q
˙
e
.Therefore, the virtual kinetic power caused by the translation can be expressed as(86)
δ
p
t
r
a
=
m
∫
0
1
δ
r
˙
T
r
¨
d
ξ
=
δ
q
˙
e
T
M
t
r
a
q
¨
e
+
F
t
r
a
,where(87)
M
t
r
a
=
m
∫
0
1
N
r
T
N
r
d
ξ
=
m
D
e
T
H
N
D
e
,
F
t
r
a
=
m
∫
0
1
N
r
T
a
r
d
ξ
=
m
D
e
T
H
N
D
˙
e
q
˙
e
,
H
N
=
∫
0
1
N
ξ
T
N
ξ
d
ξ
=
1
420
156
I
22
L
I
54
I
-
13
L
I
4
L
2
I
13
L
I
-
3
L
2
I
sym.
156
I
-
22
L
I
4
L
2
I
.The angular velocity of the cross-section is(88)
ω
=
R
θ
ω
-
,
ω
~
=
R
θ
ω
-
~
R
θ
T
.According to (70), the angular velocity in the material form and its virtual variation can be expressed as(89)
ω
-
=
N
ω
q
˙
e
,
δ
ω
-
=
N
ω
δ
q
˙
e
,where(90)
N
ω
=
W
n
T
W
t
T
W
b
T
T
.The angular acceleration is(91)
ω
˙
=
ω
˙
n
n
+
ω
˙
t
t
+
ω
˙
b
b
=
R
θ
ω
-
˙
,where(92)
ω
-
˙
=
N
ω
q
¨
e
+
a
ω
,
a
ω
=
N
˙
ω
q
˙
e
.The virtual kinetic power caused by the rotation may be expressed by(93)
δ
p
r
o
t
=
L
∫
0
1
δ
ω
T
J
ρ
ω
˙
+
ω
~
J
ρ
ω
d
ξ
=
L
∫
0
1
δ
ω
-
T
J
-
ρ
ω
-
˙
+
ω
-
~
J
-
ρ
ω
-
d
ξ
,where J
-
ρ
=
R
θ
T
J
ρ
R
θ is the centroidal inertia matrix of the cross-section with respect to the material basis. Substituting (89) and (92) into (93) yields(94)
δ
p
r
o
t
=
δ
q
˙
e
T
M
r
o
t
q
¨
e
+
F
r
o
t
,where(95)
M
r
o
t
=
L
∫
0
1
N
ω
T
J
-
ρ
N
ω
d
ξ
,
F
r
o
t
=
L
∫
0
1
N
ω
T
J
-
ρ
a
ω
+
ω
-
~
J
-
ρ
ω
-
d
ξ
.In terms of (86) and (94), the virtual power of inertial forces δ
p
i
n
e can be written as(96)
δ
p
i
n
e
=
δ
p
t
r
a
+
δ
p
r
o
t
=
δ
q
˙
e
T
M
i
n
e
q
¨
e
+
F
i
n
e
,where(97)
M
ine
=
M
tra
+
M
rot
⁡
,
F
ine
=
F
tra
+
F
rot
⁡
.


## 3.4. Governing Equation

The virtual power equation of the beam element can be written as(98)
δ
p
i
n
e
+
δ
p
i
n
t
=
δ
p
e
x
t
,in which δ
p
e
x
t
=
δ
q
˙
e
T
F
e
x
t is the external virtual power and F
e
x
t is the generalized external force. Substituting the internal virtual power (75) and the inertial virtual power (96) into the virtual power equation (98) yields(99)
δ
q
˙
e
T
M
ine
q
¨
e
+
F
ine
+
F
⌢
i
n
t
-
F
e
x
t
=
0
.By substituting the modified generalized nodal force (82) into (99), the governing equations of beam element can be expressed as(100)
M
ine
+
a
2
K
t
1
q
¨
e
+
a
1
K
t
1
+
a
2
K
˙
C
q
˙
e
+
F
ine
+
F
i
n
t
-
F
e
x
t
=
0
.Thus, the generalized acceleration is(101)
q
¨
e
=
M
⌢
-
1
F
⌢
,where the mass matrix and generalized force can be written by(102)
M
⌢
=
M
ine
+
a
2
K
t
1
,
F
⌢
=
F
e
x
t
-
F
ine
-
F
i
n
t
-
a
1
K
t
1
+
a
2
K
˙
C
q
˙
e
.The underlined items in (102) are the additional mass matrix and the additional generalized forces, which are introduced by the average strains in the expression of the nodal forces, and will vanish when the time parameter h is equal to zero. Generally, the system frequency will reduce when increasing the mass. Referring to the damping model in [28], it is found that the item a
1
K
t
1
q
˙
e is equivalent to the internal damping force. If q
˙
e represents the rigid motion only, this equivalent damping force will disappear since K
t
1
q
˙
e is equal to zero for the allowed rigid motion.The element generalized velocityq
˙
e and the generalized coordinates p
e are set as the state variables; that is,(103)
X
=
q
˙
e
p
e
,in which the relationship between the time derivative of p
e and q
˙
e becomes(104)
p
˙
e
=
Λ
q
˙
e
,
Λ
=
diag
⁡
I
,
T
φ
1
-
1
,
I
,
T
φ
2
-
1
.The governing equations, which are a set of the second order ordinary differential equations, can be written as(105)
X
˙
=
M
⌢
-
1
F
⌢
Λ
q
˙
e
.Equations (105) are the state equations, which are a set of the first order ordinary differential equations, and can be directly solved by the ordinary differential equation solvers.As mentioned in Section2, the inherent singularity problem of the rotation vector can be avoided according to the procedures as shown in Figure 6. After each successful integration step, the rotation vector of every node will be checked. Once the angle large is larger than π, the parameterization will be changed into its complement rotation vector (22). Then the next integration step can proceed.Figure 6
Flowchart for avoiding the singularity of the rotation vector.


## 4. Numerical Simulations

In this section, several numerical examples are carried out to evaluate the static and dynamic behavior of the proposed spatial Euler-Bernoulli beam element. All simulations are performed on the same PC with an Intel Core 3.0 GHz processor and 8 GB RAM. The static solution is calculated by using fsolve function in the MATLAB package. The governing dynamic equations (105) are solved by using an explicit Runge-Kutta method ode45 implemented in MATLAB and an implicit Runge-Kutta method radau5 [29] separately. In addition, the results obtained by our proposed method are compared with the solutions calculated by the existing alternative ones.


### 4.1. Cantilever Beam with an End Moment

The first problem to be considered is a cantilever beam with material properties shown in Figure7, where a concentrated end moment is imposed at the end of the beam.Figure 7
Cantilever beam with an end moment.This problem has been used by many researchers for testing nonlinear beam elements since the analytical solutions for the problem exist. The straight beam will be bent to a full circle with the momentM
=
2
π
E
I
/
L. In this problem, 10 elements are used to model the beam. The results obtained from this study are compared with the analytical solutions as shown in Figure 8. And the deformations of the cantilever beam with different bending moments are shown in Figure 9. It can be seen that the results obtained by using our proposed elements are exactly the same with the analytical solutions.Figure 8
Comparison of the displacements of the cantilever beam with an end moment.Figure 9
Large deformations of the cantilever beam with an end moment.


### 4.2. Cantilever Beam with End Point Load

The cantilever beam subjected to a concentrated end load, as shown in Figure10, is investigated in this subsection.Figure 10
Cantilever beam with an end point load.In this problem, 10 elements are used to model the beam. The obtained results are compared with the elliptic integral solutions provided by Mattiasson [30] in Figure 11 and in numeric form listed in Table 1. The comparison indicates our results are good agreement with Mattiasson’s.Table 1
Cantilever beam with an end point load.
F
L
2
/
E
I
u
/
L
v
/
L
Mattiasson
Present
Mattiasson
Present
1.0
0.05643
0.05643
0.30172
0.30172
2.0
0.16064
0.16064
0.49346
0.49347
3.0
0.25442
0.25442
0.60325
0.60327
4.0
0.32894
0.32894
0.66996
0.66999
5.0
0.38763
0.38762
0.71379
0.71383
6.0
0.43459
0.43458
0.74457
0.74462
7.0
0.47293
0.47292
0.76737
0.76743
8.0
0.50483
0.50482
0.78498
0.78505
9.0
0.53182
0.53182
0.79906
0.79913
10.0
0.55500
0.55499
0.81061
0.81070Figure 11
Displacement curve of cantilever beam with an end point load.


### 4.3. Natural Frequency of the Cantilever Beam

The cantilever beam with the material properties as shown in Figure12 is carried out to study the natural frequency of the cantilever beam by employing the proposed beam element.Figure 12
Cantilever beam.To obtain the natural frequency, the nonlinear dynamic equation (100) should be linearized at the equilibrium positions. Firstly, without considering average strains in the expression of the nodal forces which means h
=
0, the frequency equation can be given by(106)
K
t
1
-
ω
2
M
i
n
e
=
0
,where ω is the natural frequency of the cantilever beam. In this example, the beam is modeled by N elements. Then the first bending, torsion, and tension frequencies are obtained separately as shown in Table 2. It is found that the frequencies will converge very fast to the theoretical values.Table 2
Frequency of cantilever beam with different meshes.
Modes
Frequency (Hz)
1-elems
2-elems
5-elems
10-elems
20-elems
Theory
Bending
5.291
5.269
5.266
5.266
5.266
5.266
Torsion
557.328
518.511
504.440
505.961
505.571
505.441
Tension
851.382
825.971
816.855
815.500
815.159
815.046Then, the influence of introducing the average strains on the frequencies is also studied. With a small time interval parameterh, the frequency equation becomes(107)
K
t
1
-
ω
′
2
M
i
n
e
+
a
2
K
t
1
=
0
,where ω
′ is the modified frequency of the system. Equation (107) can be simplified to the equivalent form (108)
K
t
1
-
ω
′
2
1
-
a
2
ω
′
2
M
i
n
e
=
0
.By comparing (106) with (108), the relationship between ω
′ and ω can be given by(109)
ω
′
=
ω
1
+
a
2
ω
2
.Substituting a
2
=
h
2
/
6 into (109) yields(110)
ω
′
=
ω
1
+
1
/
6
ω
/
f
h
2
<
6
f
h
,where f
h
=
1
/
h is the adjustable frequency parameter associated with the given small time interval h. According to (110), the mapping between the normalized modified frequency ω
′
/
f
h and ω
/
f
h is shown in Figure 13.Figure 13
General shape of the mapping curve betweenω
′
/
f
h and ω
/
f
h.The influence of the time parameterh on the system frequencies can be summarized as follows:(a)
All the frequencies of the system decrease to varying degrees.
(b)
The frequencies with the magnitude less thanf
h decrease very slightly.
(c)
The maximum limit of the modified frequencies is6
f
h.Table3 shows the modified frequency of cantilever beam modeled by 20 elements with different time parameters.Table 3
Modified frequency of cantilever beam with different time parameterh.
Modes
Modified frequency (Hz)
h
=
0
h
=
0.0001
h
=
0.001
h
=
0.01
h
=
0.05
h
=
0.1
Bending
5.266
5.266
5.266
5.265
5.236
5.148
Torsion
505.571
505.463
495.135
220.439
48.761
24.466
Tension
815.159
814.708
773.454
234.587
48.902
24.484Obviously, from the above discussion it is found that the high frequency of the system can be filtered out by choosing an appropriate time parameterh. This property is very useful for the time integration in the stiff problem, which is proved in the following examples.


### 4.4. Right-Angle Cantilever Beam

An L-shaped cantilever beam [6] is subjected to a concentrated load applied in the direction of global z-axis at the elbow as shown in Figure 14.Figure 14
Right-angle cantilever beam.The magnitude of this load follows the pattern of a hat function, as shown in Figure14. The material properties are(111)
E
A
=
G
A
=
1
0
6
,
E
I
=
G
J
=
1
0
3
,
L
1
=
L
2
=
10
,
A
ρ
=
1
,
J
-
ρ
=
diag
⁡
20,10,10
.The beam is meshed into 10 elements to study the displacements at both the elbow and the tip. The computations are carried out by using the radau5 and the ode45 with different values of the time parameterh. The total simulation time is 30 seconds. The relative error tolerance is set to 10−4 and the absolute error tolerance is set to 10−6. Table 4 shows the total CPU time while using the different methods. And the computed response of the elbow and the tip displacements are given in Figures 15 and 16, respectively. These results correspond to the results given by Ibrahimbegović and Mikdad [31].Table 4
The comparison of CPU time for the right-angle cantilever beam.
Time parameterh
Integration algorithms
CPU time (s)
0
radau5
6279
0
ode45
—
0.01
ode45
3348
0.1
ode45
426Figure 15
Time evolution of the displacement at the elbow.Figure 16
Time evolution of the displacement at the tip.It should be pointed out that the radau5 is a stiff ODE solver, and the ode45 is only workable for solving nonstiff differential equations. With the time parameterh
=
0, the radau5 method needs 6279 seconds to accomplish the simulation while ode45 cannot finish the task since the governing equations of this example are stiff equations. When h is set to 0.01, it can be found that the ode45 solver can complete the calculation without any difficulties and spends only 3348 seconds, which is about half of that of the stiff ODE solver radau5. For h
=
0.1, the CPU time of the ode45 method is reduced to 426 seconds. From Figures 15 and 16, it can be also found that both the displacements at the elbow and at the tip become small when increasing the time parameter h. This is because the internal damping forces will become larger while increasing the equivalent damping coefficient a
1
=
h
/
2. Therefore, the displacements will become smaller in the case of no energy input.We want to emphasize that in [31] the independent interpolations are used to interpolate the incremental values of the chosen state variables for translational and rotational motion components. And the reduced integration technique is used for the tangent stiffness which is a remedy for curing the shear locking phenomena. It can be seen from Figures 15 and 16 that the results obtained by using our proposed elements are in good agreement with [31]. What is more, in this example the amplitude of vibration is the same order of the magnitude of the structure dimensions, which indicates that the present formulation can deal with these problems with the large rotations easily.


### 4.5. Flexible Beam on a Rotating Base

As shown in Figure17, a cantilever beam is fixed on a rigid base, which has a prescribed rotation φ
(
t
) with respect to the global z-axis. Relative to the moving base, the beam will have displacements at both axial and transversal directions, which are denoted by u and v, respectively.Figure 17
A flexible beam on a rotating base.This problem has been also studied in previous researches [6, 7, 32, 33] with the rotation(112)
φ
t
=
6
15
t
2
2
+
15
2
π
2
cos
⁡
2
π
t
15
-
1
0
≤
t
≤
15
s
6
t
-
45
t
>
15
swhich will also be adopted here. In this problem, 10 elements are used to study the displacements at the free end. The calculation results are plotted in Figures 18 and 19. Table 5 shows the total CPU time spent by different methods.Table 5
The comparison of CPU time for the rotating flexible beam.
Time parameterh
Integration algorithms
CPU time (s)
0
radau5
6140
0
ode45
—
0.01
ode45
3003
0.1
ode45
603Figure 18
Time evolution of the axial displacementu at the free end.Figure 19
Time evolution of the transversal displacementv at the free end.It is found that a good agreement between our simulations and reference ones [6]. With an appropriate time parameter h, the high frequencies of the system are filtered out. Consequently, the governing equation can be solved using the well-established nonstiff solver ode45.


## 4.1. Cantilever Beam with an End Moment

The first problem to be considered is a cantilever beam with material properties shown in Figure7, where a concentrated end moment is imposed at the end of the beam.Figure 7
Cantilever beam with an end moment.This problem has been used by many researchers for testing nonlinear beam elements since the analytical solutions for the problem exist. The straight beam will be bent to a full circle with the momentM
=
2
π
E
I
/
L. In this problem, 10 elements are used to model the beam. The results obtained from this study are compared with the analytical solutions as shown in Figure 8. And the deformations of the cantilever beam with different bending moments are shown in Figure 9. It can be seen that the results obtained by using our proposed elements are exactly the same with the analytical solutions.Figure 8
Comparison of the displacements of the cantilever beam with an end moment.Figure 9
Large deformations of the cantilever beam with an end moment.


## 4.2. Cantilever Beam with End Point Load

The cantilever beam subjected to a concentrated end load, as shown in Figure10, is investigated in this subsection.Figure 10
Cantilever beam with an end point load.In this problem, 10 elements are used to model the beam. The obtained results are compared with the elliptic integral solutions provided by Mattiasson [30] in Figure 11 and in numeric form listed in Table 1. The comparison indicates our results are good agreement with Mattiasson’s.Table 1
Cantilever beam with an end point load.
F
L
2
/
E
I
u
/
L
v
/
L
Mattiasson
Present
Mattiasson
Present
1.0
0.05643
0.05643
0.30172
0.30172
2.0
0.16064
0.16064
0.49346
0.49347
3.0
0.25442
0.25442
0.60325
0.60327
4.0
0.32894
0.32894
0.66996
0.66999
5.0
0.38763
0.38762
0.71379
0.71383
6.0
0.43459
0.43458
0.74457
0.74462
7.0
0.47293
0.47292
0.76737
0.76743
8.0
0.50483
0.50482
0.78498
0.78505
9.0
0.53182
0.53182
0.79906
0.79913
10.0
0.55500
0.55499
0.81061
0.81070Figure 11
Displacement curve of cantilever beam with an end point load.


## 4.3. Natural Frequency of the Cantilever Beam

The cantilever beam with the material properties as shown in Figure12 is carried out to study the natural frequency of the cantilever beam by employing the proposed beam element.Figure 12
Cantilever beam.To obtain the natural frequency, the nonlinear dynamic equation (100) should be linearized at the equilibrium positions. Firstly, without considering average strains in the expression of the nodal forces which means h
=
0, the frequency equation can be given by(106)
K
t
1
-
ω
2
M
i
n
e
=
0
,where ω is the natural frequency of the cantilever beam. In this example, the beam is modeled by N elements. Then the first bending, torsion, and tension frequencies are obtained separately as shown in Table 2. It is found that the frequencies will converge very fast to the theoretical values.Table 2
Frequency of cantilever beam with different meshes.
Modes
Frequency (Hz)
1-elems
2-elems
5-elems
10-elems
20-elems
Theory
Bending
5.291
5.269
5.266
5.266
5.266
5.266
Torsion
557.328
518.511
504.440
505.961
505.571
505.441
Tension
851.382
825.971
816.855
815.500
815.159
815.046Then, the influence of introducing the average strains on the frequencies is also studied. With a small time interval parameterh, the frequency equation becomes(107)
K
t
1
-
ω
′
2
M
i
n
e
+
a
2
K
t
1
=
0
,where ω
′ is the modified frequency of the system. Equation (107) can be simplified to the equivalent form (108)
K
t
1
-
ω
′
2
1
-
a
2
ω
′
2
M
i
n
e
=
0
.By comparing (106) with (108), the relationship between ω
′ and ω can be given by(109)
ω
′
=
ω
1
+
a
2
ω
2
.Substituting a
2
=
h
2
/
6 into (109) yields(110)
ω
′
=
ω
1
+
1
/
6
ω
/
f
h
2
<
6
f
h
,where f
h
=
1
/
h is the adjustable frequency parameter associated with the given small time interval h. According to (110), the mapping between the normalized modified frequency ω
′
/
f
h and ω
/
f
h is shown in Figure 13.Figure 13
General shape of the mapping curve betweenω
′
/
f
h and ω
/
f
h.The influence of the time parameterh on the system frequencies can be summarized as follows:(a)
All the frequencies of the system decrease to varying degrees.
(b)
The frequencies with the magnitude less thanf
h decrease very slightly.
(c)
The maximum limit of the modified frequencies is6
f
h.Table3 shows the modified frequency of cantilever beam modeled by 20 elements with different time parameters.Table 3
Modified frequency of cantilever beam with different time parameterh.
Modes
Modified frequency (Hz)
h
=
0
h
=
0.0001
h
=
0.001
h
=
0.01
h
=
0.05
h
=
0.1
Bending
5.266
5.266
5.266
5.265
5.236
5.148
Torsion
505.571
505.463
495.135
220.439
48.761
24.466
Tension
815.159
814.708
773.454
234.587
48.902
24.484Obviously, from the above discussion it is found that the high frequency of the system can be filtered out by choosing an appropriate time parameterh. This property is very useful for the time integration in the stiff problem, which is proved in the following examples.


## 4.4. Right-Angle Cantilever Beam

An L-shaped cantilever beam [6] is subjected to a concentrated load applied in the direction of global z-axis at the elbow as shown in Figure 14.Figure 14
Right-angle cantilever beam.The magnitude of this load follows the pattern of a hat function, as shown in Figure14. The material properties are(111)
E
A
=
G
A
=
1
0
6
,
E
I
=
G
J
=
1
0
3
,
L
1
=
L
2
=
10
,
A
ρ
=
1
,
J
-
ρ
=
diag
⁡
20,10,10
.The beam is meshed into 10 elements to study the displacements at both the elbow and the tip. The computations are carried out by using the radau5 and the ode45 with different values of the time parameterh. The total simulation time is 30 seconds. The relative error tolerance is set to 10−4 and the absolute error tolerance is set to 10−6. Table 4 shows the total CPU time while using the different methods. And the computed response of the elbow and the tip displacements are given in Figures 15 and 16, respectively. These results correspond to the results given by Ibrahimbegović and Mikdad [31].Table 4
The comparison of CPU time for the right-angle cantilever beam.
Time parameterh
Integration algorithms
CPU time (s)
0
radau5
6279
0
ode45
—
0.01
ode45
3348
0.1
ode45
426Figure 15
Time evolution of the displacement at the elbow.Figure 16
Time evolution of the displacement at the tip.It should be pointed out that the radau5 is a stiff ODE solver, and the ode45 is only workable for solving nonstiff differential equations. With the time parameterh
=
0, the radau5 method needs 6279 seconds to accomplish the simulation while ode45 cannot finish the task since the governing equations of this example are stiff equations. When h is set to 0.01, it can be found that the ode45 solver can complete the calculation without any difficulties and spends only 3348 seconds, which is about half of that of the stiff ODE solver radau5. For h
=
0.1, the CPU time of the ode45 method is reduced to 426 seconds. From Figures 15 and 16, it can be also found that both the displacements at the elbow and at the tip become small when increasing the time parameter h. This is because the internal damping forces will become larger while increasing the equivalent damping coefficient a
1
=
h
/
2. Therefore, the displacements will become smaller in the case of no energy input.We want to emphasize that in [31] the independent interpolations are used to interpolate the incremental values of the chosen state variables for translational and rotational motion components. And the reduced integration technique is used for the tangent stiffness which is a remedy for curing the shear locking phenomena. It can be seen from Figures 15 and 16 that the results obtained by using our proposed elements are in good agreement with [31]. What is more, in this example the amplitude of vibration is the same order of the magnitude of the structure dimensions, which indicates that the present formulation can deal with these problems with the large rotations easily.


## 4.5. Flexible Beam on a Rotating Base

As shown in Figure17, a cantilever beam is fixed on a rigid base, which has a prescribed rotation φ
(
t
) with respect to the global z-axis. Relative to the moving base, the beam will have displacements at both axial and transversal directions, which are denoted by u and v, respectively.Figure 17
A flexible beam on a rotating base.This problem has been also studied in previous researches [6, 7, 32, 33] with the rotation(112)
φ
t
=
6
15
t
2
2
+
15
2
π
2
cos
⁡
2
π
t
15
-
1
0
≤
t
≤
15
s
6
t
-
45
t
>
15
swhich will also be adopted here. In this problem, 10 elements are used to study the displacements at the free end. The calculation results are plotted in Figures 18 and 19. Table 5 shows the total CPU time spent by different methods.Table 5
The comparison of CPU time for the rotating flexible beam.
Time parameterh
Integration algorithms
CPU time (s)
0
radau5
6140
0
ode45
—
0.01
ode45
3003
0.1
ode45
603Figure 18
Time evolution of the axial displacementu at the free end.Figure 19
Time evolution of the transversal displacementv at the free end.It is found that a good agreement between our simulations and reference ones [6]. With an appropriate time parameter h, the high frequencies of the system are filtered out. Consequently, the governing equation can be solved using the well-established nonstiff solver ode45.


## 5. Conclusions

In this paper, a two-node spatial Euler-Bernoulli beam element with arbitrary rigid motion and large deformation is developed. The global displacement and rotation vectors, totally six degrees of freedom, are selected as the nodal coordinates. And the centroid position and the cross-section orientation are coupled interpolated by a special approach, which guarantees the perpendicularity between the cross-section and the deformed neutral axes. Based on this idea, a shear-locking-free beam element is proposed. By using the average strains over an appropriate time trivial to replace the transient strains in the expression of the nodal forces, the high frequencies of the system can be filtered out as discussed in the third example of Section4. This can make the governing equations of the rigid-flexible coupling system solved by using the well-established nonstiff ordinary differential equation solver. In addition, our proposed scheme can be also extended to other elements for the rigid-flexible coupling dynamic analysis.



---

*Source: 208127-2015-05-31.xml*