Eva Zupan
email: eva.zupan@fgg.uni-lj.si

Miran Saje
email: miran.saje@fgg.uni-lj.si

Dejan Zupan
email: dejan.zupan@fgg.uni-lj.si

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

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

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 [START_REF] Mortensen | Strapdown guidance error analysis[END_REF][START_REF] Phillips | Review of attitude representations used for aircraft kinematics[END_REF][START_REF] Silson | Coarse alignment of a ship's strapdown inertial attitude reference system using velocity loci[END_REF], robotics [START_REF] Chen | Parameterization and evaluation of robotic orientation workspace: a geometric treatment[END_REF][START_REF] Yun | Design, implementation, and experimental results of a quaternion-based Kalman filter for human body motion tracking[END_REF], animating rotations in computer graphics [START_REF] Faraway | Modelling orientation trajectories[END_REF][START_REF] Shoemake | Animating rotation with quaternion curves[END_REF], dynamics of rigid bodies [START_REF] Johnson | Quaternion-based rigid body rotation integration algorithms for use in particle methods[END_REF][START_REF] Udwadia | Equations of motion for general constrained systems in Lagrangian mechanics[END_REF], and flexible beams [START_REF] Lang | Multi-body dynamics simulation of geometrically exact Cosserat rods[END_REF][START_REF] Mcrobie | Simo-Vu Quoc rods using Clifford algebra[END_REF]. 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 [START_REF] Géradin | Flexible multibody dynamics : a finite element approach[END_REF][START_REF] Ghosh | Consistent quaternion interpolation for objective finite element approximation of geometrically exact beam[END_REF][START_REF] Lang | Multi-body dynamics simulation of geometrically exact Cosserat rods[END_REF][START_REF] Zupan | The quaternion-based three-dimensional beam theory[END_REF]. Our recent works in the quaternion-based beam formulations for static and dynamic analyses [START_REF] Zupan | Quaternion-based dynamics of geometrically nonlinear spatial beams using the Runge-Kutta method[END_REF][START_REF] Zupan | The quaternion-based three-dimensional beam theory[END_REF] 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 [START_REF] Zupan | Quaternion-based dynamics of geometrically nonlinear spatial beams using the Runge-Kutta method[END_REF][START_REF] Zupan | The quaternion-based three-dimensional beam theory[END_REF] are non-standard and differ considerably from the standard finite-element approaches such as [START_REF] Ghosh | Consistent quaternion interpolation for objective finite element approximation of geometrically exact beam[END_REF][START_REF] Simo | A finite strain beam formulation. The three-dimensional dynamic problem[END_REF]. 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., [START_REF] Bauchau | Interpolation of finite rotations in flexible multibody dynamics simulations[END_REF][START_REF] Bauchau | Flexible Multibody Dynamics[END_REF][START_REF] Géradin | Kinematics and Dynamics of Rigid and Flexible Mechanisms using Finite Elements and Quaternion Algebra[END_REF][START_REF] Romero | The interpolation of rotations and its finite element implementation of geometrically exact rods[END_REF][START_REF] Zupan | Integrating rotation from angular velocity[END_REF]. 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 [START_REF] Bottasso | A non-linear beam space-time finite element formulation using quaternion algebra: Interpolation of the Lagrange multipliers and the appearance of spurious modes[END_REF]. 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 [START_REF] Bottasso | A non-linear beam space-time finite element formulation using quaternion algebra: Interpolation of the Lagrange multipliers and the appearance of spurious modes[END_REF] 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 [START_REF] Ritto-Corrêa | Work-conjugacy between rotation-dependent moments and finite rotations[END_REF] 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 [START_REF] Zupan | The quaternion-based three-dimensional beam theory[END_REF] 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 [START_REF] Ritto-Corrêa | Work-conjugacy between rotation-dependent moments and finite rotations[END_REF]. 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 [START_REF] Zupan | The quaternion-based three-dimensional beam theory[END_REF]. 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 [START_REF] Zupan | The quaternion-based three-dimensional beam theory[END_REF], 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

B G = G 1 (x) , G 2 (x) , G 3 (x)
, 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., [START_REF] Argyris | An excursion into large rotation[END_REF][START_REF] Argyris | Large rotations-revisited application of Lie algebra[END_REF][START_REF] Ritto-Corrêa | On the differentiation of the Rodrigues formula and its significance for the vector-like parameterization of Reissner-Simo beam theory[END_REF]. 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:

y = x + sin ϑ (n × x) + (1 -cos ϑ)n × (n × x) = Rx, ( 1
)
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:

x = R x G and x G = R T x, (2)
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 [START_REF] Ward | Quaternions and Cayley Numbers[END_REF]. A detailed application of the quaternion representation of rotations in structural mechanics was presented in Zupan et al. [START_REF] Zupan | Integrating rotation from angular velocity[END_REF][START_REF] Zupan | The quaternion-based three-dimensional beam theory[END_REF]. 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

a • b = a 0 b 0 -a • b + b 0 a + a 0 b + a × b , for a = a 0 + a , b = b 0 + b ∈ IH , 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

B g = { g 0 , g 1 , g 2 , g 3 }, using the rule g 0 = 1+ 0 and g i = 0+ g i , i = 1, 2, 3. Let a = a 0 + a = [a 0 a 1 a 2 a 3 ] T and b = b 0 + b = [b 0 b 1 b 2 b 3 ]
T be two arbitrary quaternions in the component matrix form. Then, the quaternion multiplication can also be expressed by the matrix products as

a • b = φ L a b = ⎡ ⎢ ⎣ a 0 -a 1 -a 2 -a 3 a 1 a 0 -a 3 a 2 a 2 a 3 a 0 -a 1 a 3 -a 2 a 1 a 0 ⎤ ⎥ ⎦ ⎡ ⎢ ⎣ b 0 b 1 b 2 b 3 ⎤ ⎥ ⎦ (3)
and

