On a virtual work consistent three-dimensional Reissner–Simo beam formulation using the quaternion algebra

In the paper, we present the Reissner–Simo beam theory in which the rotations are represented by quaternions. From the generalized virtual work principle, where the unity constraint of the rotational quaternion is properly considered and the consistent energy complements of the rotational quaternions are employed, we derive the weak kinematic equations in the quaternion-based description. These equations are then employed in the extended virtual work principle to obtain the consistent governing equations of the three-dimensional beam in terms of the quaternion algebra. The quaternion moment equilibrium equation is analyzed, discussed, and interpreted. In numerical implementation, the standard Galerkin discretization is used to obtain the quaternion-based finite-element formulation. Various examples prove the suitability of the formulation.


Introduction
Three-dimensional rotations seem to be crucial and unavoidable in many engineering fields. Several mathematical models have been developed to describe the physical rotations in mathematical terms. Among them, the Euler angles, the rotation matrix, and the rotational vector seem to be most popular. In some engineering fields, an alternative description with rotational quaternions is also well recognized. The application of quaternions can be found in diverse fields such as the navigation in aerospace technologies [18,20,31], robotics [8,39], animating rotations in computer graphics [9,30], dynamics of rigid bodies [14,37], and flexible beams [15,17]. From the above-mentioned publications, it is clear that the use of the quaternion algebra brings several advantages when dealing with the rotations in three dimensions.
The rotational quaternions have only recently been studied extensively in the analysis of the threedimensional beams [10,12,15,45]. Our recent works in the quaternion-based beam formulations for static and dynamic analyses [44,45] introduced the quaternions by substituting pure quaternions for vector-like quantities, while the unity of rotational quaternions is preserved only by the update procedure. Such an approach is efficient and proves that the quaternion algebra is a suitable tool for describing the three-dimensional rotations in the beam formulations even though the formulations in [44,45] are non-standard and differ considerably from the standard finite-element approaches such as [12,32]. The choice of the parametrization of rotations has a strong impact on both the formulation of the problem and the numerical implementation, see, e.g., [5,6,11,27,43]. It should be noted that, in the field of nonlinear mechanics, we should properly consider the linearization in multiplicative configuration spaces and the effects of the chosen parametrization on the energy principles. These two important aspects of the rotational quaternion-based approach will be discussed in the present paper.
The finite-element formulation of the three-dimensional beam employing the quaternion algebra has already been considered by Bottasso [7]. His paper is mainly focused on the choice of the appropriate constraining equation for the quaternion parameters when used in the generalized virtual work principle. Our paper extends the work of Bottasso [7] in the way that the virtual work principle becomes consistent with the quaternion unknowns. The similar problem of energy-consistent moments has been discussed by Ritto-Corrêa and Camotim [26] when using various rotational parameters. They distinguish between the rotational moment, which is the energy complement to the classical variation of the rotational vector, and several other moments, belonging to different variations of other vector-like parametrizations of rotations. Their concepts will here be applied in the four-parameter quaternion description.
In the present work, the quaternion-based governing equations of the beam are derived through the energyconsistent approach. We present the consistent generalized virtual work principle which is then used for the derivation of the weak kinematic equations in the quaternion description. Next, we derive the equilibrium equations of the beam consisting of the standard equilibrium vector equation for forces and the modified equilibrium vector equation for moments. This is in a sharp contrast to the direct approach used in [45] where the standard equilibrium equations for moments were preserved. The quaternion moment as the energy complement to the rotational quaternion is defined in accordance with [26]. The present approach introduces an additional unknown of the problem by enforcing the weak form of the unity constraint to be satisfied in accord with the method of Lagrangian multipliers. We thus obtain different continuous governing equations as in [45]. The differences are carefully analyzed and the equivalence of the continuous governing equations is shown. Based on the present form of the governing equations, a numerical implementation is proposed. In contrast to [45], the numerical implementation is carried out in a standard way. The Galerkin discretization of the governing equations is made and the finite-element quaternion-based formulation is presented. A speciality of the present approach is in the treatment of the Lagrangian multiplier as the additional unknown of the discretized system, which means that the weak form of the unity constraint becomes a part equal member of the governing equations. It is important to observe that a larger number of degrees of freedom in the present approach does not increase the complexity of the numerical method. On the contrary, the finite-element formulation is rather simple and very easy to implement. Various examples prove the suitability and the efficiency of the formulation.

