1 Introduction

In this paper an extended version of the intermediate axis theorem is presented and rigorously proven. These results are used to derive a novel simulation scheme which respects the stability properties of a freely rotating body. The research is being motivated by the practical application of rockfall simulation.

Rockfall is a serious natural hazard in mountainous areas and represents a major threat to infrastructure, transportation lines, and people [7]. The identification of potential rockfall starting zones, the calculation of the rock trajectories in complex three-dimensional terrain, and the interaction of falling rocks with protection measures, such as rockfall dams, catching nets, and mountain forests, are three major components of the rockfall problem and research [27]. Methods to predict rockfall trajectories in complex three-dimensional terrain (given an unstable source area) have great practical value as they can be used to determine the danger of rockfall run-out zones as well as to dimension protection measures. Although rockfall research has been dominated in the past by simple point mass and 2D models because of their computational efficiency, 3D rockfall simulation with complex shape models is now becoming standard. Several 3D codes exist (see [16] for an extensive overview), but only few consider a 3D rock shape for the inertia properties, contact detection, and contact laws. CRSP-3D (Colorado Rockfall Simulation Program) [2] employs the distinct element method together with a penalty approach. Aspherical boulder shapes, such as a rectangular cuboid or a cylinder, are approximated in [2] by a collection of rigid spheres with interconnecting elastic elements. The rockfall simulation code STAR3-D [6, 13, 14] is a fully three-dimensional rockfall simulation code with a rigid complex shape model. However, STAR3-D is limited to parallelepiped rock shapes. A full 3D simulation technique for rockfall dynamics, taking rock shape into account and using the state-of-the-art methods of multibody dynamics and nonsmooth contact dynamics, has been developed in [16]. The rockfall simulation technique is based on the nonsmooth contact dynamics method with hard contact laws, i.e., a Moreau-type timestepping method, see [1, 15, 21]. The rock is modeled as an arbitrary convex polyhedron and the terrain model is based on a high-resolution digital elevation model. The developed numerical methods have been implemented in the code RAMMS::ROCKFALL, which is being actively used in the natural hazards research community [5].

Field observations of natural rockfall events,Footnote 1 as well as high precision measurements with instrumented experimental rocks [4], have shown that platy disk-shaped rocks have the tendency to roll and bump down the slope around their major principal axis, performing a wheel-like motion. Simulations with the present scheme implemented in RAMMS::ROCKFALL, which treats gyroscopic terms in an explicit way, fail to represent the observed wheel-like motion of platy rocks around the major principal axis. Furthermore, numerical simulations with the present scheme showed that stable rotation in free flight around the major principal axis is not possible, contrary to the intermediate axis theorem. The speed of deviation from the major principal axis can be controlled by choosing a small integrator timestep, but leads to unacceptable simulation times. The stable rotation around the major principal axis of platy disk-shaped rocks enlarges the run-out distance of rockfall trajectories. The present “explicit” scheme leads to a reduced estimation of the run-out zone of wheel-shaped rocks and therefore underestimates the rockfall danger.

The intermediate axis theorem [8, 19] is a result of the Euler equations describing the movement of a rigid body with three distinct principal moments of inertia. The theorem describes the following effect: rotation of a rigid body around its minor and major principal axes is stable, while rotation around its intermediate principal axis is unstable. These stability properties where already reported by Poinsot [22]. A rigorous proof of the stability of the minor and major principal axes of the Euler equations has been given in [11] using the energy-Casimir method, i.e., a Lyapunov function based on the Hamiltonian and a suitable Casimir function (see also [12]). The classical intermediate axis theorem, however, only involves the Euler equations for the three components of the angular velocity and not the orientation of the body.

In the course of the paper, we will see that the conservation of mechanical stability properties by numerical schemes (not to be confused with numerical stability properties) is closely related to their ability to exactly keep energy and momentum as invariants of motion. The development of energy–momentum conserving schemes is a field of active research [9, 10]. We may distinguish between projective methods and schemes which directly preserve the energy and momentum. Of the latter type, we name the ALGO-C1 algorithm of Simo and Wong [26], the symplectic Discrete Moser–Veselov Algorithm [20] (being for rigid body equations equivalent to the RATTLE algorithm) [20] and the quaternion-based scheme of Betsch and Siebert [3]. However, the application and implementation of existing energy–momentum conserving schemes in rockfall simulation are not straightforward. The scheme needs to be able to deal with a singularity-free description of orientation. The scheme needs to be able to accommodate for dissipation and, at the same time, preserve invariants of motion in the dissipation-free case. Furthermore, the need to deal with frictional unilateral constraints in rockfall simulation asks to embed the energy–momentum conserving scheme in the Moreau-type timestepping scheme, which is not a trivial task. Moreover, existing rockfall simulation codes take into account a wealth of different types of dissipation mechanisms (e.g., to take into account the effect of trees and vegetation [18]). The fully implicit coupling between position and velocity updates of most energy–momentum conserving schemes makes it difficult to judge if they can be successfully combined with these specialized dissipation and contact laws.

The aims of the current paper are threefold:

  • An extended version of intermediate axis theorem for motion in the 6-dimensional state-space, which also involves the orientation of the body, is presented and rigorously proven.

  • An energy–momentum preserving scheme which preserves the stability properties in accordance with the (extended) intermediate axis theorem is presented. The novel scheme treats the gyroscopic terms in an implicit way. It is shown how this scheme can be easily incorporated in an existing rockfall simulation code, meeting all the requirements mentioned above.

  • As demonstration of the above, we briefly study the performance of the scheme in a 3D rockfall simulation software on two test cases.

The paper is organized as follows. In Sect. 2 we describe the dynamics of a freely rotating body in state-space form using as states the three angular velocity components and an arbitrary parametrization of the orientation of the body with respect to the inertial frame. Furthermore, we explicitly describe the stationary motion of which we want to assess the stability properties. Using the method of Lyapunov functions we rigorously prove in Sect. 3 an extended version of the intermediate axis theorem in the full state-space. In Sect. 4 we use the same techniques to investigate stability properties of numerical schemes. First, we prove that the present scheme, which is fully explicit during flight phases of the rock, does not respect the intermediate axis theorem. Furthermore, we present an alternative implicit scheme which correctly describes the stability properties of a freely rotating body. Results with the novel scheme implemented in RAMMS::ROCKFALL are given in Sect. 6.

2 Dynamics of a free spinning body

In this preliminary section, we first set notation in Sects. 2.1 and 2.2 and then cast the error dynamics with respect to a stationary motion in terms of an ordinary differential equation suitable for standard Lyapunov techniques in Sect. 2.3.

2.1 Equations of motion of a free spinning body

Let \(\mathcal{V}\) be a Euclidean vector space. To describe the orientation of a rigid body, we use a body-fixed frame \(K = (\vec{\mathbf{e}}_{x}^{K},\vec{\mathbf{e}}_{y}^{K}, \vec{\mathbf{e}}_{z}^{K}) \), as well as an inertial frame \(I= (\vec{\mathbf{e}}_{x}^{I},\vec{\mathbf{e}}_{y}^{I}, \vec{\mathbf{e}}_{z}^{I})\) (Fig. 1). An arbitrary vector \(\vec{\boldsymbol{a}}\in \mathcal{V}\) can be expressed in the \(K\)-frame through the tuple \(_{K}\vec{\boldsymbol{a}}\in \mathbb{R}^{3}\), which is related to its representation \(_{I}\vec{\boldsymbol{a}}\) in the inertial frame through

$$ {}_{I}\vec{\boldsymbol{a}}= \boldsymbol{A}_{IK} {}_{K}\vec{\boldsymbol{a}}, $$
(1)

where \(\boldsymbol{A}_{IK} = \begin{pmatrix} {}_{I}\vec{\mathbf{e}}_{x}^{K} & {}_{I}\vec{\mathbf{e}}_{y}^{K} & {}_{I} \vec{\mathbf{e}}_{z}^{K} \end{pmatrix} \in \textrm{SO}(3)\) is the transformation matrix describing the orientation of the body. Likewise, we introduce \(\boldsymbol{A}_{KI} = \begin{pmatrix} {}_{K}\vec{\mathbf{e}}_{x}^{I} & {}_{K}\vec{\mathbf{e}}_{y}^{I} & {}_{K} \vec{\mathbf{e}}_{z}^{I} \end{pmatrix} = \boldsymbol{A}_{IK}^{\mathrm{T}}= \boldsymbol{A}_{IK}^{-1}\).

Fig. 1
figure 1

Coordinate frames for a free spinning body

Full size image

Let \(_{K}\vec{\boldsymbol{\Omega }}\) denote the angular velocity of the body expressed in the body-fixed frame \(K\). The spin of the body (sometimes called “angular momentum with respect to the center of gravity,” “spin angular momentum,” or simply “angular momentum”) is defined by

$$ \vec{\boldsymbol{N}}_{S} = \boldsymbol{\Theta }_{S} \,\vec{\boldsymbol{\Omega }}, $$
(2)

where \(\boldsymbol{\Theta }_{S}\) is the inertia tensor. The inertia tensor takes a constant form in the body-fixed frame, which we choose to be aligned along the principal axes of inertia, such that

$$ \boldsymbol{\Theta } := {}_{K}\boldsymbol{\Theta }_{S} = \begin{pmatrix} A & 0 & 0 \\ 0 & B & 0 \\ 0 & 0 & C \end{pmatrix} ,$$
(3)

where \(A,B,C>0\) are the principal moments of inertia around the minor, intermediate, and major principal axis, respectively. We will denote \(\boldsymbol{\Theta } := {}_{K}\boldsymbol{\Theta }_{S}\) as inertia matrix to have a short-hand notation. During flight, the rotational motion of the rock is decoupled from the translational motion and fully described by the spin invariance \(\dot{\vec{\boldsymbol{N}}}_{S}=\vec{\boldsymbol{0}}\), yielding

$$ _{K}\dot{\vec{\boldsymbol{N}}}_{S} + {}_{K}\vec{\boldsymbol{\Omega }}\times {}_{K} \vec{\boldsymbol{N}}_{S} = \vec{\boldsymbol{0}} $$
(4)

in the body-fixed frame, which leads to the Euler equations for a freely rotating rigid body

$$ \boldsymbol{\Theta }\dot{\boldsymbol{\omega }}+\boldsymbol{\omega }\times (\boldsymbol{\Theta } \boldsymbol{\omega }) = \boldsymbol{0}. $$
(5)

Herein, we used the short-hand notation \(\boldsymbol{\omega }= {}_{K}\vec{\boldsymbol{\Omega }}\) for the angular velocity expressed in the body-fixed frame.

The evolution of the orientation of the body may be derived through the use of the rigid body formula on the triad

$$ \bigl(\vec{\mathbf{e}}_{x}^{K}\bigr)^{\cdot }= \vec{\boldsymbol{\Omega }} \times \vec{\mathbf{e}}_{x}^{K} ,$$
(6)

which yields in the \(I\)-frame

$$ \begin{aligned} {}_{I}\bigl(\vec{\mathbf{e}}_{x}^{K}\bigr)^{\cdot }= \bigl({}_{I}\vec{\mathbf{e}}_{x}^{K}\bigr)^{\cdot }&= {}_{I} \vec{\boldsymbol{\Omega }}\times {}_{I}\vec{\mathbf{e}}_{x}^{K} \\ &= \widetilde{ {}_{I}\vec{\boldsymbol{\Omega }}} \,{}_{I}\vec{\mathbf{e}}_{x}^{K} \\ &= \boldsymbol{A}_{IK} \widetilde{ {}_{K}\vec{\boldsymbol{\Omega }}} \, \boldsymbol{A}_{KI} \,{}_{I}\vec{\mathbf{e}}_{x}^{K} \\ &= \boldsymbol{A}_{IK} \tilde{ \boldsymbol{\omega }} \, \boldsymbol{A}_{KI} \,{}_{I} \vec{\mathbf{e}}_{x}^{K} \end{aligned} $$
(7)