a • b = φ R b a = ⎡ ⎢ ⎣ b 0 -b 1 -b 2 -b 3 b 1 b 0 b 3 -b 2 b 2 -b 3 b 0 b 1 b 3 b 2 -b 1 b 0 ⎤ ⎥ ⎦ ⎡ ⎢ ⎣ a 0 a 1 a 2 a 3 ⎤ ⎥ ⎦. (4)
We further introduce the conjugation of a quaternion, a * , which only alters the sign of the vector part of the quaternion, a * = a 0a . In the quaternion description, Eqs. ( 1)-( 2) take the form

y = q • x • q * = φ L ( q) φ R q * x (5)
and

x = q • x G • q * = φ L ( q) φ R q * x G , ( 6
)
x G = q * • x • q = φ L q * φ R ( q)
x, [START_REF] Bottasso | A non-linear beam space-time finite element formulation using quaternion algebra: Interpolation of the Lagrange multipliers and the appearance of spurious modes[END_REF] 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:

q = cos ϑ 2 + n sin ϑ 2 = exp ϑ n 2 . ( 8
)
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,

| q| = q • q = 1 ⇐⇒ q • q = 1, (9)
which also implies

q • q * = 1 = 1 + 0 . ( 10
)
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,

a × b = 0 + a × b = a × b, ( 11
)
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

a • b = n i=1 a i • b i , a, b ∈ IR n , for a = [a i ] i=1,2,...,n and b = [b i ] i=1,2,.
..,n , the following holds true [START_REF] Rudin | Functional Analysis. 2nd edn[END_REF]:

(Fa) • b = a • F T b , ( 12
)
where F is the matrix representation of the linear operator F :

IR n → IR n . It is important to observe that φ L ( c ) T = φ L c * and φ R ( c ) T = φ R c * , ( 13
)
which follows directly from (3)-(4). After applying ( 12)-( 13) to arbitrary quaternions a, b, and c, we have

( c • a ) • b = a • c * • b , ( 14
)
( a • c ) • b = a • b • c * . ( 15
)
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 [START_REF] Zupan | The quaternion-based three-dimensional beam theory[END_REF]. 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

κ = 2 q • q * κ G = 2 q * • q , ( 16
)
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 2δ q • q * = 0 + δϑ = δ ϑ, [START_REF] Mcrobie | Simo-Vu Quoc rods using Clifford algebra[END_REF] 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

q * = -q * • q • q * . ( 18
)

The generalized virtual work principle

Following the approach of Reissner [START_REF] Reissner | On one-dimensional large-displacement finite-strain beam theory[END_REF][START_REF] Reissner | On finite deformations of space[END_REF], we assume the principle of virtual work in the form

L 0 N G (x) • δγ G (x) rel + M G (x) • (δκ G (x)) rel dx = L 0 [n (x) • δr (x) + m (x) • δϑ (x)] dx + F (x) • δr (x) + P (x) • δϑ (x) L 0 , ( 19
)
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 [START_REF] Troutman | Variational Calculus with Elementary Convexity[END_REF]. 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 [START_REF] Mukherjee | Analysis of shear locking in Timoshenko beam elements using the function space approach[END_REF]. 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 [START_REF] Mcrobie | Simo-Vu Quoc rods using Clifford algebra[END_REF], so that the generalized virtual work principle using the quaternion representation of rotation reads:

L 0 N G • δγ G rel + M G • (δκ G ) rel dx = L 0 n • δr + 2 ( m • q) • δ q dx - L 0 2λ q • δ q dx + F • δr + 2 P • q • δ q L 0 .
(20)

The kinematic equations

In the spirit of Reissner's works [START_REF] Reissner | On one-dimensional large-displacement finite-strain beam theory[END_REF][START_REF] Reissner | On finite deformations of space[END_REF], 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:

n = -N = -(R N G ) = -q • N G • q * IR 3 , ( 21
) m = -M -r × N = -(R M G ) -r × (R N G ) = -q • M G • q * -r × q • N G • q * IR 3 . ( 22
)
Here, [ . ] IR 3 indicates the vector part of a quaternion. The generalized virtual work principle given in [START_REF] Phillips | Review of attitude representations used for aircraft kinematics[END_REF] 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

x 2 x 1 N G • δγ G rel -q * • δ r • q IR 3 -2 q * • r × δ q • q * • q IR 3 dx + x 2 x 1 M G • (δκ G ) rel -2 q * • δ q • q * • q IR 3 dx + L 0 2 λ q • δ q dx + N G • q * • δ r • q IR 3 x 2 x 1 + M G • 2 q * • δ q IR 3 x 2 x 1 -F • δr + P • 2 δ q • q * IR 3 x 2 x 1 = 0. ( 23
)
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

δγ G rel = q * • δ r • q + 2 q * • r × δ q • q * • q IR 3 , ( 24
)
(δκ G ) rel = 2 q * • δ q • q * • q IR 3 , ( 25
) 0 = 2 q • δ q. ( 26
)
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 [START_REF] Zupan | The quaternion-based three-dimensional beam theory[END_REF] 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 [START_REF] Reissner | On one-dimensional large-displacement finite-strain beam theory[END_REF] who proposed a simplification of rotations, Simo [START_REF] Simo | A finite strain beam formulation. The three-dimensional dynamic problem[END_REF] 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 [START_REF] Simo | A finite strain beam formulation. The three-dimensional dynamic problem[END_REF] for the integration of the kinematic equations when expressed in terms of the rotational vector. As the details of the derivation are given in [START_REF] Zupan | The quaternion-based three-dimensional beam theory[END_REF], only the final result is displayed here:

γ G = q * • r • q + c G (27) κ G = 2 q * • q + d G . ( 28
)
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), [START_REF] Mortensen | Strapdown guidance error analysis[END_REF], and the integration by parts is used to obtain

