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

## Authors

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

The rigid body is a problem of classical mechanics that has been exhaustively studied (see, e.g., [16, 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 [40, 1, 28]). As a result, geometrical integrators for the rigid body [43, 36] are at the root of more complex numerical integration schemes for models that involve, in one way or another, rotations [41, 37, 20, 33, 7, 30].

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 [32, 6, 5].

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 . 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 [38, 27, 3, 34], 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., [8, 18]. 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 [23, 24]. 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, 10] 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  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.

## 2 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., [25, 24] for more comprehensive expositions).

### 2.1 The unit sphere

The unit sphere 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 [11, 39, 35, 31, 30]. As an embedded set on Euclidean space, it is defined as

 S2:={d∈R3∣d⋅d=1}, (1)

where the dot product is the standard one in . The tangent bundle of the unit 2-sphere is also a manifold defined as

 TS2:={(d,v), d∈S2, v∈R3, d⋅v=0}. (2)

Alternatively, tangent vectors of to a point are those defined by the expressions:

 v=w×d,withw⋅d=0 , (3)

where the product “” is the standard cross product in .

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 . As such, the covariant derivative of a smooth vector field along a vector field is the vector field that evaluated at is precisely the projection of the derivative in the direction of onto the tangent plane to . Hence, denoting as

the unit second order tensor and

the outer product between vectors, both on , this projection can be simply written as

 ∇wv:=(I−d⊗d)Dv⋅w . (4)

When is a smooth one-parameter curve on the unit sphere and its derivative, the covariant derivative of a smooth vector field in the direction of has an expression that follows from Eq. (4), that is,

 ∇d′v=(I−d⊗d)Dv⋅d′=(v∘d)′−((v∘d)′⋅d)d, (5)

which, as before, is nothing but the projection of onto the tangent space , 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 is a surjective application with a closed form expression given by the formula

 expd0(v0)=cos(|v0|)d0+sin(|v0|)v0|v0|, (6)

where must be a tangent vector on and denotes the Euclidean vector norm. The inverse of the exponential function is the logarithmic map , for which also there is a closed form expression that reads

 logd0(d)=arccos(d0⋅d)(I−d0⊗d0)d|(I−d0⊗d0)d|, (7)

with . The geodesic for with and is a great arch on the sphere obtained as the solution of the equation

 ∇d′d′=0, (8)

with the explicit form

 dG(s)=cos(|v0|s)d0+sin(|v0|s)v0|v0|. (9)

The parallel transport of along the geodesic is then given by

 w=(I−1arccos2(d0⋅dG)(logd0dG+logdGd0)⊗logd0dG)w0, (10)

and verifies

 ∇wd′=(I−d⊗d)w′=0. (11)

More details about the expressions (6) to (10) can be found in [19, 4].

Given the definitions (1) and (2) of the unit sphere and its tangent bundle, we recognize that there exists an isomorphism

 R3≅TdS2⊕span(d), (12)

for any . Given now two directors , we say that a second order tensor splits from to if it can be written in the form

 T=T⊥+T∥, (13)

where is a bijection from to with , and is a bijection from to with . The split (13) depends on the pair  but it is not indicated explicitly in order to simplify the notation.

### 2.2 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

 SO(3):={Λ∈R3×3, ΛTΛ=I, detΛ=+1}. (14)

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

 so(3):={^w∈R3×3, ^w=−^wT}. (15)

Later, it will be convenient to exploit the isomorphism that exists between three-dimensional real vectors and . To see this, consider the map such that for all , the tensor satisfies . Here, is referred to as the axial vector of

, which is a skew-symmetric tensor, and we also write

. Being a Lie group, the space of rotations has an exponential map whose closed form expression is Rodrigues’ formula

 exp[^θ]:=I+sinθθ^θ+12sin2(θ/2)(θ/2)2^θ2 , (16)

with , . The linearization of the exponential map is simplified by introducing the map that satisfies

 (17)

for every . A closed-form expression for this map, as well as more details regarding the numerical treatment of rotations can be found elsewhere [17, 36, 30].

### 2.3 The notion of twist and the non-twisting frame