and similarly for \(\vec{\mathbf{e}}_{y}^{K}\) and \(\vec{\mathbf{e}}_{z}^{K}\), where the tilde operator gives the skew-symmetric matrix \(\tilde{\boldsymbol{\omega }}\) such that \(\tilde{\boldsymbol{\omega }}\boldsymbol{c}= \boldsymbol{\omega }\times \boldsymbol{c}\) for all \(\boldsymbol{c}\in \mathbb{R}^{3}\). Hence, the time-derivative of the transformation matrix is given by

$$ \dot{\boldsymbol{A}}_{IK} = \begin{pmatrix} \bigl({}_{I}\vec{\mathbf{e}}_{x}^{K}\bigr)^{\cdot }& \bigl({}_{I} \vec{\mathbf{e}}_{y}^{K}\bigr)^{\cdot }& \bigl({}_{I}\vec{\mathbf{e}}_{z}^{K} \bigr)^{\cdot }\end{pmatrix} = \boldsymbol{A}_{IK} \tilde{ \boldsymbol{\omega }} \boldsymbol{A}_{KI} \begin{pmatrix} {}_{I}\vec{\mathbf{e}}_{x}^{K} & {}_{I}\vec{\mathbf{e}}_{y}^{K} & {}_{I} \vec{\mathbf{e}}_{z}^{K} \end{pmatrix} = \boldsymbol{A}_{IK} \tilde{ \boldsymbol{\omega }} \boldsymbol{A}_{KI} \boldsymbol{A}_{IK} ,$$
(8)

which yields the result

$$ \dot{\boldsymbol{A}}_{IK} = \boldsymbol{A}_{IK}\tilde{\boldsymbol{\omega }}. $$
(9)

The orientation of the body may be freely parameterized using, for instance, unit-quaternions \(\boldsymbol{p}\in \mathbb{R}^{4}\) or axis-angle notation (\(\boldsymbol{n}\), \(\chi \)), see [23]. The matrix differential equation (9) results in a set of ordinary differential equations for the chosen parametrization. If a quaternion \(\boldsymbol{p}= (\epsilon _{0},\boldsymbol{\epsilon }^{\mathrm{T}})^{\mathrm{T}}\) representation is used, then the rotation matrix is parametrized as

$$ \boldsymbol{A}_{IK} = \boldsymbol{I} + \frac{2}{\epsilon _{0}^{2} + \boldsymbol{\epsilon }^{\mathrm{T}}\boldsymbol{\epsilon }}( \tilde{ \boldsymbol{\epsilon }}\tilde{ \boldsymbol{\epsilon }}+\epsilon _{0} \tilde{ \boldsymbol{\epsilon }}) ,$$
(10)

resulting in

$$ \dot{\boldsymbol{p}} = \boldsymbol{F}(\boldsymbol{p}) \boldsymbol{\omega }$$
(11)

with

$$ \boldsymbol{F}(\boldsymbol{p}) = \frac{1}{2\|\boldsymbol{p}\|} \begin{pmatrix} \boldsymbol{\epsilon }^{\mathrm{T}} \\ \tilde{\boldsymbol{\epsilon }}+\epsilon _{0}\boldsymbol{I} \end{pmatrix} . $$
(12)

2.2 Invariants of motion of the Euler equations

The body has the rotational kinetic energy

$$ T = \frac{1}{2}\boldsymbol{\omega }^{\mathrm{T}}\boldsymbol{\Theta }\boldsymbol{\omega },$$
(13)

being a constant of motion as directly follows from the Euler equations (5),

$$ \frac{\mbox{d}}{\mbox{d}t}T = \boldsymbol{\omega }^{\mathrm{T}}\boldsymbol{\Theta }\dot{\boldsymbol{\omega }} = -\boldsymbol{\omega }^{\mathrm{T}}\big(\boldsymbol{\omega }\times (\boldsymbol{\Theta }\boldsymbol{\omega })\big) = 0. $$
(14)

The second invariant is the spin

$$ {}_{I}\vec{\boldsymbol{N}}_{S} = \boldsymbol{A}_{IK}\boldsymbol{\Theta }\boldsymbol{\omega }, $$
(15)

which is constant in the inertial frame, i.e., using (15) and (9),

$$ \frac{\mbox{d}}{\mbox{d}t}\left ({}_{I}\vec{\boldsymbol{N}}_{S} \right ) = \dot{\boldsymbol{A}}_{IK}\boldsymbol{\Theta }\boldsymbol{\omega }+ \boldsymbol{A}_{IK} \boldsymbol{\Theta }\dot{\boldsymbol{\omega }} = \boldsymbol{A}_{IK}\left (\boldsymbol{\omega } \times (\boldsymbol{\Theta }\boldsymbol{\omega })\right ) - \boldsymbol{A}_{IK}\left ( \boldsymbol{\omega }\times (\boldsymbol{\Theta }\boldsymbol{\omega })\right )= \boldsymbol{0}. $$
(16)

When expressed in the body-fixed frame, the spin \({}_{K}\vec{\boldsymbol{N}}_{S}\) is not constant but keeps a constant magnitude as follows from

$$ \frac{\mbox{d}}{\mbox{d}t}||{}_{K}\vec{\boldsymbol{N}}_{S}||^{2}= \frac{\mbox{d}}{\mbox{d}t}({}_{K}\vec{\boldsymbol{N}}_{S}^{\mathrm{T}}{}_{K}\vec{\boldsymbol{N}}_{S}) = 2(\boldsymbol{\Theta }\boldsymbol{\omega })^{\mathrm{T}}\boldsymbol{\Theta }\dot{\boldsymbol{\omega }} = -2(\boldsymbol{\Theta } \boldsymbol{\omega })^{\mathrm{T}}\big(\boldsymbol{\omega }\times (\boldsymbol{\Theta } \boldsymbol{\omega })\big) = \boldsymbol{0}. $$
(17)

2.3 Stationary motion

A rigid body may undergo a stationary motion for which its angular velocity is constant, i.e., \(\dot{\boldsymbol{\omega }}=\boldsymbol{0}\). We will denote such a stationary motion in state-space as \((\boldsymbol{A}_{IK\star }(t), \boldsymbol{\omega }_{\star })\). From the Euler equations (5), we infer that stationary motion is only possible if the term \(\boldsymbol{\omega }\times (\boldsymbol{\Theta }\boldsymbol{\omega })\) vanishes. Stationary motion therefore implies that \(\boldsymbol{\omega }_{\star }\) is in an eigendirection of \(\boldsymbol{\Theta }\), resulting in three stationary directions of motion \({}_{K}\vec{\mathbf{e}}_{x}^{K},{}_{K}\vec{\mathbf{e}}_{y}^{K},{}_{K} \vec{\mathbf{e}}_{z}^{K}\). Without loss of generality, let \(\boldsymbol{\omega }_{\star }= \Omega \boldsymbol{e}_{\star }\), where \(\boldsymbol{e}_{\star }= \boldsymbol{e}_{3}= \begin{bmatrix} 0 & 0 & 1 \end{bmatrix} ^{\mathrm{T}}\) agrees with \({}_{K}\vec{\mathbf{e}}_{z}^{K}\). The vector \(\boldsymbol{e}_{\star }\) may be complemented by \(\boldsymbol{e}_{1} = \begin{bmatrix} 1 & 0 & 0 \end{bmatrix} ^{\mathrm{T}}\) and \(\boldsymbol{e}_{2} = \begin{bmatrix} 0 & 1 & 0 \end{bmatrix} ^{\mathrm{T}}\) to form an orthonormal basis. The evolution of the orientation of the body during stationary motion

$$ \boldsymbol{A}_{IK\star }(t) = \boldsymbol{A}_{IK\star }(0)\mathrm{e}^{ \tilde{\boldsymbol{\omega }}_{\star }t} $$
(18)

follows from the closed form solution of the matrix differential equation (9). Without loss of generality, we set \(\boldsymbol{A}_{IK\star }(0) = \boldsymbol{I}\). From

$$ {}_{I}\vec{\mathbf{e}}_{z}^{K\star }(t) = \boldsymbol{A}_{IK\star }(t)\boldsymbol{e}_{3} = \mathrm{e}^{\tilde{\boldsymbol{\omega }}_{\star }t}\boldsymbol{e}_{\star }= \boldsymbol{e}_{\star }= {}_{I}\vec{\mathbf{e}}_{z}^{I} $$
(19)

it follows that \(\vec{\mathbf{e}}_{z}^{K\star } = \vec{\mathbf{e}}_{z}^{I}\) for all \(t\). The stationary motion \((\boldsymbol{A}_{IK\star }(t), \boldsymbol{\omega }_{\star })\) itself cannot be stable, irrespective of the principal axis which is considered, because a small error \(\Delta \Omega \) in the magnitude of the angular velocity \(\boldsymbol{\omega }= (\Omega +\Delta \Omega )\boldsymbol{e}_{\star }\) will cause \(\boldsymbol{A}_{IK}(t)\) to diverge from \(\boldsymbol{A}_{IK\star }(t)\). Hence, instead of the stationary motion, we need to study the stability of the axis of rotation, or, more precisely, of a manifold in state-space related to that. Hereto, we consider the difference between the axes of rotation \(\vec{\mathbf{e}}_{z}^{K\star } = \vec{\mathbf{e}}_{z}^{I}\) of the stationary motion and \(\vec{\mathbf{e}}_{z}^{K}\) of an arbitrary motion, which we express in the \(K\)-frame using the quantity

$$ \boldsymbol{d}(t) = {}_{K}\vec{\mathbf{e}}_{z}^{K\star } - {}_{K} \vec{\mathbf{e}}_{z}^{K} = {}_{K}\vec{\mathbf{e}}_{z}^{I} - {}_{K} \vec{\mathbf{e}}_{z}^{K} = \boldsymbol{A}_{KI}(t) \boldsymbol{e}_{\star }- \boldsymbol{e}_{\star }. $$
(20)

Furthermore, to parametrize the transformation matrix \(\boldsymbol{A}_{IK}(t)\), we also introduce the time-dependent quantities

$$\begin{gathered} \boldsymbol{h}_{1}(t) = \boldsymbol{A}_{KI}(t) \boldsymbol{e}_{1} - \boldsymbol{e}_{1}, \end{gathered}$$
(21)
$$\begin{gathered} \boldsymbol{h}_{2}(t) = \boldsymbol{A}_{KI}(t) \boldsymbol{e}_{2} - \boldsymbol{e}_{2}. \end{gathered}$$
(22)

The transformation matrix may therefore be expressed as

$$ \begin{aligned} \boldsymbol{A}_{IK}^{\mathrm{T}}(t) &= \boldsymbol{I} + \begin{bmatrix} \boldsymbol{h}_{1}(t) & \boldsymbol{h}_{2}(t) & \boldsymbol{d}(t) \end{bmatrix} \begin{bmatrix} \boldsymbol{e}_{1} & \boldsymbol{e}_{2} & \boldsymbol{e}_{\star }\end{bmatrix} ^{-1} \\ &= \boldsymbol{I} + \begin{bmatrix} \boldsymbol{h}_{1}(t) & \boldsymbol{h}_{2}(t) & \boldsymbol{d}(t) \end{bmatrix} . \end{aligned} $$
(23)

Furthermore, we introduce the quantity

$$ \boldsymbol{\alpha }(t) = \boldsymbol{\omega }(t) - \boldsymbol{\omega }_{\star }$$
(24)

to express the difference in rotation speed, which is being governed by the Euler equations (5)

$$ \boldsymbol{\Theta }\dot{\boldsymbol{\alpha }}+(\boldsymbol{\omega }_{\star }+ \boldsymbol{\alpha }) \times (\boldsymbol{\Theta }(\boldsymbol{\omega }_{\star }+ \boldsymbol{\alpha })) = \boldsymbol{0}. $$
(25)

Having introduced the quantities \(\boldsymbol{d}\) and \(\boldsymbol{\alpha }\), we can define the manifold of stationary rotation in a straightforward way by