Geometry of the three-dimensional beam
The three-dimensional beam is modeled by the position vector r (x) of the line of centroids (the beam axis) and by the family of orthonormal base vectors where G 2 and G 3 are spanning the planes of cross-sections A (x), and G 1 = G 2 × G 3 , see Fig. 1. Together with the reference point x on the beam axis, the basis defines the local coordinate system (ξ, η, ζ ). In general, the beam quantities depend on the arc-length x ∈ [0, L] of the line of centroids in the initial configuration. We further introduce a fixed orthonormal triad B g = g 1 , g 2 , g 3 which, together with the reference point O, defines the global coordinate system (X, Y, Z ). We assume that the cross-sections are rigid, but not necessarily perpendicular to the deformed beam axis. The relation between the local and the global base vectors is determined by the rotations.

Rotations and quaternions
The rotations and their parametrizations are well known and a widely studied subject in structural mechanics, see, e.g., [1,2,25]. The rotation R ϑ, n is a linear operator in IR 3 which rotates vectors for an angle ϑ about the axis of rotation described with a unit vector n . In the matrix representation and with respect to an arbitrary basis, the rotation of vector x into vector y takes the following form: where R is the matrix representation of the rotational operator R, and n is a (one-column) matrix representation of vector n . Furthermore, the rotation matrix also determines the transformation between the local and the global matrix descriptions of the same vector: where x and x G represent the description of vector x with respect to the global and the local bases, respectively, and the superscript T denotes the transpose of a matrix. Note that the rotational operator R has an equal matrix representation R in the two bases. The comprehensive presentation of the quaternion algebra can be found in textbook [38]. A detailed application of the quaternion representation of rotations in structural mechanics was presented in Zupan et al. [43,45]. The quaternion algebra (IH, +, ·, •), for IH = a = a 0 + a , a 0 ∈ IR, a ∈ IR 3 , is a four-dimensional linear space IR 4 , +, · with an additional non-commutative binary operation called the quaternion multiplication. a 0 is called the scalar part and a the vector part of quaternion a. As IR 3 is a linear sub-space in IH , it is straightforward to extend any basis in the three-dimensional Euclidean space B g = { g 1 , g 2 , g 3 } into the basis of the quaternion space T be two arbitrary quaternions in the component matrix form. Then, the quaternion multiplication can also be expressed by the matrix products as and We further introduce the conjugation of a quaternion, a * , which only alters the sign of the vector part of the quaternion, a * = a 0 − a . In the quaternion description, Eqs. (1)-(2) take the form and where q represents the rotational quaternion that belongs to the same rotational operator R as matrix R from Eq. (1). The rotational quaternion is defined as: The three-dimensional vector quantities, x and y, are in (5)-(7) replaced with pure quaternions x = 0 + x and y = 0 + y. Note that the matrix representation of the left and the right quaternion multiplication when defined by the rotational quaternion, φ L q and φ R q , satisfy the orthogonality conditions, φ S q φ T S q = φ T S q φ S q = I, for S = R, L, and thus, each represents a rotation in the four-dimensional Euclidean space IR 4 . Note also from Eq. (8) that the norm of any rotational quaternion q is equal to one, which also implies It is found suitable to introduce the cross-vector product of pure quaternions. It is a pure quaternion with its vector part equal to the standard cross-vector product, where a = 0 + a and b = 0 + b. Next, we will need some properties of the scalar product and adjoint operators to rearrange the terms in the virtual work principle. In every finite dimensional linear space IR n (n ∈ IN ) with the scalar product defined as ..,n , the following holds true [28]: where F is the matrix representation of the linear operator F : IR n → IR n . It is important to observe that which follows directly from (3)-(4). After applying (12)-(13) to arbitrary quaternions a, b, and c, we have Note that the identities (14)- (15) are essential for the rearrangement of the extended virtual work principle in terms of quaternions.

