The rotating rigid body model based on a non-twisting frame

This work proposes and investigates a new model of the rotating rigid body based on the non-twisting frame. Such a frame consists of three mutually orthogonal unit vectors whose rotation rate around one of the three axis remains zero at all times and thus, is represented by a nonholonomic restriction. Then, the corresponding Lagrange-D'Alembert equations are formulated by employing two descriptions, the first one relying on rotations and a splitting approach, and the second one relying on constrained directors. For vanishing external moments, we prove that the new model possesses conservation laws, i.e., the kinetic energy and two nonholonomic momenta that substantially differ from the holonomic momenta preserved by the standard rigid body model. Additionally, we propose a new specialization of a class of energy-momentum integration schemes that exactly preserves the kinetic energy and the nonholonomic momenta replicating the continuous counterpart. Finally, we present numerical results that show the excellent conservation properties as well as the accuracy for the time-discretized governing equations.


Introduction
The rigid body is a problem of classical mechanics that has been exhaustively studied (see, e.g., [1,2]). Its simplicity, as well as its relation with the (nonlinear) rotation group, makes of it the ideal setting to study abstract concepts of geometric mechanics, such as Poisson structures, reduction, symmetry, stability, etc. Additionally, many of the ideas that can be analyzed in the context of the rigid body can later be exploited in the study of nonlinear structural theories, such as rods and shells (as for example in [3,4,5]). As a result, geometrical integrators for the rigid body [6,7] are at the root of more complex numerical integration schemes for models that involve, in one way or another, rotations [8,9,10,11,12,13].
More specifically, the role of the rotation group is key because it is usually chosen to be the configuration space of the rigid body, when the latter has a fixed point. The rich Lie group structure of this set is responsible for much of the geometric theory of the rigid body, but it is not the only possible way to describe it. For example, the configuration of the a rigid body with a fixed point can also be described with three mutually orthogonal unit vectors. While this alternative description makes use of constraints, it has proven useful in the past for the construction of conserving numerical methods for rigid bodies, rods, and multibody systems [14,15,16].
In this article we explore a third route that can be followed to describe the kinematics of rigid bodies. This avenue is based on introducing a non-twisting or Bishop frame of motion [17]. This frame consists of three mutually orthogonal unit vectors whose rotation rate around one of the three axis remains zero at all times. Such a frame has proven useful to study the configuration of nonlinear Kirchhoff rods [18,19,20,21], but has not received enough attention in the context of the rigid body.
Formulating the equations of motion for the rigid body in the non-twisting frame demands a construction that is different from the standard one. In particular, the definition of Bishop's frame requires a constraint that is nonholonomic and does not admit a variational statement. Additionally, conservation laws take in this setting a different form when compared with the usual ones and whose identification is not straightforward, and the governing equations, i.e., the Lagrange-D'Alembert equations, are elegantly formulated by means of an splitting approach in terms of the covariant derivative on the unit sphere. Some of the conservation laws that take place under consideration of the non-twisitng frame may substantially differ from other nonholonomic cases that were investigated in the literature, e.g., [22,23]. In a pure holonomic context, some attempts to reformulate the dynamics on the unit sphere by means of advanced concepts from the differential geometry are to be found in [24,25]. However, the anisotropy of the inertial properties has been completely disregarded.
The rigid body equations with a nonholonomic constraint can be integrated in time with a conserving scheme that preserves energy and the newly identified momenta, the so-called nonholonomic momenta. This method, based on the idea of the average vector field [26,27] preserves remarkably these invariants of the motion, exactly, resulting in accurate pictures of the rigid body dynamics. The specialization of approaches based on more elaborated conservative/dissipative integration schemes like [28] is possible as well in this context, but falls outside the scope of the current work and therefore, not addressed here.
The rest of the article has the following structure. In Section 2, the fundamental concepts from the differential geometry are presented and discussed. Section 3 presents two derivations of the equations of motion for the standard rotating rigid body. The first set of equations is a split version of the well-known Euler equations that are presented within a novel geometrical framework that relies on covariant derivatives. The second one, then, produces a totally equivalent set of equations and relies on three constrained directors. Such an approach possesses a very favorable mathematical setting that will be exploited later on to derive an structure-preserving integration algorithm. In Section 4, the non-twisting condition is enforced for both fully equivalent formulations. Additionally, new conservation laws are identified in the continuous setting. Section 5 starts with the energy-momentum integration scheme for the director-based formulation and then it is modified for the nonholonomic case. The conservation properties are identified analytically for the discrete setting, replicating their continuous counterparts. In Section 6, numerical results that show the excellent performance of the new energy-momentum algorithm are presented and discussed. Finally, Section 7 closes this work with a summary.