Let indicate a curve of directors parameterized by time . Now, let us consider two other curves such that are mutually orthogonal for all . We say that the triad moves without twist if

 ˙d1⋅d2=˙d2⋅d1=0 , (18)

where the overdot refers to the derivative with respect to time. Given the initial value of the triad , there is a map transforming it to the frame that evolves without twist, and whose closed form is

 χ=d1⊗D1+d2⊗D2+d3⊗D3. (19)

The non-twisting frame has Darboux vector

 wχ=d3×˙d3, (20)

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 is said to be parallel-transported along if and only if . An consequence of this is that

 ˙d3⋅˙d1=˙d3⋅(wχ×d1)=0, (21)

and similarly

 ˙d3⋅˙d2=˙d3⋅(wχ×d2)=0. (22)

 ˙d1⋅˙d2=(wχ×d1)⋅(wχ×d2)=(d1⋅˙d3)(d2⋅˙d3), (23)

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 , with (the unit 1-sphere), and the rotated triad . Then,

 dψ1=cos(ψ)d1+sin(ψ)d2anddψ2=−sin(ψ)d1+cos(ψ)d1, (24)

and the Darboux vector of the rotated triad is

 wexp[ψ^d3]χ=d3×˙d3+˙ψd3, (25)

or, equivalently,

 wexp[ψ^d3]χ=−(dψ2⋅˙d3)dψ1+(dψ1⋅˙d3)dψ2+˙ψd3. (26)

In this last expression, we identify the twist rate and the twist angle , respectively, as the rotation velocity component on the direction and the rotation angle about this same vector (for further details, see [9, 22]). The previous calculations show that the frame — 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 without twist. In this context, a summary of recent advances and open problems are presented for instance in .

## 3 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 , the dynamics of the body that takes place in the cotangent space is separated from that one corresponding to the reciprocal normal space . In order to do this, we have employed the identifications

 R3≅Td3S2⊕span(d3)≅T∗d3S2⊕span(d3). (27)

### 3.1 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 .

Let us now study the motion of a rotating rigid body, that is, a time-parameterized curve in configuration space . The generalized velocity of the rotating rigid body belongs, for every , to the tangent bundle

 TQ:={(Λ,˙Λ),Λ∈SO(3),ΛT˙Λ∈so(3)}. (28)

The time derivative of the rotation tensor can be written as

 ˙Λ=^wΛ=ΛˆW , (29)

where and are the spatial and convected angular velocities, respectively.

Let be a fixed basis of the ambient space. Then, if , with we can use Eq. (27) to split the rotation vectors as in

 w=w⊥+w∥d3 ,W=W⊥+W∥d3 . (30)

Then, using the relations (29), we identify

 w⊥=d3×˙d3 ,W⊥=ΛTw⊥ . (31)

### 3.2 Kinetic energy and angular momentum

Let us now select a fixed Cartesian basis of denoted as , where the third vector coincides with one of the principal directions of the convected inertia tensor of the body, a symmetric, second order, positive definite tensor. Thus, this tensor splits from to and we write

 J=J⊥+J∥D3⊗D3 , (32)

where maps bijectively onto itself and satisfies .

The kinetic energy of a rigid body with a fixed point is defined as the quadratic form

 K:=12W⋅JW=12w⋅jw , (33)

where is the spatial inertia tensor, the push-forward of the convected inertia, and defined as

 j:=ΛJΛT . (34)

Let . Given the relationship between the convected and spatial inertia, it follows that the latter also splits, this time from to , and thus

 j=j⊥+j∥d3⊗d3 , (35)

where now maps bijectively onto itself and satisfies .

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:

 12w⋅jw=12W⋅JW=12W⊥⋅J⊥W⊥+12J∥W2∥=12w⊥⋅j⊥w⊥+12j∥w2∥. (36)

The angular momentum of the rotating rigid body is conjugate to the angular velocity as in

 π:=∂K∂w=jw, (37)

and we note that we can introduce a convected version of the momentum by pulling it back with the rotating tensor and defining

 Π:=ΛTπ=∂K∂W . (38)

Due to the particular structure of the inertia, the momentum can also be split, as before, as in

 π=π⊥+π∥d3 ,Π=Π⊥+Π∥d3, (39)