Derivative and variation of a rotational quaternion
The derivative of the rotational quaternion and its relation to the rotational strain vector as well as its linearization and the relation between the variation of the rotational quaternion and the variation of the rotational vector were derived and discussed in [45]. Here, we present only the results needed in the sequel. Let rotation R map the global basis B g onto the local basis B G and let q be its rotational quaternion with respect to the global basis. Then, we have where κ and κ G are descriptions of κ = 0 + κ with respect to the global and the local bases, respectively, and κ is the rotational strain vector. Analogously, we can find that where δϑ denotes the variation of the rotational vector ϑ = ϑ n represented in the global basis. The derivative of the conjugated quaternion, q * , is obtained by the differentiation of Eq. (10) and is expressed as

The generalized virtual work principle
Following the approach of Reissner [23,24], we assume the principle of virtual work in the form where N G and M G are the stress-resultant force and moment vectors of the cross-section, and δγ G (x) rel and (δκ G (x)) rel are virtual strains with respect to the local basis; n (x) and m (x) are external distributed force and moment vectors per unit of the initial length; δr and δϑ are the virtual displacement vector and the virtual rotational vector, respectively, and F (0), P (0), F (L), P (L) are the external point forces and moments at the two boundaries, x = 0 and x = L , all written with respect to the global basis. For simplicity and clearness of the notation, the dependency on x will be dropped. The quantities at the two boundaries will be denoted by an upper index, e.g., P (L) = P L . Now we will replace the three-parameter rotational vector by the four-parameter rotational quaternion using Eq. (17). In replacing the rotational vector with quaternions, we must consider the fact that the four variations of the rotational quaternion components, δq i , i = 0, 1, 2, 3, are mutually dependent; their relation can be derived through the variation of Eq. (9). This relation is accounted for by the method of Lagrangian multipliers [36]. The constraining equation (9) is first varied and then multiplied by an arbitrary unknown scalar function λ (x), independent on the primary unknowns. The product is then integrated along the length of the beam L 0 2λδ q · q dx, and added to the virtual work principle (19). The presence of quaternions in the virtual work principle extends the configuration space of the problem into the fourth dimension, which requires that some of the three-dimensional quantities are extended to four dimensions as pure quaternions, e.g., M = 0 + M. Next we recall that m · δϑ = m · δ ϑ, where m = 0 + m, and, likewise, P · δϑ = P · δ ϑ, where P = 0 + P. Equation (19) can thus be further rearranged by using Eqs. (15) and (17), so that the generalized virtual work principle using the quaternion representation of rotation reads:

The kinematic equations
In the spirit of Reissner's works [23,24], we now derive the kinematic equations from Eq. (20), which assures us that the obtained relationships between the displacements, the quaternions, and the strains will be fully consistent with the virtual work principle. The well-known local equilibrium equations of the three-dimensional beam are first rewritten in the quaternion form: Here, [ . ] IR 3 indicates the vector part of a quaternion. The generalized virtual work principle given in (20) holds true for an arbitrary part of the beam, bounded by the cross-sections at x = x 1 and x = x 2 . Considering an arbitrary integration region, [x 1 , x 2 ], in Eq. (20), inserting n and m from Eqs. (21)- (22), applying the integration by parts, collecting the terms around force N G , moment M G , and Lagrangian multiplier λ, respectively, yields Recalling that N G , M G and λ are arbitrary and mutually independent quantities and invoking the fundamental lemma of calculus of variations that the functions at the independent quantities are zero yields the relations between the virtual strains, the virtual displacements, and the virtual rotational quaternions as Equations (24)-(25) represent the linearized kinematic equations in the quaternion description and are equal to those obtained directly by converting the kinematic equations of the three-dimensional Reissner beam to the quaternion form, see [45] for details. The third kinematic equation (26) relates the rotational quaternion with its variation and implies that q and its variation, δ q, are orthogonal. Note that Eq. (26) is an exact variation of Eq. (9). Equations (24)-(26) represent a system of seven differential-algebraic equations for seven kinematic unknowns, i.e., three components of the virtual displacement vector and four components of the virtual rotational quaternion. The integration of the linearized kinematic equations becomes nontrivial in three dimensions. In contrast to Reissner [23] who proposed a simplification of rotations, Simo [32] derived the exact strong form of kinematic equations. Here, the strong kinematic equations in terms of rotational quaternions are obtained directly from the weak form (24)-(25) employing a similar procedure as introduced in [32] for the integration of the kinematic equations when expressed in terms of the rotational vector. As the details of the derivation are given in [45], only the final result is displayed here: Here c G = 0 + c G and d G = 0 + d G denote the variational constants, which are uniquely determined from the known initial configuration of the beam. Note that the variational constants are generally functions of x.
In the simplest yet a most common case, when the beam is straight in the undeformed configuration, we have c G = G 1 and d G = 0.

The equilibrium equations
We insert the above derived kinematic equations (24)- (25) into Eq. (20) and transform the equilibrium forces and moments into the global basis description. Then the terms are rearranged with the help of (14)- (15), (18), and the integration by parts is used to obtain Equation (29) represents the generalized principle of virtual work in which all components of δ r and δ q are arbitrary and mutually independent unknown functions. Collecting the terms around δr and δ q and invoking the fundamental lemma of calculus of variations give the equilibrium equations on [0, L]: and the boundary conditions: Equation (30) coincides with the standard equilibrium equation for forces, which does not depend on the rotation. In contrast, the equilibrium equation (31) for moments differs from the standard one. If we use the notation M for the standard form of the moment equilibrium equation and insert (36) into (31), we obtain Now the difference between the two expressions becomes evident. The result shows how the quaternion-based moment equilibrium equation (31) is related to the one based on the rotational vector.
Equations (30)-(31) represent a system of seven differential equations for seven unknown functions-three components of forces, three components of moments, and the Lagrangian multiplier. Equations (24)-(26) and (30)-(31) represent a set of fourteen equations of the geometrically exact beams, if described with the quaternion-based representation of rotations.

The quaternion moment and the interpretation of the Lagrangian multiplier
From (20) we can observe that the virtual work of the external distributed moments Thus, the definition of the energy-consistent quaternion moment in the sense of the work by Ritto-Corrêa and Camotim [26] is evident. Let M ϑ be the rotational moment as defined by Ritto-Corrêa and Camotim [26]. Then, the quaternion moment M q can be defined as where M ϑ = 0+ M ϑ . Please observe that the quaternion moment M q is configuration dependent, as it depends on q. It is conservative when M ϑ is conservative. It is interesting to explain the meaning of the Lagrangian multiplier λ. From Eq. (37), we can see that λ is the scalar component of the standard moment (vector) equation representing the proper extension of the moment equilibrium equation to the quaternion space. Equation (37) is now multiplied by q * using right quaternion multiplication, which yields Equation (37) is equivalent to the system of equations (39)- (40) and (9). From Eq. (39), it is obvious that λ = 0 is an exact analytical solution of Eq. (37), and it could be eliminated from the quaternion moment equilibrium equation (31) or (37) after the unity requirement of the rotational quaternion has been satisfied. In this case, Eq. (31) would be replaced by the standard moment equilibrium equation (40). This kind of approach has been followed in our previous work, see [45]. Note that Eq. (39) holds true only for a continuous system where Eq. (9) is identically satisfied for any x. In the numerical solution method, this is not the case and it is found suitable to consider λ as an additional unknown quantity, as we show next.