Relevant geometrical concepts
The governing equations of the rigid body are posed on nonlinear manifolds. We briefly summarize the essential geometrical concepts required for a complete description of the model (see, e.g., [29,25] for more comprehensive expositions).

The unit sphere
The unit sphere S 2 is a nonlinear, compact, two-dimensional manifold that often appears in the configuration spaces of solid mechanics, be it the rigid body, rods, or shells [30,31,32,33,13]. As an embedded set on Euclidean space, it is defined as where the dot product is the standard one in R 3 . The tangent bundle of the unit 2-sphere is also a manifold defined as Alternatively, tangent vectors of S 2 to a point d are those defined by the expressions: where the product "×" is the standard cross product in R 3 .
In contrast with the space of rotations, to be studied in more detail later, the unit sphere does not have a group structure, but instead it has that of a Riemannian manifold. The connection of this set can be more easily explained by considering it to be an embedded manifold in R 3 . As such, the covariant derivative of a smooth vector field v : S 2 → T S 2 along a vector field w : S 2 → T S 2 is the vector field ∇ w v ∈ T S 2 that evaluated at d is precisely the projection of the derivative Dv in the direction of w onto the tangent plane to d. Hence, denoting as I the unit second order tensor and ⊗ the outer product between vectors, both on R 3 , this projection can be simply written as When d : (a, b) → S 2 is a smooth one-parameter curve on the unit sphere and d ′ its derivative, the covariant derivative of a smooth vector field v : S 2 → T S 2 in the direction of d ′ has an expression that follows from Eq. (4), that is, which, as before, is nothing but the projection of (v • d) ′ onto the tangent space T d S 2 , and the symbol "•" stands for composition. The covariant derivative allows to compare two tangent vectors belonging to different tangent spaces. By means of the parallel transport, one vector can be transported from its tangent space to the space of the other one. Then, all the desired comparisons can be made over objects belonging to the same space. The exponential map exp : T d0 S 2 → S 2 is a surjective application with a closed form expression given by the formula where v 0 must be a tangent vector on T d0 S 2 and | · | denotes the Euclidean vector norm. The inverse of the exponential function is the logarithmic map log : S 2 → T d0 S 2 , for which also there is a closed form expression that reads is a great arch on the sphere obtained as the solution of the equation with the explicit form d G (s) = cos(|v 0 |s)d 0 + sin(|v 0 |s) v 0 |v 0 | .
The parallel transport of w 0 ∈ T d0 S 2 → w ∈ T d S 2 along the geodesic d G is then given by and verifies More details about the expressions (6) to (10) can be found in [34,35]. Given the definitions (1) and (2) of the unit sphere and its tangent bundle, we recognize that there exists an isomorphism for any d ∈ S 2 . Given now two directors d,d, we say that a second order tensor T : R 3 → R 3 splits from d tod if it can be written in the form where T ⊥ is a bijection from T d S 2 to TdS 2 with ker(T ⊥ ) = span(d), and T is a bijection from span(d) to span(d) with ker(T ) ≡ T d S 2 . The split (13) depends on the pair d,d but it is not indicated explicitly in order to simplify the notation.

The special orthogonal group
Classical descriptions of rigid body kinematics are invariably based on the definition of their configuration space as the set of proper orthogonal second order tensors, that is, the special orthogonal group, defined as This smooth manifold has a group-like structure when considered with the tensor multiplication operation, thus it is a Lie group. Its associated Lie algebra is the linear set Later, it will be convenient to exploit the isomorphism that exists between three-dimensional real vectors and so (3). To see this, consider the map(·) : R 3 → so(3) such that for all w, a ∈ R 3 , the tensorŵ ∈ so(3) satisfiesŵa = w × a. Here, w is referred to as the axial vector ofŵ, which is a skew-symmetric tensor, and we also write skew[w] =ŵ. Being a Lie group, the space of rotations has an exponential map exp : so(3) → SO(3) whose closed form expression is Rodrigues' formula exp[θ] := I + sin θ θθ + 1 2 with θ ∈ R 3 , θ = |θ|. The linearization of the exponential map is simplified by introducing the map dexp : for every θ : R → R 3 . A closed-form expression for this map, as well as more details regarding the numerical treatment of rotations can be found elsewhere [36,7,13].