N • δ r L 0 - L 0 N • δ r dx - L 0 2 r × N • q • δ q dx + 2 M • q • δ q L 0 + L 0 2 M • q • δ q dx + L 0 2 M • q • δ q dx -F • δr + 2 P • q • δ q L 0 - L 0 [n • δr + 2 ( m • q) • δ q] dx + L 0 2λ q • δ q dx = 0. ( 29
)
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]:

N + n = 0 (30) M + r × N + m -λ 1 • q = 0 (31)
and the boundary conditions:

N 0 -F 0 = 0 (32) M 0 -P 0 • q 0 = 0 (33) N L -F L = 0 (34) M L -P L • q L = 0. (35)
Equation ( 30) coincides with the standard equilibrium equation for forces, which does not depend on the rotation. In contrast, the equilibrium equation [START_REF] Silson | Coarse alignment of a ship's strapdown inertial attitude reference system using velocity loci[END_REF] for moments differs from the standard one. If we use the notation M for the standard form of the moment equilibrium equation

M = M + r × N + m = 0, (36)
and insert [START_REF] Troutman | Variational Calculus with Elementary Convexity[END_REF] into [START_REF] Silson | Coarse alignment of a ship's strapdown inertial attitude reference system using velocity loci[END_REF], we obtain

M -λ 1 • q = 0. ( 37
)
Now the difference between the two expressions becomes evident. The result shows how the quaternion-based moment equilibrium equation [START_REF] Silson | Coarse alignment of a ship's strapdown inertial attitude reference system using velocity loci[END_REF] 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 [START_REF] Phillips | Review of attitude representations used for aircraft kinematics[END_REF] we can observe that the virtual work of the external distributed moments L 0 m (x) • δϑ (x) dx now takes the equivalent form L 0 2 ( m • q)•δ q dx. Thus, the definition of the energy-consistent quaternion moment in the sense of the work by Ritto-Corrêa and Camotim [START_REF] Ritto-Corrêa | Work-conjugacy between rotation-dependent moments and finite rotations[END_REF] is evident. Let M ϑ be the rotational moment as defined by Ritto-Corrêa and Camotim [START_REF] Ritto-Corrêa | Work-conjugacy between rotation-dependent moments and finite rotations[END_REF]. Then, the quaternion moment M q can be defined as

M q = 2 M ϑ • q, ( 38
)
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

M -λ 1 = 0 or, equivalently, λ = 0 (39) M = 0. (40)
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. [START_REF] Udwadia | Equations of motion for general constrained systems in Lagrangian mechanics[END_REF], and it could be eliminated from the quaternion moment equilibrium equation [START_REF] Silson | Coarse alignment of a ship's strapdown inertial attitude reference system using velocity loci[END_REF] or [START_REF] Udwadia | Equations of motion for general constrained systems in Lagrangian mechanics[END_REF] 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 [START_REF] Zupan | Rotational invariants in finite element formulation of three-dimensional beam theories[END_REF]. This kind of approach has been followed in our previous work, see [START_REF] Zupan | The quaternion-based three-dimensional beam theory[END_REF]. 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 [START_REF] Phillips | Review of attitude representations used for aircraft kinematics[END_REF] by adding the term

L 0 δλ ( q • q -1) dx . ( 41
)
Such a generalization of the virtual work principle is fully equivalent to the original one, [START_REF] Phillips | Review of attitude representations used for aircraft kinematics[END_REF], when dealing with the continuous unknowns, see [36, pp. 129-130].

The variations of the primary unknowns are interpolated in a standard manner:

δr (x) = i P i (x) δr i , δ r i = δr (x i ) , ( 42
)
δ q (x) = i P i (x) δ q i , δ q i = δ q (x i ) , ( 43
)
δλ (x) = i P i (x) δλ i , δ λ i = δλ (x i ) , ( 44
)
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 [START_REF] Phillips | Review of attitude representations used for aircraft kinematics[END_REF], applying ( 42)-( 44) and after some rearrangement of the terms, we obtain the following discretized generalized virtual work principle:

i L 0 N P i -n P i dx -F P i L 0 • δr i + i L 0 -M • κ • q P i -2 r × N • q P i -2 ( m • q) P i + 2λ q P i + 2 M • q P i dx-2 P • q P i L 0 • δ q i + i L 0 ( q • q -1) P i dx δλ i = 0. ( 45
)
Note that the same result can also be obtained by adding the integrated unit length constraint [START_REF] Zupan | The linearized three-dimensional beam theory of naturally curved and twisted beams: the strain vectors formulation[END_REF] 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:

F (1-3) i : L 0 N P i dx -n P i dx + δ i0 F 0 -δ i L F L = 0, (46)
F (4-7) i : L 0 2 M • q P i -M • κ + 2 r × N + 2 m -2λ 1 • q P i dx + 2δ i0 P 0 • q 0 -2δ i L P L • q L = 0, (47)
F (8) i : L 0 q • q -1 P i dx = 0, (48)
for i = 0, 1, . . . , N + 1. Here,

δ i0 = 1, for i = 0 0, otherwise δ i L = 1, for i = N + 1 0, otherwise .
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

L 0 f (x) dx → I N j=0 w j f x j .
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., [START_REF] Simo | A three-dimensional finite-strain rod model. Part II: computational aspects[END_REF] and the references therein. The quaternion approach seems to inherit this inconvenience of standard discretizations. In [START_REF] Zupan | The quaternion-based three-dimensional beam theory[END_REF], the shear locking phenomenon was resolved by imposing the consistency conditions (see also [START_REF] Zupan | Rotational invariants in finite element formulation of three-dimensional beam theories[END_REF]). 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 [START_REF] Zupan | The quaternion-based three-dimensional beam theory[END_REF] is fully applicable in the present approach, yet it would demand a generalized virtual work principle augmented with two additional terms. See also [START_REF] Mukherjee | Analysis of shear locking in Timoshenko beam elements using the function space approach[END_REF][START_REF] Reddy | On locking-free shear deformable beam finite elements[END_REF] 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:

y [n+1] = y [n] + δ y.
For the displacements and the Lagrangian multiplier such an additive update is adequate:

r [n+1] = r [n] + δr, λ [n+1] = λ [n] + δλ.
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:

q [n+1] = q • q [n] , where q = exp δ q • q [n] * . ( 49
)
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:

κ [n+1] G [n+1] = κ [n] G [n] + κ G [n+1] , where κ G [n+1] = q [n+1] * • 2 q • q * • q [n+1] .
Details are given in [START_REF] Zupan | The quaternion-based three-dimensional beam theory[END_REF].

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:

N G = ⎡ ⎣ E A 1 0 0 0 G A 2 0 0 0 G A 3 ⎤ ⎦ γ G = C N γ G , ( 50
)
M G = ⎡ ⎣ G J 1 0 0 0 E J 2 0 0 0 E J 3 ⎤ ⎦ κ G = C M κ G . ( 51
)
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 G 2 and G 3 of the cross-section; J 2 and J 3 are the corresponding principal bending inertial moments of the cross-section. The strain vectors in (50)-( 51) are expressed with the displacements and the rotational quaternions by the strong form of kinematic equations ( 27)- [START_REF] Rudin | Functional Analysis. 2nd edn[END_REF]. Once ( 50)-( 51) and ( 27)-( 28) are accounted for, the system of the governing equations is complete, thus enabling us to solve the equations with respect to the unknown nodal values of displacements, rotational quaternions and Lagrangian multipliers.

The linearization of the discrete governing equations

When the rotational quaternions are used, the linearization is substantially more simple when compared to the rotational vector-based formulations. Moreover, the derivatives of variations follow directly from ( 42)-( 43), e.g.,

δr = i P i δr i and δ q = i P i δ q i .
Following the rules of variational calculus, Eq. ( 46) is linearized as

δF (1-3) i = L 0 δ N P i dx = j L 0 ∂ N ∂γ G ∂γ G ∂ r P i P j dx δr j + j L 0 ∂ N ∂γ G ∂γ G ∂ q + ∂ N ∂ q P i P j dx δ q j . ( 52
)
In the linearization of Eq. ( 47), it is suitable to introduce the notations

F m1 = -2 r × N • q, F m2 = -2 M • q, F m3 = -M • κ • q, F m4 = -2 m • q.
The linearization then takes the form:

δF (4-7) i = L 0 ∂ F m1 ∂ r P i dx δr + L 0 ∂ F m2 ∂ q P i δ q dx + L 0 ∂ F m1 ∂ q + ∂ F m3 ∂ q + ∂ F m4 ∂ q P i δ q dx + L 0 ∂ F m2 ∂ q P i δ q dx + L 0 ∂ F m3 ∂ q P i δ q dx + 2δ i0 φ L P 0 δ q 0 -2δ i L φ L P L δ q L + 2 L 0 q P i δλ dx + 2 L 0 λP i δ q dx = j L 0 ∂ F m1 ∂ r P i P j dx δr j + j 2 L 0 q P i P j dx δλ j + j L 0 ∂ F m1 ∂ q + ∂ F m3 ∂ q + ∂ F m4 ∂ q + 2λ P i P j dx + L 0 ∂ F m2 ∂ q P i P j dx + L 0 ∂ F m2 ∂ q + ∂ F m3 ∂ q P i P j dx δ q j + 2δ i0 φ L P 0 δ q 0 -2δ i L φ L P L δ q L . ( 53
)
Finally, the linearization of Eq. (48) reads

δF (8) i = 2 L 0 q • δ q P i dx = 2 j L 0 q P i P j dx • δ q j . ( 54
)
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 [START_REF] Zupan | The quaternion-based three-dimensional beam theory[END_REF]. 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 [START_REF]MATLAB, Using MATLAB[END_REF] computing environment. Finite elements with various numbers of interpolation points, N + 2, are used to investigate the influence of their number on the accuracy of numerical solutions. A number of internal points, N , will be used to describe the applied elements.

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 [START_REF] Zupan | The quaternion-based three-dimensional beam theory[END_REF], 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 [START_REF] Saje | Large deformations of thin curved plane beam of constant initial curvature[END_REF] 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

ϑ Y (x) = M Y E J 2 x.
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:

E = 2.1 × 10 4 G = 1.05 × 10 4 L = 100 A 1 = 20 A 2 = A 3 = 16 J 1 = 6.4566 J 2 = 1.6667 J 3 = 666.66.
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 [START_REF] Zupan | The quaternion-based three-dimensional beam theory[END_REF]. This might be due to a different treatment of the moment equilibrium equations.

The MacNeal and Harder test problems

In their paper from 1985 [START_REF] Macneal | A proposed standard set of problems to test finite element accuracy[END_REF], 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: where L is the length of the beam. Four load cases are considered: a unit force in all three coordinate directions resulting in extension, in-plane shear and out-of-plane shear, and a unit torsional moment resulting in twist.

A 1 = 0.02 L = 6.0 A 2 = A 3 = 0.016 J 1 = 4.8630 × 10 -5 E = 1.0 × 10 7 J 2 = 6.66667 × 10 -5 G = 3.84615 × 10 6 J 3 = 1.66667 × 10 -5 ,
The axis of the curved beam in the initial undeformed configuration has the shape of a quarter of a circle with radius R = 4.22. Only G and J 1 are changed in this beam:

G = 4.0 × 10 6 , J 1 = 4.4050 × 10 -5 .
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 [START_REF] Zupan | On "A proposed standard set of problems to test finite element accuracy": the twisted beam[END_REF]. The properties of the beam are:

A 1 = 0.3520 L = 12.0 A 2 = A 3 = 0.2816 J 1 = 0.0385 E = 2.9 × 10 7 J 2 = 0.0030 G = 1.1885 × 10 7 J 3 = 0.0355.
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 [START_REF] Macneal | A proposed standard set of problems to test finite element accuracy[END_REF] as well as to the analytical results of the linearized Reissner beam theory, see [START_REF] Zupan | The linearized three-dimensional beam theory of naturally curved and twisted beams: the strain vectors formulation[END_REF]Appendix D]. Some small differences between the theoretical [START_REF] Macneal | A proposed standard set of problems to test finite element accuracy[END_REF] and the analytical results [START_REF] Zupan | On "A proposed standard set of problems to test finite element accuracy": the twisted beam[END_REF] 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 [START_REF] Bathe | Large displacement analysis of three-dimensional beam structures[END_REF], 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 [START_REF] Bathe | Large displacement analysis of three-dimensional beam structures[END_REF]: 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 [START_REF] Zupan | The quaternion-based three-dimensional beam theory[END_REF] were observed.

Lateral buckling of a right-angle frame

This classical problem was introduced by Argyris et al. [START_REF] Argyris | On large displacements-small strain analysis of structures with rotational degrees of freedom[END_REF] 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 [START_REF] Argyris | Large rotations-revisited application of Lie algebra[END_REF] and are:

J 1 = 1.50 A 1 = 18 E = 71240 J 2 = 0.54 A 2 = 15 G = 39155.04 J 3 = 1350 A 3 = 15 L = 240.
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 [START_REF] Timoshenko | Theory of Elastic Stability[END_REF] 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 [START_REF] Zupan | The quaternion-based three-dimensional beam theory[END_REF]. 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):

δ q • q * + q • δ q * = 0, δ q * = -q * • δ q • q * = -φ L q * φ R q * δ q. ( 55
)
The linearization of the translational strain follows from [START_REF] Romero | The interpolation of rotations and its finite element implementation of geometrically exact rods[END_REF]:

δ γ G = δ q * • r • q + q * • δr • q + q * • r • δ q = φ L q * φ R ( q) δ r + φ L γ G -c G -φ R γ G -c G φ L q * δ q = ∂ γ G ∂ r δ r + ∂ γ G ∂ q δ q, ( 56
)
where

∂ γ G ∂ r = φ L q * φ R ( q) , ( 57
) ∂ γ G ∂ q = φ L γ G -c G -φ R γ G -c G φ L q * . ( 58
)
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

δγ G = ∂ γ G ∂ r IR 3×3 δr + ∂ γ G ∂ q IR 3×4 δ q, ( 59
)
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:

δγ G = ∂γ G ∂ r δr + ∂γ G ∂ q δ q, ( 60
)
where

∂γ G ∂ r = ∂ γ G ∂ r IR 3×3 and ∂γ G ∂ q = ∂ γ G ∂ q IR 3×4 .
The linearization of rotational strains from [START_REF] Rudin | Functional Analysis. 2nd edn[END_REF] gives

δ κ G = 2δ q * • q + 2 q * • δ q = 2φ L q * δ q -φ L q * φ R ( κ G ) δ q = ∂ κ G ∂ q δ q + ∂ κ G ∂ q δ q, ( 61
)
where

∂ κ G ∂ q = 2φ L q * , ∂ κ G ∂ q = -φ L q * φ R ( κ G ) .
The linearization of the stress-resultant vectors, when expressed in the local basis, follows directly from (50)-(51):

δ N G = C N δγ G , δM G = C M δκ G . ( 62
)
When vector quantities are expressed as pure quaternions, Eq. ( 62) reads

δ N G = C N δ γ G , δ M G = C M δ κ G , ( 63
)
where the hat over the matrix in IR 3×3 denotes its natural expansion into IR 4×4 :

C N = 1 0 T 0 C N , C M = 1 0 T 0 C M .
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.,

N = q • N G • q * , M = q • M G • q * .
For the resultant force vector, we have

δ N = δ q • N G • q * + q • δ N G • q + q • N G • δ q * = δ q • q * • q • N G • q * + q • C N δ γ G • q * -q • N G • q * • δ q • q * = φ L ( q) φ R q * C N δ γ G + φ R N g -φ L N g φ R q * δ q = ∂ N ∂ γ G δ γ G + ∂ N ∂ q δ q.
After considering the nature of the pure quaternions, we finally have

δ N = ∂ N ∂γ G δγ G + ∂ N ∂ q δ q, ( 64
)
where

∂ N ∂γ G = φ L ( q) φ R q * C N IR 3×3 , ∂ N ∂ q = φ R N -φ L N φ R q * IR 3×4 .
In an analogous way, we obtain

δ M = ∂ M ∂ κ G δ κ G + ∂ M ∂ q δ q, ( 65
)
where

∂ M ∂ κ G = φ L ( q) φ R q * C M , ∂ M ∂ q = φ R M -φ L M φ R q * .
Prior to the linearization of the moment equilibrium equation (47), it is useful to introduce the following notations:

F m1 = -2 r × N • q, F m2 = -2 M • q, F m3 = -M • κ • q, F m4 = -2 m • q,
and linearize them separately. The direct linearization of the above terms yields δ F m1 = -2φ R ( q) -N × δ r -2φ R ( q) r × δ N -2φ L r × N δ q = -2 φ R ( q) S (-N)

+ S r ∂ N ∂γ G ∂γ G ∂ r IR 4×3 δr -2 φ R ( q) S r ∂ N ∂ γ G ∂ γ G ∂ q + ∂ N ∂ q + φ L r × N δ q = ∂ F m1 ∂ r δr + ∂ F m1 ∂ q δ q.
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):

S (v) = ⎡ ⎢ ⎣ 1 0 0 0 0 0 -v 3 v 2 0 v 3 0 -v 1 0 -v 2 v 1 0 ⎤ ⎥ ⎦ ,
composed from components v 1 , v 2 , v 3 of a vector v. Similarly,

δ F m2 = 2φ R ( q) δ M + 2φ L M δ q = 2φ R ( q) ∂ M ∂ κ G ∂ κ G ∂ q δ q + ∂ κ G ∂ q δ q + ∂ M ∂ q δ q + 2φ L M δ q = ∂ F m2 ∂ q δ q + ∂ F m2 ∂ q δ q .
The linearization of F m3 additionally requires considering κ • q = 2 q and performing some rearrangement of terms to get

