A two-dimensional shear deformable ANCF consistent rotation-based formulation beam element

It has been demonstrated recently that the absolute nodal coordinate formulation (ANCF) can be used to develop lower-order consistent rotation-based formulations (CRBFs) that employ finite rotation parameters as nodal coordinates without the need for interpolating the rotation field. The objective of this study is to develop new planar shear deformable ANCF/CRBF beam elements and demonstrate their use. A cubic ANCF/CRBF shear deformable beam element is first developed starting with the ANCF kinematic description that employs position vector gradients as nodal coordinates. The transverse position vector gradients at the nodes are expressed in terms of finite rotation parameters, leading to a lower-dimensional beam element model that captures the shear deformation, ensures continuity of the stresses and rotations at the nodes, allows for an arbitrary large displacement, and has a kinematic description consistent with geometry methods and suitable for systematically modeling curved structures. The results, obtained using the new planar ANCF/CRBF shear deformable beam element, are compared with the original and more general ANCF shear deformable beam element. Another lower-dimension bilinear CRBF beam element which has three coordinates at each node, two translation coordinates and one rotation parameter, is also developed in this investigation. The formulations of the three finite elements, including the ANCF finite element, considered in this investigation are compared. Numerical results are presented in order to demonstrate the use of the new formulations and test their performances in the analysis of large displacements and deformations. While the ANCF/CRBF assumptions are evaluated numerically, the results obtained show, in general, a good agreement between the elements considered in this study. The results also show that the CRBF finite elements, which have nonlinear mass matrix, can be more efficient for smaller meshes. As the mesh size and the number of finite elements increase, the original higher-order ANCF finite elements, which have constant mass matrix and zero Coriolis and centrifugal forces, become more efficient. The ANCF/CRBF approach clearly demonstrates that there is no need for introducing an independent rotation field to capture the shear effect.

ing curved structures. The results, obtained using the new planar ANCF/CRBF shear deformable beam element, are compared with the original and more general ANCF shear deformable beam element. Another lower-dimension bilinear CRBF beam element which has three coordinates at each node, two translation coordinates and one rotation parameter, is also developed in this investigation. The formulations of the three finite elements, including the ANCF finite element, considered in this investigation are compared. Numerical results are presented in order to demonstrate the use of the new formulations and test their performances in the analysis of large displacements and deformations. While the ANCF/CRBF assumptions are evaluated numerically, the results obtained show, in general, a good agreement between the elements considered in this study. The results also show that the CRBF finite elements, which have nonlinear mass matrix, can be more efficient for smaller meshes. As the mesh size and the number of finite elements increase, the original higher-order ANCF finite elements, which have constant mass matrix and zero Coriolis and centrifugal forces, become more efficient. The ANCF/CRBF approach clearly demonstrates that there is no need for introducing an independent rotation field to capture the shear effect.

Introduction
In multibody system (MBS) dynamics, proper treatment of finite rotations and constraints that define mechanical joints and specified motion trajectories is necessary. MBS algorithms, in which joint constraints are defined using nonlinear algebraic equations, are designed to solve a coupled system of differentialalgebraic equations (DAEs). In this section, some important issues related to the description of the finite rotations, formulation of the stress forces, and mode of deformations are discussed. The scope of this investigation and its relevance to the important subject of the integration of computer-aided design (CAD) and analysis are also discussed in this section.