The notion of twist and the non-twisting frame
Let where the overdot refers to the derivative with respect to time. Given the initial value of the triad transforming it to the frame {d 1 , d 2 , d 3 } that evolves without twist, and whose closed form is The non-twisting frame has Darboux vector and it is related to parallel transport in the sphere. To see this relation, consider again the same non-twisting frame as before. Then, we recall that a vector field v ∈ T d3 S 2 is said to be parallel-transported along d 3 if and only if ∇ vḋ3 = 0. An consequence of this is thaṫ and similarlyḋ In addition, we have thaṫ which is merely the product of the angular velocity components and can be interpreted as a scalar curvature.
To define precisely the concept of twist, let us consider the rotation exp[ψd 3 ], with ψ : [0, T ] → S 1 (the unit 1-sphere), and the rotated triad and the Darboux vector of the rotated triad is or, equivalently, In this last expression, we identify the twist rateψ and the twist angle ψ, respectively, as the rotation velocity component on the d 3 direction and the rotation angle about this same vector (for further details, see [17,37]). The previous calculations show that the frame {d 1 , d 2 , d 3 }known as the natural or Bishop frame in the context of one-parameter curves embedded in the ambient space -is the unique one obtained by transporting {D 1 , D 2 , D 3 } without twist. In this context, a summary of recent advances and open problems are presented for instance in [38].

Standard rotating rigid body
In this section, we review the classical rotating rigid body model, which we take as the starting point for our developments. We present this model, however, in an unusual fashion. In it, the governing equations of the body appear in split form. This refers to the fact that, for a given director d 3 , the dynamics of the body that takes place in the cotangent space T * d3 S 2 is separated from that one corresponding to the reciprocal normal space N * d3 S 2 ≡ span(d 3 ). In order to do this, we have employed the identifications

Kinematic description
As customary, a rotating rigid body is defined to be a three-dimensional non-deformable body.
The state of such a body, when one of its points is fixed, can be described by a rotating frame whose orientation is given by a rotation tensor. Thus, the configuration manifold is Q ≡ SO(3).
Let us now study the motion of a rotating rigid body, that is, a time-parameterized curve in configuration space Λ : [0, T ] → Q. The generalized velocity of the rotating rigid body belongs, for every t ∈ [0, T ], to the tangent bundle The time derivative of the rotation tensor can be written aṡ where w and W are the spatial and convected angular velocities, respectively. (27) to split the rotation vectors as in Then, using the relations (29), we identify

Kinetic energy and angular momentum
Let us now select a fixed Cartesian basis of R 3 denoted as where the third vector coincides with one of the principal directions of the convected inertia tensor J : R 3 → R 3 of the body, a symmetric, second order, positive definite tensor. Thus, this tensor splits from D 3 to D 3 and we write where J ⊥ maps bijectively span(D 1 , D 2 ) onto itself and satisfies J ⊥ D 3 = 0. The kinetic energy of a rigid body with a fixed point is defined as the quadratic form where j is the spatial inertia tensor, the push-forward of the convected inertia, and defined as Let d 3 = ΛD 3 . Given the relationship between the convected and spatial inertia, it follows that the latter also splits, this time from d 3 to d 3 , and thus where now j ⊥ maps bijectively span(d 1 , d 2 ) onto itself and satisfies j ⊥ d 3 = 0.
As a consequence of the structure of the inertia tensor, the kinetic energy of a rotating rigid body can be written in either of the following equivalent ways: The angular momentum of the rotating rigid body is conjugate to the angular velocity as in and we note that we can introduce a convected version of the momentum π by pulling it back with the rotating tensor and defining Due to the particular structure of the inertia, the momentum can also be split, as before, as in with

