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

A two-node spatial beam element with the Euler-Bernoulli assumption is developed for the nonlinear dynamic analysis of slender beams undergoing arbitrary rigidmotions 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.


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][2][3][4][5][6][7][8][9][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][20][21][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 crosssection 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.

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 (g 1 , g 2 , g 3 ) is defined and the centroid line is described by the position vector r() 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 (e 1 , e 2 , e 3 ), called the cross-section basis or the material basis, is employed to describe the orientation of the crosssections.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  and crosssection base vector e  can be expressed as where R is the rotation matrix with the following properties; that is, where I is the identity matrix.

Shock and Vibration 3
The time derivative of (1) leads to In this paper, the symbol ( ̇) denotes the time derivative.For (2), we have This means that ṘR  is a skew-symmetric matrix and can be denoted as Then substituting ( 5) into (3) results in 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 arclength coordinate  leads to where the symbol ( )  denotes the derivative with respect to the arc-length coordinate , and the skew-symmetric matrix κ is given below The axial vector  corresponding to κ is defined as the curvature vector.
The curvature vector  can be expressed with respect to the material basis where Einstein notation is introduced.By substituting (9) into (7), we have where   is the permutation symbol in 3D.For (10), we have the following relationship The time derivative of (11) could be expressed as Therefore, we can obtain the following equation: Similarly the angular velocity  can be given as By substituting ( 14) into (6), finally, we get The rotation matrix R 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 by Similarly the curvature vector  can be also expressed by the arc-length derivative of the rotational vector   ; that is, where T  denotes a matrix given by According to (17), the time derivative of the rotational vector can be obtained by where is the inverse matrix of T  .According to (20), the singularity problem will occur at  = 2,  ∈  \ {0}, because 2/2 tan  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 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].

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 n() and m() per unit length of the reference line of centroids.And, in that figure, N() and M() are the stress-resultant force and resulting moment vectors over the cross-section, respectively.
The equilibrium equations of the infinitesimal element can be written as with   and J  representing the linear mass density and moment of inertia, respectively.The virtual kinetic power of the infinitesimal element is given by By substituting ( 23) into (24), we can obtain The external virtual power of the infinitesimal element can be expressed as By 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 by For 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 have Substituting ( 13) and ( 29) into the internal virtual power equation ( 27) and calculating the integral over the range [0, ] leads to where   and   are components of N and M given relative to the material basis; that is, The constitutive equation of the Euler-Bernoulli beam can be described as By substituting ( 32) into ( 30), the internal virtual power can be expressed as where  and  are the axial stiffness and the torsion stiffness, respectively. 2 and  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].

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, r 0  and R  0 represent the nodal position vector and the orientation matrix of the crosssection, 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: And the element generalized coordinates p  can be defined by The nodal position vector and the orientation matrix of the cross-section in the current configuration can be expressed as 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   = (1 +   )n  .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, where are the Hermite shape functions with  = /,  ∈ [0, 1].
The unit normal vector of the cross-section is given by 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 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 where 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 Thus 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 r().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 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 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 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.