Non-commutative finite rotations
In continuum mechanics, the rotation field is not independent from the position field, and as a consequence, the general continuum mechanics description employs only a position field. Large rotation vector formulations (LRVFs), on the other hand, employ two independent interpolations: one interpolation for the position field and the other interpolation for the rotation field [1]. Independent interpolation of the finite rotations is not consistent with the general continuum mechanics description, violates basic MBS dynamics principles, leads to a redundancy problem, and is the source of several fundamental and numerical problems as discussed in the literature [2]. For instance, finite rotations cannot be treated as vectors, and therefore, vector additions cannot be applied including interpolation which implies vector additions. While a rotation about a fixed single axis is commutative, general three-dimensional rotation is not commutative [3][4][5]. In Fig. 1a, a block is first rotated 90 • about the Y -axis followed by 90 • rotation about the Z -axis. In Fig. 1b, the same rotations in reverse order are performed, that is, the block is first rotated 90 • about the Z -axis followed by 90 • rotation about the Y -axis. It is clear from the results presented in Fig. 1 that a change in the sequence of rotations leads to different final configurations of the block, demonstrating that finite rotations in the spatial analysis are not commutative and cannot be in general added, treated as vectors, or interpolated. This comment also applies to the four Euler parameters which are quaternions and their mathematical operations are governed by the rules of quaternion algebra.

Geometrically exact formulations
It is also important to distinguish between large rotation vector formulations (LRVFs) and what is called geo- metrically exact beam formulations (GEBFs) [6][7][8][9][10]. A MBS approach is defined by the kinematic description and type of coordinates used in the formulation of the dynamic equations of motion. For example, the general continuum mechanics approach employs the matrix of position vector gradients to formulate the stress forces and the floating frame of reference (FFR) which ensures zero strains for an arbitrary rigid body displacement are geometrically exact approaches. For this reason, the authors of this paper refer to methods that interpolate finite rotations as LRVFs. With any of the MBS approaches, different methods that lead to zero strain under an arbitrary orthogonal rigid body transformation can be used in the formulation of the stress (elastic) forces. The geometrically exact beam formulation (GEBF) is one of these methods for formulating the elastic forces, and such a formulation can be used with ANCF finite elements or with LRVF finite elements as demonstrated in the literature [1,11].

Independent deformation modes and locking problems
Shear and torsion represent important independent modes of deformation that must be properly represented in the mathematical model. Their effect in large displacement and large deformation problems can be systematically captured using ANCF finite elements, which have been successfully used in many applications [11][12][13][14][15][16][17][18][19][20][21][22][23][24][25]. Capturing shear deformations does not require the use of two independent meshes and does not require interpolation of rotations. Starting with the ANCF kinematic description, a consistent rotationbased formulation (CRBF) can still be developed using one mesh that defines a unique rotation field and captures the effect of shear deformation [26]. The ANCF position vector gradients can be expressed in terms of rotation parameters. Using this coordinate transformation, a velocity transformation matrix can be developed and used to define the finite element equations of motion in terms of finite rotation parameters. In the general three-dimensional case, the ANCF/CRBF approach correctly captures the effect of shear and torsion.
Fully parameterized ANCF finite elements as well as fully parameterized conventional finite elements such as the brick and tetrahedral elements have been widely used and have proven to be very valuable in captur-ing effects that cannot be captured by gradient deficient finite elements. As evident by the large number of investigations that span more than four decades, conventional finite elements, including the brick and tetrahedral elements, can suffer from serious locking problems that need to be addressed. While locking has been an important issue in the finite element literature, it has not deterred the scientific community from exploiting the unique features of fully parameterized elements. As demonstrated in this paper, some variables can converge faster as compared to others. For example, a small number of elements is sufficient to predict displacements and deformations, while a larger number of elements is required in order to achieve convergence for some strain components including shear.
1.4 Future mechanics-based CAD/analysis systems ANCF/CRBF finite elements can also be the basis for developing new floating frame of reference (FFR) finite elements that can be converted from and to Bspline and non-uniform rational B-splines (NURBS) CAD geometry representations. Such ANCF/CRBF finite elements can be designed to have numbers and types of nodal coordinates similar to the numbers and types of nodal coordinates used in existing FFR formulations and computer algorithms. This will facilitate the development of the powerful mechanics-based CAD/analysis systems envisioned for future MBS simulations since CAD geometry can be systematically converted to FFR meshes. Therefore, further research in this direction will result in an FFR implementation that can be an important component in the integration of computer-aided design and analysis (I-CAD-A).
It is also important to point out that existing finite and infinitesimal rotation finite element formulations do not, in general, ensure the continuity of the stress field at the nodal points. The continuity of the rotation field does not imply continuity of the stress field. The ANCF/CRBF used in this investigation allows for consistently using rotations as nodal coordinates and at the same time ensuring the continuity of the stress field. This important ANCF/CRBF feature is the result of using the kinematic relationships between the ANCF position vector gradients and the rotation parameters.