The discrete governing equations
The Galerkin finite-element approach is applied to obtain the discrete equilibrium equations fully consistent with the quaternion representation of rotations. The quaternion-based formulation introduces seven scalar kinematic unknowns, three components of the displacement vector and four rotational quaternion components. To consider the Lagrangian multiplier, λ, as independent primary unknown, we further extend the generalized virtual work principle (20) by adding the term Such a generalization of the virtual work principle is fully equivalent to the original one, (20), when dealing with the continuous unknowns, see [36, pp. 129-130].
The variations of the primary unknowns are interpolated in a standard manner: where x i ∈ [0, L], i = 0, 1, 2, . . . , N + 1, with x 0 = 0 and x N +1 = L, are the discretization points and P i (x) are interpolation functions. Note that an equal order of the interpolation has also been used for the Lagrangian multiplier. Such a choice needs some comments. It is important to observe the special role of the Lagrangian multiplier in the present formulation. The derivative of λ with respect to x does not occur in any of the governing equations. Numerical studies revealed that the iteration procedure diverges if the order of interpolation for δλ is lower than for δ q. If we take the order of the interpolation for δλ higher than for δ q, the tangent stiffness matrix becomes singular. Similarly, we observe severe singularity problems, if we release the continuity condition for λ. Only for the same order of the interpolation and the full continuity conditions imposed at the nodes the method becomes numerically stable and computationally efficient. Later on it will also be shown that the present approach performs very well even for high-order interpolations. When inserting the kinematic equations (24)-(25) into the extended principle of virtual work (20), applying (42)-(44) and after some rearrangement of the terms, we obtain the following discretized generalized virtual work principle: Note that the same result can also be obtained by adding the integrated unit length constraint (41) to Eq. (29). In (45) δr i , δ q i and δλ i are independent and arbitrary unknown discrete values. From the fundamental lemma of calculus of variations, we have: for i = 0, 1, . . . , N + 1. Here, Equations (46)-(48) represent a system of 8 (N + 2) discrete algebraic equations for as many discrete unknown quantities r i , q i and λ i , i = 0, 1, . . . , N + 1.