Average Strains over a Small Time Interval and Nodal
Forces.According to (28), the normal strain within the element is given by 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 Δ.
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 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 in which the identity ṅ  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 where ṙ  = u  denotes the velocity of the element node and   is the angular velocity.The element generalized velocity is given by Substituting (48) into (47) yields Therefore, the time derivate of the average normal strain and its virtual variation can be expressed by where 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 Furthermore, 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 as The torsion component   is By substituting (40) and (42) into the two items on the righthand side of (56), respectively, we obtain where After substituting (57) into (56), we have The time derivatives of ( 55) and (60) can be expressed by in which the time derivative of the unit normal vector n is given by Therefore, the derivative of (62) with respect to arc-length coordinate  is 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 So the time derivative of  2 is where and S   is constant matrix with which the angular velocity can be denoted by   = S   q  .The time derivative of (58) is where The angular velocity  can be expressed with respect to the material basis Based on the relationship between the derivatives with respect to time  and arch-length coordinate , we have where Then the time derivative of the other base vectors t and b can be given by By substituting (63), ( 65), (67), and (72) into (61), the time derivatives of the curvature and their virtual variations can be expressed as where By substituting the normal strain   , torsion   , two bending curvatures   ,   , and their virtual time derivatives (51) and (73) into the internal virtual power equation (33), we obtain where F int is called the generalized nodal force: The time derivative of (76) can be written as where The first part K 1 is generally called the material stiffness matrix, and the second part K  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 Then the average normal strain over a small time interval ℎ is defined as where the coefficients are expressed as  1 = ℎ/2 and  2 = ℎ 2 /6.Thus the average rotational strains are given by By replacing the transient strains in (76) with average strains (80) and (81), the modified generalized nodal force can be written as where (83)

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 The 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 as where The angular velocity of the cross-section is According to (70), the angular velocity in the material form and its virtual variation can be expressed as where The angular acceleration is where The virtual kinetic power caused by the rotation may be expressed by where J  = R   J  R  is the centroidal inertia matrix of the cross-section with respect to the material basis.Substituting (89) and ( 92) into (93) yields where In terms of ( 86) and ( 94), the virtual power of inertial forces  ine can be written as where 3.4.Governing Equation.The virtual power equation of the beam element can be written as in which  ext =  q   F ext is the external virtual power and F ext is the generalized external force.Substituting the internal virtual power (75) and the inertial virtual power (96) into the virtual power equation (98) yields By substituting the modified generalized nodal force (82) into (99), the governing equations of beam element can be expressed as Thus, the generalized acceleration is where 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  1 K 1 q  is equivalent to the internal damping force.If q  represents the rigid motion only, this equivalent damping force will disappear since K 1 q  is equal to zero for the allowed rigid motion.
The element generalized velocity q  and the generalized coordinates p  are set as the state variables; that is, in which the relationship between the time derivative of p  and q  becomes The governing equations, which are a set of the second order ordinary differential equations, can be written as 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 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.

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.

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  = 2/.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.

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.

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 ℎ = 0, the frequency equation can be given by where  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 becomes where   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 by Substituting  2 = ℎ 2 /6 into (109) yields where  ℎ = 1/ℎ 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 √ 6 ℎ .
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.

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  =  = 10 6 ,  =  = 10 3 , 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 ℎ = 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 ℎ 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 ℎ = 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 ℎ.This is because the internal damping forces will become larger while increasing the equivalent damping coefficient  1 = ℎ/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 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 V, respectively.This problem has been also studied in previous researches [6,7,32,33] with the rotation (112) which 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 wellestablished nonstiff solver ode45.

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

Figure 1 :
Figure 1: Geometric configuration of the spatial beam.

Figure 2 :
Figure 2: Equilibrium of an infinitesimal element of a beam.

Figure 3 :
Figure 3: The kinematic description of the Euler-Bernoulli beam element.

Figure 4 :
Figure 4: The decomposition of beam section rotation.

Figure 6 :
Figure 6: Flowchart for avoiding the singularity of the rotation vector.

Figure 7 :
Figure 7: Cantilever beam with an end moment.

Figure 8 :
Figure 8: Comparison of the displacements of the cantilever beam with an end moment.

Figure 9 : 1 Figure 10 :
Figure 9: Large deformations of the cantilever beam with an end moment.

Figure 11 :Figure 12 :
Figure 11: Displacement curve of cantilever beam with an end point load.

Figure 16 :
Figure 16: Time evolution of the displacement at the tip.

Figure 17 :
Figure 17: A flexible beam on a rotating base.

Figure 18 :Figure 19 :
Figure 18: Time evolution of the axial displacement  at the free end.

Table 1 :
Cantilever beam with an end point load.

Table 2 :
Frequency of cantilever beam with different meshes.

Table 3 :
Modified frequency of cantilever beam with different time parameter ℎ.Figure 13: General shape of the mapping curve between   / ℎ and / ℎ . h

Table 4 :
The comparison of CPU time for the right-angle cantilever beam.