with

 π⊥=j⊥w⊥ ,Π⊥=J⊥W⊥,π∥=Π∥=j∥w∥=J∥W∥ . (40)

### 3.3 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 parametrized by the scalar and given by

 Λι(t)=exp[ιˆδθ(t)]Λ(t), (41)

where represents all arbitrary variations that satisfy

 δθ(0)=δθ(T)=0 . (42)

The curve passes through the configuration when and has tangent at this point

 ∂∂t∣∣∣ι=0Λι=ˆδθΛ. (43)

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,

 ∂∂tΛι=∂∂texp[ιˆδθ]Λ=skew[dexp[ιˆδθ]ιδ˙θ]Λ+exp[ιˆδθ]˙Λ . (44)

Then, the variation of is just

 (45)

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.

###### Theorem 3.1.

The variations of the convected angular velocities are

 δW⊥ =ΛT(I−d3⊗d3)δ˙θ, (46a) δW∥ =d3⋅δ˙θ. (46b)
###### Proof.

The convected angular velocities of the one-parameter curve of configurations are

 W⊥,ι=D3×(ΛTι˙d3,ι)andW∥,ι=d2,ι⋅˙d1,ι, (47)

where . The variation of the angular velocity perpendicular to is obtained from its definition employing some algebraic manipulations and expression (45) as follows:

 δW⊥=∂∂ι∣∣∣ι=0(D3×(ΛTι˙d3,ι))=D3×(δΛT˙ΛD3+ΛTδ˙ΛD3)=ΛT(d3×(ˆδ˙θ×d3))=ΛT(δ˙θ−(δ˙θ⋅d3)d3)=ΛT(I−d3⊗d3)δ˙θ . (48)

The variation of the angular velocity parallel to follows similar steps:

 δW∥=∂∂ι∣∣∣ι=0(d2,ι⋅˙d1,ι)=ˆδθΛD2⋅˙ΛD1+ΛD2⋅(ˆδ˙θΛ+ˆδθ˙Λ)D1=d1×d2⋅δ˙θ=d3⋅δ˙θ . (49)

### 3.4 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

 S=∫T0Kdt, (50)

with unknown fields .

###### Theorem 3.2.

The equations of motion, i.e., the Euler-Lagrange equations, for the standard rotating rigid body model in split form are:

 ∇˙d3π⊥+π∥˙d3 =0, (51a) ˙π∥+π⊥⋅˙d3 =0. (51b)

The pertaining initial conditions are:

 Λ(0)=¯Λ,w⊥(0)=¯w⊥,w∥(0)=¯w∥. (52)
###### Proof.

The theorem follows from the systematic calculation of , 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., . This is easily proven as follows

 0=˙π=(I−d3⊗d3)˙π+(d3⋅˙π)d3=(I−d3⊗d3)(˙π⊥+˙π∥d3+π∥˙d3)+d3⋅(˙π⊥+˙π∥d3+π∥˙d3)d3=∇˙d3π⊥+π∥˙d3∈T∗d3S2+(˙π∥+˙π⊥⋅d3)d3∈span(d3). (53)

###### Theorem 3.3.

The conservation laws of the rotating rigid body are:

 K =12W⋅JW=12w⋅jw=const., (54a) π =jw=ΛJW=const.. (54b)
###### 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

 δW=δθ⊥⋅mext⊥+δθ∥mext∥, (55)

and add this contribution to the variation of the action.

### 3.5 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

 Q:={q=(d1,d2,d3)∈S2×S2×S2∣di⋅dj=0,i≠j}≅SO(3), (56)

whose tangent space at the point is given by

 TqQ:={˙q=(˙d1,˙d2,˙d3)∣˙di=ω×di,ω∈R3}. (57)

Now, we start by defining the rigid body as the bounded set of points

 X=θ1D1+θ2D2+θ3D3 (58)

where are the material coordinates of the point and are three orthogonal directors, with the third one oriented in the direction of the principal axis of inertia and such that . The position of the point at time is denoted as and given by

 x(t)=φ(θ1,θ2,θ3;t)=θ1d1(t)+θ2d2(t)+θ3d3(t) (59)