$$ \mathcal{M} = \left \{ (\boldsymbol{A}_{IK}, \boldsymbol{\omega }) \in \textrm{SO}(3) \times \mathbb{R}^{3} \mid \boldsymbol{d}= \boldsymbol{0}, \boldsymbol{\alpha }= \boldsymbol{0} \right \} . $$
(26)

The dynamics of the difference \(\boldsymbol{d}(t)\) to the axis of stationary rotation is given by

$$ \begin{aligned} \dot{\boldsymbol{d}}(t) = \dot{\boldsymbol{A}}_{KI}(t) \boldsymbol{e}_{\star }&= \left ( \boldsymbol{A}_{IK}(t)\tilde{\boldsymbol{\omega }}(t)\right )^{\mathrm{T}}\boldsymbol{e}_{\star }\\ &=-\tilde{\boldsymbol{\omega }}(t)\boldsymbol{A}_{KI}(t) \boldsymbol{e}_{\star }\\ &= -\left (\tilde{\boldsymbol{\omega }}_{\star }+\tilde{\boldsymbol{\alpha }}(t) \right )\left (\boldsymbol{d}(t)+\boldsymbol{e}_{\star }\right ) \\ &= -\tilde{\boldsymbol{\omega }}_{\star }\boldsymbol{d}(t) -\tilde{\boldsymbol{\alpha }}(t) \left (\boldsymbol{d}(t)+\boldsymbol{e}_{\star }\right ), \end{aligned} $$
(27)

being only dependent on \(\boldsymbol{d}\) and \(\boldsymbol{\alpha }\). It holds that \(\dot{\boldsymbol{d}}(t) = \boldsymbol{0}\) if \(\boldsymbol{d}=\boldsymbol{\alpha }= \boldsymbol{0}\) and \(\dot{\boldsymbol{\alpha }}(t) = \boldsymbol{0}\) if \(\boldsymbol{\alpha }= \boldsymbol{0}\). Hence, the manifold ℳ of stationary rotation is invariant. The differential equations (27) and (25) can be gathered using \(\boldsymbol{y}(t) = \begin{pmatrix} \boldsymbol{d}^{\mathrm{T}}& \boldsymbol{\alpha }^{\mathrm{T}}\end{pmatrix} ^{\mathrm{T}}\) in the system of ordinary differential equations

$$ \dot{\boldsymbol{y}}(t) = \boldsymbol{f}(\boldsymbol{y}(t)), $$
(28)

which is time-autonomous. The stability of the axis of rotation now bears down to the stability of the invariant manifold ℳ, i.e., the stability of the equilibrium \(\boldsymbol{y}^{\star }= \boldsymbol{0}\) of system (28). The stability properties of stationary rotation in the 6-dimensional state-space \(\boldsymbol{y}\) will be referred to as the extended intermediate axis theorem.

3 Extended intermediate axis theorem

3.1 Stability of rotation around the major principal axis

We consider motion in the vicinity of the stationary rotation \(\boldsymbol{\omega }_{\star }= \Omega \boldsymbol{e}_{\star }\), \(\Omega >0\), and \(C\geq \max (A,B)\) such that \(\boldsymbol{e}_{\star }\) is the major principal axis of inertia. Herein, the \(K\star \)-frame is the body-fixed frame of stationary motion, whereas we will reserve the \(K\)-frame for the body-fixed frame of an arbitrary motion in the vicinity of the stationary motion. The frames are related through \(\boldsymbol{A}_{KK\star } = \boldsymbol{A}_{KI}\boldsymbol{A}_{IK\star } = \boldsymbol{A}_{IK}^{\mathrm{T}}\boldsymbol{A}_{IK\star }\) and it therefore holds that

$$ \boldsymbol{A}_{KK\star }\boldsymbol{e}_{\star }= \boldsymbol{A}_{KI}\boldsymbol{A}_{IK\star }\boldsymbol{e}_{\star }= \boldsymbol{A}_{KI}\boldsymbol{e}_{\star }= \boldsymbol{d}+ \boldsymbol{e}_{\star }. $$
(29)

The following theorem may be viewed as an extended intermediate axis theorem, as it not only proves that the angular velocity \(\boldsymbol{\omega }\) remains close to \(\boldsymbol{\omega }_{\star }\) (i.e., the classical intermediate axis theorem), but also proves that \(\boldsymbol{d}\) remains small, i.e., the orientation of the axis of rotation is stable.

Theorem 3.1

Stability of rotation around the major principal axis

Let \(C\geq \max (A,B)\) be such that \(\boldsymbol{e}_{\star }\) is the major principal axis of inertia. The invariant manifold

$$ \mathcal{M} = \left \{ (\boldsymbol{A}_{IK}, \boldsymbol{\omega }) \in \textrm{SO}(3) \times \mathbb{R}^{3} \mid \boldsymbol{d}= \boldsymbol{0}, \boldsymbol{\alpha }= \boldsymbol{0} \right \} $$

is Lyapunov stable.

Proof

In order to set up a Lyapunov function for the equilibrium \(\boldsymbol{y}^{\star }= \boldsymbol{0}\) of system (28), we consider the function

$$ \bar{V}(\boldsymbol{d},\boldsymbol{\alpha }) = \| {}_{I}\vec{\boldsymbol{N}}_{S} - {}_{I} \vec{\boldsymbol{N}}_{S\star }\|^{2} = \| {}_{K}\vec{\boldsymbol{N}}_{S} - {}_{K} \vec{\boldsymbol{N}}_{S\star }\|^{2}. $$
(30)

The spin \(\vec{\boldsymbol{N}}_{S\star }\) can be easily expressed in the \(K\star \)-frame as \({}_{K\star }\vec{\boldsymbol{N}}_{S\star } = \boldsymbol{\Theta } \boldsymbol{\omega }_{\star }= C\Omega \boldsymbol{e}_{\star }\), which can be cast in the \(K\)-frame through

$$ {}_{K}\vec{\boldsymbol{N}}_{S\star } = \boldsymbol{A}_{KK\star } \,{}_{K\star } \vec{\boldsymbol{N}}_{S\star } = C\Omega \, \boldsymbol{A}_{KK\star }\boldsymbol{e}_{\star }= C \Omega (\boldsymbol{d}+ \boldsymbol{e}_{\star }). $$
(31)

Using \({}_{K}\vec{\boldsymbol{N}}_{S} = \boldsymbol{\Theta }\boldsymbol{\omega }= \boldsymbol{\Theta }( \boldsymbol{\alpha }+\boldsymbol{\omega }_{\star }) = \boldsymbol{\Theta }\boldsymbol{\alpha }+ C \Omega \boldsymbol{e}_{\star }\), we arrive at

$$ \bar{V}(\boldsymbol{d},\boldsymbol{\alpha }) = \| \boldsymbol{\Theta }\boldsymbol{\alpha }- C \Omega \boldsymbol{d}\|^{2} = (A\alpha _{x} - C\Omega d_{x})^{2} + (B\alpha _{y} - C\Omega d_{y})^{2} + (C\alpha _{z} - C\Omega d_{z})^{2}. $$
(32)

Furthermore, we introduce the function

$$ \hat{V}(\boldsymbol{\alpha }) = 2C\;T(\boldsymbol{\omega })-\|\vec{\boldsymbol{N}}_{S}( \boldsymbol{\omega })\|^{2} + \frac{1}{\Omega ^{2}}(2\;T(\boldsymbol{\omega })-2\;T( \boldsymbol{\omega }_{\star }))^{2} \;, $$
(33)

which can be expressed as

$$ \begin{aligned} \hat{V}(\boldsymbol{\alpha }) &= C\boldsymbol{\omega }^{\mathrm{T}}\boldsymbol{\Theta } \boldsymbol{\omega }-\boldsymbol{\omega }^{\mathrm{T}}\boldsymbol{\Theta }^{2}\boldsymbol{\omega }+ \frac{1}{\Omega ^{2}}(\boldsymbol{\omega }^{\mathrm{T}}\boldsymbol{\Theta } \boldsymbol{\omega }-\boldsymbol{\omega }_{\star }^{\mathrm{T}}\boldsymbol{\Theta } \boldsymbol{\omega }_{\star })^{2} \\ &=(\boldsymbol{\alpha }+ \boldsymbol{\omega }_{\star })^{\mathrm{T}}(C\boldsymbol{\Theta } - \boldsymbol{\Theta }^{2}) (\boldsymbol{\alpha }+ \boldsymbol{\omega }_{\star }) + \frac{1}{\Omega ^{2}}\left ((\boldsymbol{\alpha }+ \boldsymbol{\omega }_{\star })^{\mathrm{T}}\boldsymbol{\Theta }(\boldsymbol{\alpha }+ \boldsymbol{\omega }_{\star })- \boldsymbol{\omega }_{\star }^{\mathrm{T}}\boldsymbol{\Theta }\boldsymbol{\omega }_{\star }\right )^{2} \\ &= \boldsymbol{\alpha }^{\mathrm{T}}(C\boldsymbol{\Theta } - \boldsymbol{\Theta }^{2}) \boldsymbol{\alpha }+ \frac{1}{\Omega ^{2}}\left (\boldsymbol{\alpha }^{\mathrm{T}} \boldsymbol{\Theta }\boldsymbol{\alpha }+ 2C\boldsymbol{\omega }_{\star }^{\mathrm{T}} \boldsymbol{\alpha }\right )^{2} \\ &=A(C-A)\alpha _{x}^{2}+B(C-B)\alpha _{y}^{2}+\frac{1}{\Omega ^{2}}(A \alpha _{x}^{2}+B\alpha _{y}^{2}+C\alpha _{z}^{2}+2C\Omega \alpha _{z})^{2}. \end{aligned} $$
(34)

To prove stability of the trivial equilibrium of (28) for \(C\geq \max (A,B)\), we consider the Lyapunov candidate function

$$ V(\boldsymbol{y}) = \bar{V}(\boldsymbol{d},\boldsymbol{\alpha }) + \hat{V}(\boldsymbol{\alpha }), $$
(35)

which purely consists of the invariants of motion \(\vec{\boldsymbol{N}}_{S} = \) const. and \(T(\boldsymbol{\omega }) = \) const. and is therefore constant, i.e., it holds that \(\dot{V} = 0\) along solutions of the system. The (local) positive definiteness of \(V\) remains to be investigated in the following.

First, we prove that the function \(\hat{V}(\boldsymbol{\alpha })\) is a locally positive definite function of its argument \(\boldsymbol{\alpha }\). We infer that \(\hat{V}(\boldsymbol{\alpha })\geq 0\) for arbitrary \(\boldsymbol{\alpha }\in \mathbb{R}^{3}\) as it is a sum of squares with positive coefficients for \(C\geq \max (A,B)\). Moreover, the points where \(\hat{V}(\boldsymbol{\alpha })= 0\) are characterized by \(\alpha _{x} = 0\), \(\alpha _{y}=0\) and \(\alpha _{z}(\alpha _{z}+2\Omega ) = 0\). This implies that \(\hat{V}\) vanishes at the origin and at the point \(\boldsymbol{\alpha }= (0\ 0\ -2\Omega )^{\mathrm{T}}\) and is strictly positive for all other \(\boldsymbol{\alpha }\in \mathbb{R}^{3}\), proving that \(\hat{V}(\boldsymbol{\alpha })\) is locally positive definite.

As \(V\) is the sum of \(\bar{V}(\boldsymbol{d},\boldsymbol{\alpha })\geq 0\) and \(\hat{V}(\boldsymbol{\alpha })\geq 0\), it may only vanish if \(\bar{V}(\boldsymbol{d},\boldsymbol{\alpha })\) and \(\hat{V}(\boldsymbol{\alpha })\) vanish simultaneously. From the local positive definiteness of \(\hat{V}(\boldsymbol{\alpha })\) it is clear that, in the neighborhood of the origin, \(V\) may only vanish for \(\boldsymbol{\alpha }= \boldsymbol{0}\).