δ F m3 = -2φ R q δ M -2φ L M δ q = 2φ R q ∂ M ∂ κ G ∂ κ G ∂ q δ q + ∂ κ G ∂ q δ q + ∂ M ∂ q δ q + 2φ L M δ q = ∂ F m3 ∂ q δ q + ∂ F m3 ∂ q δ q .
Finally, the linearization of F m4 gives δ F m4 = -φ L ( m) δ q = ∂ F m4 ∂ q δ q.

Fig. 1

 1 Fig. 1 Model of the three-dimensional beam

Fig. 2

 2 Fig. 2 Normalized vertical tip displacement versus G L 2 /Eh 2

Fig. 4 Fig. 5

 45 Fig.[START_REF] Bathe | Large displacement analysis of three-dimensional beam structures[END_REF] The MacNeal and Harder test problems

Fig. 6

 6 Fig.[START_REF] Bauchau | Flexible Multibody Dynamics[END_REF] Lateral buckling of a right-angle frame

Table 1

 1 Rotations of the cantilever under an in-plane moment

		n e = 1	n e = 2	n e = 10	n e = 50
	N = 1	2.825955	2.854684	2.857139	2.857143
	N = 4	2.855784	2.857135	2.857143
	N = 6	2.857142	2.857143
	N = 9	2.857143
	Exact nonlinear [29]				2.857143
	N = number of internal points, n e = number of elements

Table 2

 2 Present results and analytical solutions of the MacNeal and Harder test problems

	Initial beam type	Type of solution	Displacement in the direction of the load
			Extension (F 1 )	In-plane shear (F 2 )	Out-of-plane shear (F 3 )	Twist(M 1 )
	Straight	6 straight elements	3.0 × 10 -5	0.10810	0.43210	0.032079
		Theoretical [16]	3 .0 × 10 -5	0.1081	0.4321	0.03208
		Analytical	3.0 × 10 -5	0.10810	0.43210	0.032079
	Curved	6 straight elements	-	0.087349	0.50222	-
		40 straight elements	-	0.088576	0.50608	-
		6 curved elements	-	0.088577	0.50615	-
		Theoretical [16]	-	0 .08734	0.5022	-
		Analytical	-	0.088604	0.50617	-
	Twisted	6 straight elements	-	0.0054290	0.0017494	-
		Theoretical [16]	-	0 .005422	0.001754	-
		Analytical [45]	-	0 .0054293	0.0017496	-

Table 3

 3 Free-end position of the 45 • bend cantilever

	Formulation	r X	r Y	r Z
	Present, N = 0	16.5541	46.6026	53.5528
	Present, N = 1	15.2860	47.4156	53.4733
	Present, N = 3	15.2862	47.4159	53.4725
	Present, N = 7	15.2862	47.4159	53.4725
	Present, N = 15	15.2861	47.4159	53.4725
	Zupan [45], N = 7 Zupan [45], N = 15 a	15.74 15.74	47.15 47.15	53.43 53.43
	Bathe and Bolourchi [4]	15.9	47.2	53.4
	Simo [32]	15.79	47.23	53.37
	Number of elements = 8, N = number of internal points a Problems with convergence

Table 4

 4 Critical moment M c of the simply supported right-angle frame

	Element type	n e = 1	n e = 2	n e = 10
	Present, N = 1	±686.085	±626.876	±622.219
	Present, N = 2	±645.430	±624.547	±622.215
	Present, N = 4	±622.282	±622.211	±622.210
	Jelenić and Saje [13], N = 4	±622.2
	Nour-Omid and Rankin [21], mixed			±626.7
	Simo [32], N = 2			±626
	Zupan and Saje [40], N = 4	±622.40	±622.22
	Zupan [45], N = 5	±622.19	±622.22	±622.22
	Analytical solution ±622.210
	n e = number of elements, N = number of internal points

Acknowledgments This work was supported by the Slovenian Research Agency through the research program Z2-4158. The support is gratefully acknowledged.


Hongjuan Ran

Xianwu Luo

Dezhong Wang

Numerical Simulation of Instability Characteristic in Pump turbines

Keywords:

In classical pumps, there is single positive slope part on the curve of pump performance, named head-drop phenomenon. However, a second head-drop phenomenon is found in our experiments in the pump mode of some medium pump turbines. The flow rate of this second head-drop phenomenon is 0.68 -0.71 times of the design flow rate. As shown in our experimental results, it is significantly influenced by the complex vortices in the impeller inlet domain. Using SST k-ɷ turbulence model, numerical simulations are conducted to reveal its mechanism. After comparing numerical simulation results with the experimental results, it can be concluded that this second head-drop is closely related with the complex vortices in the impeller inlet. The complex vortices, due to the transmission flow between separated flow and back flow, is responsible for this second head-drop phenomenon.

INTRODUCTION

Head-drop phenomonon occurs widely in pumps and compressors, and it influences the system's security and reliability. Our paper studies a scaled model of some pump turbine, which is a key component in the hydraulic storage power plant, whose stability influences the flexibility and stability of grid net. The whole passage of pump turbine includes volute, diffuser, runner and draft tube, which are the same parts with that of general pumps. However, the angles of the diffuser component is adjustable in pump turbines. Hence, there are two head-drops on one performance curve, due to the adjustable angles of the diffuser part, whereas single head-drop is found on one curve of general pumps with the fixed angles of the diffuser part. Lots of researchers study the head-drop phenomenon with experimental and numerical methodology. Some conclude that rotating stall (near vane less domain between the diffure and the impeller) plays a key role in the head-drop phenomenon in pumps [1, 2, 9, 10, 11, 12,and 13], while other researchers state that the head-drops are closely tied with the flow separation near the runner shroud domain in other pumps [3,4 and 5]. Though lot of experiments study and CFD research are conducted, the mechanism of the head-drop are not clear. Based on the conclusion of our research, two different head-drops are observed on one curve of pump turbine through experimental study and numerical study [START_REF] Ran | Study on Double positive slopes phenomena in pump turbines[END_REF][START_REF] Ran Hongjuan | Experimental Study of the Pressure Fluctuations in a Pump Turbine at Large Partial Flow Conditions[END_REF]. This paper investigate the second head-drop phenomonon by methods of numerical simulations. Numerical results are compared with the experimental results. The results show that this second head-drop is caused by the complex vortices in the runner inlet domain. The complex vortices, due to transmission flow between seperated flow and back flow, includes high vortex flow pattern in the impeller inlet. Its detailed mechanism is partly revealed in this paper.