Numerical implementation
In the present numerical implementation, the Lagrangian polynomials are used for the interpolation. The integrals are evaluated numerically. We use a classical numerical integration where an integral is replaced by a weighted sum The values of the weights, w q , and the locations of the points, x q , are dependent on their number I N + 1 and on the chosen numerical integration method. Because of its high accuracy, the Gaussian integration is implemented. A typical property of standard displacement-based beam finite elements of a low order is the possible occurrence of shear locking. To avoid shear locking in such elements, several solutions have been proposed. One of the widely used is the reduced integration, see, e.g., [33] and the references therein. The quaternion approach seems to inherit this inconvenience of standard discretizations. In [45], the shear locking phenomenon was resolved by imposing the consistency conditions (see also [40]). For simplicity, we will here employ a standard approach which uses the reduced integration (I N = N ) for the terms consisting of the extensional and shear-strain components for the linear (N = 0) and the quadratic (N = 1) elements. For higher-order elements (N ≥ 2), a full integration (I N = N + 1) is used in all the terms. Note that the solution employed in [45] is fully applicable in the present approach, yet it would demand a generalized virtual work principle augmented with two additional terms. See also [19,22] for further details on the shear locking phenomenon.
The solution of the system of nonlinear algebraic equations (46)-(48) is obtained by the Newton method. The method requires the linearization of the discrete nonlinear equations and an update of the vector of unknowns with their corrections δ y = (δr i , δ q i , δλ i ) must be made in each iteration. The corrections are obtained by the solution of the system of linear algebraic equations K [n] δ y = − f [n] , for n = 0, 1, 2, . . . , obtained by the linearization of Eqs. (46)-(48). Here, K [n] is the global tangent matrix, and f [n] is the residual vector in iteration n. In linear vector spaces, new, improved solutions are obtained by summing the current iterative value and its correction: For the displacements and the Lagrangian multiplier such an additive update is adequate: This kind of update is not admissible for the rotational quantities. The quaternion representation of rotations conserves the multiplicative nature of rotations, so that the resultant rotational quaternion must be obtained by the quaternion product of the two subsequent rotational quaternions. Having the correction of the rotational quaternion δ q, we must first obtain the multiplicative incremental rotational quaternion q that, subsequently, is combined with the current rotational quaternion q [n] . From Eq. (17), it follows that 2δ q • q * is a more suitable measure for the variation of rotations than only δ q. The term 2δ q • q * is also analogous to the standard definitions of rotational strain in statics and angular velocity in dynamics. The advantage of considering the more sophisticated measure for the variation of rotations is that its natural projection from the tangent space onto the set of rotational quaternions follows directly from Eq. (8).
The consistent update of the rotational quaternions hence takes the following form: The update procedures for the derivatives of the rotational quaternions need to consider that only the variation of rotational quaternion has been interpolated. Thus, we first obtain q by the differentiation of q with respect to x employing Eq. (49) and then express the incremental strain using Eq. (28). The update of the rotational strains can then be carried out by taking into account their additive nature. When written with respect to the suitably chosen bases, it takes the following form: Details are given in [45].
In the present paper, we limit our studies to the linear elastic material of the beam. A simple linear relation between stress resultants and strain vectors is thus assumed in the following form: E and G denote the elastic and the shear moduli of material; A 1 is the cross-sectional area; J 1 is the torsional inertial moment of the cross-section; A 2 and A 3 are the effective shear areas in the principal inertial directions Following the rules of variational calculus, Eq. (46) is linearized as In the linearization of Eq. (47), it is suitable to introduce the notations The linearization then takes the form: Finally, the linearization of Eq. (48) reads δF (8) The details regarding the partial derivatives in the above equations are presented in the Appendix. From (52)-(54), we can observe that the linearization is relatively simple, which results in its easier implementation in a computer code and a higher computational efficiency when compared to our previous quaternion-based formulation [45]. By studying the sparseness of the tangent stiffness matrix, we found that the additional unknown only slightly enhances the complexity of the linearized system of equations, while at the same time the component blocks of the tangent stiffness matrix are easier to evaluate. The sparseness of these additional blocks due to taking the Lagrangian multiplier as an additional unknown is also related to the singularity problems at the structural level, if we allow λ to be discontinuous over the boundary. By satisfying the continuity conditions for λ and by taking the same order of interpolation for all the interpolated unknowns, we completely avoid any such numerical problems.

Numerical examples
In order to demonstrate validity, efficiency and accuracy of the proposed numerical approach, we now discuss several standard numerical examples. The examples also show the capability of the present formulation to consider the beams at arbitrary deformed configuration. All the calculations were performed in the Matlab [34] computing environment. Note that the present algorithm exactly preserves the unity of rotational quaternions at the interpolation and integration points.

The shear locking test
In the present quaternion-based formulation, we alleviate shear locking in low-order elements in a standard way by using the reduced integration, see Sect. 5. For the sake of demonstration, we study a straight cantilever beam under a vertical tip load F = 1 at the free end. The geometric and material properties of the cantilever are: E = 10 7 , G = 10 13 , L = 1, t = 0.1, t is the width of the rectangular cross-section. Note that shear modulus G was taken exceedingly large as the problem of shear locking becomes severe by increasing G. On the other hand, in conventional displacementbased finite elements, locking progresses rapidly as the height of the cross-section decreases. The height, h, of the beam was therefore varied from 0.1 (small) to 10 (large). Our numerical results for the tip displacement u Z are compared with the exact solution u Z ,ref in the form of quotient u Z /u Z ,ref . The variation of u Z /u Z ,ref is shown with respect to the structural parameter G L 2 /Eh 2 in the logarithmic scale in Fig. 2.
We can observe that the linear elements (N = 0) are prone to shear locking indeed and that locking is efficiently resolved by using the reduced integration. The quadratic elements (N = 1) do not seem to suffer severe locking. Nonetheless, the use of the reduced integration is advisable, because it gives more accurate results. The higher-order elements (N ≥ 2) enable us to use the full integration with no harm.
Note that the origin of the shear locking problems demonstrated here is the use of the standard Galerkin finite-element approach and are not imposed by the use of quaternions. This was indirectly confirmed in [45], where the non-standard, yet a quaternion-based finite-element approach is presented which is locking-free.