Now we consider \(\bar{V}(\boldsymbol{d},\boldsymbol{\alpha })\) and note that \(\bar{V}(\boldsymbol{d},\boldsymbol{0}) = C^{2}\Omega ^{2}\|\boldsymbol{d}\|^{2}\) can only vanish if \(\boldsymbol{d}= \boldsymbol{0}\). This proves local positive definiteness of \(V\) and, thereby, that rotation around the major principal axis is stable. □

Remark 1

In [11] the stability of stationary rotation around the major principal axis of the 3-dimensional Euler equations (5), or equivalently of the error dynamics (25), was proven in a coordinate-free way using a Lyapunov function constructed using the sum of the Hamiltonian \(T(\boldsymbol{\omega })\) and a suitable Casimir function \(\Phi (\frac{1}{2}\|\vec{\boldsymbol{N}}_{S}(\boldsymbol{\omega })\|^{2})\). The result of [11] would have given an alternative choice for the function \(\hat{V}(\boldsymbol{\alpha })\) in Theorem 3.1, being constant and locally positive definite in \(\boldsymbol{\alpha }\). Here, the aim is to derive a numerical scheme and we have refrained from a coordinate-free description. The choice (34) for the function \(\hat{V}(\boldsymbol{\alpha })\) involves less powers of the components of \(\boldsymbol{\alpha }\) than a construction by using a Casimir function.

3.2 Stability of rotation around the minor principal axis

We now consider the stability of stationary rotation around the minor principal axis by setting again \(\boldsymbol{\omega }_{\star }= \Omega \boldsymbol{e}_{\star }\), \(\Omega >0\) with \(\boldsymbol{e}_{\star }=\boldsymbol{e}_{3}\), but assuming \(C\leq \min (A,B)\). The proof of the stability of stationary rotation around the minor principal axis is completely analogous to the proof for the major principal axis.

Theorem 3.2

Stability of rotation around the minor principal axis

Let \(C\leq \min (A,B)\) be such that \(\boldsymbol{e}_{\star }\) is the minor principal axis of inertia. The invariant manifold

$$ \mathcal{M} = \left \{ (\boldsymbol{A}_{IK}, \boldsymbol{\omega }) \in \textrm{SO}(3) \times \mathbb{R}^{3} \mid \boldsymbol{d}= \boldsymbol{0}, \boldsymbol{\alpha }= \boldsymbol{0} \right \} $$

is Lyapunov stable.

Proof

We consider again the Lyapunov candidate function of the form

$$ V(\boldsymbol{y}) = \bar{V}(\boldsymbol{d},\boldsymbol{\alpha }) + \hat{V}(\boldsymbol{\alpha }), $$
(36)

where the function

$$ \bar{V}(\boldsymbol{d},\boldsymbol{\alpha }) = \| \boldsymbol{\Theta }\boldsymbol{\alpha }- C \Omega \boldsymbol{d}\|^{2} = (A\alpha _{x} - C\Omega d_{x})^{2} + (B\alpha _{y} - C\Omega d_{y})^{2} + (C\alpha _{z} - C\Omega d_{z})^{2} $$
(37)

is defined as before, but \(\hat{V}(\boldsymbol{\alpha })\) is chosen as

$$ \begin{aligned} \hat{V}(\boldsymbol{\alpha }) &= -C\boldsymbol{\omega }^{\mathrm{T}}\boldsymbol{\Theta } \boldsymbol{\omega }+\boldsymbol{\omega }^{\mathrm{T}}\boldsymbol{\Theta }^{2}\boldsymbol{\omega }+ \frac{1}{\Omega ^{2}}(\boldsymbol{\omega }^{\mathrm{T}}\boldsymbol{\Theta } \boldsymbol{\omega }-\boldsymbol{\omega }_{\star }^{\mathrm{T}}\boldsymbol{\Theta } \boldsymbol{\omega }_{\star })^{2} \\ &=A(A-C)\alpha _{x}^{2}+B(B-C)\alpha _{y}^{2}+\frac{1}{\Omega ^{2}}(A \alpha _{x}^{2}+B\alpha _{y}^{2}+C\alpha _{z}^{2}+2C\Omega \alpha _{z})^{2}. \end{aligned} $$
(38)

The Lyapunov function is constant along solutions and is locally positive definite for \(C\leq \min (A,B)\), proving stability of rotation around the minor principal axis. □

3.3 Instability of rotation around the intermediate axis

Within Lyapunov stability theory, instability of equilibria is usually proven using Chetayev’s instability theorem, which involves a continuously differentiable function \(V\) for which \(\dot{V}>0\) in a wedge \(V>0\) near the origin. The construction of such a function can therefore not be based on invariants of the system. For the instability proof, we therefore resort to an eigenvalue analysis. The angular velocity dynamics (25) is decoupled from the orientation dynamics (27). It therefore suffices to prove instability of the equilibrium \(\boldsymbol{\alpha }= \boldsymbol{0}\) of (25), similar to the classical intermediate axis theorem, of which we briefly show the main steps.

Rearranging (25) leads to

$$ \begin{aligned} \dot{\boldsymbol{\alpha }} &= -\boldsymbol{\Theta }^{-1}\big((\boldsymbol{\omega }_{\star }+ \boldsymbol{\alpha })\times (\boldsymbol{\Theta }\boldsymbol{\omega }_{\star }+ \boldsymbol{\Theta }\boldsymbol{\alpha })\big) \\ &= -\boldsymbol{\Theta }^{-1}\big(\boldsymbol{\omega }_{\star }\times (\boldsymbol{\Theta } \boldsymbol{\omega }_{\star }) +\boldsymbol{\omega }_{\star }\times (\boldsymbol{\Theta } \boldsymbol{\alpha })+\boldsymbol{\alpha }\times (\boldsymbol{\Theta }\boldsymbol{\omega }_{\star }+ \boldsymbol{\Theta }\boldsymbol{\alpha })\big) \\ &= -\boldsymbol{\Theta }^{-1}\big(\boldsymbol{\omega }_{\star }\times (\boldsymbol{\Theta } \boldsymbol{\alpha })+\boldsymbol{\alpha }\times (\boldsymbol{\Theta }\boldsymbol{\omega }_{\star }) + \boldsymbol{\alpha }\times (\boldsymbol{\Theta }\boldsymbol{\alpha })\big). \end{aligned} $$
(39)

Neglecting terms of second-order leads to the linearized dynamics

$$ \begin{aligned} \dot{\boldsymbol{\alpha }} = -\boldsymbol{\Theta }^{-1}\big(\boldsymbol{\omega }_{\star }\times (\boldsymbol{\Theta }\boldsymbol{\alpha })+\boldsymbol{\alpha }\times ( \boldsymbol{\Theta }\boldsymbol{\omega }_{\star })\big) =: -\boldsymbol{\Theta }^{-1}\boldsymbol{A} \boldsymbol{\alpha }. \end{aligned} $$
(40)

The newly defined matrix \(\boldsymbol{A}\) can be obtained from

$$ \boldsymbol{\omega }_{\star }\times (\boldsymbol{\Theta }\boldsymbol{\alpha })+\boldsymbol{\alpha } \times (\boldsymbol{\Theta }\boldsymbol{\omega }_{\star }) = \boldsymbol{A}\boldsymbol{\alpha }$$
(41)

and, therefore, it holds that

$$ \boldsymbol{A} = \boldsymbol{\omega }_{\star }\boldsymbol{\Theta }- \widetilde{\boldsymbol{\Theta }\boldsymbol{\omega }_{\star }}= \begin{pmatrix} 0 & \Omega (C-B) & 0 \\ \Omega (A-C) & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix} . $$
(42)

The eigenvalues of the error dynamics, i.e., the roots \(\lambda \) of the characteristic polynomial

$$ \begin{aligned} \det (\lambda \boldsymbol{I}-(-\boldsymbol{\Theta }^{-1}\boldsymbol{A})) &= \det \begin{pmatrix} \lambda & \frac{\Omega }{A}(C-B) & 0 \\ \frac{\Omega }{B}(A-C) & \lambda & 0 \\ 0 & 0 & \lambda \end{pmatrix} \\ &=\lambda \left (\lambda ^{2}-\frac{\Omega ^{2}}{AB}(A-C)(C-B)\right ) \end{aligned} $$
(43)

are given by

$$ \lambda _{1} = 0,\quad \lambda _{2} = \sqrt{\frac{\Omega ^{2}}{AB}(A-C)(C-B)}, \quad \lambda _{3} = -\sqrt{\frac{\Omega ^{2}}{AB}(A-C)(C-B)}\;. $$
(44)

We see that the eigenvalue analysis is inconclusive for rotation around the major principal axis (\(C\geq \max (A,B)\)) and minor principal axis \(C\leq \min (A,B)\) as all eigenvalues have a vanishing real part. For rotation around the intermediate axis (\(A< C< B\) or \(A>C>B\)), all eigenvalues are real with \(\lambda _{2}>0\) and \(\lambda _{3}<0\). Stationary rotation around the intermediate axis is therefore a saddle point and unstable.

4 Numerical schemes for a free rotating body

We first briefly discuss the fully explicit scheme (or pseudoimplicit scheme) which is currently implemented in RAMMS::ROCKFALL in Sect. 4.1 and subsequently present a novel scheme in Sect. 4.2. A numerical comparison is given in Sect. 4.3.

4.1 Fully explicit scheme

Let \(\boldsymbol{\omega }^{k}\) denote the angular velocity of the rock at time instant \(t^{k}\). The present scheme calculates the angular velocity \(\boldsymbol{\omega }^{k+1}\) at time instant \(t^{k+1} = t^{k}+\Delta t\) (in the absence of contact with the terrain) through

$$ \boldsymbol{\Theta }(\boldsymbol{\omega }^{k+1}-\boldsymbol{\omega }^{k})+ \frac{\Delta t}{2}\boldsymbol{G}(\boldsymbol{\omega }^{k})(\boldsymbol{\omega }^{k+1}+ \boldsymbol{\omega }^{k}) = 0, $$
(45)

where \(\boldsymbol{G}(\boldsymbol{\omega })=\boldsymbol{\Theta }\tilde{\boldsymbol{\omega }}+ \tilde{\boldsymbol{\omega }}\boldsymbol{\Theta }\). The resulting explicit velocity update is

$$ \begin{aligned} \boldsymbol{\omega }^{k+1} &= \left ( \boldsymbol{\Theta } + \frac{\Delta t}{2}\boldsymbol{G}(\boldsymbol{\omega }^{k})\right )^{-1}\left ( \boldsymbol{\Theta } - \frac{\Delta t}{2}\boldsymbol{G}(\boldsymbol{\omega }^{k})\right ) \boldsymbol{\omega }^{k} \\ &= \boldsymbol{\omega }^{k} -\Delta t\left ( \boldsymbol{\Theta } + \frac{\Delta t}{2}\boldsymbol{G}(\boldsymbol{\omega }^{k})\right )^{-1} \boldsymbol{G}( \boldsymbol{\omega }^{k}) \boldsymbol{\omega }^{k} , \end{aligned} $$
(46)

where the last step is the result of a telescopic expansion. The rotation matrix is parametrized using a quaternion \(\boldsymbol{p}^{k}\) with a midpoint update rule

$$ \boldsymbol{p}_{\textrm{pre}}^{k+1} = \boldsymbol{p}^{k+\frac{1}{2}} + \frac{\Delta t}{2} \boldsymbol{F}(\boldsymbol{p}^{k+\frac{1}{2}})\boldsymbol{\omega }^{k+1}, \qquad \boldsymbol{p}^{k+1} = \frac{\boldsymbol{p}_{\textrm{pre}}^{k+1}}{\|\boldsymbol{p}_{\textrm{pre}}^{k+1}\|}, $$
(47)

where \(\boldsymbol{p}^{k+\frac{1}{2}} = \boldsymbol{p}^{k} + \frac{\Delta t}{2}\boldsymbol{F}( \boldsymbol{p}^{k})\boldsymbol{\omega }^{k}\).

The rationale behind the explicit scheme is that it preserves the kinetic energy. The change in kinetic energy over the time step is