Variations of the motion rates
The governing equations of the rigid body will be obtained using Hamilton's principle of stationary action, using calculus of variations. We gather next some results that will prove necessary for the computation of the functional derivatives and, later, for the linearization of the model. To introduce these concepts, let us consider a curve of configurations Λ ι (t) parametrized by the scalar ι and given by where δθ : [0, T ] → R 3 represents all arbitrary variations that satisfy The curve Λ ι passes through the configuration Λ when ι = 0 and has tangent at this point For future reference, let us calculate the variation of the derivativeΛ. To do so, let us first define the temporal derivative of the perturbed rotation, that is, Then, the variation ofΛ is just With the previous results at hand, we can now proceed to calculate the variations of the convected angular velocities, as summarized in the following theorem.
Proof. The convected angular velocities of the one-parameter curve of configurations Λ ι are where d i,ι = Λ ι D i . The variation of the angular velocity perpendicular to D 3 is obtained from its definition employing some algebraic manipulations and expression (45) as follows: The variation of the angular velocity parallel to D 3 follows similar steps:

Governing equations and invariants
Here, we derive the governing equations of the rotating rigid body model and the concomitant conservation laws. Hamilton's principle of stationary action states that the governing equations are the Euler-Lagrange equations of the action functional with unknown fields (Λ,Λ) ∈ T Q.
Theorem 3.2. The equations of motion, i.e., the Euler-Lagrange equations, for the standard rotating rigid body model in split form are: The pertaining initial conditions are: Proof. The theorem follows from the systematic calculation of δS, the variation of the action, based on the variation of the convected angular velocities of Eq. (46), thus we omit a detailed derivation. Eq. (51), which is presented in its split form, is equivalent to Euler's equations which state that the spatial angular momentum is preserved, i.e.,π = 0. This is easily proven as follows Theorem 3.3. The conservation laws of the rotating rigid body are: Proof. This is an standard result and thus, we omit further details.
Remark. To include external moments acting on the standard rotating rigid body it is necessary to calculate the associated virtual work as follows and add this contribution to the variation of the action.

Model equations based on directors
Here, we present an alternative set of governing equations for the rotating rigid body model that will be used later to formulate a structure preserving algorithm. For this purpose, let us define the following configuration space whose tangent space at the point q is given by Now, we start by defining the rigid body as the bounded set B 0 ⊂ R 3 of points where (θ 1 , θ 2 , θ 3 ) are the material coordinates of the point and {D i } 3 i=1 are three orthogonal directors, with the third one oriented in the direction of the principal axis of inertia and such that D 3 = D 1 × D 2 . The position of the point X at time t ∈ [0, T ] is denoted as x(t) ∈ R 3 and given by with (d 1 , d 2 , d 3 ) = q ∈ Q, for all t. On this basis, there must be a rotation tensor Λ(t) = d i (t) ⊗ D i , where we have employed the sum convention for repeated indices, such that d i (t) = Λ(t)D i . The material velocity of the particle X is the vectorẋ(t) ∈ R 3 that can be written aṡ with (ḋ 1 ,ḋ 2 ,ḋ 3 ) =q ∈ T q Q representing three director velocity vectors.
To construct the dynamic equations of the model, assume the body B 0 has a density ρ 0 per unit of reference volume and hence its total kinetic energy, or Lagrangian, can be formulated as To employ Hamilton's principle of stationary action, but restricting the body directors to remain orthonormal at all time, we define the constrained action where K is given by Eq. (61), λ ∈ R 3 is a vector of Lagrange multipliers, and h is of the form such that h(d 1 , d 2 , d 3 ) = 0 expresses the directors' orthonormality.
Theorem 3.4. The alternative equations of motion, i.e., the Euler-Lagrange equations, for the standard rotating rigid body model are: The generalized momenta (π 1 , π 2 , π 3 ) = p ∈ T * q Q are defined as where Euler's inertia coefficients are for i and j from 1 to 3. In addition, the splitting of the inertia tensor implies J 13 = J 23 = 0 and H i ∈ L(T di S 2 , R n ) stands for ∂h ∂di .
The pertaining initial conditions are: Proof. The theorem follows from the systematic calculation of δS.
The reparametrized equations presented above are totally equivalent to the commonly used equations for the standard rotating rigid body. Consequently, the conservation laws described previously apply directly to this equivalent model. For a in-depth discussion on this subject, the reader may consult [39].
Remark. As before, to include external moments acting on the standard rotating rigid body, the following additional terms need to be added to the variation of the action 4 Rotating rigid body based on the non-twisting frame In this section, we introduce the nonholonomic rotating rigid body, which incorporates the nonintegrable constraint that is necessary to set the non-twisting frame according to Eq. (18). This is a non-variational model, since it cannot be derived directly from a variational principle. For this purpose, we modify Eq. (51) to account the non-integrable condition W = w = 0 according to the usual nonholonomic approach. We also introduce the concomitant conservation laws. Additionally, we present an alternative formulation that relies on constrained directors, whose particular mathematical structure enables the application of structure-preserving integration schemes.

