Abstract
Based on geometrically exact beam theory, a hybrid interpolation is proposed for geometric nonlinear spatial Euler-Bernoulli beam elements. First, the Hermitian interpolation of the beam centerline was used for calculating nodal curvatures for two ends. Then, internal curvatures of the beam were interpolated with a second interpolation. At this point, C1 continuity was satisfied and nodal strain measures could be consistently derived from nodal displacement and rotation parameters. The explicit expression of nodal force without integration, as a function of global parameters, was founded by using the hybrid interpolation. Furthermore, the proposed beam element can be degenerated into linear beam element under the condition of small deformation. Objectivity of strain measures and patch tests are also discussed. Finally, four numerical examples are discussed to prove the validity and effectivity of the proposed beam element.
1. Introduction
In engineering practice, a majority of three-dimensional beams can be considered slender beams and can thus be regarded as Euler-Bernoulli beams. Much research in recent years has been focused on modeling geometrical nonlinearity for 3D Euler-Bernoulli beams subjected to large deformations. Because spatial large rotations are physically nonadditive, an improper discretization of spatial rotations may lead to nonobjectivity of strain measures [1]. Different from the linear beam theory, large rotations of beam cross sections cannot be totally determined by the beam centerline. However, displacement and rotation of the cross section of a slender beam must satisfy the Kirchhoff constraints. Therefore, there is a critical need to develop an efficient interpolation method to solve these problems.
Geometrically exact beam theory, developed by Reissner [2, 3] and Simo [4, 5], has been widely adopted for modeling of geometrical nonlinear beams [6–11]. For a geometrically exact beam, the configuration of a spatial beam can be described by the position vector of the beam centerline and an orthogonal matrix, which specifies rotation of a cross section. A curvature vector is used to describe the rotational strain measure, which is anisotropic and proportional to the stress resultants. Relationship between rotation and rotational strain was found to be consistent with the virtual power principle and is valid for configurations with arbitrary large displacements and rotations. Due to independent interpolations of displacement and rotation, use of geometrically exact beam theory may lead to shear-locking problem for slender beams and nonobjective interpolation of strain measure [1]. In order to avoid these problems, many researchers presented different formulations of geometrically exact beam elements. The core technology in developing these formulations has been on the discretization of rotational degrees of freedom or the incremental rotation [6, 11–13]. In spite of the efforts, a problem still exists with accurate determination of an analytical expression of elastic forces using the displacement-based geometrically exact beam theory.
To avoid the interpolation of rotation parameters, an absolute nodal coordinate formulation (ANCF) was proposed by Shabana [14–18]. Three polynomials were used as the assumed displacement field in the ANCF-based beam element. In ANCF, 12 variables were used for each node of a 3D beam element, which includes the position vector and 9 slopes. Therefore, dimension of governing equations is significantly larger in comparison with traditional beam elements. Based on the interpolation of beam curvatures, a new beam element was proposed by Zupan [7, 19], in which displacement and rotational vector were not interpolated at all. In this formulation, strain measure-interpolation based element was demonstrated to preserve the objectivity, namely, invariant with rigid body motions, and path independence of strain measures. However, use of the beam element based on the curvature interpolation was not convenient when used in calculating nodal displacement and rotation, since they cannot be analytically integrated by strain measures.
In this paper, a hybrid interpolation scheme for geometrically exact Euler-Bernoulli beams is proposed. Global nodal displacement and the rotational vector are the basic unknown variables for the formulation. Interpolation of the beam centerline was used to determine nodal value and its derivatives of the curvatures. These new parameters were also interpolated. Using combined interpolation of the nodal displacement and strain measures, an analytical expression for nodal forces was formulated and objectivity and path independence of strain measures were satisfied.
The outline of this work is listed as follows. In Section 2, basic theory of geometrically exact Euler-Bernoulli beams with large deformation is briefly reviewed. In Section 3, based on the principle of virtual power, weak form of the dynamic equilibrium equations for 3D Euler-Bernoulli beams was deduced using the novel method proposed here, which is concise and convenient for large rotations. Next, hybrid interpolations of nodal displacements and curvatures are introduced in Section 4. In Section 5, governing equations for the discrete finite beam element model are given. Property tests and applications of the beam element are presented in Sections 6 and 7, respectively.
2. General Theory of Geometrically Exact Euler-Bernoulli Beams
In this section, geometrically exact Euler-Bernoulli beam theory is reviewed, including kinematics of large deformation beams, strain measures, stress resultants, and constitutive relations.
2.1. Kinematics and Strain Measures
The initial and deformed configurations of a spatial Euler-Bernoulli beam in the inertial reference frame () are shown in Figure 1. The main assumption of 3D Euler-Bernoulli beam theory is that all cross sections keep rigid and perpendicular to the tangent vector of the centerline during deformation. The configuration can be fully described by the following.
(1) The position vector of beam centerline, where is the arc-length coordinate of the beam centerline at the initial configuration and is the initial length of the beam.
(2) The cross-sectional reference frame : the orientation of a cross section is denoted by the right-handed orthonormal triads attached to the cross section centroid. With the Kirchhoff constraint, the normal vector of the cross section is parallel to the tangent vector of the beam centerline , which yieldswhere a symbol with a superscript “” denotes its derivative with respect to the arc-length coordinate . is the axial strain. The second and third axial vectors and are in the directions of two principal axes of inertia of the cross section. Since , it is obvious that the tangent vector of the beam centerline does not have any components in the directions and ; that is,The equation means that the beam is shear-free. In consequence, orientation of a cross section is in some sense determined by the beam centerline and the position vector derivatives seeing (1) and (2). That is to say, the so-called Kirchhoff constraint is carried out accurately in a strong manner.
Under the definition of the local coordinate system of the beam cross section, the rotational matrix of the cross section with respect to the inertial one is given by
The cross-sectional rotation is parameterized by rotational vector . Hence, the rotational matrix can be expressed as here with . In this paper, a skew-symmetric tensor with symbol “” denotes the cross product of two vectors; that is, . The skew-symmetric tensor is given by
The current cross-sectional orientation matrix in (4) could be decomposed into the original orientation of the cross section and the rotation of the local coordinate system; that is,where is the rotation of the cross section. The time derivative of base vectors of the cross-sectional reference frame can be obtained by the cross product of its angular velocity and the base vectors:where the superscript dot denotes the derivative with respect to time and the angular velocity of the cross section can be expressed byhere withThe skew-symmetric tensor of the angular velocity can be expressed in terms of the time derivative of the rotational matrix referring to (7):
Similarly, the derivative of base vector with respect to the arc-length coordinate could be given by the cross product of the curvature and base vectors:Curvature is used as the strain measure in Reissner’s work [2, 3]. Simo [4, 5] has also discussed it in detail in geometrically exact beam theory. The skew-symmetric tensor of the curvature vector can be formulated by the derivative of rotational matrix with respect to the arc-length parameter referring to (11):Furthermore, the curvature vector is also represented by the derivative of rotational vector with respect to the arc-length coordinate:Angular velocity and curvature vector can be expressed in the cross-sectional coordinate system as follows:where denotes the cross-sectional torsional curvature and and are the cross-sectional bending curvatures in the two perpendicular directions.
The axial strain can be obtained by the elongation referring to (1):
2.2. Stress Resultants and Constitutive Relations
The stress-resultant force and moment of a cross section can be expressed in the cross-sectional coordinate system as follows:
The constitutive relation between the strain measure and the resultant force and moment of stress for the beam of linear elastic material is given by [3, 4]where , , , are the axial tensile, torsional, and principal bending stiffness, respectively. In some sense, when beams are regarded as one-dimensional structures, their constitutive relations are anisotropy. This feature must be taken into account in the modeling of beam elements.
3. Formulation of 3D Geometrically Exact Euler-Bernoulli Beam
In most formulations of linear finite elements, the principle of virtual work is adopted. However, if rotations are large, it is difficult to obtain expression of the rotational virtual work. Due to the nonadditive characteristic of 3D rotations, product of the virtual rotational vector and moment do not constitute virtual work. There are kinds of parameters describing large rotation. Different parameters have different meanings. The principle of virtual power has a natural advantage in finite element formulation of large rotations since rational virtual power is the product of the virtual angular velocity and moment. Virtual angular velocity is definitely conjugate to internal moment, and angular velocity has clear physical meanings. In this work, principle of virtual power is adopted to deduce the finite element formulation of 3D Euler-Bernoulli beams with large rotations.
A beam can be regarded as a series of rigid bodies under the assumption of rigid cross section. In such case, current arc-length coordinate of the beam centerline is a function of time and the initial arc-length coordinate . Therefore, all strains, stresses, and internal forces can be regarded as functions of these two variables. Thus, axial strain is given by
An infinitesimal beam is shown in Figure 2, where and are the stress-resultant force and moment, respectively. In order to simplify the derivation, a reasonable assumption could be made that there is no external force on the infinitesimal beam. The strong form of the dynamic equilibrium equations of the infinitesimal beam can be expressed aswhere is the cross-sectional moment of inertia and is the mass density per unit length before deformation.
The virtual power produced by applied forces of the infinitesimal beam can be formulated asThe virtual power of inertia forces of the infinitesimal beam is given byAccording to the principle of virtual power, the virtual power of internal forces can be given by the difference of the virtual powers of applied forces and inertia forces, which yields
From another point of view, the virtual power of inertia forces also can be obtained through the integration of (22) along the beam length:
The components in (14) could be given byThe first-order derivatives of (14) with respect to arc-length coordinate and time , respectively, areAs arc-length and time are independent, the derivation orders of base vectors can be exchanged. Therefore, the following relationships can be obtained referring to (24):Considering (28) and (27), (26) could be rewritten as Taking the time derivative of (1), one yieldsWith (29)-(30) and (23), the virtual power of internal forces can be expressed as
All parameters in the virtual power equations are defined in the cross-sectional coordinate system, which do not contain rigid body motions, and the beam element formulated using the present method cannot produce self-strains, even for large deformations.
4. Discretization
It can be concluded from Section 3 that the internal virtual power of 3D Euler-Bernoulli beams is directly dependent on the strain measure rather than the displacement and rotation. Unlike the rotational vector, the strain measure is additive. It cannot cause nonobjective strain measure [7] by the strain interpolation instead of the displacement or rotational interpolation. Nodal coordinates which are global nodal displacement and the rotational vector are the basic unknown variables. Based on a Hermitian interpolation of the beam centerline, nodal value and its derivatives of the bending curvature are consistently determined which are afterwards reinterpolated again by Hermitian polynomials. A constant axial strain field as well as a constant torsion field is constructed from the nodal degrees of freedom based on the geometrical approximations.
4.1. Nodal Coordinates
Basis vectors of cross sections at the two ends are and , respectively. The rotational vectors on the two ends are and , respectively. A beam is divided into a series of beam elements. The left and right cross-sectional orientations of a node are not the same when two elements are joined at the node in an intersection angle different from 180°. The left and right original orientations of the cross sections of a node may be different in a nonstraight beam. However, the rotations of these two cross sections are the same. Nodal displacements , and rotational vectors , are adopted as the basic unknown functions of a two-node beam element, namely, nodal coordinates:The corresponding generalized velocities include the nodal velocities , and the angular velocities , :
4.2. Interpolation of the Beam Centerline
The relative displacement between two nodes of an element after deformation can be expressed as
If the beam element is small enough, the centerline after deformation could be regarded as a circular arc, as shown in Figure 3. The angle between the two base vectors and is denoted by ; that is, . The current element length is related to the chord length , which yieldswhere . In an element, the angular cannot be greater than . Thus, is definitely positive and approaches 1 when tends to 0. This conclusion conforms to the reality.
To fulfill C1 continuity at the element boundaries, the beam centerline of an element is assumed by the Hermite interpolation. Due to the situation that beams undergo large displacements and rotations, three components of the position vector on the beam centroid line need to be treated in the same way:where , are nodal position vectors and , are unit tangent vectors of beam centerline at the two nodes. The influence of axial deformation has been considered in the expression of the current element length , as shown in (38). The Hermite shape functions are
Equation (36) for the interpolation of the beam centerline is of great help for calculating the nodal bending curvatures.
4.3. Resolving Nodal Strain Measures
As a kind of slender structure, the beam axial tensile and torsional strains are small. The axial strain could be assumed as a constant along the element centerline; that is,
Based on Euler-Bernoulli beam theory, the cross-sectional normal vector can be derived by taking the derivative of (36) with respect to . The first- and second-order derivatives of with respect to arc-length are
The bending curvatures are expressed by cross-sectional base vectors and their derivatives with respect to arc-length referring to (11):Referring to (11) and (25), the derivatives of bending curvatures in (40) with respect to arc-length are given as follows:Taking the derivatives of shape functions, the derivatives of base vector on the beam ends can be obtained referring to (39):Nodal curvature parameters in two perpendicular bending directions are resolved by substituting (42) with (40)-(41):Nodal curvature parameters of a beam element are rewritten asWith the given process, nodal curvature parameters and their derivatives with respect to arc-length parameter are expressed by nodal coordinates.
4.4. Discretization of Curvature Components in Space
In the original theory of geometrically exact beam, the analytical expression of centerline curvatures is almost impossible to be calculated by the nodal rotational parameters. Thus, a numerical integration technique is implemented. Unfortunately, the numerical integration needs to be taken for all elements at every step. In this work, the internal centerline curvatures are given directly by Hermitian interpolation, which makes the curvature expression explicit and analytically integrable. This property can benefit the calculation of explicit nodal forces.
The torsional curvature could be assumed as a constant in the beam element. Basis vector has three components in the left end cross-sectional frame after deformation, as shown in Figure 4:where and .
The relative torsional angle between the left and right end sections of a beam element can be calculated by the arcsine function according to the geometrical relationship:Torsional curvature in the element can be approximated as the average:
Nodal bending curvatures and their arc-length derivatives could be consistently derived from the displacement field in Sections 4.2 and 4.3. It affords a possibility to employ cubic Hermite polynomials for bending curvatures and :
Cubic Hermite interpolation for bending curvatures has third-order accuracy, which can make the interpolation curve comparatively smooth. With the above description, four independent strain measures in the cross-sectional reference frame are discretized in space.
4.5. Time Derivatives of Strain Measures
Rates of the strains can also be formulated by nodal parameters. The time derivative of the current beam element length in (35) is expressed aswhereThe time derivative of the axial strain can be given referring to (49):The time derivative of the torsional curvature equation (47) is given in the following equation:where . Taking the time derivative of nodal curvatures (43), the rates of the nodal curvature parameters can be given as follows:
5. Discrete Governing Equations
In this section, the discrete governing equations are formulated based on the dynamic equilibrium equations introduced in Section 3. Virtual power of a beam element is determined by the sum of the virtual power of internal, inertia, and external forces. These forces are introduced separately.
5.1. Elastic Nodal Forces
Using the interpolations of the strain measures (38) and (47)-(48), the internal virtual power of an element is calculated by integrating (31) along the beam domain as follows:where the coefficient matrix isBy substituting the virtual variations of the strain rates in (51)–(53) into (54), the virtual power produced by internal forces of the element can be rewritten in the following form:where , , , and are corresponding nodal forces and moments, which are power-conjugated to , , , and , respectively. Expression of the left end nodal force isThe right end nodal force isThe left end nodal moment of the element isThe right end nodal moment of the element iswhere the nodal parameters are defined in matrix form as
Examining these expressions of nodal forces, it is easy to verify thatIt can be proved that nodal forces described in (62) constitute an equilibrium force system. Hence, the presented beam element in this work is applicable to dynamic analysis of flexible multibody systems. Furthermore, the stiffness matrix of the proposed beam element is symmetric. Nodal forces were derived in the current configuration taking into consideration the coupled effects of the axial tensile, torsional, and bending deformations. Finally, the proposed beam element allows the explicit expression of nodal forces to be determined without using the numerical integration. This has a distinct advantage in the numerical solution of governing equations.
5.2. Mass Matrix
The displacement and rotation of the beam element are not directly discretized in space. Hence, the mass matrix and inertia force of a beam element are given by independent interpolations in this section. The position vector of an arbitrary point in a rigid cross section of Euler-Bernoulli beam could be expressed as where is the position vector of the centroid and is the relative position vector with respect to the centroid. The virtual power of the inertia forces of the beam element can be written as Using global nodal coordinates, the equation of flexible multibody system dynamics contains no centripetal force or Coriolis force.
If the acceleration caused by the axial extension is ignored, the velocity and acceleration of the beam centerline are given referring to (36):The change of the cross-sectional angular velocity in a beam element is usually small. Hence, a linear interpolation is employed; that is, the angular velocity of the beam cross section is expressed in the cross-sectional frame as follows:where and are angular velocities of the nodes. Angular velocities in inertia coordinate system areThe inertia tensor of an arbitrary cross section in its cross-sectional frame can be expressed as . The virtual power of inertia forces in (64) can be rewritten aswhere , is the mass matrix of the beam element, and is the inertia force:
5.3. The Governing Equations
Based on the principle of virtual power, the virtual power produced by external forces is equal to the sums of the virtual powers of the internal and inertia forces:where is the generalized external force. The governing equation of the beam element can be written as The state variables of the governing equation are the nodal coordinates in (32) and the general velocity in (33):The governing equation can be rewritten as a set of first-order differential equations, which can be solved by the well-developed numerical methods for ordinary differential equations:
6. Tests of the Proposed Beam Element
As mentioned in Section 5.1, an explicit expression of nodal forces without using the numerical integration can be given by the proposed large deformation beam element. Therefore, a linear beam element can be evolved. Because of the additivity of strain measures, the path independence of numerical solutions has been proved by Zupan and Saje [7]. In addition, the objectivity of strain measures and the convergence of the beam element are tested.
6.1. Degenerate Linear Beam Element
With the assumption of small displacements, the equilibrium equation is formulated in the initial configuration. Simplifications for the basis vectors of cross-sectional frames on the beam ends are given by where the frame is the initial cross-sectional reference frame of a straight uniform cross-sectional element. Hence, . Linearizing the relative displacement in (34), one obtainsAs and , the following approximations could be made for the intermediate parameters:Neglecting the order of strains, the linearized nodal forces of the beam element are expressed as follows:
Under small deformation, the linearized nodal forces still satisfy the equilibrium equations (62). The parameters , (, ), , and are corresponding bending, torsional, and axial strains, respectively. The proposed element can degenerate to classical linear beam element.
6.2. Objectivity of Strain Measures
In order to test the objectivity of strain measures, the rigid body motions, including the translation and the rotation , are superimposed onto a deformed beam element. Since is a unit matrix, the increment of intermediate variables and strain measures are zeros referring to (47), (38), and (43): The discrete bending curvatures given by (48) are also unaffected by the rigid body motions:It is proved that the strain measures in material reference frame are invariant under the rigid body motions.
6.3. Convergence and Patch Test of the Beam Element
In this section, the convergence of the proposed beam element is tested. The proposed beam formulation is displacement-based, considering its primary variables. Since the interpolation curve of the beam centerline has C1 continuity, the finite element solution converges to the exact solution when the length of a beam element tends to zero.
Instead of displacement and rotation parameters, strain measures are the discretized parameters in the proposed beam element. Hence, patch test is needed to verify the convergence of the element. In this section, patch tests of the three kinds of constant strains are given below.
Without loss of generality, an arbitrary finite element mesh for a beam is given in Figure 5. The node connects the elements and . and are lengths of the elements. The material properties are uniform. Loads and displacements of these nodes are chosen to make sure that strains of the elements are constant. No body force is applied on the elements and no external force is acted on the node . The element will be convergent if the elastic force on the node satisfies the equilibrium equation.
(1) The axial strain in elements is a constant, and other strain measures are all zeros. A pair of equal and opposite forces along -axis is applied on nodes and . Displacements and rotations of the nodes should be chosen from the following modes: where with the corresponding subscript are displacements of the nodes, the first subscript of basis vectors is the number of node, and the second subscript is the number of base vector. Parameters () are constants. The curvature parameters are zeros in the aforementioned displacement modes. By substituting (80) into the expression of nodal forces, the node forces on node have only one nonzero component which is induced by the axial strains:Thus, the following equilibrium equation of node is deduced:
(2) The torsional curvature along the elements is constant; that is, , and other strain measures are all zeros. A pair of equal and opposite torques along -axis is applied on nodes and . Displacements and rotations of the nodes are chosen from the following modes:where is the angle between two base vectors and . The base vectors , , and are along one straight line, which yields . Therefore,The equilibrium equations of node in (82) are satisfied.
(3) The bending curvature along the elements is constant; that is, , and other strain measures are all zero. A pair of equal and opposite moments along -axis is acted on nodes and . Displacements and rotations of the nodes are chosen from the following modes:The derivatives of the curvature components with respect to arc-length coordinate tend to zero when the angles , are infinitesimal. The nonzero parameters for the pure bending have the following relationships:By substituting (86) into node forces formulas, it obtainsThe equilibrium equations of node in (82) are fulfilled.
The equilibrium equations of node verify that the proposed element passes the patch test when the element size becomes infinitesimal. Hence, the proposed element is convergent.
The Hermitian interpolation of the beam centerline is suitable for linear beams. However, it gives lower precision for large curvatures. The projection of the vector on the base vector for Euler-Bernoulli beams is theoretically zero. However, the projection at the beam end node equals zero only when their included angles are infinitesimal. More elements need to be employed for beams with large deformation to derive the real curvatures.
7. Numerical Examples
Four numerical examples are examined to establish the validity, accuracy, and the convergence performance of the proposed beam element. First three examples deal with static analysis of a cantilever beam, an initially deformed beam, and a spatial beam structure. The last example is a dynamic problem of a pendulum with the stiffening effect taken into consideration.
The procedure is written in Matlab computing environment. The fsolve function in Optimization Toolbox of Matlab was employed to solve static problems. TolFun was set to 10−18 as the size of the latest change in sum of squares of function values and gradient of the sum of squares. TolX was set to 10−18, which is the lower bound of change in a step size. Iterations were terminated when the last step was smaller than TolFun or TolX. The ode15s function was employed to solve dynamic equations. The initial step was set as 10−4. Error control properties were chosen as relative error tolerance RelTol = 10−4 and absolute error tolerance AbsTol = 10−4.
Example 1 (a cantilever beam subjected to an end moment). An initially straight cantilever beam is subjected to a concentrated moment at the free end, as shown in Figure 6. The cantilever beam undergoing large deformation, which can be described as of pure bending, was investigated by many researchers [20–22], because an analytical solution exists in this example. Here, curvature radius is related to the applied moment and the bending stiffness .
When satisfies the equation , the beam is curved into a complete circle and rotational angle of the end section is . The cantilever beam is divided into elements. The rotational angle of the end section is given in Table 1. Furthermore, the errors with respect to the analytical solution are also given in Table 1, which shows that errors reduce when number of elements increases. Figure 7 shows deformed shapes of the cantilever beam when subjected to different end moments . In this particular case, the cantilever beam is divided into 10 elements. Figure 8 shows results for dimensionless vertical and horizontal tip displacements as a function of dimensional load , which are also compared with the analytical results. It can be seen from the comparison that numerical and analytical results are in good agreement.
Example 2 (cantilever beam initially curved to ). A cantilever beam initially 45-degree bending subjected to a concentrated end load is shown in Figure 9. This example was first proposed by Bathe and Bolourchi [23] and has been employed by numerous researchers [5, 10, 11, 21, 24].
The bending beam lies in the plane and has an average radius of 100 m, as shown in Figure 9. Initial tip position is (70.71 m, 29.29 m, 0). A vertical concentrated load N was applied at the tip along the direction. The beam cross section was a square with an area of 1 m2. The material was linear elastic. Elastic modulus and Poisson’s ratio were Pa and , respectively.
Initially, the structure can be regarded as a curved beam with constant curvature, which maintains equilibrium under load and displacement boundary conditions. Eight equal straight beam elements were used to approximate the curved beam. An increasing magnitude of concentrated load on the beam end tip was applied in 40 load step increments. Figure 10 shows the evolution of tip displacement as a function of the applied load modulus. This result was compared to results reported by Bathe and Bolourchi [23], in which beam shear deformations were neglected. In Bathe’s results, the Newton-Cotes formula of order and Gauss integration of order were employed for the isoparametric elements. Our results coincide with results reported by Bathe, as shown in Figure 8. Relative error of the numerical solution was found to increase with increase in displacement amplitude. Displacement magnitude of the beam end tip under 600 N force was determined as 73.5308 m in this work whereas the reference value is 73.022 m in Bathe’s work, as shown in Table 2. Relative error was 0.7%. In the proposed beam element, assumptions about displacement model have been given up. The proposed beam element avoids overstiff conditions and tends to be flexible.
Example 3 (spatial beam structure under coupling deformation). A spatial beam structure was anchored on the support, as shown in Figure 11. Three straight beams of 1 m were connected perpendicular to each other. Cross section was a square with an area of 0.01 m2. Two forces were loaded on the tip along the and directions. Elastic modulus and Poisson’s ratio were Pa and , respectively.
Every straight beam was divided into 4 elements. 80 load steps were used to observe load-deflection curves. As shown in Figure 12, load-deflection curves of the beam tip were compared to results from absolute nodal coordinate formulation (ANCF) [25]. Force-displacement curves were found to agree with results of ANCF. Displacement error compared to results in [25] was found to reduce with more elements in every straight beam, as shown in Table 3.
Example 4 (flexible pendulum subjected to spin-up motion). This example considers a clockwise spinning motion of a flexible pendulum, as shown in Figure 13. The pendulum was attached to a shaft orthogonal to the beam axis and driven by a prescribed rotation. Left end of the beam was on the origin of the inertia coordinate system. Angular velocity is given by This is a popular example to test geometric stiffening effect and has been studied by many researchers [20, 26, 27]. This spin-up model is usually motivated by a spacecraft antenna and helicopter blade. The beam has a length of m. Cross-sectional area is m2 and the second-area moment of inertia is m4. Elastic modulus used was Pa. Mass density was kg/m3. First three natural frequencies were determined as 2.9 Hz, 18.2 Hz, and 51.1 Hz. The pendulum reached a steady-state angular velocity after the spin-up time s. For this simulation, gravity was not considered.
The pendulum was divided into 16-beam element. Steady-state angular velocity of rad/s was used for simulation. Transversal deformation of the tip is given in Figure 14(a). A comparison of the deformation curve agrees well with the results of [27]. In the first 15 seconds, angular acceleration due to the inertia forces caused bending deformation of the flexible pendulum. When angular velocity was in steady state, centrifugal force straightened the flexible beam. Steady-state angular velocity rad/s was also tested. Transversal deformation of the tip is given in Figure 14(b). Transversal deformation decreases under higher rotational speed, which can be observed from the results shown in Figures 14(a) and 14(b). Geometric stiffening effects could be reflected in this dynamic simulation.
(a) rad/s
(b) rad/s
8. Conclusions
In this paper, a new 3D geometrically exact Euler-Bernoulli beam element was proposed based on a hybrid interpolation of the beam centerlines and curvatures. The interpolation of the beam centerlines determined nodal curvatures, where the C1 continuity was satisfied, and nodal strain measures were consistently formulated. Analytical expression of nodal forces can be obtained using this method, which has an advantage in numerical simulation. While deformation is small, the proposed beam element can degenerate into classical linear beam element. Furthermore, objectivity and path independence of strain measures were satisfied by the cubic Hermitian interpolation of the beam curvatures. Through four numerical examples, validity of the proposed beam element was proved by the static analysis of 3D Euler-Bernoulli elements with large deformation and dynamic analysis of flexible multibody systems.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
The authors are grateful for the National Natural Science Foundation of China (11372057). Their sincere thanks are due to Gang Wang for assistance with the revision of this paper.