A cantilever beam subject to a free-end moment
A straight cantilever beam is subjected to the in-plane moment M Y = 1000 at its free end (Fig. 3). The deformed shape is characterized by large rotations and large displacements both in the X and Z direction. The analytical solution [29] shows that the beam deforms into a part of a circle. The in-plane rotation of the beam is thus a linear function of x: Fig. 3 The cantilever under free-end in-plane moment 2.857143 Exact nonlinear [29] 2.857143 N = number of internal points, n e = number of elements This simple test is unpretentious for the rotational vector-based formulations but very significant when using rotational quaternions. In contrast to the rotational vector-based formulations, where the first-order interpolation describes the actual rotations exactly, the rotational quaternion-based approach employs the trigonometric functions of the angle of rotation. This example indicates that higher-order elements may sometimes be needed when using the rotational quaternions, which is also confirmed in the convergence study presented in Table 1.
The following geometric and material properties of the cantilever are taken into account: We can observe from Table 1 that a single high-order element gives the result accurate in all significant digits. The accuracy of low-order elements can be improved by increasing the density of the finite-element mesh. Clearly, the use of a small number of elements in conjunction with a high order of interpolation is computationally more efficient. Note that for this problem, the present formulation shows better convergence, especially for higher magnitudes of the applied moment, when compared to the quaternion-based elements, presented in [45]. This might be due to a different treatment of the moment equilibrium equations.

The MacNeal and Harder test problems
In their paper from 1985 [16], MacNeal and Harder proposed a standard set of problems to test the finiteelement accuracy of beam, shell, and solid elements. Three of them are applicable in the spatial beams: (i) the straight cantilever beam; (ii) the curved beam, and (iii) the twisted beam. The data for the straight beam are:  The twisted beam has a straight axis, while the cross-sections are linearly rotated from an angle ϑ = 0 at the clamped edge to an angle ϑ = π/2 at the free end, see also discussion in [42]. The properties of the beam are: For both the curved and the twisted beam, we consider two load cases: the unit force in the in-plane and in the out-of-plane direction, see Fig. 4. The results obtained by the cubic elements are presented in Table 2 and compared to the theoretical results of MacNeal and Harder [16] as well as to the analytical results of the linearized Reissner beam theory, see [41,Appendix D]. Some small differences between the theoretical [16] and the analytical results [42] can be observed which are due to the different theoretical assumptions of the two beam theories. It is important to observe that the results using the present elements agree well with the analytical results. The mesh of six cubic elements (N = 2) with an appropriate description of the initial curved geometry gives results accurate to at least four significant digits. The importance of a proper description of the initial geometry for the curved beam is demonstrated in Table 2. When the beam is modeled with six straight elements, the numerical results agree with the analytical ones only in two significant digits; a mesh of 40 straight elements is needed for an accuracy, obtained with a mesh of six initially curved elements.
The present quaternion-based elements are thus capable of well describing the beam behavior with high accuracy even for an initially curved and/or twisted geometry of the three-dimensional beams.

Bending of 45 • cantilever
This problem, first presented by Bathe and Bolourchi [4], represents one of the classical finite-element tests since it includes all possible natural modes of deformation of a beam: bending, shear, extension, and torsion. The data are taken from [4]: the circular arch with radius 100 is located in the horizontal plane and clamped at one end (Fig. 5). The cross-section is taken to be a unit square. The arch is subjected to a vertical load F = 600 at the free end. The elastic and the shear moduli of material are E = 10 7 and G = E/2.
In Table 3, we compare the results for the position vector of the free end of the cantilever. Eight straight elements of various order were used to obtain the results of the present formulation. The present results agree well with the results of other authors. An excellent performance and numerical stability of our method can be observed even for very high order of interpolation and a rapid convergence with respect to the order is evident. In contrast, for N = 15, some convergence problems with formulation [45] were observed.