with , for all . On this basis, there must be a rotation tensor , where we have employed the sum convention for repeated indices, such that . The material velocity of the particle is the vector that can be written as

 ˙x(t)=˙φ(θ1,θ2,θ3;t)=θ1˙d1(t)+θ2˙d2(t)+θ3˙d3(t) (60)

with representing three director velocity vectors.

To construct the dynamic equations of the model, assume the body has a density per unit of reference volume and hence its total kinetic energy, or Lagrangian, can be formulated as

 K=∫B0ρ02|˙x|2dB0 . (61)

To employ Hamilton’s principle of stationary action, but restricting the body directors to remain orthonormal at all time, we define the constrained action

 S=∫T0(K−h(d1,d2,d3)⋅λ)dt , (62)

where is given by Eq. (61), is a vector of Lagrange multipliers, and is of the form

 h(d1,d2,d3)=⎛⎜⎝d2⋅d3d1⋅d3d1⋅d2⎞⎟⎠, (63)

such that 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:

 ˙π1(d1,d2,d3)+H1(d1,d2,d3)Tλ =0, (64a) ˙π2(d1,d2,d3)+H2(d1,d2,d3)Tλ =0, (64b) ˙π3(d1,d2,d3)+H3(d1,d2,d3)Tλ =0, (64c) h(d1,d2,d3) =0. (64d)

The generalized momenta are defined as

 πi=Ji1˙d1+Ji2˙d2+Ji3˙d3, (65)

where Euler’s inertia coefficients are

 Jij=Jij=∫B0ϱ0θiθjdB0, (66)

for and from to . In addition, the splitting of the inertia tensor implies and stands for .

The pertaining initial conditions are:

 d1(0)=¯d1,d2(0)=¯d2,d3(0)=¯d3,w⊥(0)=¯d3×˙¯d3,w∥(0)=¯w∥ (67)
###### Proof.

The theorem follows from the systematic calculation of . ∎

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 .

###### 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

 δW=δW1+δW2+δW3withδWi=12(di×δdi)⋅(mext⊥+mext∥d3). (68)

## 4 Rotating rigid body based on the non-twisting frame

In this section, we introduce the nonholonomic rotating rigid body, which incorporates the non-integrable 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 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.

### 4.1 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:

 ∇˙d3π⊥+π∥˙d3 =0, (69a) ˙π∥+π⊥⋅˙d3+μ =0, (69b) w⋅d3 =0. (69c)

The pertaining initial conditions are:

 Λ(0)=¯Λ,w⊥(0)=¯w⊥,w∥(0)=0. (70)

Moreover, Eq. (69) can be rewritten as

 ∇˙d3π⊥ =0, (71a) π⊥⋅˙d3+μ =0, (71b) w∥ =0, (71c)

where must satisfy the parallel transport along the curve .

###### Proof.

The first part follows from the inclusion of the force associated to the presence of the nonholonomic restriction given by

 g=w⋅d3=0 (72)

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

 δWnh=μ∂(w⋅d3)∂w⋅δθ=δθ⋅(μd3)=δθ∥⋅(μd3) (73)

where denotes the corresponding Lagrange multiplier. The second part follows from noticing that implies . ∎

###### Theorem 4.2.

The conservation laws of the rotating rigid body based on the non-twisting frame are:

 K =12W⋅JW=12w⋅jw=const., (74a) Π1 =D1⋅Π=d1⋅π=% const., (74b) Π2 =D2⋅Π=d2⋅π=% const.. (74c)
###### Proof.

To prove the conservation of kinetic energy, let us consider the following equilibrium statement

 δθ⊥⋅(∇˙d3π⊥+π∥˙d3)+δθ∥⋅((π⊥⋅˙d3+μ)d3)+δμ(w⋅d3)=0, (75)

where , and are admissible variations. Now by choosing , and , we have that

 0=w⊥⋅(∇˙d3π⊥+π∥˙d3)+w∥⋅((π⊥⋅˙d3+μ)d3)=w⋅˙π=W⋅˙Π=∂∂t(12W⋅JW)=˙K, (76)

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

 π⊥=(d1⋅π)d1+(d2⋅π)d2=Π1d1+Π2d2. (77)

Now by introducing the former expression into the first statement of Eq. (71), we have that

 0=∇˙d3π⊥=∇˙d3Π1d1+Π2d2