Scope and contributions of this investigation
In this investigation, the first CRBF elements are developed, and their use is demonstrated using numerical examples. A distinction is made between the ANCF/CRBF finite elements and other CRBF finite elements. ANCF/CRBF ensures the continuity of the rotation and stress fields, while CRBF ensures only the rotation field as the result of using lower order of interpolation. The ANCF kinematic description of 12-coordinate planar beam element is first used [12]. Crucial in the development of the ANCF/CRBF finite elements is the use of the ANCF position gradient vectors which cannot be zero vectors. Displacement gradient vectors, on the other hand, can be zero vectors. In the planar formulation presented in this paper, the ANCF transverse position gradient vector at the node is expressed in terms of a finite rotation parameter, leading to an element with a lower dimension. In the case of shear non-deformable element, one can show that this element can be reduced to an element which has three or four coordinates per node, depending on whether or not extensibility is accounted for. Another CRBF shear deformable element which has three coordinates per node based on bilinear interpolation is also proposed in this investigation. This bilinear element does not ensure the continuity of the longitudinal gradient vector, and therefore, it will be referred to as CRBF finite element only. The results obtained using the elements considered in this investigation are compared, and the comparative study shows that as the mesh size increases, the original ANCF finite elements can become more efficient due to the fact that these elements lead to a constant mass matrix and zero Coriolis and centrifugal forces, while the ANCF/CRBF finite elements lead to a nonlinear mass matrix. Nonetheless, the new elements proposed in this investigation can be effectively and efficiently used in static applications, in dynamics problems which have small to moderate size meshes, and can be the foundation for a new FFR implementation that is consistent with future mechanics-based I-CAD-A.