Lateral buckling of a right-angle frame
This classical problem was introduced by Argyris et al. [3] and was later studied by many authors. The simply supported thin right-angle frame is analyzed to check the ability of the present formulation to consider the torsional-bending coupling properly. At the two supported ends, the frame is subjected to the in-plane moments in the opposite directions, see Fig. 6. The geometrical and material parameters are taken from [2] and are: Because the slenderness of the rectangular cross-section is extreme, the frame buckles out of the plane. This phenomenon is numerically detected by the onset of singularity of the tangent stiffness matrix. The accuracy of the solution is dependent on the density of the mesh used, see Table 4. Due to the symmetry of the problem only a half of the structure was modeled. Our results are compared to the analytical solution by Timoshenko and Gere [35] and to the numerical results of other authors. From Table 4, we can observe that the present approach gives the critical moment accurate in all significant digits for a mesh of only two elements with four internal points. This is also more accurate with respect to the results of other formulations including the quaternion-based approach in [45]. The advantage of the present approach for this particular problem might be due to a convenient form of the quaternion moment, which results in the boundary point moments being the direct members of the tangent stiffness matrix, see Eq. (53).

Conclusion
We presented a fully consistent standard formulation of the Reissner three-dimensional beam theory based on the rotational quaternions. We strictly use the quaternion algebra in the derivation of the governing equations. The quaternion moment and the consistent quaternion-based moment equilibrium equation have been derived in the continuous and the discrete forms. The numerical implementation is carried out in detail and some numerical examples have been presented. The proposed method proves to be numerically robust, accurate, and computationally efficient. These important properties stem from the present kinematically consistent approach, in which the algebraic constraint for the rotational quaternions is properly considered. Without much additional effort, the present model is very suitable for the implementation in existing finite-element computer environments.

Appendix: The details about the linearization
Here, we derive the linearization of several rotation-dependent terms which will serve well in the linearization process of the governing equations. Special care is taken in considering the nature of the varied quantities. Note that the vector quantities need to be extended to pure quaternions whenever the rotations are applied. The variation of a pure quaternion is properly reduced to three dimensions whenever needed. The linearization of the conjugated rotational quaternion follows directly from Eq. (9): The linearization of the translational strain follows from (27): where Note that ∂ γ G ∂ r and ∂ γ G ∂ q are 4 × 4 matrices. Displacements and strains are pure quaternions, and thus, their variations δ r and δ γ are pure quaternions as well. The trivial scalar components of pure quaternions will be eliminated from expressions. Considering the nature of the pure quaternions, we introduce the following combined vector-quaternion notation to express the variations where [.] IR 3×3 denotes the sub-matrix obtained by removing the first column and the first row of the original matrix, and [.] IR 3×4 denotes the matrix obtained by removing only the first row. For the sake of simplicity, the following notation will be used instead: where The linearization of rotational strains from (28) gives where The linearization of the stress-resultant vectors, when expressed in the local basis, follows directly from (50)-(51): When vector quantities are expressed as pure quaternions, Eq. (62) reads where the hat over the matrix in IR 3×3 denotes its natural expansion into IR 4×4 : The stress resultants as they appear in the equilibrium equations are expressed with respect to the fixed basis, and thus, we also need to linearize the coordinate transformation, which is, in the present paper, expressed using quaternions (see Eq. 6), e.g., For the resultant force vector, we have After considering the nature of the pure quaternions, we finally have where In an analogous way, we obtain where Prior to the linearization of the moment equilibrium equation (47), it is useful to introduce the following notations: and linearize them separately. The direct linearization of the above terms yields S is a natural extension of the skew-symmetric matrix associated with the cross vector product into four dimensions and defines the cross-vector product of the pure quaternions as defined in Eq. (11): The linearization of F m3 additionally requires considering κ • q = 2 q and performing some rearrangement of terms to get Finally, the linearization of F m4 gives