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.

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 [110].

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 [1922]. 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 Section 2 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 Figure 1, where a fixed global frame with orthonormal base vectors is defined and the centroid line is described by the position vector with representing the arc-length coordinate of the reference configuration.

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 , called the cross-section basis or the material basis, is employed to describe the orientation of the cross-sections. The base vectors and are directed along the principal axes of inertia of the cross-section, and is the normal vector of the cross-section; that is, . The relationship between the global base vector and cross-section base vector can be expressed aswhere is the rotation matrix with the following properties; that is,where is the identity matrix.

The time derivative of (1) leads toIn this paper, the symbol denotes the time derivative. For (2), we haveThis means that is a skew-symmetric matrix and can be denoted asThen substituting (5) into (3) results inwhere 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 leads towhere the symbol denotes the derivative with respect to the arc-length coordinate , and the skew-symmetric matrix is given belowThe axial vector corresponding to is defined as the curvature vector.

The curvature vector can be expressed with respect to the material basiswhere Einstein notation is introduced. By substituting (9) into (7), we havewhere is the permutation symbol in 3D. For (10), we have the following relationshipThe time derivative of (11) could be expressed asTherefore, we can obtain the following equation:Similarly the angular velocity can be given asBy substituting (14) into (6), finally, we getThe rotation matrix could be written as a function of rotational vector ; that is,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 bySimilarly the curvature vector can be also expressed by the arc-length derivative of the rotational vector ; that is,where denotes a matrix given byAccording to (17), the time derivative of the rotational vector can be obtained bywhereis the inverse matrix of . According to (20), the singularity problem will occur at , because 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 ; that is, ; then its complement rotation vector is defined by It can be seen that the magnitude of the new rotation vector is , which means . Note that the rotation vector and its complement represent the same rotation; that is, . 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 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 Figure 2, where the beam is subjected to the external distributed forces and moment vectors and per unit length of the reference line of centroids. And, in that figure, and are the stress-resultant force and resulting moment vectors over the cross-section, respectively.

The equilibrium equations of the infinitesimal element can be written aswith and representing the linear mass density and moment of inertia, respectively. The virtual kinetic power of the infinitesimal element is given byBy substituting (23) into (24), we can obtainThe external virtual power of the infinitesimal element can be expressed asBy using the principle of virtual power, which states that at any time , 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 byFor the Euler-Bernoulli beam, the cross-section remains perpendicular to the tangent of the centerline during deformation; that is,where is the normal strain, and we also haveSubstituting (13) and (29) into the internal virtual power equation (27) and calculating the integral over the range leads towhere and are components of and given relative to the material basis; that is,The constitutive equation of the Euler-Bernoulli beam can be described asBy substituting (32) into (30), the internal virtual power can be expressed aswhere and are the axial stiffness and the torsion stiffness, respectively. and 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 length and unchanged cross-section is considered as depicted in Figure 3.

In the initial configuration, and represent the nodal position vector and the orientation matrix of the cross-section, respectively. For simplicity, the cross-section basis is denoted by , which means that , , . The global displacement and rotation vectors are chosen as the nodal coordinates:And the element generalized coordinates can be defined byThe nodal position vector and the orientation matrix of the cross-section in the current configuration can be expressed asFor 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 . 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, . 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 is obtained by employing the Hermite interpolation; that is,whereare the Hermite shape functions with .

The unit normal vector of the cross-section is given by The other two orthogonal unit vectors of the material basis and can be determined by the successive application of the following two rotation vectors to the left end section basis .

Firstly, an intermediate reference frame arrived when rotating an angle with respect to the direction , which is illustrated in Figure 4. The corresponding rotation vector is expressed aswhereThen, the final material basis is achieved by rotating the by an angle with respect to the vector . The corresponding rotation vector isThus the orientation matrix of the cross-section can be expressed by

The above analysis shows that the first rotation vector can be determined by the interpolation of the position vector . Hence, we just need to interpolate the second rotation vector within the element. Moreover, the direction of rotation vector is , and the angle can be interpolated linearly bywere and the relative torsion angle of the right end section can be obtained according to the following steps:(a)The rotation vector can be calculated according to . The extraction of from 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 .

It should be mentioned that (40) cannot provide a unique solution if the first rotation angle is equal to , which means . 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 byIt 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 .

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 byFor 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 in which the identity is used in the simplification. According to the interpolation of the global position vector (37), the velocity of the centroid can be expressed as where denotes the velocity of the element node and is the angular velocity. The element generalized velocity is given bySubstituting (48) into (47) yieldsTherefore, the time derivate of the average normal strain and its virtual variation can be expressed bywhereBased 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 byFurthermore, the components of curvature with respect to the material basis are given by The two bending curvature components and , in terms of (11), are written asThe torsion component isBy substituting (40) and (42) into the two items on the right-hand side of (56), respectively, we obtainwhereAfter substituting (57) into (56), we haveThe time derivatives of (55) and (60) can be expressed byin which the time derivative of the unit normal vector is given byTherefore, the derivative of (62) with respect to arc-length coordinate isAccording to the decomposition of the cross-section rotation (43), the relationship of the angular velocities at two ends of the element can be expressed asSo the time derivative of iswhereand is constant matrix with which the angular velocity can be denoted by .