METHODS

Numerical simulation with commercial software CFX is conducted on the pump turbine model, based on the study of its experiments' results. The proper turbulence model is important for predicting the unsteady flow at the off-design conditions for the rotating machine [START_REF] Miyake | Full Simulation of a Flow in a Single Stage Axial-Rotor in Rotating Stall[END_REF] [14] [START_REF] Yexiang | Experimental and Numerical Analysis of Blade Channel Vortices in a Francis Turbine Runner[END_REF]. Numerical simulation comparison among five different rotating models and three different turbulence models are carried out, and their results are compared with experimental results individually, the Shear Stress Transport k-ω turbulence model is demonstrated to simulate this phenomenon more precisely than other models. Shear Stress Transport k-ω turbulence model is adopted in our paper. y + range is from 30 to 300, and the main area of fluid flow is 100. Since SST model accounts for the transport of the turbulent shear stress, it is more proper to simulate the flow separation with the consideration of adverse pressure gradients. It captures viscous sublyaer between the solid surface and fluid flow. In our case, based on our analysis of experiments results, the head-drop phenomena is related with separation flow near the impeller inlet domain, especially the cavitation vortices on the shroud surface, and these vortices are supposed to be related with viscous sublayer. Hence, SST k-ω is adopted.

Study model and boundary set up

The case study is a model of some pump turbine with 9 blades and 20 wicket gates, as shown in Figure 1. Experiments are carried out at test rig of swiss federal institute of technology in Lausanne according to IEC 60193 standard.

Grid independence check

Numerical simulation domain includes the whole passages of the pump turbine. Mesh number independent certifies as follows. Pump performance of three different meshes are simulated, and Table 1 shows the different mesh numbers. Table .2 displays the losses in different components. With the consideration of the computer sources and time consumption, mesh 2 is adopted in our simulation.

Results comparison

As shown in Figure 2, pump performance are carried out in test rig in Swiss Federal institute of Techonolgy in Lausanne. In Figure 2, there are two curves and each one is with two head-drop parts. These two curves are due to the hysterisis characteristic of head-drop phenoma. The lowest head of the second head-drop part is shown with the black arrow. This paper simulates this valley point to study the characteristic of the head-drop. Test H/HBEP Simulation /HBEP ε(%)

1.04 1.10 5.50

Table .3 demonstrates the head comparisons between numerical simulation results and experimental results. In the results of the numerical simulation, we choose their average head of different time steps as their compared head. Table .3 shows their error is 5.5%. Their error is caused by the complex flow pattern in the head drop phenomnon, including cavitation, noise and unstable flow. Two curves are with similar tendancy. Their differences are due to complex flow pattern with cavitation and noise. whole passage of pump. Figure 4 shows the work (or loss) of the different components varies with different runner steps. These components include a runner, the wicket gates, and the whole passage of the pump. The vertical coordinate represents the pressure, normalized by the total pressure of the operating point of the best efficiency. With the runner rotating, the pump head vibrates, and the work of the runner and the loss of the wicket gates verifies synchronically. They both cause head-drop phenomonon. It is concluded that the head (or loss) vibrations of the runner and the wicket gates both take responsible for the second head-drop phenomonon. This paper studies the flow pattern of the runner passages, since this flow pattern are also oberved by the experimets. As shown in Figure 5, the different performances between 4.9r and 5.7r (the time steps shown in Figure 4) are displayed. The work of the runner at 4.9r is higher than that in 5.73r (as shown in Figure 4). Briefly, the complex flow near runner inlet varies with the impeller rotating, and so do the runner performances, such as pressure and mass flow. Hence, the complex vortices induce the vibration of the impeller work, which further causes the head-drop.

Results Discussion

Components head (or loss) versus runner revolution

Unsteady flow pattern near the runner inlet domain

Complex vortices near the runner inlet domain

(a) Streamlines in several spans near the shroud domain

Figure 1 .

 1 Figure 1. Numerical simulation model with the whole passage of some pump turbine The uniform flow velocity profile is set at the inlet of calculation domain according to the mass flow. The average pressure at the domain outlet is set. Standard wall function is applied for the near wall, and non-slip condition is set for those solid walls. The simulated flow rate of the valley point is 0.694 Q/QBEP with wicket gate opening 26°, as described in Figure 2. The time step is set up as one degree of the impeller revolution.1.2 Grid independence checkNumerical simulation domain includes the whole passages of the pump turbine. Mesh number independent certifies as follows.

Figure 2 .

 2 Figure 2. Pump performance curves with two head-drop phenomena with 26 o of wicket gates ( 'Inc' denoted increasing flow rate from part load to full load during experimental procedure; and 'Dec' denoted decreasing flow rate from full load to part load)

Figure 3 .

 3 Figure 3. Pressure fluctuation comparison between numerical results and experimental results on the draft tube.

Figure 3

 3 Figure 3 shows pressure flucutation comparison on the draft tube near the runner inet. The x-axis represents the runner revolution, and the y-axis represents the relative amplitude of pressure fluctuation to the pump head of best efficiency point. Two curves are with similar tendancy. Their differences are