$$ \begin{aligned} T(\boldsymbol{\omega }^{k+1})-T(\boldsymbol{\omega }^{k}) &= \frac{1}{2}( \boldsymbol{\omega }^{k+1})^{\mathrm{T}}\boldsymbol{\Theta }\boldsymbol{\omega }^{k+1}- \frac{1}{2}(\boldsymbol{\omega }^{k})^{\mathrm{T}}\boldsymbol{\Theta }\boldsymbol{\omega }^{k} \\ &= \frac{1}{2}(\boldsymbol{\omega }^{k+1} + \boldsymbol{\omega }^{k})^{\mathrm{T}} \boldsymbol{\Theta }(\boldsymbol{\omega }^{k+1}-\boldsymbol{\omega }^{k}). \end{aligned} $$
(48)

Substitution of the explicit scheme (45) yields

$$ T(\boldsymbol{\omega }^{k+1})-T(\boldsymbol{\omega }^{k}) = - \frac{\Delta t}{4}( \boldsymbol{\omega }^{k+1} + \boldsymbol{\omega }^{k})^{\mathrm{T}}\boldsymbol{G}( \boldsymbol{\omega }^{k})(\boldsymbol{\omega }^{k+1}+\boldsymbol{\omega }^{k}) = 0 $$
(49)

due to the skew-symmetry of \(\boldsymbol{G}(\boldsymbol{\omega })\) which shows that the kinetic energy is conserved. A stationary solution \(\boldsymbol{\omega }_{\star }=\boldsymbol{\omega }^{k}=\boldsymbol{\omega }^{k+1}\) of the scheme respects

$$ \boldsymbol{G}(\boldsymbol{\omega }_{\star })\boldsymbol{\omega }_{\star }= \boldsymbol{\omega }_{\star }\times (\boldsymbol{\Theta }\boldsymbol{\omega }_{\star })=0, $$
(50)

and corresponds to a stationary rotation of the Euler equations around principal axes.

Let \(\boldsymbol{\alpha }^{k}\) be the perturbation of the angular velocity with respect to the stationary rotation \(\boldsymbol{\omega }_{\star }= \Omega \boldsymbol{e}_{\star }\) around the principal axis \(\boldsymbol{e}_{\star }= \boldsymbol{e}_{3}\). The perturbation dynamics is obtained using \(\boldsymbol{\omega }^{k} = \boldsymbol{\omega }_{\star }+ \boldsymbol{\alpha }^{k}\) as

$$ \boldsymbol{\Theta }(\boldsymbol{\alpha }^{k+1}-\boldsymbol{\alpha }^{k})+ \frac{\Delta t}{2}\boldsymbol{G}(\boldsymbol{\omega }_{\star })(\boldsymbol{\alpha }^{k+1}+ \boldsymbol{\alpha }^{k}) +\frac{\Delta t}{2}\boldsymbol{G}(\boldsymbol{\alpha }^{k})( \boldsymbol{\alpha }^{k+1}+\boldsymbol{\alpha }^{k}) + \Delta t \boldsymbol{G}( \boldsymbol{\alpha }^{k})\boldsymbol{\omega }_{\star }= 0, $$
(51)

where the linearity of \(\boldsymbol{G}\) and (50) have been used. For small perturbations we can neglect higher order terms in \(\boldsymbol{\alpha }^{k}\) giving the linearized perturbation dynamics

$$ \begin{aligned} \boldsymbol{0} &=\boldsymbol{\Theta }(\boldsymbol{\alpha }^{k+1}-\boldsymbol{\alpha }^{k})+ \frac{\Delta t}{2}\boldsymbol{G}(\boldsymbol{\omega }_{\star })(\boldsymbol{\alpha }^{k+1}+ \boldsymbol{\alpha }^{k}) + \Delta t \boldsymbol{G}(\boldsymbol{\alpha }^{k})\boldsymbol{\omega }_{\star }\\ &= \left [\boldsymbol{\Theta }+\frac{\Delta t}{2}\boldsymbol{G}(\boldsymbol{\omega }_{\star })\right ](\boldsymbol{\alpha }^{k+1}-\boldsymbol{\alpha }^{k})+\Delta t\big( \boldsymbol{G}(\boldsymbol{\alpha }^{k})\boldsymbol{\omega }_{\star }+\boldsymbol{G}(\boldsymbol{\omega }_{\star })\boldsymbol{\alpha }^{k}\big). \end{aligned} $$
(52)

The linearized perturbation dynamics can be solved for \(\boldsymbol{\alpha }^{k+1}\) explicitly. Hereto, the second term is reformulated as

$$ \begin{aligned} \boldsymbol{G}(\boldsymbol{\alpha }^{k})\boldsymbol{\omega }_{\star }+\boldsymbol{G}(\boldsymbol{\omega }_{\star })\boldsymbol{\alpha }^{k} &= \boldsymbol{\Theta }\tilde{\boldsymbol{\alpha }}^{k} \boldsymbol{\omega }_{\star }+\tilde{\boldsymbol{\alpha }}^{k}\boldsymbol{\Theta } \boldsymbol{\omega }_{\star }+\boldsymbol{\Theta }\tilde{\boldsymbol{\omega }}_{\star }\boldsymbol{\alpha }^{k}+\tilde{\boldsymbol{\omega }}_{\star }\boldsymbol{\Theta } \boldsymbol{\alpha }^{k} \\ &= \tilde{\boldsymbol{\alpha }}^{k}\boldsymbol{\Theta }\boldsymbol{\omega }_{\star }+ \tilde{\boldsymbol{\omega }}_{\star }\boldsymbol{\Theta }\boldsymbol{\alpha }^{k} \\ &= \boldsymbol{\alpha }^{k}\times (\boldsymbol{\Theta }\boldsymbol{\omega }_{\star })+ \boldsymbol{\omega }_{\star }\times (\boldsymbol{\Theta }\boldsymbol{\alpha }^{k}) \\ &=: \boldsymbol{A}\boldsymbol{\alpha }^{k}, \end{aligned} $$
(53)

where \(\boldsymbol{A} = -\widetilde{(\boldsymbol{\Theta }\boldsymbol{\omega }_{\star })} + \tilde{ \boldsymbol{\omega }}_{\star }\boldsymbol{\Theta }\). Next, the matrix \(\hat{\boldsymbol{\Theta }} = \boldsymbol{\Theta }+\frac{\Delta t}{2}\boldsymbol{G}( \boldsymbol{\omega }_{\star })\) in (52) has the inverse

$$ \hat{\boldsymbol{\Theta }}^{-1} = \frac{1}{\det \hat{\boldsymbol{\Theta }}} \begin{pmatrix} BC & \frac{\Delta t}{2}\Omega C(A+B) & 0 \\ -\frac{\Delta t}{2}\Omega C(A+B) & AC & 0 \\ 0 & 0 & \frac{\det \boldsymbol{B}}{C} \end{pmatrix} , $$
(54)

where the determinant of \(\hat{\boldsymbol{\Theta }}\) is given by \(\det \hat{\boldsymbol{\Theta }} = ABC +\frac{\Delta t^{2}}{4}C(A+B)^{2}\). Hence, the linearized perturbation dynamics can be given in the explicit form

$$ \boldsymbol{\alpha }^{k+1} = \left (\boldsymbol{I}-\Delta t\hat{\boldsymbol{\Theta }}^{-1} \boldsymbol{A}\right )\boldsymbol{\alpha }^{k} = \boldsymbol{D}\boldsymbol{\alpha }^{k}, $$
(55)

in which the matrix \(\boldsymbol{D}\) has the nonzero components

$$ \begin{aligned} D_{11} &= 1-\frac{\Delta t^{2}}{2\det \hat{\boldsymbol{\Theta }}}\Omega ^{2}C(A+B)(A-C), & D_{12} &= -\frac{\Delta t}{\det \hat{\boldsymbol{\Theta }}}\Omega AC(A-C), \\ D_{21} &= -\frac{\Delta t}{\det \hat{\boldsymbol{\Theta }}}\Omega BC(C-B), & D_{22}, &= 1+\frac{\Delta t^{2}}{2\det \hat{\boldsymbol{\Theta }}}\Omega ^{2}C(A+B)(C-B), \\ D_{33} &= 1. & \end{aligned} $$
(56)

The stability of rotation around the major principal axis is determined through the eigenvalues of \(\boldsymbol{D}\), being the roots of the characteristic polynomial

$$ \begin{aligned} \det (\lambda \boldsymbol{I}-\boldsymbol{D}) &= (\lambda -1)\left [(\lambda -D_{11})( \lambda -D_{22})-D_{12}D_{21})\right ] \\ &= (\lambda -1)\left [\lambda ^{2}-(D_{11}+D_{22})\lambda +(D_{11}D_{22}-D_{12}D_{21}) \right ] \end{aligned} $$
(57)

given by

$$ \lambda _{1/2} = \frac{b\pm \sqrt{b^{2}-4c}}{2}\quad \text{and}\quad \lambda _{3} = 1. $$
(58)

Herein, the parameters \(b = D_{11}+D_{22}\) and \(c = D_{11}D_{22}-D_{12}D_{21}\) can be calculated as

$$ \begin{aligned} b &= 2- \frac{\Delta t^{2}\Omega ^{2}}{2\det \hat{\boldsymbol{\Theta }}}C(A+B)(A+B-2C), \\ c &= 1+\frac{\Delta t^{2}\Omega ^{2}}{2\det \hat{\boldsymbol{\Theta }}}C \left ( \left (A+B\right )\left (A+B-2C\right )-2\left (A-C\right ) \left (C-B\right )\right ). \end{aligned} $$

The eigenvalues \(\lambda _{1/2}\) in (58) take the form

$$ \lambda _{1/2} = 1-\frac{d}{2}\pm \frac{\sqrt{d^{2}+4e}}{2} $$
(59)

with

$$ d = \frac{\Delta t^{2}}{2\det \hat{\boldsymbol{\Theta }}}\Omega ^{2}C(A+B)(A+B-2C), \qquad e = \frac{\Delta t^{2}}{\det \hat{\boldsymbol{\Theta }}}\Omega ^{2}C(A-C)(C-B). $$
(60)

For rotation around the major principal axis, it holds that \(C\geq \max (A,B)\), from which follows that \(d<0\). Therefore, at least one of the eigenvalues has a magnitude larger than unity which proves instability of rotation around the major principal axis, contrary to the intermediate axis theorem. Hence, the explicit scheme cannot correctly represent the stability properties of a freely rotating body.

4.2 A stability preserving implicit scheme

Here, we propose an alternative scheme for the rotational motion of a free body which preserves the stability properties of the principal axes of rotation in accordance with the extended intermediate axis theorem.

As a side remark, one needs to keep in mind that the aim of the paper is the improved prediction of rockfall. There exists of course a wealth of general purpose schemes which are of higher order, e.g. (explicit or implicit) Runge–Kutta schemes, BDF schemes, etc. However, any scheme that is not energy–momentum preserving will perform poorly on the dynamics of a freely rotating rigid body as it cannot respect the stability properties. Moreover, for rockfall simulation a scheme is needed which can deal with a singularity-free description of orientation (e.g., using unit quaternions). Of course, as explained in the introduction, there already exist energy–momentum preserving schemes for rigid body motion (e.g., ALGO-C1), having much better properties than the scheme presented here. However, for rockfall simulation a few constraints need to be met:

  1. 1.

    The scheme needs to be energy–momentum preserving;

  2. 2.

    It needs to use a singularity-free description of the orientation;

  3. 3.

    It needs to fit within the Moreau timestepping scheme, together with specialized dissipation and contact laws.

The last requirement rules out the use of existing dedicated schemes for rigid body motion as for these schemes there is a close mutual coupling between the angular velocity and orientation dynamics. Moreover, the use of higher-order approximations within the Moreau timestepping scheme makes little sense as, due to the fact that it has to deal with collisions, the Moreau scheme is intrinsically of low order. We therefore present a very simple alternative scheme in this section and will show in the remainder of the paper that this leads to an improved rockfall simulation code.