ANCF/CRBF finite elements
A simple rigid body motion cannot be described using kinematic equations that are linear in the rotation parameters. Therefore, the position equations of a finite element whose geometry is invariant under rigid body transformation cannot be linear functions of the rotation parameters. Nonetheless, the element velocity equations can be linear functions in the time derivatives of the rotation parameters. This fact will be used in this paper to develop a consistent formulation that employs rotations as nodal coordinates without violating basic principles of dynamics and continuum mechanics. In this section, the equations required for developing the planar ANCF/CRBF finite beam elements are briefly discussed. These equations can be considered as a special case of the three-dimensional equations recently introduced [26]. The assumed displacement field of planar ANCF finite elements can be written as r ( where r is the global position vector of an arbitrary point on the element, S is the element shape function matrix, e is the vector of nodal coordinates, x = x y T is the vector of the element spatial coordinates, and t is time. At a given node of a planar fully parameterized ANCF finite element, the vector of nodal coordinates can be defined using the position and gradient coordinate vectors r, r x , and r y , where r x = ∂r/∂ x and r y = ∂r/∂ y. In the ANCF description, the position gradient vectors r x and r y , which are not in general orthogonal unit vectors, define the strain components at the nodes. In the planar analysis, an ANCF/CRBF finite element can be developed by considering the gradient vector r y as unit vector expressed in terms of a rotation θ (t) as r y = − sin θ cos θ T . It follows thaṫ In this equation, I is the 2 × 2 identity matrix, and a = − cos θ sin θ T . By using the preceding equation, the time derivative of the vector of element nodal coordinates can be written in terms of the time derivatives of the position coordinates and rotation parameters asė = Bṗ, where p is the vector of nodal coordinates that include position coordinates, longitudinal gradient vectors, and rotation parameters, and B is the element velocity transformation matrix that depends nonlinearly on the independent rotation parameters [26]. More details on the formulation of the velocity transformation matrix B will be provided in the following sections.
The ANCF finite element equations of motion can be written as Më = Q, where M is the constant symmetric ANCF mass matrix and Q is the vector of nodal forces including the elastic forces [22]. Using the velocity transformationė = Bṗ, one can writeë = Bp+Ḃṗ. Substituting this acceleration equation into the element equation of motion Më = Q and pre-multiplying by the transpose of the velocity transformation matrix B, one obtains the ANCF/CRBF element equations associated with the new set of coordinates that include finite rotation coordinates asMp =Q, whereM = B T MB andQ = B T Q − MḂṗ . It is clear that while the ANCF mass matrix M is constant, the ANCF/CRBF mass matrixM is nonlinear, and as a consequence, quadratic velocity inertia forces appear explicitly in the equations of motion.
In order to describe stress-free initial curved geometry, the conditions that define the position gradient vectors in terms of the rotation parameters can be developed after the mesh assembly. In this case, ANCF finite elements can be used first to develop the curved stressfree geometry using the position vector gradients. This is particularly important in the integration of CAD and analysis since advantage can be taken of the position vector gradients in defining the accurate geometry of the system components.

ANCF/CRBF shear deformable beam
In this section, a new ANCF/CRBF planar beam element is developed starting with the assumed displacement field of the planar ANCF shear deformable beam element [12]. The new ANCF/CRBF finite beam element correctly describes an arbitrary rigid body motion including arbitrary finite rotations, ensures continuity of the rotation and stresses at the nodal points, and captures shear deformation.

Displacement field
The ANCF assumed displacement field r (x, t) = S (x) e (t) is used as the starting point, where r is the global position vector of an arbitrary point on the element, S is the element shape function matrix, and e is the vector of nodal coordinates that includes position vector gradients as previously explained, x = x y T is the vector of the element spatial coordi-nates, and t is time. In the case of the two-dimensional shear deformable beam element, the displacement field r is defined in the global coordinate system as where a i and b i are the polynomial coefficients, and x and y are the spatial coordinates defined in a beam coordinate system. The spatial coordinate x is chosen to be along the beam axis (0 ≤ x ≤ l), where l is the element length. The vector of the element nodal coordinates e can be defined using the position and gradient coordinates r, r x , and r y , where r x = ∂r/∂ x and r y = ∂r/∂ y as where the functions s i = s i (ξ, η) are defined as and ξ = x/l, η = y/l, and l is the element length.

Finite rotations and velocity transformation
In order to obtain the ANCF/CRBF shear deformable beam element using the ANCF description, the transverse gradient vectors at the nodes are assumed to be unit vectors. This can be accomplished using the following definitions: where θ 1 and θ 2 are the rotation angles at the nodes. The preceding equations, upon differentiation with respect to time, lead to the linear velocity relationṙ 1 y = a 1θ1 , andṙ 2 y = a 2θ2 , where a 1 = − cos θ 1 sin θ 1 T , and a 2 = − cos θ 2 sin θ 2 T . Using these kinematic relationships, the time derivatives of the vector of element nodal coordinates can be written as In this equation, I is the 2 × 2 identity matrix and The velocity transformation matrix B can be used to write the element equations in terms of the new set of nodal coordinates p asMp =Q, as previously explained. The element equations of motionMp =Q can be solved for the acceleration vectorp which can be integrated to determine the coordinate vector p which includes the rotation coordinates and the longitudinal position vector gradients. The rotation coordinates can be used to determine the transverse position gradient vectors at the nodes. The use of the procedure described in this section reduces the number of nodal coordinates from 6 coordinates to 5 coordinates at the node. The resulting element ensures the continuity of the rotations and stresses at the nodal points. The resulting element mass matrixM = B T MB can be calculated using the mass matrix M of the original ANCF shear deformable beam element given explicitly in [12] and the velocity transformation matrix B defined in the preceding equation.

Formulation of the stress forces
Having determined the rotations at the nodes, the longitudinal and transverse position gradient vectors can then be used in a geometrically exact approach to evaluate the element elastic forces that enter into the formulation of the equations of motion. Other elastic force formulations based on classical beam theories can also be used with the proposed ANCF/CRBF finite element.
In this investigation, the general continuum mechanics approach for formulating the elastic forces is used. In this approach, the matrix of position vector gradient J = r x r y is used to define the Green-Lagrange strain tensor ε, which in turn is used to define the second Piola-Kirchhoff stress tensor σ using a linear Hookean material model. The virtual work of the elastic forces δW s = − V σ : δεdV , where V is the element volume, is used to define the generalized elastic forces associated with coordinates of the elements including the finite rotation coordinates at the node. It is important to point out that the same approach for formulating the elastic forces is used when the bilinear CRBF finite element is introduced in a latter section of this paper.

Comparison with the analytical solution
In order to compare the results obtained using the new higher-order ANCF/CRBF element introduced in this section against the analytical solution, a simple cantilever beam problem is considered. The elastic line approach is used instead of the general continuum mechanics approach for formulating the elastic forces in order to get closer to the assumptions used in the simplified analytical approach in which the deflection of the tip point is defined as δ = P/ 3E I /l 3 , where P is the load, E is the modulus of elasticity, I is the second moment of area, and l is the length of the beam. Since the simplified approach is based on a linear theory, the properties used are selected to ensure small deformation. The length of the beam is assumed to be 1 m, the load P is assumed to be 10 N, the modulus of elasticity E and Poisson's ratio are assumed to be 2 × 10 11 N/m 2 and 0.3, respectively, and the rectangular cross section is assumed to have length and height of 0.01 m. Neglecting the effect of gravity, the analytical solution is δ = 1.9908 × 10 −2 and the solution predicted using 12 ANCF/CRBF elements is δ = 1.990 × 10 −2 , showing a relatively good agreement despite the fundamental formulation differences.

Discussion
For the finite elements presented in this section and the following section, the transverse gradient vector is interpolated linearly, that is, Since a 1 and a 2 remain unit vectors, it follows that the square of the norm of the transverse gradient vector r y within the finite element is given by r T y r y = 1 + 2 (1 − ξ ) ξ a T 1 a 2 − 1 . This equation shows that the norm of the transverse vector r y does not remain equal to one, and the maximum deviation from a unit vector occurs at ξ = 0.5. Simple calculations can show that this deviation from unity remains very small as long as the relative rotation of the unit vector a 2 with respect to the unit vector a 1 does not exceed 30 • [26].
A shear-non-deformable beam element that has 3 coordinates per node can be also developed using the ANCF assumed displacement field presented in this section. In this case, the two gradient vectors at the node must be expressed in terms of the rotation parameter. To this end, one can use the following two conditions at node k: r k x = cos θ k sin θ k T , r k y = − sin θ k cos θ k T . Applying these two conditions at each node leads to ANCF/CRBF finite element that has three coordinates per node, two translations and one rotation. In this case, since no extensibility parameter is used with the longitudinal gradient vector r x , the use of the general continuum mechanics approach to formulate the elastic forces will lead to zero stresses at the nodal points. This is consistent with the Euler-Bernoulli beam theory in which the centerline of the beam is assumed to be stress free. An ANCF/CRBF shear non-deformable finite element can be used to shed light on the assumptions used in Euler-Bernoulli beam theory and the role of the position vector gradients and interpolation in defining the classical beam theory linear stress distribution along the cross section. That is, classical beam formulations can be systematically obtained from the general continuum mechanics strain and stress formulations using the ANCF/CRBF finite elements.
In the case of shear deformable beam elements, and when higher order of interpolation is used, the vector p can include the gradient vectors r x at the nodes. If a lower-order bilinear element is used, the inclusion of this position gradient vector is not necessary [1]. Using a lower order of interpolation, another new CRBF shear deformable finite element that has three coordinates per node can be developed as described in the following section. This element, however, will not be referred to as ANCF finite element because of the low order of interpolation.

Lower-order bilinear CRBF finite element
If the continuity at the element nodes is not required for r x , lower-order polynomials in x can be used to reduce the number of coordinates to 3 coordinates per node. Linear polynomials were used for both the position and rotation interpolations in the literature [1]. In such a case of lower-order element, r x is not considered as a coordinate vector, and one can only use r and r y as nodal coordinates. The displacement field r of the bilinear element considered in this section is defined in the global coordinate system as In this case of lower-order shear deformable beam element, the vector of nodal coordinates is defined as Using the assumed displacement field and the element nodal coordinates, the element shape function can be defined as where the functions s i = s i (ξ, η) are defined as and ξ = x/l, η = y/l, and l is the element length.
Using the velocity transformation matrix B, the finite element equations can be obtained as previously described asMp =Q. This equation can be integrated numerically to determine the coordinate vector p which includes the rotation parameters. Knowing the rotations, the transverse gradient vectors at the nodes can be determined and used in the formulation of the elastic forces.
The bilinear element discussed in this section ensures the continuity of the rotation field. It does not, however, ensure the continuity of the longitudinal gradient vector r x at the nodal points. The use of g X Y Fig. 2 Free falling flexible pendulum the low order of interpolation from the outset does not allow imposing continuity on r and r x at the same time. This element also reduces to a truss element at the element centerline, and its use in accurately modeling curved geometry may require the use of very fine mesh. Nonetheless, one can develop a higher-order element that has three coordinates per node (two translations and one rotation) by using internal nodes.

Numerical examples
In this section, a numerical example that has results reported in the literature is considered in order to evaluate the performance of the two proposed CRBF shear deformable finite beam elements introduced in this investigation. The example considered, shown in Fig. 2, is the free falling of a flexible pendulum under its own weight [12]. As shown in the figure, the beam is connected to the ground by a pin joint at one end and is assumed to be initially horizontal. The beam has length of 1.2 m, cross-sectional area of 0.0016 m 2 , second moment area of 8.533 × 10 −6 m 4 , a mass density of 5540 kg/m 3 , Poisson's ratio of 0.3, and a modulus of elasticity of 0.70 × 10 6 Pa. The gravity constant is assumed to be 9.81 m/s 2 . Three element types are considered in the analysis of the falling pendulum. The first type is the planar shear deformable ANCF finite element [12], the second type is the higher-order ANCF/CRBF element, and the third type is the lowerorder bilinear CRBF element. The simulations of the beam with the three element types are performed using different numbers of elements in order to check the convergence of the solutions of different models. Since different variables (displacements and strain variables for example) converge at different rates when different elements are used, a number of elements that achieves convergence for the variable examined is used. For this reason, some results are reported using 12 elements and some other results are reported using 100 elements. The CPU times for the simulations performed are shown in Table 1. The CPU times are reported for the models with numbers of elements that achieve convergence. For the third element type (lower order), because the continuity at the element nodes is not required for r x , and lower-order polynomials in x reduce the number of coordinates to 3 coordinates per node only, significant saving in computer time was achieved as compared to the first two element types in the case of a small mesh size. It was observed that the new element is about five times faster than the general ANCF element when the mesh size is not very large. As the mesh size and number of elements increase, the original ANCF finite element becomes efficient since the mass matrix is constant. This fact is particularly important when developing complex MBS models that may include thousands of finite elements. The CPU times reported in Table 1 are obtained using a MATLAB code and sequential PC computations. Figure 3 shows the position of the tip point of the beam using 6 and 12 higher-order ANCF/CRBF elements when the gravity constant is equal to 9.81 m/s 2 . It is clear from the results presented in this figure that there is a good agreement between the two models. The  results demonstrate that the solution converges with small number of elements. Figure 4 shows a comparison between the vertical displacements of the beam tip point obtained using the higher-order ANCF/CRBF element and the ANCF finite element previously presented in the literature [12]. The results presented in this figure show a good agreement between the two different models. Figure 5 compares the magnitudes of the transverse position gradient vector r y at the midpoint of the last element as predicted by the higher-order ANCF/CRBF element model and the ANCF finite element model. In the case of the ANCF/CRBF model, because of imposing the unit vector condition on the gradient vector r y at the nodes, the magnitude of r y at the nodes is equal to 1, but it is not in general equal to 1 inside the element. From the results presented in Fig. 5, one can observe that the magnitude of r y at the midpoint of Element 100 as predicted by the ANCF/CRBF model is very close to 1, and therefore, there is no significant violation of this condition as the result of the interpolation. Figure 6 shows the shear angle at the midpoint of the last element as predicted by the ANCF/CRBF and ANCF finite element models. The results presented in Fig. 6 show that there is again a good agreement between the solutions obtained using the two different models. Figure 7 shows the shear strain as predicted by the higher-order ANCF/CRBF and ANCF finite element models. The results presented in Fig. 7 show again a good agreement.

Bilinear CRBF element
In this section, the performance of the CRBF finite element developed using a lower-order bilinear polynomial is examined. This element will not be referred to as an ANCF finite element because of the low interpolation order. Figure 8 shows the position of the tip point of the beam using 6, 12, and 24 lower-order bilinear CRBF elements when the gravity constant is equal to 9.81 m/s 2 . It is clear from the results presented in this figure that there is a good agreement between 12-element and 24-element solutions. These results demonstrate that this model requires more elements to converge than the other two higher-order elements, and the solution for this element converges with 12 elements. Figure 9 shows the vertical displacement of the beam tip point as predicted by the bilinear CRBF element model and the higher-order ANCF/CRBF element model. The results presented in Fig. 9 show a good agreement between the two models. Figure 10 shows the magnitude of r y at the midpoint of the last element (Element 12) as predicted using the bilinear CRBF element and higher-order ANCF/CRBF element models. In these two models, the magnitude of r y at the nodes is equal to 1, but it is not necessarily equal to 1 inside the finite element. The results show that, in general, the magnitude of r y at the midpoint of Element 12 is close to 1, demonstrating that there is no significant violation of the constraints inside the element as the result of the interpolation. Figure 11 shows the shear angle at the midpoint of the last element (Element 100) predicted using the bilinear CRBF element and higher-order ANCF/CRBF element models. The results presented in this figure show a good agreement between the solutions obtained using the two CRBF models. Figure 12 shows a comparison of the shear strains at the midpoint of the last element obtained using the two CRBF models. Figure 12 shows that the variation trends of the two curves are similar, and the biggest difference is about 0.0005. In order to explain the results of the two models, one can examine the magnitude and orientation of the position gradient vector r x predicted by the two models. Figure 13 shows the magnitude of r x predicted by the two models, while Fig. 14 shows the orientation of r x . Figures 13 and 14 show a good agreement between the two models. Figure 15 shows the orientation of r y at the midpoint of Element 100. The results presented in this figure show also a good agreement between the solutions obtained using the two models. Using the ANCF/CRBF approach, large rotations due to deformation can still be described using a relatively small number of elements. The use of an algebraic equation to enforce the condition r y = 1 at the nodes does not lead to significant violation of this condition inside the finite element. While this condition leads to the elimination of a degree of freedom, the shear is not necessarily constant inside the element since r x can change depending on the order of interpolation in x, and consequently, the shear strain (1/2) r x · r y can vary from one point to another within the element. The ANCF/CRBF approach clearly demonstrates that there is no need for introducing an independent rotation field to capture the shear effect. Such an CRBF approach can lead to reasonably accurate results without violating basic dynamics and continuum mechanics principles. The analysis presented in this paper also shows that the interpolation of the gradient vectors is not equivalent to the interpolation of rotations. It is important also to note that the gradient vectors are not interpolated independently.

Summary and conclusions
In the case of a simple finite rotation, the rigid body kinematic equations cannot be written as linear functions of the rotation parameters. The nonlinearity of the position equations is necessary in order to ensure that the geometry is invariant under rigid body coordinate transformation. The velocity equations, on the other hand, are linear in the derivatives. This fact is used in this investigation to develop new shear deformable finite elements that employ rotations as nodal coordinates. The approach used in developing these new finite elements does not violate basic principles of dynamics and/or continuum mechanics and observes the non-commutativity nature of the finite rotations. The use of this approach clearly demonstrates that there is no need for the interpolation of finite rotations in order to capture shear deformations; finite rotations are in general non-commutative, and interpolation implies vector addition and commutativity. As also pointed out in this paper, it is important to distinguish between the MBS approach used and the method of formulating the elastic forces. Geometrically exact beam formulations (GEBFs) are methods for formulating the elastic forces and can be used with any approach that is defined by the kinematic description and the coordinates used; for instance, the general continuum mechanics approach for formulating the elastic forces is a geometrically exact formulation. For this reason, the authors refer to approaches that interpolate finite rotations as large rotation vector formulations (LRVFs). As it is known, a finite rotation about a fixed axis is commutative, an example of which is the planar case considered in this paper. Nonetheless, the approach presented in this paper can be generalized to three-dimensional finite elements as will be reported by the authors in future publications. The use of fixed axis of rotation, as it is the case in many rotorcraft and wind turbine examples reported in the literature, is not sufficient to justify using an approach as the basis for developing general MBS algorithms that are based on non-incremental non-commutative rotation procedures. In this investigation, an attempt is made to shed light on some of the fundamental issues discussed above by using a consistent rotation-based formulation (CRBF) to develop new planar finite elements that employ rotation parameters as nodal coordinates. A planar ANCF/CRBF shear deformable beam element is developed for the large rotation and large deformation analysis. The formulation defines a unique rotation field, employs one interpolation, captures shear deformations, does not suffer from the redundancy problem, allows for systematically describing curved geometry, and leads to elastic force definitions that eliminate highfrequency modes associated with the deformation of the element cross section. Another lower-order bilinear CRBF finite element was also developed, and the results obtained using the two new CRBF finite elements are compared with results published in the literature obtained using a planar ANCF shear deformable beam element that employs position vector gradients as nodal coordinates. The numerical comparative study presented in this paper demonstrates that the kinematic constraints imposed on the transverse gradient vectors at the nodes are not significantly violated within the finite elements. The results also show that as the mesh size and number of elements increase, the original ANCF finite element that employ position vector gradients as nodal coordinates can be more efficient since such an element leads to a constant mass matrix and zero Coriolis and centrifugal forces. This conclusion is particularly important in MBS applications that require large meshes that may include thousands of finite elements.
Future investigations will be focused on developing three-dimensional beam, plate, and shell elements using the ANCF/CRBF approach. While ANCF finite elements are more general and more powerful for MBS implementation and can be much more efficient when large meshes are considered, ANCF/CRBF finite elements can be of great value in developing new beam, plate, and shell models for static applications, lowerdimension dynamic problems, and new floating frame of reference (FFR) finite elements that can be converted from and to B-spline and non-uniform rational B-splines (NURBS) CAD geometry. ANCF/CRBF finite elements can be developed to have numbers and types of nodal coordinates similar to the numbers and types of nodal coordinates used by the finite elements employed in existing FFR formulations and algorithms. This, as previously mentioned, will facilitate the development of mechanics-based CAD/analysis systems envisioned for future MBS simulations.