The time derivative of (58) iswhereThe angular velocity can be expressed with respect to the material basisBased on the relationship between the derivatives with respect to time and arch-length coordinate , we havewhereThen the time derivative of the other base vectors and can be given byBy substituting (63), (65), (67), and (72) into (61), the time derivatives of the curvature and their virtual variations can be expressed aswhereBy substituting the normal strain , torsion , two bending curvatures , , and their virtual time derivatives (51) and (73) into the internal virtual power equation (33), we obtainwhere is called the generalized nodal force:The time derivative of (76) can be written aswhereThe first part is generally called the material stiffness matrix, and the second part 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 byThen the average normal strain over a small time interval is defined aswhere the coefficients are expressed as and . Thus the average rotational strains are given byBy replacing the transient strains in (76) with average strains (80) and (81), the modified generalized nodal force can be written aswhere

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 asThe accelerated velocity of the centroid can be obtained by the time derivation of (48):Therefore, the virtual kinetic power caused by the translation can be expressed aswhereThe angular velocity of the cross-section isAccording to (70), the angular velocity in the material form and its virtual variation can be expressed aswhereThe angular acceleration iswhereThe virtual kinetic power caused by the rotation may be expressed bywhere is the centroidal inertia matrix of the cross-section with respect to the material basis. Substituting (89) and (92) into (93) yieldswhereIn terms of (86) and (94), the virtual power of inertial forces can be written aswhere

3.4. Governing Equation

The virtual power equation of the beam element can be written asin which is the external virtual power and is the generalized external force. Substituting the internal virtual power (75) and the inertial virtual power (96) into the virtual power equation (98) yieldsBy substituting the modified generalized nodal force (82) into (99), the governing equations of beam element can be expressed asThus, the generalized acceleration iswhere the mass matrix and generalized force can be written by

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 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 is equivalent to the internal damping force. If represents the rigid motion only, this equivalent damping force will disappear since is equal to zero for the allowed rigid motion.

The element generalized velocity and the generalized coordinates are set as the state variables; that is,in which the relationship between the time derivative of and becomesThe governing equations, which are a set of the second order ordinary differential equations, can be written asEquations (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 Section 2, 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.

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 Figure 7, where a concentrated end moment is imposed at the end of the beam.

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

4.2. Cantilever Beam with End Point Load

The cantilever beam subjected to a concentrated end load, as shown in Figure 10, is investigated in this subsection.

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.

4.3. Natural Frequency of the Cantilever Beam

The cantilever beam with the material properties as shown in Figure 12 is carried out to study the natural frequency of the cantilever beam by employing the proposed beam element.

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 , the frequency equation can be given bywhere is the natural frequency of the cantilever beam. In this example, the beam is modeled by 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.

Then, the influence of introducing the average strains on the frequencies is also studied. With a small time interval parameter , the frequency equation becomeswhere is the modified frequency of the system. Equation (107) can be simplified to the equivalent form By comparing (106) with (108), the relationship between and can be given bySubstituting into (109) yieldswhere is the adjustable frequency parameter associated with the given small time interval . According to (110), the mapping between the normalized modified frequency and is shown in Figure 13.

The influence of the time parameter 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 than decrease very slightly.(c)The maximum limit of the modified frequencies is .

Table 3 shows the modified frequency of cantilever beam modeled by 20 elements with different time parameters.

Obviously, from the above discussion it is found that the high frequency of the system can be filtered out by choosing an appropriate time parameter . 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 -axis at the elbow as shown in Figure 14.

The magnitude of this load follows the pattern of a hat function, as shown in Figure 14. The material properties are

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 parameter . 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].

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 parameter , 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 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 , 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 . This is because the internal damping forces will become larger while increasing the equivalent damping coefficient . 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 Figure 17, a cantilever beam is fixed on a rigid base, which has a prescribed rotation with respect to the global -axis. Relative to the moving base, the beam will have displacements at both axial and transversal directions, which are denoted by and , respectively.

This problem has been also studied in previous researches [6, 7, 32, 33] with the rotationwhich 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.

It is found that a good agreement between our simulations and reference ones [6]. With an appropriate time parameter , 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 Section 4. 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.

Conflict of Interests

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

Acknowledgment

The authors are grateful for the National Science Foundation of China (11372057).