The alternative scheme consists of two consecutive parts:

  1. 1.

    As update rule for the angular velocity, we use the implicit midpoint scheme proposed by [24];

  2. 2.

    As update rule for the orientation parametrization, we propose a novel scheme which preserves the spin.

The implicit midpoint scheme for the angular velocity calculates \(\boldsymbol{\omega }^{k+1}\) at time instant \(t^{k+1} = t^{k}+\Delta t\) during a flight phase as

$$ \boldsymbol{\Theta }(\boldsymbol{\omega }^{k+1}-\boldsymbol{\omega }^{k})+\Delta t \boldsymbol{\omega }^{k+\frac{1}{2}}\times (\boldsymbol{\Theta }\boldsymbol{\omega }^{k+ \frac{1}{2}}) = 0\quad \mathrm{with}\quad \boldsymbol{\omega }^{k+ \frac{1}{2}} = \frac{1}{2}(\boldsymbol{\omega }^{k}+\boldsymbol{\omega }^{k+1}). $$
(61)

Substitution of the scheme (61) in the kinetic energy expression (48) yields

$$ T(\boldsymbol{\omega }^{k+1})-T(\boldsymbol{\omega }^{k}) = -\frac{\Delta t}{2}( \boldsymbol{\omega }^{k+\frac{1}{2}})^{\mathrm{T}}\left [\boldsymbol{\omega }^{k+ \frac{1}{2}}\times (\boldsymbol{\Theta }\boldsymbol{\omega }^{k+\frac{1}{2}}) \right ] = 0, $$
(62)

showing that the kinetic energy is preserved by the implicit scheme. The magnitude of the spin \(\|\vec{\boldsymbol{N}}_{S}\| = \|{}_{K}\vec{\boldsymbol{N}}_{S}\| = \| \boldsymbol{\Theta }\boldsymbol{\omega }\|\), defined by (15), is only dependent on the angular velocity. The implicit scheme for the angular velocity also conserves the magnitude of the spin as follows from

$$ \begin{aligned} ||\vec{\boldsymbol{N}}_{S}^{k+1}||^{2}-||\vec{\boldsymbol{N}}_{S}^{k}||^{2}&\stackrel{\phantom{\text{(61)}}}{=}( \boldsymbol{\Theta }\boldsymbol{\omega }^{k+1})^{\mathrm{T}}\boldsymbol{\Theta } \boldsymbol{\omega }^{k+1}-(\boldsymbol{\Theta }\boldsymbol{\omega }^{k})^{\mathrm{T}} \boldsymbol{\Theta }\boldsymbol{\omega }^{k} \\ &\stackrel{\phantom{\text{(61)}}}{=} 2(\boldsymbol{\Theta } \boldsymbol{\omega }^{k+\frac{1}{2}})^{\mathrm{T}}\boldsymbol{\Theta }(\boldsymbol{\omega }^{k+1}- \boldsymbol{\omega }^{k}) \\ &\stackrel{\text{(61)}}{=}-2\Delta t(\boldsymbol{\Theta } \boldsymbol{\omega }^{k+\frac{1}{2}})^{\mathrm{T}}\left [\boldsymbol{\omega }^{k+ \frac{1}{2}}\times (\boldsymbol{\Theta }\boldsymbol{\omega }^{k+\frac{1}{2}}) \right ] = 0. \end{aligned} $$
(63)

We now propose an update rule for the orientation parametrization. The update of the orientation \(\boldsymbol{A}_{IK}^{k+1}\) is chosen such that the spin remains constant, i.e.,

$$ {}_{I}\vec{\boldsymbol{N}}_{S}^{k+1} = {}_{I}\vec{\boldsymbol{N}}_{S}^{k}, $$
(64)

and such that the kinematic equation is correctly approximated in the sense that

$$ \lim _{\Delta t\downarrow 0} \frac{\boldsymbol{A}_{IK}^{k+1} - \boldsymbol{A}_{IK}^{k}}{\Delta t} = \dot{\boldsymbol{A}}_{IK}(t^{k}), $$
(65)

where \(\dot{\boldsymbol{A}}_{IK} = \boldsymbol{A}_{IK}(t)\tilde{\boldsymbol{\omega }}\). We choose an update of the form

$$ \boldsymbol{A}_{IK}^{k+1} = \boldsymbol{A}_{IK}^{k}\mathrm{e}^{ \tilde{\boldsymbol{\omega }}^{k+\frac{1}{2}}\Delta t} \boldsymbol{B}(\boldsymbol{\omega }^{k}, \boldsymbol{\omega }^{k+1}), $$
(66)

where the matrix \(\boldsymbol{B}\) depends on \(\boldsymbol{\omega }^{k}\) and \(\boldsymbol{\omega }^{k+1}\) and needs to fulfill \(\boldsymbol{B}^{\mathrm{T}}\boldsymbol{B} = \boldsymbol{I}\) to ensure that \((\boldsymbol{A}_{IK}^{k+1})^{\mathrm{T}}\boldsymbol{A}_{IK}^{k+1} = \boldsymbol{I}\). Furthermore, we demand that \(\boldsymbol{B}(\boldsymbol{x},\boldsymbol{x}) = \boldsymbol{I}\) for all \(\boldsymbol{x}\) such that

$$ \lim _{\Delta t\downarrow 0} \boldsymbol{B}(\boldsymbol{\omega }^{k},\boldsymbol{\omega }^{k+1}) = \boldsymbol{B}(\boldsymbol{\omega }^{k},\boldsymbol{\omega }^{k}) = \boldsymbol{I}, $$
(67)

from which follows the kinematic consistency (65).

To ensure the conservation of spin (64), we demand

$$ \boldsymbol{A}_{IK}^{k+1}\boldsymbol{\Theta }\boldsymbol{\omega }^{k+1} = \boldsymbol{A}_{IK}^{k} \boldsymbol{\Theta }\boldsymbol{\omega }^{k} $$
(68)

which, after substitution of the update rule, gives

$$ \boldsymbol{A}_{IK}^{k}\mathrm{e}^{\tilde{\boldsymbol{\omega }}^{k+\frac{1}{2}} \Delta t} \boldsymbol{B}(\boldsymbol{\omega }^{k},\boldsymbol{\omega }^{k+1}) \boldsymbol{\Theta }\boldsymbol{\omega }^{k+1} = \boldsymbol{A}_{IK}^{k}\boldsymbol{\Theta } \boldsymbol{\omega }^{k}. $$
(69)

The matrix \(\boldsymbol{B}\) will now be chosen such that

$$ \mathrm{e}^{\tilde{\boldsymbol{\omega }}^{k+\frac{1}{2}}\Delta t} \boldsymbol{B}( \boldsymbol{\omega }^{k},\boldsymbol{\omega }^{k+1}) \boldsymbol{\Theta }\boldsymbol{\omega }^{k+1} = \boldsymbol{\Theta }\boldsymbol{\omega }^{k} $$
(70)

and therefore

$$ \boldsymbol{B}(\boldsymbol{\omega }^{k},\boldsymbol{\omega }^{k+1}) \boldsymbol{\Theta } \boldsymbol{\omega }^{k+1} = \mathrm{e}^{-\tilde{\boldsymbol{\omega }}^{k+ \frac{1}{2}}\Delta t} \boldsymbol{\Theta }\boldsymbol{\omega }^{k} ,$$
(71)

from which we see that it is indeed only dependent on \(\boldsymbol{\omega }^{k}\) and \(\boldsymbol{\omega }^{k+1}\). The matrix \(\boldsymbol{B}\) is a pure rotation as it is orthogonal, i.e., \(\boldsymbol{B}^{\mathrm{T}} \boldsymbol{B} = \boldsymbol{I}\), and has a determinant \(+1\).

We recall that, by using Rodrigues formula, every pure rotation \(\boldsymbol{R}\) around the unit vector \(\boldsymbol{k}\) with rotation angle \(\theta \) can be represented as

$$ \boldsymbol{R} = \boldsymbol{I} + \sin \theta \tilde{\boldsymbol{k}} + (1-\cos \theta ) \tilde{\boldsymbol{k}}^{2}. $$
(72)

If \(\boldsymbol{k}\) and \(\theta \) are chosen as

$$ \boldsymbol{k}= \frac{\boldsymbol{b}\times \boldsymbol{c}}{\|\boldsymbol{b}\times \boldsymbol{c}\|}, \quad \cos \theta = \frac{\boldsymbol{b}\cdot \boldsymbol{c}}{\|\boldsymbol{b}\| \|\boldsymbol{c}\|} $$
(73)

then it holds that \(\boldsymbol{R}\boldsymbol{b}= \boldsymbol{c}\) if \(\|\boldsymbol{b}\| =\|\boldsymbol{c}\|\). Furthermore, it holds that \(\boldsymbol{R} \rightarrow \boldsymbol{I}\) for \(\boldsymbol{b}\rightarrow \boldsymbol{c}\).

Hence, we choose

$$ \boldsymbol{B}(\boldsymbol{\omega }^{k},\boldsymbol{\omega }^{k+1}) = \boldsymbol{I} + \sin \theta \tilde{\boldsymbol{k}} + (1-\cos \theta )\tilde{\boldsymbol{k}}^{2} $$
(74)

with

$$ \boldsymbol{k}= \frac{\boldsymbol{\Theta }\boldsymbol{\omega }^{k+1}\times (\mathrm{e}^{-\tilde{\boldsymbol{\omega }}^{k+\frac{1}{2}}\Delta t} \boldsymbol{\Theta }\boldsymbol{\omega }^{k}) }{\|\boldsymbol{\Theta }\boldsymbol{\omega }^{k+1}\times (\mathrm{e}^{-\tilde{\boldsymbol{\omega }}^{k+\frac{1}{2}}\Delta t} \boldsymbol{\Theta }\boldsymbol{\omega }^{k}) \|}, \quad \cos \theta = \frac{\boldsymbol{\Theta }\boldsymbol{\omega }^{k+1}\cdot (\mathrm{e}^{-\tilde{\boldsymbol{\omega }}^{k+\frac{1}{2}}\Delta t} \boldsymbol{\Theta }\boldsymbol{\omega }^{k}) }{\|\boldsymbol{\Theta }\boldsymbol{\omega }^{k+1}\| \|\boldsymbol{\Theta }\boldsymbol{\omega }^{k}\|}. $$
(75)

To prove that the proposed scheme has the desired stability properties of the principal axes of inertia, we use the Lyapunov functions \(V(\boldsymbol{y}^{k})\) for the major and minor principal axes as presented before. As the proposed scheme conserves the kinetic energy and the spin by construction it holds that \(V(\boldsymbol{y}^{k+1}) = V(\boldsymbol{y}^{k})\), whereas positive definiteness has already been shown. The proposed scheme therefore preserves the stability properties of the principal axes of inertia in accordance with the extended intermediate axis theorem.

4.3 Numerical comparison