Governing equations and invariants
Theorem 4.1. The nonholonomic equations of motion, i.e., Lagrange-D'Alembert equations, for the rotating body model based on the non-twisting frame are: The pertaining initial conditions are: Moreover, Eq. (69) can be rewritten as where π ⊥ ∈ T * d3 S 2 must satisfy the parallel transport along the curve d 3 ∈ S 2 . Proof. The first part follows from the inclusion of the force associated to the presence of the nonholonomic restriction given by which ensures that the rotating frame renders no twist at all. The virtual work performed by the force associated to the presence of this nonholonomic restriction can be computed as where µ ∈ R denotes the corresponding Lagrange multiplier. The second part follows from noticing that w = w · d 3 = 0 implies π = 0.
Theorem 4.2. The conservation laws of the rotating rigid body based on the non-twisting frame are: Proof. To prove the conservation of kinetic energy, let us consider the following equilibrium statement where δθ ⊥ ∈ T d3 S 2 , δθ ∈ N d3 S 2 and δµ are admissible variations. Now by choosing δθ ⊥ = w ⊥ , δθ = w and δµ = 0, we have that which shows that the kinetic energy is preserved by the motion.
To prove the conservation of the first and second components of the material angular momentum, let us consider the fact that Now by introducing the former expression into the first statement of Eq. (71), we have that in which the parallel transport of d 1 and d 2 , both in T d3 S 2 , has been accounted for. This shows that the first and second components of the angular momentum are preserved by the motion.
Remark. To include external moments acting on the rotating rigid body based on the nontwisting frame, it is necessary to compute the associated virtual work as follows

Alternative governing equations
Here, we present an alternative formulation for the rotating rigid body based on the non-twisting frame that relies on constrained directors. The extension of the standard rotating rigid body model to the one relying on the non-twisting frame requires the introduction of the constraint given by in which a ∈ [0, 1] is a parameter that can be freely chosen for convenience. This will be used later on for the proof of the conservation properties of the specialized structure preserving algorithm.
The pertaining initial conditions are: Proof. The theorem follows from the computation of the virtual work associated to the presence of the nonholonomic restriction, namely where µ ∈ R denotes the corresponding Lagrange multiplier.
Remark. To include external moments acting on the rotating rigid body based on the nontwisting frame, it is necessary to compute the associated virtual work as follows

Structure-preserving time integration
A fundamental aspect to produce acceptable numerical results in the context of nonlinear systems is the preservation of mechanical invariants whenever possible, e.g., first integrals of motion. These conservation properties ensure that beyond the approximation errors, the computed solution remains consistent with respect to the underlying physical essence. Here then, we chose the family of integration methods that is derived by direct discretization of the equations of motion.

Basic energy-momentum algorithm
Next, we describe the application of the energy-momentum integration algorithm [6,40] to the "standard rotating rigid body" case. For this purpose, the following nomenclature is necessary: While q is the vector of generalized coordinates, p collects the generalized momenta and Q ext contains the generalized external loads, if present. The discrete version of Eq. (64) can be expressed at time n + 1 2 as where ·, · stands for the dual pairing.
A key point to achieve the desired preservation properties, is to define the momentum terms by using the midpoint rule, i.e., where q n andq n are known from the previous step, q n+1 areq n+1 are unknown, andq n+1 is computed as 2 h (q n+1 − q n ) −q n once q n+1 has been determined by means of an iterative procedure, typically the Newtown-Raphson method.
The mass matrix takes the form and J ij for i and j running from 1 to 3 being defined above. This very simple construction satisfies, only for the standard rigid body case, the preservation of linear and angular momenta in combination with the kinetic energy in absence of external loads.
The discrete version of the Jacobian matrix of the constraints can be computed with the average vector field [41,42] as where q(ξ) is defined as 1 2 (1 − ξ)q n + 1 2 (1 + ξ)q n+1 for ξ ∈ [−1, +1]. The algorithmic Jacobian matrix defined in this way satifies for any admissible solution the discrete version of the hidden constraints, i.e., H d (q n , q n+1 )(q n+1 − q n ) = 0. (90) Theorem 5.1. The discrete conservation laws of the energy-momentum integration algorithm specialized to the standard rotating rigid body are: Proof. This is an standard result and thus, we omit further details.

