Abstract
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.
Similar content being viewed by others
Avoid common mistakes on your manuscript.
1 Introduction
The rigid body is a problem of classical mechanics that has been exhaustively studied (see, e.g., Goldstein 1980; Arnold 1989). 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. 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 (for example in Simo et al. 1988; Antman 1995; Mielke and Holmes 1988). As a result, geometric integrators for the rigid body (Simo and Wong 1991; Romero 2008) are at the root of more complex numerical integration schemes for models that involve, in one way or another, rotations (Simo and Vu-Quoc 1986; Lewis and Simo 1994; Sansour and Bednarczyk 1995; Jelenić and Crisfield 1998; Romero and Armero 2002b; Betsch and Steinmann 2003; Romero and Arnold 2017).
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 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 (Romero and Armero 2002a; Betsch and Steinmann 2002; Betsch and Leyendecker 2006).
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 (Bishop 1975). 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 (Shi and Hearst 1994; McMillen and Goriely 2002; Audoly et al. 2007; Bergou et al. 2008; Giusteri and Fried 2018; Romero and Gebhardt 2020), 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. General geometric formulations for nonholonomic systems have been, in the last 2 or 3 decades, subject of a very active research. Its modern treatment can be found in seminal works such as Koiller (1992), Bates and Śniatycki (1993), de León and de Diego (1996), Bloch et al. (1996) and de León and de Diego (1997). Moreover, its extension toward a general nonholonomic field theory is investigated in Vankerschaver et al. (2005). In the nonholonomic setting, conservation laws take a different form when compared with the usual ones and their identification is not straightforward. For the rigid body model based on the non-twisting frame, 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-twisting frame may substantially differ from other nonholonomic cases that were investigated in the literature, e.g., (Betsch 2006; Hedrih 2019). 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 Lee et al. (2009, (2018). However, the anisotropy of the inertial properties has been completely disregarded.
There already exists a plethora of references dealing with the geometric integration of nonholonomic systems based on discrete Lagrangian mechanics (e.g., Cortés and Martínez 2001; de León et al. 2004). These works are based directly on a discrete version of the Lagrange–D’Alembert principle and result in integrators that inherit the geometric structure of the problem, for example, preserving the structure of the evolution of the symplectic form and the nonholonomic momentum along trajectories.
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 (McLachlan et al. 1999; Celledoni et al. 2012) preserves remarkably these invariants of the motion, exactly, resulting in accurate pictures of the rigid body dynamics. A comprehensive description of such ideas in the context of general nonholonomic dynamics is to be found in Celledoni et al. (2019). The specialization of approaches based on more elaborated conservative/dissipative integration schemes like Gebhardt et al. (2020) is possible as well in this context, but falls outside the scope of the current work and therefore, not addressed here. Lastly, be aware that the proposed model should not be understood as an alternative formulation to the well-known standard rigid body model but as a different mechanical problem that possesses interesting properties, which can find applications in fields like multibody systems, n-body problems on manifolds, computer graphics and ballistics among others.
The rest of the article has the following structure. In Sect. 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 geometric 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 favorable mathematical setting that will be exploited later on to derive an structure-preserving integration algorithm. In Sect. 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 Sect. 6, numerical results that show the excellent performance of the new energy–momentum algorithm are presented and discussed. Finally, Sect. 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., Marsden and Ratiu 1994; Lee et al. 2018 for more comprehensive expositions).
2.1 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 (Eisenberg and Guy 1979; Simo et al. 1989; Romero 2004; Romero et al. 2014; Romero and Arnold 2017). As an embedded set on Euclidean space, it is defined as
where the dot product is the standard one in \(\mathbb {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 \(\varvec{d}\) are those defined by the expressions:
where the product “\(\times \)” is the standard cross-product in \(\mathbb {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 \(\mathbb {R}^3\). As such, the covariant derivative of a smooth vector field \(\varvec{v}{:}\,S^2\rightarrow TS^2\) along a vector field \(\varvec{w}{:}\,S^2\rightarrow TS^2\) is the vector field \(\nabla _{\varvec{w}}\varvec{v}\in TS^2\) that evaluated at \(\varvec{d}\) is precisely the projection of the derivative \(D\varvec{v}\) in the direction of \(\varvec{w}\) onto the tangent plane to \(\varvec{d}\). Hence, denoting as \(\varvec{I}\) the unit second-order tensor and \(\otimes \) the outer product between vectors, both on \(\mathbb {R}^3\), this projection can be simply written as
When \(\varvec{d}{:}\,(a,b)\rightarrow S^2\) is a smooth one-parameter curve on the unit sphere and \(\varvec{d}'\) its derivative, the covariant derivative of a smooth vector field \(\varvec{v}{:}\,S^2\rightarrow TS^2\) in the direction of \(\varvec{d}'\) has an expression that follows from Eq. (4), that is,
which, as before, is nothing but the projection of \((\varvec{v}\circ \varvec{d})'\) onto the tangent space \(T_{\varvec{d}}S^2\), and the symbol “\(\circ \)” 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_{\varvec{d}_0}S^2 \rightarrow S^2\) is a surjective application with a closed-form expression given by the formula
where \(\varvec{v}_0\) must be a tangent vector on \(T_{\varvec{d}_0}S^2\) and \(|\cdot |\) denotes the Euclidean vector norm. The inverse of the exponential function is the logarithmic map \(\log {:}\,S^2 \rightarrow T_{\varvec{d}_0}S^2\), for which also there is a closed-form expression that reads
with \(\varvec{d}_0\ne \varvec{d}\). The geodesic \(\varvec{d}_\mathrm {G}(s)\) for \(s\in [0,1]\) with \(\varvec{d}_\mathrm {G}(0)=\varvec{d}_0\) and \(\varvec{d}'_\mathrm {G}(0)=\varvec{v}_0\) is a great arch on the sphere obtained as the solution of the equation
with the explicit form
The parallel transport of \(\varvec{w}_0\in T_{\varvec{d}_0}S^2 \mapsto \varvec{w} \in T_{\varvec{d}}S^2\) along the geodesic \(\varvec{d}_\mathrm {G}\) is then given by
and verifies
More details about expressions (6)–(10) can be found in Hosseini and Uschmajew (2017) and Bergmann et al. (2018).
Given definitions (1) and (2) of the unit sphere and its tangent bundle, we recognize that there exists an isomorphism
for any \(\varvec{d}\in S^2\). Given now two directors \(\varvec{d},\tilde{\varvec{d}}\), we say that a second-order tensor \(\varvec{T}{:}\,\mathbb {R}^3\rightarrow \mathbb {R}^3\) splits from \(\varvec{d}\) to \(\tilde{\varvec{d}}\) if it can be written in the form
where \({\varvec{T}}_\perp \) is a bijection from \(T_{\varvec{d}} S^2\) to \(T_{\tilde{\varvec{d}}}S^2\) with \(\ker ({\varvec{T}}_\perp ) = \hbox {span}(\varvec{d})\), and \({\varvec{T}}_\parallel \) is a bijection from \(\hbox {span}(\varvec{d})\) to \(\hbox {span}(\tilde{\varvec{d}})\) with \(\ker ({\varvec{T}}_\parallel ) \equiv T_{\varvec{d}}S^2\). Split (13) depends on the pair \(\varvec{d},\tilde{\varvec{d}}\) 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
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 \(\hat{(\cdot )}{:}\,\mathbb {R}^3\rightarrow so(3)\) such that for all \(\varvec{w},\varvec{a}\in \mathbb {R}^3\), the tensor \(\hat{\varvec{w}}\in so(3)\) satisfies \(\hat{\varvec{w}}\varvec{a} = \varvec{w}\times \varvec{a}\). Here, \(\varvec{w}\) is referred to as the axial vector of \(\hat{\varvec{w}}\), which is a skew-symmetric tensor, and we also write . Being a Lie group, the space of rotations has an exponential map \(\exp {:}\,so(3)\rightarrow SO(3)\) whose closed-form expression is Rodrigues’ formula
with \(\varvec{\theta }\in \mathbb {R}^3, \theta =|\varvec{\theta }|\). The linearization of the exponential map is simplified by introducing the map \(\mathrm {dexp}{:}\,so(3)\rightarrow \mathbb {R}^{3\times 3}\) that satisfies
for every \(\varvec{\theta }{:}\,\mathbb {R}\rightarrow \mathbb {R}^3\). A closed-form expression for this map, as well as more details regarding the numerical treatment of rotations, can be found elsewhere (Hairer et al. 2006; Romero 2008; Romero and Arnold 2017).
2.3 The Notion of Twist and the Non-twisting Frame
Let \(\varvec{d}_3{:}\,\left[ 0,T\right] \rightarrow S^2\) indicate a curve of directors parameterized by time \(t\in [0,T]\). Now, let us consider two other curves \(\varvec{d}_1,\varvec{d}_2{:}\,[0,T]\rightarrow S^2\) such that \(\{\varvec{d}_1(t), \varvec{d}_2(t), \varvec{d}_3(t)\}\) are mutually orthogonal for all \(t\in [0,T]\). We say that the triad \(\{\varvec{d}_1,\varvec{d}_2,\varvec{d}_3\}\) moves without twist if
where the overdot refers to the derivative with respect to time. Given the initial value of the triad \(\{\varvec{d}_1(0), \varvec{d}_2(0), \varvec{d}_3(0) \} = \{ \varvec{D}_1,\varvec{D}_2, \varvec{D}_3\}\), there is a map \(\varvec{\chi }{:}\,[0,T]\rightarrow SO(3)\) transforming it to the frame \(\{ \varvec{d}_1, \varvec{d}_2, \varvec{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 \(\varvec{v} \in T_{\varvec{d}_3}S^2\) is said to be parallel-transported along \(\varvec{d}_3\) if and only if \(\nabla _{\dot{\varvec{d}}_3}{\varvec{v}}=\varvec{0}\). An consequence of this is that
and similarly
In addition, we have that
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 [\psi \,\,\hat{\varvec{d}}_3]\), with \(\psi {:}\,[0,T]\rightarrow S^1\) (the unit 1-sphere), and the rotated triad \(\{\varvec{d}_1^\psi ,\varvec{d}_2^\psi ,\varvec{d}_3\} = \exp [\psi \,\,\hat{\varvec{d}}_3] \{\varvec{d}_1,\varvec{d}_2,\varvec{d}_3\}\). Then,
and the Darboux vector of the rotated triad is
or, equivalently,
In this last expression, we identify the twist rate \(\dot{\psi }\) and the twist angle \(\psi \), respectively, as the rotation velocity component on the \(\varvec{d}_3\) direction and the rotation angle about this same vector (for further details, see Bishop 1975; Langer and Singer 1996). The previous calculations show that the frame \(\{\varvec{d}_1,\varvec{d}_2,\varvec{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 \(\{\varvec{D}_1, \varvec{D}_2,\varvec{D}_3\}\) without twist. In this context, a summary of recent advances and open problems is presented for instance in Farouki (2016).
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 \(\varvec{d}_3\), the dynamics of the body that takes place in the cotangent space \(T^*_{\varvec{d}_3}S^2\) is separated from that one corresponding to the reciprocal normal space \(N^*_{\varvec{d}_3}S^2 \equiv \text {span}(\varvec{d}_3)\). In order to do this, we have employed the identifications
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 \(Q\equiv SO(3)\).
Let us now study the motion of a rotating rigid body, that is, a time-parameterized curve in configuration space \(\varvec{\varLambda }{:}\, [0,T] \rightarrow Q\). The generalized velocity of the rotating rigid body belongs, for every \(t\in [0,T]\), to the tangent bundle
The time derivative of the rotation tensor can be written as
where \(\varvec{w}\) and \(\varvec{W}\) are the spatial and convected angular velocities, respectively.
Let \(\{\varvec{E}_1,\varvec{E}_2, \varvec{E}_3\}\) be a fixed basis of the ambient space. Then, if \(\varvec{d}_i=\varvec{\varLambda } \varvec{E}_i\), with \(i=1,2,3\) we can use Eq. (27) to split the rotation vectors as in
Then, using relations (29), we identify
3.2 Kinetic Energy and Angular Momentum
Let us now select a fixed Cartesian basis of \(\mathbb {R}^3\) denoted as \(\{ \varvec{D}_i \}_{i=1}^3\), where the third vector coincides with one of the principal directions of the convected inertia tensor \(\varvec{J}{:}\,\mathbb {R}^3\rightarrow \mathbb {R}^3\) of the body, a symmetric, second-order, positive definite tensor. Thus, this tensor splits from \(\varvec{D}_3\) to \(\varvec{D}_3\) and we write
where \(\varvec{J}_\perp \) maps bijectively \(\mathrm {span}(\varvec{D}_1,\varvec{D}_2)\) onto itself and satisfies \(\varvec{J}_\perp \varvec{D}_3 = \varvec{0}\).
The kinetic energy of a rigid body with a fixed point is defined as the quadratic form
where \(\varvec{j}\) is the spatial inertia tensor, the push-forward of the convected inertia, and defined as
Let \(\varvec{d}_3=\varvec{\varLambda }\varvec{D}_3\). Given the relationship between the convected and spatial inertia, it follows that the latter also splits, this time from \(\varvec{d}_3\) to \(\varvec{d}_3\), and thus
where now \(\varvec{j}_\perp \) maps bijectively \(\mathrm {span}(\varvec{d}_1,\varvec{d}_2)\) onto itself and satisfies \(\varvec{j}_\perp \varvec{d}_3 = \varvec{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 \(\varvec{\pi }\) 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
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 \(\varvec{\varLambda }_\iota (t)\) parameterized by the scalar \(\iota \) and given by
where \(\delta \varvec{\theta }{:}\,[0,T]\rightarrow \mathbb {R}^3\) represents any arbitrary variation that satisfies
The curve \(\varvec{\varLambda }_\iota \) passes through the configuration \(\varvec{\varLambda }\) when \(\iota =0\) and has tangent at this point
For future reference, let us calculate the variation of the derivative \(\dot{\varvec{\varLambda }}\). To do so, let us first define the temporal derivative of the perturbed rotation, that is,
Then, the variation of \(\dot{\varvec{\varLambda }}\) 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.
Theorem 1
The variations of the convected angular velocities \((\varvec{W}_\perp , W_\parallel )\) are
Proof
The convected angular velocities of the one-parameter curve of configurations \(\varvec{\varLambda }_\iota \) are
where \(\varvec{d}_{i,\iota } = \varvec{\varLambda }_\iota \varvec{D}_i\). The variation of the angular velocity perpendicular to \(\varvec{D}_3\) is obtained from its definition employing some algebraic manipulations and expression (45) as follows:
The variation of the angular velocity parallel to \(\varvec{D}_3\) follows similar steps:
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
with unknown fields \((\varvec{\varLambda }, \dot{\varvec{\varLambda }})\in TQ\).
Theorem 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 \(\delta S\), the variation of the action, based on the variation of the convected angular velocities of Eq. (1); thus, we omit a detailed derivation. Equation (2), which is presented in its split form, is equivalent to Euler’s equations which state that the spatial angular momentum is preserved, i.e., \(\dot{\varvec{\pi }}=\varvec{0}\). This is easily proven as follows
Theorem 3
The conservation laws of the rotating rigid body are:
Proof
This is an standard result, and thus, we omit further details.
Remark 1
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.
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
whose tangent space at the point \(\varvec{q}\) is given by
Now, we start by defining the rigid body as the bounded set \(\mathcal{B}_0\subset \mathbb {R}^3\) of points
where \((\theta ^1,\theta ^2,\theta ^3)\) are the material coordinates of the point and \(\{ \varvec{D}_i\}_{i=1}^3\) are three orthogonal directors, with the third one oriented in the direction of the principal axis of inertia and such that \(\varvec{D}_3=\varvec{D}_1\times \varvec{D}_2\). The position of the point \(\varvec{X}\) at time \(t\in [0,T]\) is denoted as \(\varvec{x}(t)\in \mathbb {R}^3\) and given by
with \((\varvec{d}_1,\varvec{d}_2,\varvec{d}_3) = \varvec{q} \in Q\), for all t. On this basis, there must be a rotation tensor \(\varvec{\varLambda }(t)= \varvec{d}_i(t)\otimes \varvec{D}_i\), where we have employed the sum convention for repeated indices, such that \(\varvec{d}_i(t) = \varvec{\varLambda }(t) \varvec{D}_i\). The material velocity of the particle \(\varvec{X}\) is the vector \(\dot{\varvec{x}}(t)\in \mathbb {R}^3\) that can be written as
with \((\dot{\varvec{d}}_1,\dot{\varvec{d}}_2,\dot{\varvec{d}}_3)=\dot{\varvec{q}}\in T_{\varvec{q}}Q\) representing three director velocity vectors.
To construct the dynamic equations of the model, assume the body \(\mathcal{B}_0\) has a density \(\rho _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), \(\varvec{\lambda }\in \mathbb {R}^3\) is a vector of Lagrange multipliers, and \(\varvec{h}\) is of the form
such that \(\varvec{h}(\varvec{d}_1,\varvec{d}_2,\varvec{d}_3)=\varvec{0}\) expresses the directors’ orthonormality.
Theorem 4
The alternative equations of motion, i.e., the Euler–Lagrange equations, for the standard rotating rigid body model are:
The generalized momenta \((\varvec{\pi }^1, \varvec{\pi }^2, \varvec{\pi }^3) = \varvec{p}\in T_{\varvec{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 \(\mathscr {J}^{13}=\mathscr {J}^{23} = 0\) and \(\varvec{H}_i \in L(T_{\varvec{d}_i}S^2, \mathbb {R}^n)\) stands for \(\frac{\partial \varvec{h}}{\partial \varvec{d}_i}\).
The pertaining initial conditions are:
Proof
The theorem follows from the systematic calculation of \(\delta S\).
The reparameterized 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 (Romero and Armero 2002a).
Remark 2
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 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. (2) to account the non-integrable condition \(W_\parallel =w_\parallel =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.
4.1 Governing Equations and Invariants
Theorem 5
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. (5) can be rewritten as
where \(\varvec{\pi }_\perp \in T^{*}_{\varvec{d}_3}S^2\) must satisfy the parallel transport along the curve \(\varvec{d}_3\in S^2\).
Proof
The first part follows from the inclusion of the force associated with 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 with the presence of this nonholonomic restriction can be computed as
where \(\mu \in \mathbb {R}\) denotes the corresponding Lagrange multiplier. The second part follows from noticing that \(w_\parallel =\varvec{w}\cdot \varvec{d}_3=0\) implies \(\pi _\parallel =0\).
Theorem 6
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 \(\delta \varvec{\theta }_\perp \in T_{\varvec{d}_3}S^2, \delta \varvec{\theta }_\parallel \in N_{\varvec{d}_3}S^2\) and \(\delta \mu \) are admissible variations. Now by choosing \(\delta \varvec{\theta }_\perp = \varvec{w}_\perp , \delta \varvec{\theta }_\parallel = \varvec{w}_\parallel \) and \(\delta \mu = 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. (5), we have that
in which the parallel transport of \(\varvec{d}_1\) and \(\varvec{d}_2\), both in \(T_{\varvec{d}_3}S^2\), has been accounted for. This shows that the first and second components of the material angular momentum are preserved by the motion.
Remark 3
To include external moments acting on the rotating rigid body based on the non-twisting frame, it is necessary to compute the associated virtual work as follows
4.2 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\in [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.
Theorem 7
The alternative nonholonomic equations of motion, i.e., the Lagrange–D’Alembert equations, for the rotating rigid body model based on the non-twisting frame are:
The pertaining initial conditions are:
Proof
The theorem follows from the computation of the virtual work associated with the presence of the nonholonomic restriction, namely
where \(\mu \in \mathbb {R}\) denotes the corresponding Lagrange multiplier.
Remark 4
To include external moments acting on the rotating rigid body based on the non-twisting frame, it is necessary to compute the associated virtual work as follows
5 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.
5.1 Basic Energy–Momentum Algorithm
Next, we describe the application of the energy–momentum integration algorithm (Simo and Wong 1991; Simo and Tarnow 1992) to the “standard rotating rigid body” case. For this purpose, the following nomenclature is necessary:
While \(\varvec{q}\) is the vector of generalized coordinates, \(\varvec{p}\) collects the generalized momenta and \(\varvec{Q}^\text {ext}\) contains the generalized external loads, if present. The discrete version of Eq. (4) can be expressed at time \(n+\frac{1}{2}\) as
where \(\left\langle \cdot ,\cdot \right\rangle \) 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 \(\varvec{q}_n\) and \(\dot{\varvec{q}}_n\) are known from the previous step, \(\varvec{q}_{n+1}\) are \(\dot{\varvec{q}}_{n+1}\) are unknown, and \(\dot{\varvec{q}}_{n+1}\) is computed as \(\frac{2}{h}(\varvec{q}_{n+1}-\varvec{q}_n)-\dot{\varvec{q}}_n\) once \(\varvec{q}_{n+1}\) has been determined by means of an iterative procedure, typically the Newtown–Raphson method.
The mass matrix takes the form
and \(\mathscr {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 (Gebhardt et al. 2019a, b) as
where \(\varvec{q}(\xi )\) is defined as \(\frac{1}{2}(1-\xi )\varvec{q}_{n}+\frac{1}{2}(1+\xi )\varvec{q}_{n+1}\) for \(\xi \in [-1,+1]\). The algorithmic Jacobian matrix defined in this way satisfies for any admissible solution the discrete version of the hidden constraints, i.e.,
Theorem 8
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.
5.2 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 \(\varvec{q}(\xi )\) is defined as before. In this way, the nonholonomic constraint is identically satisfied at the midpoint, i.e.,
Theorem 9
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., \(\varPi ^1\), we need to consider the following discrete variation
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
which can be easily shown by considering that
since the angular velocity \(\varvec{w}_{n+\frac{1}{2}}\) has the form \(\varvec{d}_{3,n+\frac{1}{2}}\times \dot{\varvec{d}}_{3,n+\frac{1}{2}}\) due to the satisfaction of the non-twisting condition, see Sect. 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., \(\varPi ^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.
6 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]\,\,\text {s}\) with a time step size of \(h = 0.005\,\,\text {s}\) and relative tolerance \(10^{-10}\). Additionally, we present a brief comparison between the results obtained with the proposed algorithm and those results obtained with a well-established nonholonomic integrator.
6.1 Case 1: Response to Nonzero Initial Conditions
For this first case we consider
and
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 \(\omega _\parallel = 0\,\,\mathrm {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 \(\varvec{d}_1\) and \(\varvec{d}_2\) in \(T_{\varvec{d}_3}S^2\) are being parallel transported along the time-parameterized solution curve \(\varvec{d}_3\).
Figure 2 presents the time history for the kinetic energy and the second quotient of precision as defined in “Appendix A.” 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 \(\varvec{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.
6.2 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 \(\varvec{m}^{\mathrm {ext}}_\perp \) 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 varies 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 \(\varvec{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.
6.3 Case 3: Response to Nonzero Initial Conditions and a Vanishing Load
For this last case we consider
and
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 varies 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 \(\varvec{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.
6.4 Comparison with the Discrete Lagrange–D’Alembert Algorithm
Here we present, for the rigid body model based on the non-twisting frame, a brief comparison between the results obtained with the specialized energy–momentum method proposed and those results obtained with a well-established nonholonomic integrator (Cortés and Martínez 2001). We note that the aim of this section is not to provide a detailed comparison between numerical methods, but rather to put the proposed ideas into a suitable context. Such a detailed comparative analysis falls outside the present scope.
For the particular case of a holonomically (at most quadratically) constrained mechanical system with a constant mass matrix and without a potential function, the variational integrator of second order based on the discrete Euler–Lagrange algorithm is able to exactly preserve the Hamiltonian and energy functions at time \(t_i\) and \(t_{i+\frac{1}{2}}\), respectively. Moreover, for the unforced case, the Hamiltonian function practically matches the exact energy of the system. The Hamiltonian and the energy have similar, but different, values. This behavior is comprehensively reported in Leyendecker et al. (2008). In presence of a potential function, even the simplest linear one, this property is lost. For those cases, the Hamiltonian and energy functions oscillate about an average value that is close to the exact energy. Typically, the amplitude of such oscillations tends to zero as the time step tends to zero as well. Notwithstanding, the linear and angular momenta are always exactly preserved. The non-variational discrete Lagrange–D’Alembert algorithm can be considered as a natural extension of the discrete Euler–Lagrange algorithm to the realm of nonholonomic systems. This well-established nonholonomic integrator is able to preserve the nonholonomic momenta associated with horizontal symmetries and under the same conditions, i.e., a mechanical system with a constant mass matrix and without a potential function, it seems to preserve the Hamiltonian and the energy functions as well. Nonetheless, a formal proof of this particular behavior also falls outside the current scope. “Appendix B” summarizes the discrete Lagrange–D’Alembert algorithm for nonholonomic systems and provides the formulas used in this subsection to compute the responses for the three cases investigated.
In the first case, which studies the response to nonzero initial conditions, Table 7 shows the stationary values of the motion invariants. For sake of brevity, we introduce the following nomenclature: EM stands for energy–momentum and DLA stands for discrete Lagrange–D’Alembert. We can observe that the numerical values of the first and second components of the material angular momentum are indistinguishable. In addition, the value of the Hamiltonian is identical for both methods. The energy computed with the EM algorithm matches the value of the Hamiltonian, but the energy computed with the DLA algorithm differs a little as expected. Figure 10 shows the first component of the third director. The responses are similar when plotted for the time interval [0, 5] s. Nevertheless, if we take a closer look, for instance in the time interval [4, 5] s, the responses can be distinguished from each other, with the EM showing slower oscillations. A similar behavior is observed for the second and third components of the third director, see Figs. 11 and 12, respectively. The two methods converge quadratically and take in average 3–4 Newton–Raphson iterations to satisfy the relative tolerance criterion chosen. Moreover, we computed with both methods their second quotients of precision which evaluate to 4, approximately, certifying the implementation’s correctness.
For the second case, the one analyzing the response to a vanishing load, Table 8 shows the stationary values of the motion invariants. This time, all the invariants computed with the DLA algorithm, even if the agreement is excellent, differ from the invariants computed with the EM algorithm. The comparison of the components of the third director shows a very similar behavior to the one observed in the first case, and therefore, these additional results are omitted.
Finally, for the third case that corresponds to the solution for nonzero initial conditions and a vanishing load, Table 9 shows the stationary values of the motion invariants. Once again, the invariants computed with the DLA algorithm, even if the agreement with those calculated with the EM is very good, differ from the latter. As before, the comparison of the components of the third director shows a very similar behavior as observed in the first and second cases. Thus, these results are omitted as well.
In summary, we can claim that, for the investigated cases, the results obtained with both methods are in excellent agreement with each other. Moreover, the computational costs of carrying out the numerical simulations are practically the same for both approaches. Therefore, we can state that our approach is very competitive and produces correct pictures of the underlying nonholonomic system. These observations are very encouraging and motivate further investigations along this direction.
7 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 analysis done for the proposed model, an implicit, second-order accurate, energy and momentum conserving algorithm is presented, which discretizes in time the rigid body, nonholonomic 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.
References
Antman, S.S.: Nonlinear Problems of Elasticity. Springer, New York (1995)
Arnold, V.I.: Mathematical Methods of Classical Mechanics. Springer, Berlin (1989)
Audoly, B., Clauvelin, N., Neukirch, S.: Elastic knots. Phys. Rev. Lett. 99, 137–4 (2007)
Bates, L., Śniatycki, J.: Nonholonomic reduction. Rep. Math. Phys. 32, 99–115 (1993)
Bergmann, R., Fitschen, J.H., Persch, J., Steidl, G.: Priors with coupled first and second order differences for manifold-valued image processing. J. Math. Imaging Vis. 60, 1459–1481 (2018)
Bergou, M., Wardetzky, M., Robinson, S., Audoly, B., Grinspun, E.: Discrete elastic rods. In: ACM Transactions on Graphics, vol. 27 (2008)
Betsch, P.: Energy-consistent numerical integration of mechanical systems with mixed holonomic and nonholonomic constraints. Comput. Methods Appl. Mech. Eng. 195, 7020–7035 (2006)
Betsch, P., Leyendecker, S.: The discrete null space method for the energy consistent integration of constrained mechanical systems: Part II. Multibody dynamics. Int. J. Numer. Methods Eng. 67, 499–552 (2006)
Betsch, P., Steinmann, P.: Frame-indifferent beam finite elements based upon the geometrically exact beam theory. Int. J. Numer. Methods Eng. 54, 1775–1788 (2002)
Betsch, P., Steinmann, P.: Constrained dynamics of geometrically exact beams. Comput. Mech. 31, 49–59 (2003)
Bishop, R.L.: There is more than one way to frame a curve. Am. Math. Mon. 82, 246–251 (1975)
Bloch, A.M., Krishnaprasad, P.S., Marsden, J.E., Murray, R.M.: Nonholonomic mechanical systems with symmetry. Arch. Ration. Mech. Anal. 136, 21–99 (1996)
Celledoni, E., Grimm, V., McLachlan, R.I., McLaren, D.I., O’Neale, D., Owren, B., Quispel, G.R.W.: Preserving energy resp. dissipation in numerical PDEs using the ‘Average Vector Field’ method. J. Comput. Phys. 231, 6770–6789 (2012)
Celledoni, E., Farré Puiggalí, M., Høiseth, E., de Diego, D.M.: Energy-preserving integrators applied to nonholonomic systems. J. Nonlinear Sci. 29, 1523–1562 (2019)
Cortés, J., Martínez, S.: Non-holonomic integrators. Nonlinearity 14, 1365–1392 (2001)
de León, M., de Diego, D.M.: On the geometry of non-holonomic Lagrangian systems. J. Math. Phys. 37, 3389–3414 (1996)
de León, M., de Diego, D.M.: A constraint algorithm for singular Lagrangians subjected to nonholonomic constraints. J. Math. Phys. 38, 3055–3062 (1997)
de León, de Diego, D.M., Santamaría-Merino, A.: Geometric numerical integration of nonholonomic systems and optimal control problems. Eur. J. Control 10(5), 515–521 (2004)
Eisenberg, M., Guy, R.: A proof of the hairy ball theorem. Am. Math. Mon. 86, 571–574 (1979)
Farouki, R.T.: Rational rotation-minimizing frames-recent advances and open problems. Appl. Math. Comput. 272, 80–91 (2016)
Gebhardt, C.G., Hofmeister, B., Hente, C., Rolfes, R.: Nonlinear dynamics of slender structures: a new object-oriented framework. Comput. Mech. 63, 219–252 (2019a)
Gebhardt, C.G., Steinbach, M.C., Rolfes, R.: Understanding the nonlinear dynamics of beam structures: a principal geodesic analysis approach. Thin-Walled Struct. 140, 357–372 (2019b)
Gebhardt, C.G., Romero, I., Rolfes, R.: A new conservative/dissipative time integration scheme for nonlinear mechanical systems. Comput. Mech. 65, 405–427 (2020)
Giusteri, G., Fried, E.: Importance and effectiveness of representing the shapes of cosserat rods and framed curves as paths in the special euclidean algebra. J. Elast. 132, 43–65 (2018)
Goldstein, H.: Classical Mechanics, 2nd edn. Addison-Wesley, Boston (1980)
Hairer, E., Lubich, C., Wanner, G.: Geometric Numerical Integration, 2nd edn. Springer, Berlin (2006)
Hedrih, K.R.: Rolling heavy ball over the sphere in real Rn3 space. Nonlinear Dyn. 97, 63–82 (2019)
Hosseini, S., Uschmajew, A.: A Riemannian gradient sampling algorithm for nonsmooth optimization on manifolds. SIAM J. Optim. 27, 173–189 (2017)
Jelenić, G., Crisfield, M.A.: Interpolation of rotational variables in nonlinear dynamics of 3d beams. Int. J. Numer. Methods Eng. 43, 1193–1222 (1998)
Koiller, J.: Reduction of some classical non-holonomic systems with symmetry. Arch. Ration. Mech. Anal. 118(2), 113–148 (1992)
Kreiss, H.O., Ortiz, O.E.: Introduction to Numerical Methods for Time Dependent Differential Equations. Wiley, London (2014)
Langer, J., Singer, D.: Lagrangian aspects of the Kirchhoff elastic rod. SIAM Rev. 38, 605–618 (1996)
Lee, T., Leok, M., McClamroch, N.H.: Lagrangian mechanics and variational integrators on two-spheres. Int. J. Numer. Methods Eng. 79, 1147–1174 (2009)
Lee, T., Leok, M., McClamroch, N.H.: Global Formulations of Lagrangian and Hamiltonian Dynamics on Manifolds. Springer, Berlin (2018)
Lewis, D., Simo, J.C.: Conserving algorithms for the dynamics of hamiltonian systems on lie groups. J. Nonlinear Sci. 4, 253–299 (1994)
Leyendecker, S., Marsden, J.E., Ortiz, M.: Variational integrators for constrained dynamical systems. Z. Angew. Math. Mech. 88, 677–708 (2008)
Marsden, J.E., Ratiu, T.S.: Introduction to Mechanics and Symmetry, 1st edn. Springer, New York (1994)
McLachlan, R.I., Quispel, G.R.W., Robidoux, N.: Geometric integration using discrete gradients. Philos. Trans. R. Soc. Lond. Ser. A Math. Phys. Sci. Eng. 357, 1021–1045 (1999)
McMillen, T., Goriely, A.: Tendril perversion in intrinsically curved rods. J. Nonlinear Sci. 12, 241–281 (2002)
Mielke, A., Holmes, P.: Spatially complex equilibria of buckled rods. Arch. Ration. Mech. Anal. 101, 319–348 (1988)
Romero, I.: The interpolation of rotations and its application to finite element models of geometrically exact rods. Comput. Mech. 34, 121–133 (2004)
Romero, I.: Formulation and performance of variational integrators for rotating bodies. Comput. Mech. 42, 825–836 (2008)
Romero, I., Armero, F.: An objective finite element approximation of the kinematics of geometrically exact rods and its use in the formulation of an energy-momentum conserving scheme in dynamics. Int. J. Numer. Methods Eng. 54, 1683–1716 (2002)
Romero, I., Armero, F.: Numerical integration of the stiff dynamics of geometrically exact shells: an energy-dissipative momentum-conserving scheme. Int. J. Numer. Methods Eng. 54, 1043–1086 (2002)
Romero, I., Arnold, M.: Computing with rotations: algorithms and applications. In: Stein, E., de Borst, R., Hughes, T.J.R. (eds.) Encyclopedia of Computational Mechanics. Wiley, London (2017)
Romero, I., Gebhardt, C.G.: Variational principles for nonlinear Kirchhoff rods. Acta Mech. 231, 625–647 (2020)
Romero, I., Urrecha, M., Cyron, C.J.: A torsion-free nonlinear beam model. Int. J. Non-Linear Mech. 58, 1–10 (2014)
Sansour, C., Bednarczyk, H.: The Cosserat surface as a shell model, theory and finite-element formulation. Comput. Methods Appl. Mech. Eng. 120, 1–32 (1995)
Shi, Y., Hearst, J.E.: The Kirchhoff elastic rod, the nonlinear Schrödinger equation, and DNA supercoiling. J. Chem. Phys. 101, 5186–5200 (1994)
Simo, J.C., Fox, D.D.: On a stress resultant geometrically exact shell model. I. Formulation and optimal parametrization. Comput. Methods Appl. Mech. Eng. 72, 267–304 (1989)
Simo, J.C., Tarnow, N.: The discrete energy-momentum method. Conserving algorithms for nonlinear elastodynamics. Z. Angew. Math. Phys. 43, 757–792 (1992)
Simo, J.C., Vu-Quoc, L.: A three-dimensional finite-strain rod model. Part II: Computational aspects. Comput. Methods Appl. Mech. Eng. 58, 79–116 (1986)
Simo, J.C., Wong, K.K.: Unconditionally stable algorithms for rigid body dynamics that exactly preserve energy and momentum. Int. J. Numer. Methods Eng. 31, 19–52 (1991)
Simo, J.C., Marsden, J.E., Krishnaprasad, P.S.: The Hamiltonian structure of nonlinear elasticity: the material and convective representations of solids, rods, and plates. Arch. Ration. Mech. Anal. 104, 125–183 (1988)
Vankerschaver, J., Cantrijn, F., de León, J.M., de Diego, D.M.: Geometric aspects of nonholonomic field theories. Rep. Math. Phys. 56, 387–411 (2005)
Acknowledgements
Open Access funding provided by Projekt DEAL.
Author information
Authors and Affiliations
Corresponding author
Additional information
Communicated by Jorge Cortes.
Publisher's Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Appendices
A Second quotient of precision
To investigate the correctness of the integration scheme, we can employ the second quotient of precision (Kreiss and Ortiz 2014). Any numerical solution of an initial value problem can be written as
with \(\varvec{\xi }(t)\) representing the exact solution of the initial value problem under consideration and \(\varvec{\psi }_{i}\) for \(i=1,\ldots ,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
For the numerator, we have that
meanwhile, for the denominator, 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_{\text {II}}(t)\approx 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 \(\Vert \varvec{\psi }_{n}(t)\Vert \) 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.
B Discrete Lagrange–D’Alembert algorithm (nonholonomic integrator)
In absence of a potential energy function, the discrete Lagrange–D’Alembert equations for the forced case can be written as:
For this, we are required to introduce the following definitions:
The first discrete Lagrange–D’Alembert equation, the discrete balance equation, establishes then for the unforced case a matching that leads to the existence of a unique momentum \(\varvec{p}_i\) at time instant \(t_i\). In addition, we are also required to consider a new definition for \(\varvec{G}_d\), i.e.,
which substantially differs from the definition introduced in the context of the energy–momentum method proposed in this paper. Given the case that horizontal symmetries do exist, it is well known that this nonholonomic integrator of second order exactly preserves the associated nonholonomic momenta, namely
For more details the reader may consult (Cortés and Martínez 2001). Finally, the discrete energy at time instant \(t_{i+\frac{1}{2}}\) can be defined in terms of subsequent configurations as
In turn, the discrete Hamiltonian at time instant \(t_i\) can be determined in terms of the momenta by means of
Both the discrete energy and the discrete Hamiltonian tend to the same value when the time steps size tends to zero (Leyendecker et al. 2008).
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
About this article
Cite this article
Gebhardt, C.G., Romero, I. The Rotating Rigid Body Model Based on a Non-twisting Frame. J Nonlinear Sci 30, 3199–3233 (2020). https://doi.org/10.1007/s00332-020-09648-3
Received:
Accepted:
Published:
Issue Date:
DOI: https://doi.org/10.1007/s00332-020-09648-3
Keywords
- Rotating rigid body model
- Non-twisting frame
- Nonholonomic system
- Conservation laws
- Structure preserving integration