We compare the explicit and the implicit scheme on a numerical example. We consider a cuboid of mass \(m = 1\) kg with length \(a=3\) m, width \(b=2\) m, and height \(c=1\) m in the \(\vec{\mathbf{e}}_{x}^{K}\), \(\vec{\mathbf{e}}_{y}^{K}\) and \(\vec{\mathbf{e}}_{z}^{K}\) direction, respectively. The principal moments of inertia are therefore \(A = \frac{m}{12}(b^{2}+c^{2})\), \(B = \frac{m}{12}(a^{2}+c^{2})\), and \(C = \frac{m}{12}(a^{2}+b^{2})\) such that \(A< B< C\). Rotation in the neighborhood of stationary rotation \(\boldsymbol{\omega }_{\star }= \Omega \boldsymbol{e}_{\star }\) around the major principal axis \(\boldsymbol{e}_{\star }= \boldsymbol{e}_{3}\) is considered. As initial conditions we choose \(\boldsymbol{A}_{IK}(0) = \boldsymbol{I}\) and \(\boldsymbol{\omega }= \begin{pmatrix} \omega _{x} & \omega _{y} & \omega _{z} \end{pmatrix} ^{\mathrm{T}}= \begin{pmatrix} 10^{-3} & 10^{-3} & 10 \end{pmatrix} ^{\mathrm{T}}\) rad/s. We simulate 20 s using a timestep of \(\Delta t = 0.01\) s using the explicit and implicit scheme. The results of both schemes are shown in Figs. 2 and 3. The body initially rotates in the vicinity of the major principal axis with angular velocity \(\omega _{z} = \Omega = 10\) rad/s, which is stable as follows from the extended intermediate axis theorem. However, in the numerical solution of the explicit scheme, the body deviates from stationary rotation around the major principal axis (approximately at \(t=10\) s) and tends to stable rotation around the minor principal axis with angular speed \(\omega _{x} = -\sqrt{\frac{C}{A}}\Omega \). If a smaller timestep is taken in the explicit scheme, then the change of axis will be slower and will take place at a later point in time. This also follows from the fact that the eigenvalues \(\lambda _{1,2}\) in (59) tend to 1 if the timestep is reduced. The solution of the implicit scheme remains very close to the major principal axis, both in angular velocity and in the difference vector \(\boldsymbol{d}(t) = {}_{K}\vec{\mathbf{e}}_{z}^{K\star } - {}_{K} \vec{\mathbf{e}}_{z}^{K}\) and is therefore much more accurate.

Fig. 2
figure 2

Angular velocities for initial rotation around the major principal axis. The explicit scheme does not correctly represent the stability of stationary rotation around the major principal axis and converges to rotation around the minor principal axis. The implicit scheme correctly reflects the stability of the major principal axis. (Color figure online)

Full size image
Fig. 3
figure 3

Difference vector \(\boldsymbol{d}(t) = {}_{K}\vec{\mathbf{e}}_{z}^{K\star } - {}_{K} \vec{\mathbf{e}}_{z}^{K}\) for initial rotation around the major principal axis, showing that the explicit scheme cannot represent the stability of the major principal axis. The orientation of the body changes drastically. The implicit scheme is able to keep the orientation difference small. (Color figure online)

Full size image

In Fig. 4, the results of both schemes are compared for an initial rotation around the unstable intermediate axis principal axis. An initial angular velocity of \(\boldsymbol{\omega }= \begin{pmatrix} \omega _{x} & \omega _{y} & \omega _{z} \end{pmatrix} ^{\mathrm{T}}= \begin{pmatrix} 10^{-3} & 10 & 10^{-3} \end{pmatrix} ^{\mathrm{T}}\) rad/s is chosen. For the explicit scheme, the rotation around the \({}_{K}\vec{\mathbf{e}}_{y}^{K}\)-axis is unstable, and the motion asymptotically tends towards rotation around the minor principal axis. The implicit scheme correctly predicts that the motion will move on the heteroclinic orbit and will keep on swapping back and forth from rotation around the \({}_{K}\vec{\mathbf{e}}_{y}^{K}\)-axis in positive direction to rotation around this axis in negative direction.

Fig. 4
figure 4

Angular velocities for initial rotation around the unstable intermediate principal axis. The explicit scheme converges to stationary rotation around the minor principal axis, whereas the implicit scheme reveals the back-and-forth swapping of rotation around the intermediate axis in positive and negative direction. (Color figure online)

Full size image

5 Implementation in a rockfall simulation code

The translational and rotational equations of motion and impact equations for rockfall simulation are described in detail in [16], together with the modeling of rock shape, contact interactions, dissipation mechanisms, and terrain. The explicit and implicit scheme are discussed in Sect. 4 for the rotational motion of a free body. Here we will explain how the implicit scheme can be implemented in RAMMS::ROCKFALL taking into account dissipation and contact interaction. First, we will briefly explain the current implementation of the explicit scheme.

The motion of the rock is described through a set of generalized coordinates

$$ \boldsymbol{q}= \begin{pmatrix} \boldsymbol{r} \\ \boldsymbol{p}\end{pmatrix} , $$
(76)

where \(\boldsymbol{r}: = {}_{I}\vec{\boldsymbol{r}}_{OS}\) is the position of the center of mass and \(\boldsymbol{p}\) is the quaternion which encodes the transformation matrix \(\boldsymbol{A}_{IK}\). The equations of motion are expressed in the generalized velocities

$$ \boldsymbol{u}= \begin{pmatrix} \boldsymbol{v} \\ \boldsymbol{\omega }\end{pmatrix} , $$
(77)

where \(\boldsymbol{v}:= {}_{I} \vec{\boldsymbol{v}}_{S} = {}_{I}\dot{\vec{\boldsymbol{r}}}_{OS}\) is the velocity of the center of mass, i.e., the translational motion is expressed in the inertial frame and the rotational motion in the body-fixed frame.

Various kinds of dissipation (e.g., viscous drag in layer close to the terrain) and gravity are taken into account through an external force \(\boldsymbol{f}_{\textrm{ext}}(\boldsymbol{q},\boldsymbol{u}) := {}_{I}\vec{\boldsymbol{F}}\) and torque \(\boldsymbol{M}_{\textrm{ext}}(\boldsymbol{q},\boldsymbol{u}) := {}_{K}\vec{\boldsymbol{M}}_{S}\) with respect to the center of mass \(S\). The velocity dependence of the external forces and torques stems from drag terms which are assumed to be linear [16]

$$ \begin{aligned} \boldsymbol{f}_{\textrm{ext}}(\boldsymbol{q},\boldsymbol{u}) &= - c_{\textrm{drag}} \boldsymbol{v}+ \boldsymbol{f}_{\textrm{grav}}, \\ \boldsymbol{M}_{\textrm{ext}}(\boldsymbol{q},\boldsymbol{u}) &= - \boldsymbol{C}_{\textrm{drag}} \boldsymbol{\omega }, \end{aligned} $$
(78)

where the drag coefficients \(c_{\textrm{drag}}\) and \(\boldsymbol{C}_{\textrm{drag}}\) are only active in a drag layer above the terrain and therefore depend on the position \(\boldsymbol{r}\) of the rock. Contact interactions with the terrain lead to nonimpulsive contact forces and impulsive contact forces which are governed by frictional contact and impact laws. Velocity–impulse-based integration techniques, such as Moreau’s timestepping scheme, are a very robust way to deal with multibody systems with frictional unilateral constraints, forming an active field of research within the Nonsmooth Dynamics community. The implementation in RAMMS::ROCKFALL may be seen as a variant of the classical Moreau timestepping scheme. Hence, the contact forces and impulsive forces are gathered in a discrete percussion \(\boldsymbol{P}\) which enters into the translational and rotational part of the equations of motion through the generalized force directions \(\boldsymbol{W}_{v}(\boldsymbol{q})\) and \(\boldsymbol{W}_{\omega }(\boldsymbol{q})\). Specially tailored contact laws are used for the discrete percussion \(\boldsymbol{P}\) to correctly describe the physics of rock–terrain interaction. In particular, a dynamic slippage dependent friction law is employed to describe the scarring phenomenon, i.e., upon collision with the terrain the rock tends to slide along the slope pushing soil in front of it which increases the friction and finally lets it tumble over. These highly specialized contact laws have been developed and extensively reported in [16].

5.1 Implementation of the explicit scheme

In the style of Moreau’s timestepping scheme, the explicit scheme first calculates a midpoint

$$ \begin{aligned} \boldsymbol{r}^{k+\frac{1}{2}} &= \boldsymbol{r}^{k} + \frac{\Delta t}{2}\boldsymbol{v}^{k}, \\ \boldsymbol{p}^{k+\frac{1}{2}} &= \frac{\boldsymbol{p}_{\textrm{pre}}^{k+\frac{1}{2}}}{\|\boldsymbol{p}_{\textrm{pre}}^{k+\frac{1}{2}}\|}, \qquad \boldsymbol{p}_{\textrm{pre}}^{k+\frac{1}{2}} = \boldsymbol{p}^{k} + \frac{\Delta t}{2} \boldsymbol{F}(\boldsymbol{p}^{k})\boldsymbol{\omega }^{k}, \end{aligned} $$
(79)

which is used to approximate the generalized force directions \(\boldsymbol{W}_{v}^{k+\frac{1}{2}} = \boldsymbol{W}_{v}(\boldsymbol{q}^{k+\frac{1}{2}})\) and \(\boldsymbol{W}_{\omega }^{k+\frac{1}{2}} = \boldsymbol{W}_{\omega }(\boldsymbol{q}^{k+ \frac{1}{2}})\).

Subsequently, the generalized velocities \(\boldsymbol{u}^{k+1}\) are obtained by solving the balance laws on velocity–impulse level

$$ \begin{aligned} & m(\boldsymbol{v}^{k+1}-\boldsymbol{v}^{k}) = \boldsymbol{f}_{\textrm{ext}}^{k+ \frac{1}{2}} \Delta t + \boldsymbol{W}_{v}^{k+\frac{1}{2}}\boldsymbol{P}^{k+1}\;, \\ &\boldsymbol{\Theta }(\boldsymbol{\omega }^{k+1}-\boldsymbol{\omega }^{k})+ \frac{\Delta t}{2}\boldsymbol{G}(\boldsymbol{\omega }^{k})(\boldsymbol{\omega }^{k+1}+ \boldsymbol{\omega }^{k}) = \boldsymbol{M}_{\textrm{ext}}^{k+\frac{1}{2}}\Delta t+ \boldsymbol{W}_{\omega }^{k+\frac{1}{2}}\boldsymbol{P}^{k+1}\;, \end{aligned} $$
(80)

together with set-valued combined contact–impact laws for the discrete percussions \(\boldsymbol{P}^{k+1}\), which implicitly depend on \(\boldsymbol{u}^{k+1}\). Also the external forces \(\boldsymbol{f}_{\textrm{ext}}^{k+\frac{1}{2}} = \boldsymbol{f}_{\textrm{ext}}( \boldsymbol{q}^{k+\frac{1}{2}},\boldsymbol{u}^{k+\frac{1}{2}})\) and external torques \(\boldsymbol{M}_{\textrm{ext}}^{k+\frac{1}{2}} = \boldsymbol{M}_{\textrm{ext}}( \boldsymbol{q}^{k+\frac{1}{2}},\boldsymbol{u}^{k+\frac{1}{2}})\) depend through \(\boldsymbol{u}^{k+\frac{1}{2}} = \frac{1}{2}( \boldsymbol{u}^{k} + \boldsymbol{u}^{k+1}) \) on \(\boldsymbol{u}^{k+1}\). The stepping equations (80), together with the contact–impact laws, form a set of implicit equations for the unknowns \((\boldsymbol{u}^{k+1},\boldsymbol{P}^{k+1})\). The linearity of the drag terms in the external forces and torques, together with the explicit character of scheme, allows us to express \(\boldsymbol{u}^{k+1}\) directly in the unknown percussion \(\boldsymbol{P}^{k+1}\), i.e.,

$$ \boldsymbol{v}^{k+1} = \left ( m + \frac{\Delta t}{2}c_{\textrm{drag}}\right )^{-1} \left ((m-\frac{\Delta t}{2} c_{\textrm{drag}})\boldsymbol{v}^{k} + \boldsymbol{f}_{\textrm{grav}}\Delta t+ \boldsymbol{W}_{v}^{k+\frac{1}{2}}\boldsymbol{P}^{k+1}\right ), $$
(81)
$$ \boldsymbol{\omega }^{k+1} = \left ( \boldsymbol{\Theta } + \frac{\Delta t}{2} \left (\boldsymbol{G}(\boldsymbol{\omega }^{k}) + \boldsymbol{C}_{\textrm{drag}}\right ) \right )^{-1} \left (\left (\boldsymbol{\Theta } - \frac{\Delta t}{2}( \boldsymbol{G}(\boldsymbol{\omega }^{k}) + \boldsymbol{C}_{\textrm{drag}})\right ) \boldsymbol{\omega }^{k} + \boldsymbol{W}_{\omega }^{k+\frac{1}{2}}\boldsymbol{P}^{k+1} \right ), $$
(82)