Figure 4 .

 4 Figure 4. Unsteady numerical simulation results: the distribution of work (or loss) in different passages: guide vanes, runner andwhole passage of pump. Figure4shows the work (or loss) of the different components varies with different runner steps. These components include a runner, the wicket gates, and the whole passage of the pump. The vertical coordinate represents the pressure, normalized by the total pressure of the operating point of the best efficiency. With the runner rotating, the pump head vibrates, and the work of the runner and the loss of the wicket gates verifies synchronically. They both cause head-drop phenomonon. It is concluded that the head (or loss) vibrations of the runner and the wicket gates both take responsible for the second head-drop phenomonon. This paper studies the flow pattern of the runner passages, since this flow pattern are also oberved by the experimets.3.2 Unsteady flow pattern near the runner inlet domain

  (a) The span value (runner inlet l=0.23; runner outlet l=0.50) and the streamline value (shroud: span =1; hub: span=0) (b). Mass flow contour (c). Pressure contour (d). Vortices contour Figure.5 Contours' comparisons between 4.90r and 5.73r (4.90r stands fo r the higher value of pump head, and 5.73r denoted the lower value of pump head, as shown in Figure 4) near the shroud domain of the runner

Figure 5 (

 5 a) shows the streamline value and the span value corresponded to the position of the runner. The span value of 0 is on the hub, while the span value of 1 is on the shroud. The streamline value from 0.25 to 0.5 corresponded to the runner from the inlet to the outlet. As shown in Figure5(b), Figure 5 (c) and Figure 5 (d), The horizontal coordinate shows the span value, and the vertical coordinate represents the streamline value, and these curves in these figures stand for the contour distributions of the mass flow, the pressure and the vortices individually. As the mass flow distributions of 4.90r and 5.73 r shown in Figure 5 (b), the backflow domain is represented by contour value of -0.01 and no separation domain is denoted by the contour value of 0.01, and contour value of 0 represents the totally blockage flow domain. The seperation domain (contour value less than 0.01 value) is larger in 5.73r than that in 4.90r. So does the totally blocked domain. However, the backflow domain (contour value is larger than 0.01 value) is larger in 4.9r than that in 5.73r. Shortly, the separation flow and blockage flow and complex vortices among them cause the diminish of the runner work in 5.73r. It's must be noted that the backflow seems not harmful for the runner work. The contour of Figure 5 (c) shows the pressure distribution comparison between 4.90r and 5.73r. From span value 0.90 to span value 0.97, The pressure value is less in 5.73 r than that in 4.90r during this range of span, corresponding to the dramatically unsteady flow pattern in Figure 5 (b).

Figure 5 (

 5 d) shows the Vortices distribution comparison between 4.90r and 5.73r. Near the runner leading edge (stream value from 0.25 to 0.28), the vortices domain of the same value is larger in 5.73 than that in 4.90r, corresponding to the unsteady flow pattern in Figure 5(b) and pressure distribution in Figure 5(c). The change of the complex vortices is significantly important for the head-drop phenomonon.

  Velocity near the shroud domain (span=0.96) Fig.6several flow patterns in runner passages near the inlet domain, And the time step is 4.9 runner revolution step (as shown in Figure4).

Figure. 6

 6 Figure. 6 shows different flow patterns occured in the runner inlet domain. They include the well flow, the separated flow, the complex vortices and the back flow.As shown in Figure6(a), there are several different flow patterns with the span value 0.90 -0.97. As shown in Figure6(a) with the span value 0.9, some vortex occur in the middle of the passage, near the suction surface of the impeller. The flow pattern are also similar with the span value 0.6 -0.9; At the span value 0.96, the flow becomes complex, which includes three parts: A part of the flow in one passage flows to the adjcent passage, and they are totally blocked, coressponding to Figure6(b). A part of the flow developes new vortex because they collapes with the leading edge of the baldes. Another part of these seperated flow generates a circulating flow, which flows around the circle band near the impeller inlet (as shown in Figure 6 (b)). and these complicated flow changes significantly with the different time steps (as shown in Figure 5). The flow is back flow near the shroud domain, and it is from the runner inlet to the draft tube cone, as shown in third image of Figure. 6 (a). Briefly, the complex vortices, with high value Vortices and high performance of instability, (as shown in Figure 6 (a)) are caused by the transition flow pattern from the initial separation flow to the backflow. It is concluded that the transition flow induces the second head-drop phenomonon. 4 Conclusions This paper study the head-drop phenomonon of a pump turbine. The flow rate of the head-drop is from 0.68 QBEP to 0.71 QBEP. Its valley point is 0.68 QBEP, which is as the simulated point. Pump head comparisons between numerical results and experimental results are carried out, and their error is 5%. The pressure fluctuation near the impeller inlet are also compared between numerical results and experimental results, and they are with similar tendency. The errors between numerical and experimental results due to the complex flow pattern with cavitation and noise, yet the numerical simulation doesn't consider the influence of cavitation and noise. The results of numerical simulation demonstrate that the transition flow in the runner inlet domain causes the head-drop phenomenon. The transition flow are transmission flow from seperated flow to back flow near the runner inlet domain. The transition flow is very complex. It includes the leakage flow from one passage to the adjecent passage, the blocked flow near the blade inlet, the high vortex flow near the blades' leading edges, the circulating flow near the runner inlet domain and little back flow near the shroud domain.

Table .

 .

				1 three different meshes
	items	Guide vane	runner	Stay vane	cone	Draft tube	Volute	total
	Mesh 1	19943	60243	17037 155700 259638 172773 685334
	mesh 2	132714 358420 188623 410172 653700 360396 2104025
	mesh 3	738888 1426414 826608 1814089 653700 1443028 6902727
			Table.2 Losses in different components
	cases	CPrunner	CPvolume	CPdiffuser	CPexcept
							volute
	mesh1	5735	-29		330.15	-333.27
	mesh 2	5740	-78.21	327.29	-346.68
	mesh 3	5662	-119.8	328.8	-334.29

Table . 3 Pump performance comparison

 .

Acknowledgments

The author really appreciates the help from the lab for hydraulic machine in EFPL, the Chinese natural science funding, Industrial funding provided by HEC, China.