Specialized energy-momentum algorithm
The energy-momentum integration algorithm can be further specialized to the nonholonomic case, where the discrete governing equations are: Once again, we can use the average vector field to compute that arises from the nonholonomic constraint, where q(ξ) is defined as before. In this way, the nonholonomic constraint is identically satisfied at the midpoint, i.e., Theorem 5.2. The discrete conservation laws of the energy-momentum integration algorithm specialized to the rotating rigid body based on the non-twisting frame are: Proof. To prove the conservation of kinetic energy, let us consider the following discrete variation By inserting the previous discrete variation in Eq. (92), we get For the first component of the angular momentum, i.e., Π 1 , we need to consider the following discrete variation (δq n+ 1 2 , δλ n+1 ) = (δd 1,n+ 1 2 , δd 2,n+ 1 2 , δd 3, and let a be equal to 1. By inserting the previous discrete variation in Eq. (92), we get where By using Taylor's approximations, we have that and Then insomuch asḋ 3,n+1 · π 2 n +ḋ 3,n · π 2 n+1 = 2ḋ 3,n+ 1 2 · π 2 which can be easily shown by considering thaṫ since the angular velocity w n+ 1 2 has the form d 3,n+ 1 2 ×ḋ 3,n+ 1 2 due to the satisfaction of the non-twisting condition, see Subsection 2.3.
In this way, the first term of Eq. (99) becomes and with same reasoning, the second term of Eq. (99) turns to be By replacing the two previous expressions in (99), we have that which is true since Finally, for the second component of the angular momentum, i.e., Π 2 , we need to consider the following discrete variation and let a be equal to zero. Then, the rest of the proof follows as before.

Numerical results
In this section, we present numerical results of the motion of a rotating rigid body based on the non-twisting frame with (non-physical) inertia Next, we evaluate the qualitative properties of the proposed numerical setting in a reduced picture. For the first case, we consider the dynamic response to an initial condition different from the trivial one. For the second case, we consider the dynamic response to a vanishing load. The third and last case is a combination of both, i.e., initial condition different from the trivial one and a vanishing load. All the three cases were numerically solved in the time interval [0, 5] s with a time step size of h = 0.005 s and relative tolerance 10 −10 .

Case 1 -response to nonzero initial conditions
For this first case we consider and m ext ⊥ (t) = 0 Kg m 2 /s 2 .
(113) Figure 1 presents the time history for the spatial and material components of the angular momentum. On the left, we can observe that the components of the spatial angular momentum (SAM) oscillate with constant amplitude and frequency and therefore, they are not constant as in the case of the standard rotating rigid body. On the right we can observe that the components of the material angular momentum (MAM) are identically preserved. While the first and second components are constant and different from zero, the third one is zero as expected from the imposition of the nonholonomic restriction ω = 0 rad/s. As shown before for the analytical setting as well as for the numerical setting, this non-intuitive behavior results from the fact that the dynamics of the system is not taking place in the environment space, but on the 2-sphere. Therefore, this behavior is truly native on the 2-sphere since the directors d 1 and d 2 in T d3 S 2 are being parallel transported along the time-parameterized solution curve d 3 . Figure 2 presents the time history for the kinetic energy and the second quotient of precision as defined in the Appendix. On the left we can observe that the kinetic energy is identically preserved as expected. On the right we see that the second quotient of precision is approximately 4, see also Table 1, which means that the integrator is really achieving second-order accuracy. Figure 3 shows the trajectory followed by d 3 , which as expected takes place on a plane that separates the sphere into two half spheres. Such trajectory minimizes locally the distance on S 2 and thus, this is geodesic. Finally, and to summarize the excellent performance of the numerical setting, Table 2 presents the stationary values for the motion invariants, i.e., the first and second components of the material angular momentum and kinetic energy.

Case 2 -response to a vanishing load
For this second case we consider and where Figure 4 presents the time history for the spatial and material components of the angular momentum, where the applied material load m ext ⊥ is active only during the first second of simulation. On the left figure, we observe that the components of the spatial angular momentum vary starting from zero since the rotating rigid body is initially at rest. After the load vanishes, the components of the spatial angular momentum oscillate with constant amplitude and frequency and therefore, they are not constant, but indicate a steady state. On the right figure we can observe that the components of the material angular momentum also vary from zero, except the third one that remains always equal to zero. After the material load vanishes, the components of the material angular momenta are identically preserved. Once again, the first and second components are constant and different from zero. Figure 5 presents the time history for the kinetic energy and the second quotient of precision. On the left, we can observe that the kinetic energy vary during the first second, where the applied material load is active. After this vanishes, the kinetic energy is identically preserved. On the right figure we confirm again that the second quotient of precision is approximately 4, see also Table 3, which means that the integrator is second-order accurate even during the time in which the applied material load is active. Figure 6 shows the trajectory followed by d 3 , which due to the fixed relation among components of the applied material load takes place on a plane that separates the sphere in two half spheres. Such trajectory minimizes locally the distance on S 2 and thus, this is geodesic as well. Table 4 presents the stationary values for the motion invariants.

Case 3 -response to nonzero initial conditions and a vanishing load
with f (t) defined as in Eq. (116). Figure 7 presents the time history for the spatial and material components of the angular momentum. On the left, we can observe that the components of the spatial angular momentum vary starting from the values corresponding to the initial condition adopted. After the material load vanishes, the components of the spatial angular momentum oscillate with constant amplitude and frequency indicating a steady state. On the right, we observe that the components of the material angular momentum also vary from the values corresponding to the initial condition adopted, except the third one that remains always equal to zero. After the material load vanishes the components of the material angular momenta are identically preserved as in the previous cases. Figure 8 presents the time history for the kinetic energy and the second quotient of precision. On the left we can observe that the kinetic energy vary during the first second, where the applied material load is active. After this vanishes, the kinetic energy is identically preserved. To the right we see that the second quotient of precision is approximately 4, see also Table 5. Figure 9 shows the trajectory followed by d 3 . During the first second, the trajectory does not render a distance minimizing curve on S 2 . This is due to the combination of initial condition adopted and the load applied that produces a change in the direction of the axis of rotation. After the material load vanishes, the trajectory describes a circle of radius 1. Table 6 provides the stationary values for the motion invariants.

Summary
This article describes the governing equations of the rotating rigid body in a nonholonomic context, and discusses their relation with other, well-known, equivalent models based on rotations and orthonormal vectors. The equations obtained are non-variational and possess first invariants of motion. Some of them, i.e., the nonholonomic momenta (first and second components of the material angular momentum), are neither evident from the standard descriptions nor intuitive. To the best of our knowledge, there is no work in the literature that reports similar observations and thus, it represents a main innovation of the current work.
Complementing the rigorous mathematical analyis done for the proposed model, an implicit, second order accurate, energy and momentum conserving algorithm is presented, which discretizes in time the rigid body, nonholomic equations. Such a time integration scheme preserves exactly the energy and nonholonomic momenta and thus, this represents also a main innovation of the current work. Finally, simple examples, which make use of all elements of the approach proposed, are provided and confirm the excellent conservation properties and second order accuracy of the new scheme.

A Second quotient of precision
To investigate the correctness of the integration scheme, we can employ the second quotient of precision [43]. Any numerical solution of an initial value problem can be written as with ξ(t) representing the exact solution of the initial value problem under consideration and ψ i for i = 1, . . . , n representing smooth functions that only depend on the time parameter t. h stands for the time step and k is a positive integer number that enables defining finer solutions computed from the original resolution.
Let us introduce the second quotient of precision given by Q II (t) = ξ(t, h, 1) − ξ(t, h, 2) ξ(t, h, 2) − ξ(t, h, 4) . (120) For the numerator, we have that If we have time steps that are sufficiently small, it can be showed that the second quotient of precision is very close to take the value 2 n , i.e., In this work, we consider an integration algorithm whose accuracy is warranted to be of second order and thus, Q II (t) ≈ 4. Be aware that for the computation of precision quotients, h needs to be chosen small enough. Such a selection strongly depends on the case considered. Moreover, when ψ n (t) is very small, the test may fail even if correctly implemented. Therefore, we recommend to play with initial conditions and time resolutions for achieving correct pictures.