which can be used to reduce the set of equations to the unknowns \(\boldsymbol{P}^{k+1}\). In fact, the addition of the contact percussions to the explicit scheme renders it a semiimplicit scheme. Finally, the positions are updated through

$$ \begin{aligned} \boldsymbol{r}^{k+1} &= \boldsymbol{r}^{k+\frac{1}{2}} + \frac{\Delta t}{2} \boldsymbol{v}^{k+1}, \\ \boldsymbol{p}^{k+1} &= \frac{\boldsymbol{p}_{\textrm{pre}}^{k+1}}{\|\boldsymbol{p}_{\textrm{pre}}^{k+1}\|}, \qquad \boldsymbol{p}_{\textrm{pre}}^{k+1} = \boldsymbol{p}^{k+\frac{1}{2}} + \frac{\Delta t}{2} \boldsymbol{F}(\boldsymbol{p}^{k+\frac{1}{2}})\boldsymbol{\omega }^{k+1}. \end{aligned} $$
(83)

5.2 Implementation of the implicit scheme

We now discuss how to embed the implicit scheme for rotational motion within the Moreau-type timestepping method in a way which is as little invasive as possible by using a two-step method.

First, similar to the explicit method a midpoint is calculated as in (79). In a first implicit step, a free angular velocity \(\boldsymbol{\omega }_{\textrm{free}}^{k+1}\) is calculated through the implicit equation

$$ \boldsymbol{\Theta }(\boldsymbol{\omega }_{\textrm{free}}^{k+1}-\boldsymbol{\omega }^{k})+ \Delta t\boldsymbol{\omega }_{\textrm{free}}^{k+\frac{1}{2}}\times ( \boldsymbol{\Theta }\boldsymbol{\omega }_{\textrm{free}}^{k+\frac{1}{2}}) = 0\quad \text{with}\quad \boldsymbol{\omega }_{\textrm{free}}^{k+\frac{1}{2}} = \frac{1}{2}(\boldsymbol{\omega }^{k}+\boldsymbol{\omega }_{\textrm{free}}^{k+1}), $$
(84)

which does not account for the external torque and contact interaction. In a second implicit step, the final generalized velocity \(\boldsymbol{u}^{k+1}\) are obtained by solving the implicit equations

$$ \begin{aligned} m(\boldsymbol{v}^{k+1}-\boldsymbol{v}^{k}) &= \boldsymbol{f}_{\textrm{ext}}^{k+ \frac{1}{2}} \Delta t + \boldsymbol{W}_{v}^{k+\frac{1}{2}}\boldsymbol{P}^{k+1}, \\ \boldsymbol{\Theta }(\boldsymbol{\omega }^{k+1} - \boldsymbol{\omega }_{\textrm{free}}^{k+1}) &= \boldsymbol{M}_{\textrm{ext}}^{k+\frac{1}{2}}\Delta t+ \boldsymbol{W}_{\omega }^{k+ \frac{1}{2}}\boldsymbol{P}^{k+1}, \end{aligned} $$
(85)

together with the contact–impact laws. In the case of linear drag, similar to (81) and (82), equations (85) can be solved for \(\boldsymbol{u}^{k+1}\), allowing us to reduce the equations to the unknowns \(\boldsymbol{P}^{k+1}\). The final velocity \(\boldsymbol{v}^{k+1}\) is used to obtain the position update

$$ \boldsymbol{r}^{k+1} = \boldsymbol{r}^{k+\frac{1}{2}} + \frac{\Delta t}{2} \boldsymbol{v}^{k+1}, $$
(86)

whereas we cannot use \(\boldsymbol{\omega }^{k+1}\) for the orientation update (66) because the spin is not conserved in the second implicit step. Instead, the free velocity \(\boldsymbol{\omega }_{\textrm{free}}^{k+1}\) is used for the orientation update

$$ \boldsymbol{A}_{IK}^{k+1} = \boldsymbol{A}_{IK}^{k}\mathrm{e}^{ \tilde{\boldsymbol{\omega }}_{\textrm{free}}^{k+\frac{1}{2}}\Delta t} \boldsymbol{B}( \boldsymbol{\omega }^{k},\boldsymbol{\omega }_{\textrm{free}}^{k+1}), $$
(87)

where \(\boldsymbol{A}_{IK}^{k}\) is calculated from \(\boldsymbol{p}^{k}\) using (10) and \(\boldsymbol{B}\) follows from (74). Finally, \(\boldsymbol{p}^{k+1}\) is retrieved from \(\boldsymbol{A}_{IK}^{k+1}\) using, e.g., Shepperd’s method [25].

6 Simulation results with RAMMS::ROCKFALL

A realistic modeling of platy shaped rocks, which descend downslope in a wheel-like motion, is of paramount importance in rockfall simulations. The following example simulations show the implications of the implicit integration scheme for accurate simulation results, especially for rocks that rotate around their major principal axis. Comparisons between the newly proposed implicit scheme versus the old explicit scheme are visualized for RAMMS::ROCKFALL simulations. All simulations were performed with a so-called EOTA221-platy shaped rock (see Table 1) which was also used in rockfall experiments with instrumented rocks at Chant Sura in Davos, Switzerland [4]. The rock mass is 780 kg, homogeneously distributed in a rigid body with dimensions of 0.93, 0.93, and 0.47 m along the three principal axes (\(x\)-, \(y\)-, and \(z\)-axis of the \(K\)-frame).

Table 1 Comparison of the computational time for the explicit and the implicit schemes (1000-rock simulation)
Full size table

6.1 Idealized inclined slope

The first group of simulations was performed on an idealized digital elevation model (DEM) with an inclined slope of 40∘. In total 1000 simulations were performed. The initial rock orientation was the only varying input parameter. The same soil conditions as well as the same time step (0.002 s) was employed for both the explicit and the implicit scheme [17]. Figure 5 visualizes the maximal kinetic rock energy of the simulation ensemble of the thousand trajectories using the statistic mode of RAMMS::ROCKFALL for the explicit scheme and for the new implicit scheme (Figs. 5a and 5b, respectively). The maximal spatial spread of the explicit scheme is delineated with a red line for comparison purposes (Fig. 5b). The implicit scheme leads to a narrower distribution in the inclined section but longer and more curved run-out behavior in the flat run-out area in comparison to the explicit scheme. Figures 5c and 5d show the altitude of the center of mass of the rock as a function of the projected distance, being the arclength of the horizontal projection of the trajectory. Furthermore, Figs. 5c and 5d compare the development of rock rotations along its principal axes for a single simulation with the same initial orientations. In the explicit scheme, the stability of the major principal axis (\(z\)-axis in blue) is not well pronounced in the beginning, and is completely lost halfway down the slope. The implicit scheme on the other hand, stabilizes very fast along the major principal axis, and remains in this configuration after reaching the horizontal run-out zone. The longer run-out in the implicitscheme is due to the stable rotation along the major principal axis (wheel-motion). The new integrator is thus able to simulate the effective lateral spreading behavior during the run-out section with larger accuracy whereas the explicit integrator is producing a broad trajectory area coverage induced by the erratic tumbling motion. Table 2 compares the corresponding mean and standard deviation values obtained for both numerical schemes. Rotations and kinetic energy values are larger for the implicit scheme, whereas the jump heights are smaller and the velocities are only slightly smaller. The stable motion along the principal axis in the implicit scheme favors higher rotations and lower jump heights.

Fig. 5
figure 5

EOTA221 rock rolling down an idealized inclined slope: (a) the explicit scheme, (b) the implicit scheme (1000 simulations). Figures (c) and (d) show angular velocities of all three axes (\(x\)-axis = green, \(y\)-axis = red, \(z\)-axis = blue) plotted during the descent for an exemplary trajectory with identical starting orientation (brown = terrain altitude, black = COG altitude). A clear stabilization around the major principal axis (blue) is visible for the implicit scheme (d), whereas erratic tumbling motion prevails for the major part of the descent in the explicit scheme (c). (Color figure online)

Full size image
Table 2 Comparison of key data values of 1000 rocks on an idealized slope for the explicit and the implicit schemes
Full size table

6.2 Chant Sura

The second group of simulations was performed using the DEM of the Chant Sura rockfall test site located near Davos, Switzerland [4]. Again, 1000 simulations were carried out, with identical soil parameters and the same time step (0.002 s) for both the explicit and the implicit schemes. The implicit scheme gives a similar extent in the acceleration zone for the rocks, see the first row in Fig. 6, whereas the lateral spreading and the run-out distance in the run-out zone are much longer when compared with the explicit scheme. The maximal spatial spread of the explicit scheme is indicated with a red line in Fig. 6b for comparison purposes. The implicit scheme captures the fast and stable rotation along the major principal axis (\(z\)-axis in blue) of the flat rock during the entire trajectory while the explicit scheme cannot reproduce the rock’s stable rotation around its major principal axis throughout the entire traveling path. Table 3 lists the mean and standard deviation values for both numerical schemes. Equivalent trends as for the inclined plane are distinguishable, although the differences in kinetic energy are more obvious; the differences for the jump heights are less pronounced. Again, the stable motion along the principal axis in the implicit scheme favors higher rotations and velocities (and thus kinetic energies) and lower jump heights.

Fig. 6
figure 6

EOTA221 rock rolling down an actual slope at Chant Sura (Davos, Switzerland): (left) the explicit scheme, (right) the implicit scheme. For the comparison the same starting orientations were used for all rocks. Figures (c) and (d) show angular velocities of all three axes (\(x\)-axis = green, \(y\)-axis = red, \(z\)-axis = blue) plotted during the descent for an exemplary trajectory with identical starting orientation (brown = terrain altitude, black = COG altitude). The implicit scheme leads to wheel-like motion and predicts a longer run-out distance. (Color figure online)

Full size image
Table 3 Comparison of key data values of 1000 rocks at Chant Sura for the explicit and the implicit schemes
Full size table

Finally, Table 1 displays the computational time recorded for the 1000-rock simulations for the explicit and the implicit schemes. The implicit scheme is slower by a factor of around 2 in comparison to the explicit scheme, which is expected as the updating of rock rotations requires more computational effort. However, the former does respect the stability properties of rotating 3D objects, which is a key improvement for RAMMS::ROCKFALL. It is anticipated that for extremely long rolling phases of platy rocks the differences between the explicit and the implicit schemes could be even larger.

7 Conclusions

In this paper an extended intermediate axis theorem has been presented as a 6-dimensional extension of the classical stability results for the Euler equations. The stability of the major and minor principal axes of the 6-dimensional system has been proven using the method of Lyapunov functions. Using the same Lyapunov functions, we have given a rigorous proof that the implicit scheme conserves the energy and spin and, thereby, respects the extended intermediate axis theorem. Here, we have chosen an implicit scheme which uses the implicit midpoint rule for the angular velocity dynamics. It is evident that all energy–momentum conserving schemes respect the extended intermediate axis theorem. Existing energy–momentum conserving schemes may have better accuracy than the scheme we proposed here, but the proposed scheme can be easily incorporated in existing rockfall simulation codes. Further research needs to be conducted to investigate how more sophisticated energy–momentum conserving schemes can be combined with the specialized dissipation and contact laws of rockfall simulation using Moreau-type timestepping methods.

The computational cost per time-step is larger for the implicit scheme than the explicit scheme, as Newton iterations are needed to solve the implicit equations. However, numerical simulations show that the implicit scheme is far more accurate as it respects both the energy conservation and the invariance of the spin.

Numerical simulations with RAMMS::ROCKFALL using the newly developed implicit scheme show that for the downward motion of a platy rock on an actual slope the rotation around the major principal axis is stable, even in the presence of intermediate collisions and contact phases with the slope as is also observed in field experiments. Furthermore, the simulations show that the proper simulation of wheel-like motion is essential for the correct statistical prediction of the run-out length. Future work will focus on the validation of simulation results against the real rockfall events and field experiments with instrumented rocks.