Trajectories of a ball moving inside a spherical cavity using first integrals of the governing nonlinear system

Abstract

Analytical study of ball vibration absorber behavior is presented in the paper. The dynamics of trajectories of a heavy ball moving without slipping inside a spherical cavity are analyzed. Following our previous work, where a similar system was investigated through various numerical simulations, research of the dynamic properties of a sphere moving in a spherical cavity was carried out by methods of analytical dynamics. The strategy of analytical investigation enabled definition of a set of special and limit cases which designate individual domains of regular trajectories. In order to avoid any mutual interaction between the domains along a particular trajectory movement, energy dissipation at the contact of the ball and the cavity has been ignored, as has any kinematic excitation due to cavity movement. A governing system was derived using the Lagrangian formalism and complemented by appropriate non-holonomic constraints of the Pfaff type. The three first integrals are defined, enabling the evaluation of trajectory types with respect to system parameters, the initial amount of total energy, the angular momentum of the ball and its initial spin velocity. The neighborhoods of the limit trajectories and their dynamic stability are assessed. Limit and transition special cases are investigated along with their individual elements. The analytical means of investigation enabled the performance of broad parametric studies. Good agreement was found when comparing the results achieved by the analytical procedures in this paper with those obtained by means of numerical simulations, as they followed from the Lagrangian approach and the Appell–Gibbs function presented in previous papers.

1 Introduction

The motivation for this work comes from the results of two recent papers by the authors which examined the behavior of a ball-type absorber using numerical simulations [36, 37]. In these engineering-oriented articles, a number of particular phenomena were identified, the causes of which had to be estimated from numerical results alone. However, these phenomena can have a major impact on the stability of the whole structure including the absorber. The presented analysis based on the first integrals, which is made possible by the assumption of a fixed cavity, enables a detailed analysis of the movement of the sphere inside the cavity. Therefore, this work does not intend to primarily examine the design details and effectiveness of the ball-type absorber. Instead, it focuses on the theoretical explanation and interpretation of phenomena observed during numerical and physical experiments, and on a clear separation of individual cases corresponding to the natural properties of the system.

The family of passive tuned mass vibration absorbers is well established in the engineering literature; see the exhaustive review paper [13] which reflects the situation until 2017. Conventional passive absorbers are of the pendulum type, mostly based on the auto-parametric principle; see, for instance, [16, 30, 31, 45]. The basic principle of ball-type absorbers comes from the rolling movement of a metallic ball of radius r inside a metallic spherical cavity of radius $R>r$ (Fig. 1), possibly with a rubber lining. Such a simple device is practically maintenance free; this is a topic which is gaining increasing popularity [10]. Its vertical dimension can be relatively small, and thus, it can be used in cases where a pendulum absorber is inapplicable due to lack of vertical space or difficult maintenance. Not surprisingly, this type of absorber is growing in popularity in connection with wind turbines, e.g., [8].

Fig. 1

Ball vibration absorber in a dynamic testing laboratory

The first papers dealing with the theory and practical aspects of ball absorbers were published during the past decades, see [41, 42], and are based on engineering approaches. However, to the best of the authors’ knowledge, the first papers that dealt with this issue from the perspective of rational (analytical) dynamics were published only some years ago; see [38, 39]. They presented a basic nonlinear model in 2D together with its numerical evaluation, and a report on its practical application, including some results of long-term in situ measurements. The 2D approach is satisfactory when the absorber is limited to acting in a specific direction, e.g., in bridges, systems of multiple one-directionally acting elements, etc. However, structures like masts or towers require equipment which works simultaneously in both horizontal directions, and hence, a 3D model must be considered.

Modeling of a homogeneous sphere rolling on a perfectly rough surface has a long tradition in classical mechanics. Several 3D approaches are available which respect the spatial and strongly nonlinear character of the system. The system is non-holonomic with linear constraints in the first derivatives with respect to time. The classical setting of several particular cases, including the one that this paper is concerned with, is considered by Routh [43]. He and other authors of popular monographs use the general Lagrangian methodology [4, 32, 40], which is by far the most popular approach in classical mechanics; it offers many advantages which are actively applied in practice and research, including the rolling sphere problems [17]. For example, a study regarding the rolling of a ball over a spherical surface based on homogeneous Lagrange equations has recently been published [18]. There the author analyzes the free and undamped movement of a ball in the vertical plane using phase portraits for different initial conditions together with generalizations regarding vibro-impact dynamics.

The approach based on the Appell–Gibbs function, even though not frequently used, has proven very effective in subsequent numerical simulations. For details and some specific attributes of this approach, see, e.g., papers [11, 14, 26, 28, 46, 48]. Recently, the rolling of a ball over a curved surface was dealt with on an abstract basis using the Lie group theoretical methods, e.g., [19]. Borisov et al. [7] use a similar abstract methodology and present a more complete analysis of the Routh solution for the solid of revolution. They derived new integrals for the ball rolling on non-symmetrical surfaces of the second order. Jurdjevic and Zimmerman [20] extended the problem to a hyperbolic analogue in which the spheres are replaced by the hyperboloids, and rolling is taken in an isometric sense in either Euclidean or Riemann geometry. A numerical analysis of the ball rolling in a spherical recess, also based on the Appell–Gibbs function, was studied numerically by Legeza [24].

Another derived topic with high publication activity regards the rolling of a so-called Chaplygin sphere—a dynamically non-symmetric non-homogeneous sphere—which represents one of the best known integrable systems of classical non-holonomic mechanics [5], either with or without considering a spin [22]. There is no doubt that devices based on similar effects find their use in absorbing unwanted vibrations. The usage of non-homogeneous spheres, hemispheres or semi-elliptic spheres would allow the absorber to be fine-tuned for a precisely limited nonlinear damping effect or multidirectional damping [27, 28]. These topics are already popular in the engineering literature, however, still only partially treated analytically. For example, a prospective tuned vibration absorber based on the nested ball principle has been vaguely described in a patent proposal [33]. The theoretical analysis of this setup is partially addressed in an abstract way by Borisov [6], and a case with a semi-elliptic cavity was described by Legeza [25]. It is also worth noting that the overview [13] does not mention any text concerning the dynamic stability analysis of structures equipped with vibration absorbers. This is despite the fact that the subject is mentioned among the topics that require considerable attention.

The above-mentioned abstract solutions offer valuable tools of investigation. From the engineering perspective, however, the actual trajectories that can be encountered in a particular vibration absorbing device are interesting because they determine its efficiency. The authors addressed two possible strategies for this purpose: (1) the Lagrangian formalism in 2D [38] and (2) the Appell–Gibbs function in 3D; see [36, 37]. In each case, a governing differential system was composed, and the solution itself was conducted numerically with subsequent analysis of the extensive data set. Important physical properties of the “ball—spherical cavity” system were discovered in these papers, and many particular trajectory types were identified numerically, however, without proper analytical explanation. This indicates that the problem should also be investigated at an analytical level using energy, momentum and angular momentum balance principles related to relevant first integrals. This strategy makes a qualitative investigation of the system behavior possible at least in the setting of the homogeneous problem. A similar technique was also adopted for investigating a nonlinear spherical pendulum; see [35]. In general, such an approach can significantly improve insight into the internal character of the system and increase the possibilities for practical applications. It enables a systematic identification of limit trajectories and the definition of relevant categories. The existence of limit trajectories depends on the system parameters and the initial conditions setting. The possibility of delimiting the domains of parameters can significantly contribute to an analysis of a system’s reliability and its lifetime period. Moreover, the obtained results can serve in constructing effective forms of the Lyapunov function intended for particular cases; see [9, 23] and other resources in the domain of dynamic stability.

Finally, it is worth mentioning that many of the algebraic manipulations required to derive or verify formulas in this paper were done using the Wolfram Mathematica Package [47].

The paper is organized as follows. First, after this introduction, the governing system and three of the first integrals based on the Lagrangian approach are derived. A particular characteristic function is then introduced as a basic tool for further classification of the response trajectories. Next, the particular “separation circle” trajectory is introduced, which separates two main trajectory groups. Sections 4 and 5 describe settings with and without consideration of the initial spin of the ball. A particular type of an almost planar rocking is then analyzed in detail in Sect. 6. Finally, the last section concludes.

2 Governing system and first integrals

2.1 Lagrangian system and first integrals

The rolling without slipping of a ball on a surface is a non-holonomic problem because constraints relating to the mutual movement of a ball and a surface include velocity components. When putting together an expression for the kinetic and potential energies T, V, and external forces $\mathbf{Q}=[Q_1,\ldots ,Q_n]$, the relevant Lagrangian equations should be written as follows:

$$\begin{aligned} \frac{{ \text{ d }}}{{ \text{ d }}t}\left( \frac{\partial T}{\partial \dot{q}_j}\right) - \frac{\partial T}{\partial q_j}+\frac{\partial V}{\partial q_j} = Q_j+ \sum \limits _{m=1}^l \lambda _m\cdot B_{mj}, \end{aligned}$$

(1)

$j=1,\ldots ,n$, and with non-holonomic constraints:

$$\begin{aligned} \sum \limits _{j=1}^n B_{mj}\cdot \dot{q}_j+ B_m = 0. \qquad m=1,\ldots ,l, \end{aligned}$$

(2)

where $q_j\;(j=1,\ldots ,n)$ are generalized coordinates. Symbols $B_{mj}$ and $B_m$ are generally functions of $q_j$. Explicit time can be usually omitted as constraints are considered scleronomous as a rule if external kinematic excitation is absent and only initial conditions provide energy to the system. Provided that kinematic excitation works, the respective parameters $B_m\ne 0$ and should then be considered as functions of time. Symbols $\lambda _m\;(m=1,\ldots ,l)$ are Lagrangian multipliers which are used to add non-holonomic conditions to the Hamiltonian functional. Therefore, the system Eqs. (1), (2) includes $n+l$ unknowns $q_j,\lambda _m$; see, for instance, the popular monographs [2, 3, 12, 15]. The non-holonomic constraints Eq. (2) are formulated as linear functions of velocities $\dot{q}_j$; this has been shown to be satisfactory concerning the problems considered. For more details and generalization, see [4, 44].

In order to arithmetize the mathematical model, we need to introduce three adequate coordinate systems. In accordance with Fig. 2, the fixed Cartesian coordinates (x, y, z) are obvious. Their origin is in the “Southern Pole of the Cavity” (SPC) denoted A. The position of the ball center is described in standard spherical coordinates with their origin being in the center of the cavity and $\alpha ,\gamma$ denoting the polar and azimuthal angles, respectively. The origin of the moving coordinates is located in the center of the moving ball, so that it lies on the concentric sphere with radius $\varrho =R-r$. Moving axis p follows a tangent of the concentric sphere meridian in the vertical plane ($x',z$), axis q is always horizontal, and axis n has the direction of the upward directed normal at the contact of both bodies; see Fig. 2. Rotation of the ball with respect to the moving coordinates is represented by the Euler angles $\varphi ,\theta ,\psi$. Components of the angular velocity vector ${\varvec{\omega }}=[\omega _p,\omega _q,\omega _n]^T$ in moving coordinates are positive as corresponds to the usual convention of the right hand rule [40]. The velocities of the ball center with respect to global coordinates are $\mathbf{v}=[v_p,v_q,v_n]^T$.

In further text, only the components of vectors $\mathbf{v},{\varvec{\omega }}$ will be used. Nevertheless, their relation to the velocity components $\dot{q}_j$ used in Eqs. (1), (2) is obvious.

Fig. 2

Outline of coordinate systems; left: axonometric view; right: plane $\xi z$ view—along $\gamma$ orientation

The basic formulae for kinetic and potential energies with respect to moving coordinates read

$$\begin{aligned} T&= \frac{1}{2}m\left( v_p^2+v_q^2+v_n^2+\frac{2}{5}r^2\left( \omega _p^2+\omega _q^2+\omega _n^2\right) \right) , \end{aligned}$$

(3a)

$$\begin{aligned} V&= mg\varrho (1-\cos \alpha ),\quad (\varrho =R-r), \end{aligned}$$

(3b)

where m, g represent the mass of the ball and gravitational acceleration, respectively. In order to relate the angular velocity vector ${\varvec{\omega }}$ and the velocity components expressed in Euler angles $\varphi ,\theta ,\psi$, we write

$$\begin{aligned} \omega _p&=-{\dot{\varphi }}\sin \theta \cos \psi +{\dot{\theta }}\sin \psi , \end{aligned}$$

(4a)

$$\begin{aligned} \omega _q&= {\dot{\varphi }}\sin \theta \sin \psi +{\dot{\theta }}\cos \psi , \end{aligned}$$

(4b)

$$\begin{aligned} \omega _n&= {\dot{\varphi }}\cos \theta +{\dot{\psi }}. \end{aligned}$$

(4c)

Because the spherical cavity has a constant curvature $1/\varrho$ at every point regardless of direction, the relation between angles $\alpha ,\gamma$ and velocities $v_p,v_q$ can be expressed simply, see Fig. 2,

$$\begin{aligned} v_p\;=\;\varrho {\dot{\alpha }},\quad v_q\;=\;\varrho {\dot{\gamma }}\sin \alpha , \quad v_n\;=\;0, \end{aligned}$$

(5)

where the expression $v_n=0$ represents one of the three contact constraints. For the same reasons, the Pfaff contact conditions of perfect rolling without slipping can be easily expressed. They specify the relations between angular velocities $\omega _p,\omega _q,\omega _n$ and angles $\alpha ,\gamma$ or displacements $v_p,v_q$. With respect to Eq. (5), the contact conditions can be reformulated as follows:

$$\begin{aligned} \begin{aligned} v_p-r\omega _q=0&,&r\omega _q - \varrho {\dot{\alpha }}=0&, \\ v_q+r\omega _p=0&,&\ \Longrightarrow \quad&r\omega _p + \varrho {\dot{\gamma }}\sin \alpha =0&, \\ v_n =0&,&v_n =0&. \end{aligned} \end{aligned}$$

(6)

Taking into account Eqs. (5) and (6), the expressions for energies from Eq. (3) can be rewritten in the form

$$\begin{aligned} T&= \frac{1}{2}m\left( \frac{7}{5}\varrho ^2({\dot{\alpha }}^2+{\dot{\gamma }}^2\sin ^2\alpha ) + \frac{2}{5}r^2 ({\dot{\psi }}+{\dot{\gamma }}\cos \alpha )^2\right) , \end{aligned}$$

(7a)

$$\begin{aligned} V&= mg\varrho (1-\cos \alpha ). \end{aligned}$$

(7b)

With reference to the problem definition, no external excitation or energy dissipation is assumed. Hence, the internal energy introduced to the system by non-homogeneous initial conditions has the form:

$$\begin{aligned} \begin{aligned} E_0=T+V=&\frac{1}{2}m\left( \frac{7}{5}\varrho ^2({\dot{\alpha }}^2+{\dot{\gamma }}^2\sin ^2\alpha ) + \frac{2}{5}r^2 \omega _n^2\right) \\&+ mg\varrho (1-\cos \alpha ). \end{aligned} \end{aligned}$$

(8)

Here, the spin of the moving sphere is given in the form $\omega _n={\dot{\psi }}+{\dot{\gamma }}\cos \alpha$ with respect to Eq. (4c) and the geometric properties of the cavity.

Every conservative Lagrangian system (in the sense of energy balance) possesses at least one first integral, which can be considered a multiple of the total energy of the system. Therefore, with respect to Eq. (8), it can be formulated as follows:

$$\begin{aligned} {\dot{\alpha }}^2+{\dot{\gamma }}^2\sin ^2\alpha + \mu \omega _n^2+2\omega _0^2 (1-\cos \alpha )\;=\;E, \end{aligned}$$

(9)

where the following notations were adopted:

$$\begin{aligned} \mu =\frac{2r^2}{7\varrho ^2},\quad \omega _0^2=\frac{5g}{7\varrho }, \quad E=\frac{10E_0}{7m\varrho ^2}. \end{aligned}$$

(10)

The dimensionless constant $\mu$ reflects the ratio between the ball radius and the distance between the centers of the ball and the cavity; it is, in fact, closely related to the ratio of both radii. $\omega _0^2$ denotes the squared natural frequency of the rocking of the homogeneous ball in the spherical cavity, and E relates to the internal energy of the system.

Due to the transparent structure of the problem, although it is non-holonomic, it can be rewritten in three degrees of freedom only, and no procedure via Lagrangian multipliers has to be applied. Moreover, since no explicit external excitation is present, both terms on the right-hand side of Eq. (1) identically vanish. Therefore, after some modifications, we obtain three Lagrangian equations for the three unknowns $\alpha$, $\gamma$ and $\omega _n$:

$$\begin{aligned}&\alpha {:}\, \ddot{\alpha }-({\dot{\gamma }}^2\cos \alpha -\mu \omega _n{\dot{\gamma }}-\omega _0^2)\sin \alpha = 0, \end{aligned}$$

(11a)

$$\begin{aligned}&\gamma {:}\, \frac{{ \text{ d }}}{{ \text{ d }}t}\left( {\dot{\gamma }}\sin ^2\alpha + \mu \omega _n\cos \alpha \right) = 0, \end{aligned}$$

(11b)

$$\begin{aligned}&\psi {:}\, \frac{{ \text{ d }}}{{ \text{ d }}t}\left( \omega _n\right) = 0. \end{aligned}$$

(11c)

It is obvious that Eqs. (11b,c) are the first integrals, because $\gamma$ and $\omega _n$ are cyclic coordinates. Therefore, we can write

$$\begin{aligned}&\gamma {:}\,\, {\dot{\gamma }}\sin ^2\alpha + \mu \omega _n\cos \alpha = H,\quad H= \frac{5H_0}{7m\varrho ^2}, \end{aligned}$$

(12a)

$$\begin{aligned}&\omega _n{:}\,\, \omega _n = S. \end{aligned}$$

(12b)

Parameters $E_0$ and $H_0$ represent the energy and angular momentum of the system, respectively; they are given by Eqs. (8,12a) or (9, 12a), where the initial conditions $\alpha ,{\dot{\alpha }},{\dot{\gamma }},\omega _n$ are substituted.

Equations (12) are the second and third first integrals of the system. The first one represents the conservation of the system’s angular momentum with respect to axis z at the constant level H, while Eq. (12b) shows the spin velocity (proportional to S) of the ball with respect to the normal n. We can see that the spin velocity is constant throughout the whole period of the system movement, although it interacts with the other angular velocity components of the ball through relations Eq. (4). Note that this very special character of the spin is a direct consequence of the spherical shape of both the cavity and the ball. Any other combination of two bodies in contact would lead to a variable spin velocity. The general form of Eq. (11c) would encompass an additional term dependent on differences in the principal curvatures of the two bodies. This difference is zero for spherical surfaces.

2.2 Characteristic equation

The three first integrals Eqs. (9, 12) represent certain invariants and an alternative description of the system behavior; they provide a possibility for a transparent introduction of energy and movement through initial parameters. Thus, they enable the investigation of particular states of the system, an analytical formulation of various characteristics of trajectories and, consequently, much greater insight into the nature of system parameters than one based solely on a numerical integration of the original differential system.

Let $\delta$ denote the height of the ball above the bottom of the cavity, i.e., the z coordinate of the center of the sphere:

$$\begin{aligned} \delta = 1-\cos \alpha \quad \Rightarrow \quad {\dot{\alpha }}={\dot{\delta }}/\sin \alpha , \end{aligned}$$

(13)

If ${\dot{\gamma }}$ is eliminated from Eq. (9) using Eq. (12a), one obtains after some manipulation

$$\begin{aligned} {\dot{\delta }}^2&=\; (E-\mu \omega _n^2-2\omega _0^2\delta )(2\delta -\delta ^2)-(H-\mu \omega _n(1-\delta ))^2 \nonumber \\&=\ f(\delta ) \end{aligned}$$

(14a)

$$\begin{aligned} {\dot{\gamma }}&=\; \displaystyle {\frac{H-\mu \omega _n(1-\delta )}{2\delta -\delta ^2}}. \end{aligned}$$

(14b)

In further text, we term $f(\delta )$ defined by Eq. (14) the “Characteristic Function” (CF), and $f(\delta )=0\;$ the “Characteristic Equation” (CE). Symbols E and H in Eq. (14) represent the measure of energy introduced into the ball at the instant $t=0$ by means of the initial conditions: ${\delta _c}$—the initial height of the ball, ${{\dot{\gamma }}_c}$—“Initial Horizontal Velocity” (IHV), and ${\omega _{n}}$—“Initial Spin Velocity” (ISV); see Eqs. (9, 12a). Based on these initial conditions, the energy in the system may be quantified using the following relations:

$$\begin{aligned} \begin{aligned} E&= {\dot{\gamma }}_0^2\delta _{c}(2-\delta _{c})+2\omega _0^2\delta _{c}+\mu \omega _n^2, \\ H&= {\dot{\gamma }}_0\delta _{c}(2-\delta _{c})+\mu \omega _n(1-\delta _c). \end{aligned} \end{aligned}$$

(15)

Function $f(\delta )$ is a cubic parabola which attains positive or negative values based on system parameters $\mu ,\omega _0$, state variables $\alpha , {\dot{\alpha }},{\dot{\gamma }},\omega _n$ together with their initial values $\alpha _c,\dot{\alpha _c},{\dot{\gamma }}_0,\omega _{n0}$ (hidden in E, H) and with respect to independent variable $\delta$. The general form of the CF is obvious in Fig. 3. The interval $\delta \in (0,2)$ spans the whole diameter 2R of the cavity from the SPC to the NPC (NPC—“Northern Pole of the Cavity”).

Obviously, it holds that

$$\begin{aligned} \begin{aligned} f(-\infty )&< 0,\\ f(0)&=-(H-\mu \omega _n)^2,\ f(2)=-(H+\mu \omega _n)^2, \\ f(1)&=E-H^2-2\omega _0^2-\mu \omega _n^2,\\ f(\infty )&>0. \end{aligned}\nonumber \\ \end{aligned}$$

(16)

Fig. 3

General shape of characteristic function $f(\delta )$ and of the area delimiting the active spherical strip in the interval $\delta \in (\delta _1,\delta _2)$

The cubic polynomial Eq. (14a) is physically meaningful only for values where $f(\delta )>0$, as ${\dot{\delta }}$ is considered to be real. The first derivative of $f(\delta )$ may be formally written as a general quadratic polynomial:

$$\begin{aligned} \frac{ { \text{ d }}f(\delta )}{{ \text{ d }}\delta }&=A\delta ^2-2B\delta +C, \nonumber \\ A&=6\omega _0^2,\quad B=4\omega _0^2+E-\mu \omega _n^2(1-\mu ),\nonumber \\ C&=2\left( E-H\mu \omega _n -\mu \omega _n^2(1-\mu )\right) . \end{aligned}$$

(17)

Quadratic equation Eq. (17) has two real roots because its discriminant is always positive:

$$\begin{aligned}&B^2-A\cdot C>0\;\Rightarrow \;\nonumber \\&\left( E-\mu \omega _n^2(1-\mu )-2\omega _0^2\right) ^2+12\omega _0^2\left( \omega _0^2+\mu \omega _nH\right) >0. \end{aligned}$$

(18)

The inequality Eq. (18) is fulfilled trivially for $\omega _n H>0$. Otherwise, introduction of initial conditions $\delta _c,{\dot{\gamma }}_0$ and $\omega _n$ into E, H, defined by Eqs. (9, 12a), implies that Eq. (18) is valid for $0\le \delta _c<2$. Therefore, the cubic parabola Eq. (14a) has two extremes. The analogous procedure confirms that $C>0$. Thus,

$$\begin{aligned} B^2>B^2-AC, \end{aligned}$$

(19)

and both the extremes are situated on the positive semi-axis: $0\le \delta _{e1}\le \delta _{e2}$; the first one is positive, and the second is negative:

$$\begin{aligned} f(\delta _{e1})\ge 0,\quad f(\delta _{e2})\le 0. \end{aligned}$$

(20)

Summarizing the above contemplation, one can conclude that the CE Eq. (14a) has three real roots satisfying the following conditions:

$$\begin{aligned} 0\le \delta _1\le \delta _2<2<\delta _3. \end{aligned}$$

(21)

The first two roots are physically meaningful, as they delimit an interval on axis $\delta$ where ${\dot{\delta }}^2\ge 0$, which is a necessary condition for the energy accumulated in the system to be real. For geometrical reasons, the values $\delta>\delta _3>2$ do not represent a physically meaningful state, although ${\dot{\delta }}^2\ge 0$ there as well. Note that zero and the coinciding roots can occur. As we will see later, they represent physically important cases.

The values used in the majority of the plots are as follows: $R=1$ m, $r=1/4R$, $M=1$ kg, $J=2/5 M r^2$ kg m$^2$, $g=9.81\,\hbox {m}\,\hbox {s}^{-2}$. If not stated otherwise, the initial position used in the example is given as $\gamma =0$ and $\delta _c=1-1/\sqrt{2}\doteq 0.3$, which corresponds to the polar angle $\alpha _c=\pi /4$; in the Cartesian coordinates, it is $(x,y,z)=(1/\sqrt{2},0,1-1/\sqrt{2})$.

3 The separation circle and its neighborhood

3.1 Definition and relevance of the Separation Circle

Previous studies by the authors [36, 37] have shown that the most important separation limit between trajectory types (or groups) that start at a certain point is a trajectory running at constant angular velocity $\varGamma$ along a parallel of the cavity. This trajectory will be called the “Separation Circle” (SC), because it separates qualitatively different trajectories. The SC can be characterized as follows: The ball is pulled up along a meridian of the cavity to a certain level $\delta _c$ to the “Starting Point of the Trajectory” (SPT), and subsequently, a horizontal impulse is applied to it. The intensity of the impulse is so high that it sets the ball onto a horizontal circular trajectory at the vertical level $\delta _c$. This condition can be used reversely to determine the necessary IHV ${\dot{\gamma }}_c=\varGamma$ or $v_{qc}$, the values of which represent the constant angular or tangential velocities relevant to movement along the SC. The initial spin velocity $\omega _n$ is considered zero in the first step. The spin velocity can also be set such that it compensates for ${\dot{\gamma }}_c$ when different from $\varGamma$ in order to maintain the SC trajectory.

Let us evaluate the notion of the SC from the viewpoint of the CE discussed above. Figure 4a shows parabola $f(\delta )$ for a given height $\delta _c$ and three different values of the IHV. Intervals where $f(\delta )\ge 0$ for $\delta \in \langle 0,2\rangle$ represent possible heights at which the trajectory of the ball in the cavity can occur. In the case of the SC trajectory, i.e., for IHV ${\dot{\gamma }}_c=\varGamma$, one double root $\delta _1=\delta _2=\delta _c$ of $f(\delta )=0$ occurs; this case provides an active area of zero width, see the dotted curve in Fig. 4. The dependence of the width of the active area on the initial horizontal velocity ${\dot{\gamma }}_c$ is illustrated in Fig. 4b.

For ${\dot{\gamma }}_c<\varGamma$, the active area spans between $\delta _{1}$ and $\delta _{2}$ which coincide with $\delta _{c}$. An initial velocity higher than $\varGamma$ leads to a trajectory within the spherical strip above the $\delta _{c}=\delta _{1}$ boundary and goes up to $\delta _{2}$. In such a case, it can occur that $\delta _{2}>1$, which means that the ball passes into the upper hemisphere of the cavity. The limit case is reached when the IHV approaches an infinite value. Then, $\delta _{1},\delta _{2}$ are symmetrically distributed with respect to $\delta =1$. The upper boundary of the active strip is represented by the SC mirrored in the upper hemisphere. The trajectory becomes planar again although the plane is slanted passing the SPT and the center of the cavity.

The graph in Fig. 4b also shows the effect of nonzero initial spin $\omega _n$. The dashed curve shows the case when $\omega _n<0$, and the dot-dashed curve corresponds to positive $\omega _n$. Each of the curves also has a second branch above $\delta =2$ that have no physical meaning.

Fig. 4

Active area: a solid curve —trajectory below the SC: $\delta \in (\delta _{1,b},\delta _{2,b}=\delta _c)$; solid curve —trajectory above the SC: $\delta \in (\delta _{1,a}=\delta _c,\delta _{2,b})$; dashed curve —transition case representing the SC—no active area due to coincidence $\delta _{c}=\delta _1=\delta _2$. b Schema of the dependence of $\delta$ on the initial horizontal angular velocity ${\dot{\gamma }}_c/\varGamma$ for three values of $\omega _n$: solid curve : $\omega _n=0$, dashed curve : $\omega _n<0$, dot-dashed curve : $\omega _n>0$. Verticals – correspond to curves – in plot a)

The classification strategy based on the SC was intuitively adopted by the authors in their earlier studies, e.g.,   [36, 37]. Actually it appears that this classification well describes all possible trajectories of a ball rolling inside a spherical cavity and starting from a certain point. Whatever the orientation and intensity of the initial impulse are, the movement of the ball takes place within a uniquely defined spherical strip delimited by the two lower roots $\delta _1,\delta _2$ of the CE, Eq. (14), which are related to the energy contained in the system.

3.2 Position of the separation circle

Let us inspect Eq. (11a). Provided the ball follows the SC at a height given by initial condition $\alpha _c<\pi /2$, i.e., $\delta _c<1$, its angular velocity is constant, ${\dot{\gamma }}_0=\varGamma$, and vertical acceleration $\ddot{\alpha }$ vanishes, so that we can write

$$\begin{aligned} \varGamma ^2\cos \alpha _c -\mu \omega _n\varGamma -\omega _0^2 = 0, \quad 0<\alpha _c<\pi /2. \end{aligned}$$

(22)

This quadratic equation has two real roots:

$$\begin{aligned} \varGamma= & {} \frac{\mu \omega _n \pm \sqrt{\mu ^2\omega _n^2+4\omega _0^2\cos \alpha _c}}{2\cos \alpha _c} \;\Rightarrow \; \nonumber \\ v_{qc}= & {} \frac{1}{2}\varrho \cdot { \text{ tg }}\alpha _c \left( \mu \omega _n \pm \sqrt{\mu ^2\omega _n^2+4\omega _0^2\cos \alpha _c}\right) .\nonumber \\ \end{aligned}$$

(23)

When zero spin velocity is assumed, $\omega _n=0$, these relations simplify greatly:

$$\begin{aligned} \varGamma _0= \pm \frac{\omega _0}{\sqrt{\cos \alpha _c}}= \pm \frac{\omega _0}{\sqrt{1-\delta _c}}. \end{aligned}$$

(24)

The IHV $v_{qc}$ or angular velocity $\varGamma$ satisfying Eq. (23) produces the SC for the given polar angle $\alpha _c$. The roots represent two opposite directions with respect to the trajectory initial point. Their ratio depends on the sign of the spin velocity $\omega _n$. For zero initial spin, the image is symmetrical in the horizontal plane with respect to the initial point of the trajectory. The schematic plots in Fig. 5 demonstrate the dependence of velocity $v_{qc}$ on both the initial height (given by parameter $\alpha _c$) and the ratio of the sphere to the cavity. Graph (a): fixed height $\alpha _c= \pi /4$, plot (b): fixed ratio r/R. In this latter case, the dependence starts from zero for the SPC ($\alpha _c=0$) and tends to infinity for the “Equator of the Cavity” (EQC) ($\alpha _c=\pi /2$). The bold-black curves in plots (a) and (b) represent the spin-free state $(\omega _n=0)$, while the color curves show the influence of the initial spin. The relation between the velocity $v_{qc}$ and the initial spin $\omega _n$, which maintains the trajectory in the SC, may be deduced from Eq. (22) and is illustrated in picture (c). A decrease in horizontal velocity is compensated by a negative spin, whereas an increase in horizontal velocity implies a proportional increase in positive spin. One-sided limits of $\omega _n$ for $v_{qc}\rightarrow 0^+$ and $v_{qc}\rightarrow \infty$ exist and equal $\pm \infty$.

Fig. 5

Initial horizontal velocity producing the trajectory of the SC: a fixed initial height $\alpha _c= \pi /4$, varying ratio r/R, b fixed ratio r/R, varying polar angle $\alpha _c$; color curves in a, b correspond to various spin values, c compensation spin. (Color figure online)

3.3 Dynamic stability of the separation circle

We now examine the neighborhood near the SC. We revisit Eq. (11a) and also the first integral Eq. (12a). The vertical position of the ball on the SC is given by the angle $\alpha _c$, and the horizontal angular velocity $\varGamma$ is determined by Eq. (23). These are the nominal values which are subjected to small perturbations generated by dispersion in the initial conditions setting. Thus, we reformulate the initial conditions as

$$\begin{aligned} {\dot{\gamma }}_0\approx {\varGamma }+\eta ,\quad \alpha \approx \alpha _c+\zeta , \end{aligned}$$

(25)

where $\zeta$ and $\eta$ are small values.

Substituting perturbed initial conditions Eq. (25) into Eq. (11a), a linearized equation for $\zeta (t)$ can be deduced. Disregarding the higher-order terms $\zeta ^2,\eta ^2,\eta \cdot \zeta$, one obtains

$$\begin{aligned} \begin{aligned}&\ddot{\alpha }_c-\left( \varGamma ^2\cos \alpha _c-\mu \omega _n\varGamma \!-\omega _0^2\right) \sin \alpha _c \\&\quad +\,\ddot{\zeta }-\left( \left( \varGamma ^2\cos \alpha _c-\mu \omega _n\varGamma \! -\omega _0^2\right) \cos \alpha _c \right. \\&\quad \left. -\,\varGamma ^2\sin ^2\alpha _c+\right) \zeta \\&\quad -\,\left( 2\varGamma \cos \alpha _c-\mu \omega _n\right) \sin \alpha _c \cdot \eta = 0. \end{aligned} \end{aligned}$$

(26)

The first line of Eq. (26) vanishes due to Eq. (11a), and the coefficient of the term with $\cos \alpha _c$ on the second line disappears because of Eq. (22). Hence, it holds that

$$\begin{aligned}&\ddot{\zeta } + \varGamma ^2\sin ^2\alpha _c\cdot \zeta \nonumber \\&\quad +\,(\mu \omega _n-2\varGamma \cos \alpha _c)\sin \alpha _c \cdot \eta = 0. \end{aligned}$$

(27)

This equation is solvable for zero initial conditions in a closed form:

$$\begin{aligned} \zeta = \frac{2 \eta \left( 2 \varGamma \cos \alpha _c -\mu \omega _n\right) }{\varGamma ^2\sin \alpha _c}\sin ^2\!\left( t\cdot \frac{1}{2}\varGamma \sin \alpha _c\right) .\nonumber \\ \end{aligned}$$

(28)

At the level of the first approximation, supposing that a small increase in initial tangential velocity is considered, it is obvious that the trajectory lies above the SC within the narrow spherical strip $\delta \in (\delta _{c}=\delta _1, \delta _2)$, where $\delta _1<\delta _2$. Similarly, decreasing the IHV, we get a trajectory below the SC within the limits $0<\delta _1<\delta _2=\delta _{c}$. The width of the strip in both cases is $|\eta |/\varOmega _v$.

The explicit form of perturbation Eq. (28) represents a form of harmonic ripples which oscillate above or below the separation circle. The period of perturbations is

$$\begin{aligned} T_{\mathrm{per}}=\frac{4\pi }{\varGamma \sin \alpha _c}. \end{aligned}$$

(29)

One loop around the SC takes $T_{\mathrm{loop}}\approx 2\pi /(\varGamma +\eta )$. Therefore, the number of small waves during one loop is

$$\begin{aligned} N_{\mathrm{per}}=\frac{2(\varGamma +\eta )}{\varGamma \sin \alpha _c}\approx \frac{2}{\sin \alpha _c}, \end{aligned}$$

(30)

which, in general, is not a rational number, and the trajectory alternating above or below the SC does not pass the SPT.

Fig. 6

Example of the trajectory of the SC type without influence of the ISV $(\omega _n=0)$: a time history, b top view, c axonometric demonstration

We should be aware that the estimates Eq. (25) and also the results Eq. (28) are applicable if $0<\alpha _c<\pi /2$. Indeed, for small $\alpha _c$, the values $\alpha _c$ and $\zeta$ are commeasurable as are the values $\varGamma ,\eta$, and classification according to the powers of a small parameter becomes invalid. Of course, the perturbations $\eta ,\zeta$ should also remain small in order for linear approximation to be justified.

Finally, one can conclude that the SC (perhaps with the exception of the SPC neighborhood) is dynamically stable, and a small initial perturbation does not cause a receding of this trajectory from the original SC.

Let us show an example of a time history of the trajectory corresponding to the SC; see Fig. 6. The time history shows a harmonic process in both horizontal coordinates and a constant value in the vertical coordinate, plot (a).

4 Trajectories not influenced by initial spin of the ball

4.1 Trajectories above the SC

4.1.1 Roots of the CE

Let us examine the case where $\delta _c=\delta _1\le \delta _2$ and $\omega _n=0$, i.e., the ball is rolling above the SC and no spin is assumed, ${\dot{\gamma }}_0>\varGamma _0$; see also curve (a) in Fig. 4.

Setting the position for the SPT on a meridian of the cavity at a certain level characterized by polar angle $0<\alpha _c<\pi /2$ or equivalently by height $\delta _1\in (0,1)$ (in the lower hemisphere of the cavity) in fact means that the lowest root $\delta _1=\delta _{c}$ of the CE is fixed. Considering that one root is known, we can rewrite Eq. (14a) in the form of a partial decomposition with respect to the root factors:

$$\begin{aligned} f(\delta )&=(\delta -\delta _{c})(K\delta ^2+2L\delta +M)= 0, \nonumber \\ K&=2\omega _0^2, \quad L=-\left( 2\omega _0^2+\frac{1}{2}{\dot{\gamma }}_0^2\delta _{c}(2-\delta _{c})\right) , \nonumber \\ M&={\dot{\gamma }}_0^2\delta _{c}\left( 2-\delta _{c}\right) ^2, \qquad {\dot{\gamma }}_0>\varGamma _0. \end{aligned}$$

(31)

The detailed form of coefficients K, L, M is derived from the full CE, Eq. (14a), where a zero ISV was substituted. The energy E, Eq. (9), and the angular momentum H, Eq. (12a), contained in Eq. (14a), were included in the initial parameter values for the SPT.

The SC is regarded as the lower boundary of the spherical strip on the cavity surface within which a particular trajectory runs. The upper boundary is then the lower of the remaining roots $\delta _2,\delta _3$. They can be calculated from the quadratic equation, which is outlined in Eq. (31):

$$\begin{aligned} \delta _{2,3}&= \frac{1}{K}\left( -L\pm \sqrt{L^2-K\cdot M}\right) , \end{aligned}$$

(32a)

$$\begin{aligned} D&=L^2-K\cdot M= \end{aligned}$$

(32b)

$$\begin{aligned}&= \left( 2\omega _0^2+\frac{1}{2}{\dot{\gamma }}_0^2\delta _{c}(2-\delta _{c})\right) ^2 -2\omega _0^2{\dot{\gamma }}_0^2\delta _{c}(2-\delta _{c})^2 \nonumber \\&= \left( 2\omega _0^2-\frac{1}{2}{\dot{\gamma }}_0^2\delta _{c}(2-\delta _{c})\right) ^2 +2\omega _0^2{\dot{\gamma }}_0^2\delta _{c}^2(2-\delta _{c})\;>\;0. \end{aligned}$$

(32c)

The discriminant D is always positive, Eq. (32c), which fact confirms that Eq. (14a) has three real roots. Furthermore, $0<D<L^2$ and, consequently, all roots are positive and fulfill conditions $0<\delta _{c}=\delta _1<\delta _2<2$ and $\delta _3>2$. Of course, root $\delta _3$ is geometrically out of scope, and it is no longer considered; see Sect. 2.2.

4.1.2 Circumferential periodicity of trajectories

In order to outline the basic character of a trajectory occurring between boundaries $\delta _{c}=\delta _1, \delta _2$, we inspect the expression for the circumferential velocity ${\dot{\gamma }}$ following from Eq. (14b). In general, it can be assumed that the trajectory is periodical in the vertical direction, where the ratio of the period to the SC lengths is not a rational number. Therefore, one round along the SC will not contain a whole number of periods. The angular momentum H is always positive, see Eq. (12a), and the second term of the numerator vanishes since $\omega _n=0$. The denominator in this fraction is also positive, because $\delta (2-\delta )>0$ for $\delta \in (0,2)$. Therefore, ${\dot{\gamma }}>0$ regardless of the settings for the initial conditions at the SPT. Moreover, it can be supposed that the variability of ${\dot{\gamma }}$ during one period for given initial settings will not be dramatic. Indeed, it is obvious that the angular momentum H is constant during one period. Consequently, the horizontal angular velocity at the point where the trajectory touches the upper boundary of the strip, $\delta _2$, follows from Eq. (12), where $\omega _n=0$, and it holds that

$$\begin{aligned} \dot{\gamma _2}\sin ^2\alpha _2\;=\;{\dot{\gamma }}_0\sin ^2\alpha _c, \end{aligned}$$

(33)

where $\alpha _c,{\dot{\gamma }}_0$ or $\alpha _2,{\dot{\gamma }}_2$ are values of the respective parameters at the initial point $(\delta _c=\delta _1)$ or at the touching point on the upper boundary $(\delta _2)$. Hence, it can be written:

$$\begin{aligned} {\dot{\gamma }}_{2T}&={\dot{\gamma }}_0\;\left( \frac{\sin \alpha _c}{\sin \alpha _2}\right) ^2 ={\dot{\gamma }}_0\;\frac{\left( \delta _c-2\right) \delta _c}{\left( \delta _2-2\right) \delta _2} \end{aligned}$$

(34a)

$$\begin{aligned} {\dot{\gamma }}_{2T}&\approx {\dot{\gamma }}_0\left( 1-2\varDelta _{\alpha }\cot \alpha _c\right) \approx {\dot{\gamma }}_0\left( 1-\frac{2\left( \delta _c-1\right) \varDelta _{\delta }}{\left( \delta _c-2\right) \delta _c}\right) \end{aligned}$$

(34b)

where Eq. (34b) is valid for small values of the difference $\varDelta _\alpha =\alpha _2-\alpha _c$ or $\varDelta _{\delta }=\delta _2-\delta _c$.

The horizontal velocity is slightly lower at the upper apex due to an increase in potential energy. At the same time, some qualitative confirmation of this fact follows from the constant angular momentum H and $\delta$, which changes more or less monotonously. This implies that no singular points emerge during one vertical period, and that the angular velocity along the trajectory is mildly variable. Thus, the trajectory has the shape of a simple periodic curve without any turnabout points reversing velocity ${\dot{\gamma }}$.

4.1.3 Vertical periodicity of trajectories

Here, we assess the length and duration of one vertical period of the trajectory. The IHV leads to a trajectory that is symmetrical with respect to both points where it touches the lower boundary ($\delta _c=\delta _1$, where the IHV is applied) and the upper boundary ($\delta _2$); see Fig. 7. Thus, it is sufficient to examine only the first half of the period for $\delta$ increasing between $\delta _1$ and $\delta _2$. We revisit both relations in Eq. (14). Assuming $\omega _n=0$, the following differential system can be written:

$$\begin{aligned} {\dot{\delta }}^2&= (E-2\omega _0^2\delta )(2\delta -\delta ^2)-H^2, \end{aligned}$$

(35a)

$$\begin{aligned} {\dot{\gamma }}&= \frac{H}{2\delta -\delta ^2}. \end{aligned}$$

(35b)

The first equation is ${\dot{\gamma }}$ independent, and, therefore, it can be solved as the first step. The pair $\pm \delta$ could be put into the second equation to obtain $\gamma (t)$ by means of integration between $\delta _c$ and $\delta _2$. However, it comes to light that both points on the strip boundaries $\delta _c,\delta _2$ are singular and represent bifurcation points. In addition, the relevant Jacobi matrix is also singular and, consequently, does not enable us to predict the principal directions in the point neighborhood. Three solutions start from the SPT (including the constant $\delta _c$), and all of them have the zero derivative. Therefore, neither an analytical nor a numerical stable solution can be deduced from these points. This difficulty can be overcome by differentiating Eq. (35a) with respect to time. After reducing by ${\dot{\delta }}$, the modified system reads

$$\begin{aligned} \ddot{\delta } =&\; E(1-\delta )-\omega _0^2\delta (4-3\delta ), \end{aligned}$$

(36a)

$$\begin{aligned} {\dot{\gamma }} =&\; \frac{H}{2\delta -\delta ^2}, \end{aligned}$$

(36b)

where $E= {\dot{\gamma }}_0^2\delta _{c}(2-\delta _{c})+2\omega _0^2\delta _{c},\ H\;=\; {\dot{\gamma }}_0\delta _{c}(2-\delta _{c})$, cf. Eq. (15).

Fig. 7

Shape of the trajectory above the SC for various IHV, $\alpha _c=\pi /4$, $\delta _c={1-1/\sqrt{2}\approx }0.3$

Fig. 8

Trajectory above the SC (no spin applied), a spatial period dependent on the IHV $({\dot{\gamma }}_0)$; b upper boundary of the strip $\delta _2$ as a function of ${\dot{\gamma }}_0$ ($\alpha _c=\pi /4$, $\delta _c={1-1/\sqrt{2}\approx }0.3$)

The denominator in Eq. (36b) is always positive because $\delta \in (\delta _c,\delta _2)$. Then, Eq. (36) can be solved for initial conditions $\gamma (0)=0$, $\delta (0)=\delta _c$, ${\dot{\delta }}(0)=0$ sequentially putting partial results of Eq. (36a) into Eq. (36b). Finally, the time is eliminated, and one obtains $\delta$ as a function of $\gamma$. Some samples of the trajectory are plotted in Fig. 7, in which a set of five cases were ${\dot{\gamma }}_0>\varGamma _0$ for the selected $\delta _c={1-1/\sqrt{2}\approx }0.3$ and the relevant $\varGamma _0=3.63$ is demonstrated.

Notice the shape of the high curves in this figure. They correspond to the dominating first term in the right-hand side of Eq. (36a) when ${\dot{\gamma }}_0$ is large. Omitting the second term when ${\dot{\gamma }}_0\gg \varGamma _0$, we obtain a linear equation of harmonic oscillation with an adequate constant right-hand side. The resulting equation for ${\dot{\gamma }}_0\rightarrow \infty$ is analytically resolvable; its complicated resulting expression $\delta (\gamma )$ is $2\pi$ periodic with the extremal values $\delta _c$ and $2-\delta _c$ at $\gamma =2k\pi$ and $2(k+1)\pi ,\, k\in {\mathbb {Z}}$, respectively. The shape of the limiting function is indicated by the dashed black curve segments. This limit case of ${\dot{\gamma }}_0\rightarrow \infty$ will be discussed separately in Sect. 4.2. On the other hand, $\gamma$ changes in steps from 0 to $\pi$ whenever $\delta (t)$ passes the SPC for ${\dot{\gamma }}_0=0$, which makes $\gamma _T=\pi$.

In a general case, for given initial height $\delta _c$, the vertical component of the trajectory $\delta (\gamma )$ increases its angular period and amplitude as the initial angular velocity ${\dot{\gamma }}_0$ increases. This can be observed in Fig. 8; the blue curves above the SPT, ${\dot{\gamma }}_0>\varGamma _0$ are regarded in this section. For a fixed $\varGamma _0$ and increasing IHV, the angular period ${\gamma }_T$ slowly rises and approaches the full angle for ${\dot{\gamma }}_0\rightarrow \infty$; see picture (a). Picture (b) demonstrates the rising of the strip’s upper level $\delta _2$ with increasing IHV. The parameter $\delta _2$ approaches the horizontal asymptote $2-\delta _c$, symmetrically placed with respect to the EQC, in the particular case from 0.3 to 1.7. When $\delta _c$ changes its value, the dependences shown in Fig. 8 remain valid, but the strip $(\delta _c,2-\delta _c)$ decreases/increases its width.

Fig. 9

Example of the trajectory above the SC with no ISV $(\omega _n=0)$: a time history, b top view, c axonometric demonstration

4.1.4 Discussion

Using the above analytical results, we can outline some detailed trajectory properties with respect to IHV and the height of the SPT above the SPC. A typical trajectory shape is plotted in Fig. 9. Comparing Figs. 6 and 9, we can see that the trajectory time history is still a simple periodic curve synchronous in both horizontal coordinates with a slight modulation. This modulation depends on both values of $\delta _c$ and ${\dot{\gamma }}_0$. It relates to the difference between frequencies of $\gamma (t)$ and $\delta (t)$ and vanishes when ${\dot{\gamma }}_0=\varGamma _0$ or ${\dot{\gamma }}_0\rightarrow \infty$.

The frequency of vertical displacement $u_{cz}$ in Fig. 9a is related to that of horizontal displacements via the angular period ${\gamma }_T$ which depends on the $\delta _1=\delta _{c}$ and $\delta _2$ boundaries of the spherical strip. No visible influence of higher harmonics is observed, although a very light quasi-periodic character can be noticed, as it follows from the geometric character of the system. The shape in the top view is helical and resembles the form of a prolate-type hypotrochoid close to a hypotrochoid form of prolate types. As no sharp apexes or loop multiple points are detected, we can conclude that no basic cycloid or curtate trochoid is approached.

4.2 High initial horizontal velocity

4.2.1 General shape of the trajectory

Let us discuss the trajectories above the SC that emerge when the IHV distinctly exceeds the velocity $\varGamma _0$ or $v_{qc}$, and simultaneously, no ISV is applied, $\omega _{n}=0$. The CF for an infinite IHV has a form corresponding to Fig. 3, which means that zero points $\delta _1,\delta _2$ are symmetrically distributed with respect to the point $\delta =1$. A lower IHV shifts $\delta _2$ to the left, and the position of $\delta _1$ remains the same.

The trajectories are again concentrated within the spherical strip delimited by roots $0<\delta _1=\delta _{c}<\delta _2<2$ of the characteristic equation Eq. (14a). We will inspect its evolution when ${\dot{\gamma }}_0\gg \varGamma _0$. While it is always true that $\delta _1<1$, the upper boundary, being given by $\delta _2$, can enter the upper hemisphere of the cavity reaching a value in the interval $1<\delta _2<2$. The IHV corresponding to the transition case $\delta _2=1$ follows from equation Eq. (14a), where the two lowest roots $\delta _1,\delta _2$ are considered as known. After some manipulation, one can write

$$\begin{aligned} \begin{aligned} {\dot{\gamma }}_0^2&=\displaystyle {\frac{8\omega _0^2\cos \alpha _c}{\sin ^22\alpha _c}},\\&\quad 1-\cos \alpha _c=\delta _1{=\delta _c},\quad 0<\alpha _c<\pi /2, \end{aligned} \end{aligned}$$

(37)

where both the boundary values of $\alpha _c$ leading to an infinite IHV are obviously not admissible, as could be expected.

Increasing the IHV beyond all limits, the upper limit of the strip approaches a theoretical maximum:

$$\begin{aligned} \delta _2=2-\delta _1. \end{aligned}$$

(38)

It is an asymptotic position, which is monotonously approached as the initial velocity rises to infinity. In this theoretical state, the trajectory becomes planar, having a circular form with the diameter 2R. This plane is inclined, being determined by the horizontal tangent at the SPT and by the cavity center; see Fig. 11a.

4.2.2 Rotation of the osculating plane

In the case when ${\dot{\gamma }}_0$ is high but finite, $0<\varGamma _0\ll |{\dot{\gamma }}_0|<\infty$, the root $\delta _2$ is adequately lower:

$$\begin{aligned} \delta _2=2-\delta _1-\epsilon ,\quad 0<\epsilon \ll 1. \end{aligned}$$

(39)

Trajectories maintain their spatial character, but it is worthwhile to define an osculating plane at each point. The osculating plane makes sense from a physical perspective if it enables us to characterize the basic position of the trajectory using simpler elements than in the general case for lower IHV. In the case we are discussing, it can be assumed that the spiral does not differ much from a planar shape, and that its projection into the osculating plane at each point represents a good approximation. In other words, we can define an affine space which meets a sub-manifold at a point in such a way as to have a second order of contact at that point. The osculating plane passes the initial point $\delta _1$ and a point in the upper hemisphere at the upper boundary of the strip $\delta _2-\epsilon$, which is slightly below the level of the limit case Eq. (38). Therefore, its inclination is slightly lower than that corresponding to the limit case for ${\dot{\gamma }}_0\rightarrow |\infty |$. The osculating plane rotates around the vertical axis of the cavity with rotational speed $\varOmega _v$, which decreases as ${\dot{\gamma }}_0$ rises, and vanishes for an infinite IHV $({\dot{\gamma }}_0)$.

The velocity $\varOmega _v$ can be estimated based on the position of the contact point of the trajectory and the upper boundary $\delta _2$ evaluated at the conclusion of a single period. Although an explicit formula cannot be expressed, the process can be carried out observing the outline in Fig. 10. A sketch of the osculating plane behavior is graphically demonstrated in Fig. 11 for infinite IHV—picture (a), and finite IHV—pictures (b) and (c), showing position of this plane after the 1st and 3rd half-period, respectively. (A negative movement sense has been selected for graphical reasons.)

For a high but finite ${\dot{\gamma }}_0$, a ball starting at the SPT covers a distance of $2\pi -\eta$ along a nearly planar trajectory in the advancing osculating plane, where $\eta \rightarrow 0$ when ${\dot{\gamma }}_0\rightarrow |\infty |$; see Fig. 10. Let us again consider Eq. (36), in particular equation (a) and the expression for E. Assuming that velocity ${\dot{\gamma }}_0$ is high, the second term of the right-hand side in Eq. (36a) can be neglected and also enables E to be reduced. Then, it can be approximately written:

$$\begin{aligned} \ddot{\delta }+E_\gamma \delta \;=\;E, \end{aligned}$$

(40)

where coefficient $E_\gamma$ can be determined due to the equality of values at the lower and upper limits, $E_\gamma =E_{{0}}={\dot{\gamma }}_0^2\delta _c(2-\delta _c)$ and $E_\gamma =E_{\eta }=({\dot{\gamma }}_0-{\dot{\eta }})^2\delta _c(2-\delta _c)$, respectively. For initial conditions $\delta (0)=\delta _c, {\dot{\delta }}(0)=0$, the solution to equation Eq. (40) has the form

$$\begin{aligned} \delta =1+(\delta _c-1)\cos \left( t\sqrt{E_\gamma }\right) . \end{aligned}$$

(41)

The length of the period in time is

$$\begin{aligned} T_{{0}}^t=2\pi {E_{{0}}^{-\frac{1}{2}}}\quad \text{ or } \quad T_{\eta }^t=2\pi {E_{\eta }^{\frac{1}{2}}}, \end{aligned}$$

(42)

and, therefore, for the rotational speed of the osculation plane around the z axis, we can approximately write

$$\begin{aligned} \varOmega _v\;\approx \;1-({\dot{\gamma }}_0-\eta )/{\dot{\gamma }}_0. \end{aligned}$$

(43)

Since $\eta \rightarrow 0$ when ${\dot{\gamma }}_0\rightarrow |\infty |$, $\lim _{{\dot{\gamma }}_0\rightarrow |\infty |}\varOmega _v=0$ and $\delta _2\rightarrow 2-\delta _c$ as it corresponds with Eq. (38). These results can be confirmed intuitively by examining both of the graphs in Fig. 8. Indeed, the length of the period approaches the horizontal asymptote on the $2\pi$ level more or less exponentially. Because the tangential velocity along the trajectory is approximately constant, ${\dot{\gamma }}(t)\approx {\dot{\gamma }}_0+{\dot{\eta }}$, the relation between $\gamma (t)$ and t can be expressed simply as $\gamma (t)\approx ({\dot{\gamma }}_0+{\dot{\eta }}) t$. Then, it holds approximately that

$$\begin{aligned} \displaystyle { \varOmega _v\;\approx \;-\frac{ 2T_{{0}}-2T_{\eta }}{2T_{\eta }}\cdot {\dot{\gamma }}_0 \;\approx \;\frac{-\exp {(-\kappa {\dot{\gamma }}_0)}}{2T_{\eta }}\cdot {\dot{\gamma }}_0}, \end{aligned}$$

(44)

where $\kappa ,[s]$ is a positive constant. Employing L’Hospital’s rule, it is obvious that $\lim _{{\dot{\gamma }}_0\rightarrow |\infty |}\varOmega _v=0$ holds as before.

Fig. 10

Outline of the trajectory layer for high and infinite IHV

Fig. 11

Osculating plane (light gray) of a trajectory for high IHV (ISV is not applied): a infinite IHV; b finite IHV:—3rd half-period; c finite IHV:—5th half-period; $\gamma _{b1},\gamma _{b2}$: starting or finishing points of one full period

4.2.3 Discussion

Let us point out some properties of the trajectory discussed above, such as features of a curve passing through a layer of thickness $\epsilon$; see Eq. (39) and Fig. 10. An analysis using the small parameter approach cannot be done directly because the solution for ${\dot{\gamma }}_0\rightarrow |\infty |$, which serves as a zero approximation, does not exist in an explicit form for several reasons (infinite energy, infinite angular momentum, indefinite derivatives with respect to time, etc.). Despite this fact, we have seen that all trajectories are stable, whatever their parameter setting is. Consequently, the existence of the zero approximation can be assumed in an implicit meaning of a certain limit. Therefore, analogously with Eq. (25), we are entitled to write

$$\begin{aligned} {\dot{\gamma }}\approx {\dot{\gamma }}_0-{\dot{\eta }},\qquad \alpha \approx \alpha _2+\zeta , \end{aligned}$$

(45)

where $\alpha _2,{\dot{\gamma }}_0$ are relevant to the planar trajectory incident with the inclined plane, see Fig. 11a, and ${\dot{\eta }},\zeta$ are small unidirectional deviations. They enable us to define high values of ${\dot{\gamma }},\alpha$ but still finite values of the initial approximations $\alpha _2,{\dot{\gamma }}_0$.

We recall Eq. (11). By the way, let us note that despite the fact that the ISV is not included in this section $(\omega _n=0)$, we can see that the proportion of the spin energy at high velocity ${\dot{\gamma }}$ is negligible anyway. Thus, putting approximations Eq. (45) into Eqs. (11a,b), and comparing terms that involve the same powers of ${\dot{\eta }},\zeta$, we can write

$$\begin{aligned}&{\dot{\eta }}^0{:}\ \ddot{\alpha }_2- \left( {\dot{\gamma }}_0^2 \cos \alpha _2-\omega _0^2\right) \sin \alpha _2 =\,0, \end{aligned}$$

(46a)

$$\begin{aligned}&{\dot{\gamma }}_0 \sin ^2\alpha _2=\, H, \end{aligned}$$

(46b)

$$\begin{aligned}&{\dot{\eta }}^1{:}\ \ \ddot{\zeta } +{\dot{\gamma }}_0 \sin 2\alpha _2\cdot {\dot{\eta }} \nonumber \\&\quad -\, \left( {\dot{\gamma }}_0^2 \cos 2\alpha _2 - \omega _0^2 \cos \alpha _2 \right) \cdot \zeta =0, \end{aligned}$$

(46c)

$$\begin{aligned}&\sin ^2\alpha _2\cdot {\dot{\eta }} -{\dot{\gamma }}_0\sin 2\alpha _2 \cdot \zeta =\,0. \end{aligned}$$

(46d)

Equations (46a,b) are implicitly fulfilled and approximately represent a circular trajectory in the inclined osculating plane. Equations (46c,d) represent a rough approximation of the trajectory behavior, provided a high IHV in the sense of Eq. (45) is applied. The system Eqs. (46c,d) is linear because $\alpha _2,{\dot{\gamma }}_0$ are known parameters. Variable ${\dot{\eta }}$ can be eliminated, so it holds that

$$\begin{aligned} \ddot{\zeta }+\left( {\dot{\gamma }}_0^2\left( 2+\cos 2\alpha _2\right) +\omega _0^2\cos \alpha _2\right) \cdot \zeta =0. \end{aligned}$$

(47)

It is obvious that for high values of ${\dot{\gamma }}_0$, the term $\omega _0^2\cos \alpha _2$ is negligible. The remaining coefficient is always positive. In a ratio to the length of the circular trajectory, it represents a parameter analogous with Eq. (30). It is proportional to the velocity $\varOmega _v$ of the osculating plane rotation around the z axis, which approaches zero for ${\dot{\gamma }}_2\rightarrow |\infty |$. This result is identical with the one obtained above in this section.

Fig. 12

Example of the trajectory above the SC with no ISV $(\omega _n=0)$ for a high IHV: a time history, b top view, c axonometric demonstration

To demonstrate the character of the upper root $\delta _3$ during the IHV limitation to $\pm \infty$, let us assess its value respecting Eq. (39), and the fact that the trajectory is running in a thin layer following the scheme in Fig. 10. The process is limited to within this domain, and, consequently, it can be linearized in the framework of this layer. Making use of these facts and a factored form of the polynomial, we can reformulate the CE, Eq. (14a), because two roots $\delta _c,\delta _2$ are known:

$$\begin{aligned} \begin{aligned}&(\delta -\delta _c)\left( \delta -(2-\delta _c-\epsilon )\right) (K\delta +L)=0, \\&K=2\omega _0^2,\quad L=({\dot{\gamma }}_0^2\delta _c(2-\delta _c)+2\omega _0^2(\delta _c+\epsilon )), \end{aligned} \end{aligned}$$

(48)

which results in the third root:

$$\begin{aligned} \delta _3\;=\;\displaystyle {\frac{1}{2\omega _0^2} \left( {\dot{\gamma }}_0^2\delta _c(2-\delta _c)+2\omega _0^2(\delta _c-\epsilon ) \right) .} \end{aligned}$$

(49)

This formula demonstrates that $2<\delta _3\rightarrow \infty$ for ${\dot{\gamma }}_0^2\rightarrow \pm \infty$.

For completeness, let us demonstrate an example of a trajectory based on a high IHV (Fig. 12). The time history due to the high and nearly constant tangential velocity of the ball is homogeneous at all coordinates which interact rather on a geometrical basis. The circular character of the trajectory is obvious in pictures (b) and (c).

4.3 Trajectories below the SC

4.3.1 Roots of the CF

The CF has a form corresponding to curve (b) in Fig. 4. In general, $\delta _1$ can descend to zero, as we can easily deduce from Eq. (50) or from the original characteristic equation [Eq. (14)], with the vanishing absolute term. This case, together with the neighborhood of this value $(0\le \delta _1<\epsilon )$, will be discussed in a separate section (Sect. 6). For these initial settings, the spin-free and the spin-considered cases of initial settings intermingle, and, hence, it is worthwhile to discuss the two of them together.

The spherical strip in which the trajectory for IHV ${\dot{\gamma }}_0<\varGamma _0$ emerges is below the SC. The strip is limited by the SC, which forms its upper boundary, $\delta _{c}=\delta _2$, and by the lower boundary $\delta _1$, which is given by the quadratic term in Eq. (31). The basic analysis is similar to that which has been performed in the beginning of Sect. 4.1. The only difference is that the roots are ordered as follows: $0<\delta _1\le \delta _{c}=\delta _2<1<\delta _3$. The roots $\delta _{1,3}$ may be symbolically written as in Eq. (32):

$$\begin{aligned} \delta _{1,3}= \frac{1}{K}\left( -L\pm \sqrt{L^2-K\cdot M}\right) , \end{aligned}$$

(50)

where the discriminant has the same form as in Eq. (32b) for $\delta _c=\delta _2$.

In order to determine velocity ${\dot{\gamma }}$ at the tangent point at the $\delta _1$ boundary, we refer to Eq. (34), where $\alpha _1$ is to be substituted instead of $\alpha _2$:

$$\begin{aligned} {\dot{\gamma }}_{1T}={\dot{\gamma }}_0\;\left( \frac{\sin \alpha _c}{\sin \alpha _1}\right) ^2 ={\dot{\gamma }}_0\;\frac{\left( \delta _c-2\right) \delta _c}{\left( \delta _1-2\right) \delta _1}. \end{aligned}$$

(51)

It is obvious that ${\dot{\gamma }}_{1T}>{\dot{\gamma }}_0$ due to a lower potential energy at the $\delta _1$ level. Inspecting Eq. (14), where $\omega _n=0$ is substituted, we can see that ${\dot{\gamma }}$ is a simple continuous function of $\delta$, which is integrable on any interval $\delta \in (a,b)$ with $0<a,b<2$. This implies that no singular points emerge within one period, and the velocity along the trajectory is mildly variable. The trajectory has the shape of a simple curve without any multiple or turnabout points reversing its velocity.

4.3.2 Vertical periodicity of trajectories

As for the length of a single vertical period, a similar deduction can be made like in Sect. 4.1. However, care should be taken when the IHV approaches zero and the system Eqs. (36) becomes unstable or discontinuous in the neighborhood of $\gamma =\pi /2$. For details, see Sect. 6. Nevertheless, in a common case, like in Sect. 4.1, we can assume once again that one period consists of two identical halves symmetrically distributed around the tangent point on the lower boundary $\delta _1$, and, therefore, it is sufficient to examine only one half of the period.

Hence, the differential system Eq. (36) can also be used here, except that the zero initial conditions are formulated for the upper strip boundary, $\delta (0)=\delta _c, {\dot{\delta }}(0)=0$ and the position of the lower boundary $\delta _1$ should either be evaluated using Eq. (50e) or during solution of system Eq. (36).

Fig. 13

Shape of the trajectory below the SC for various IHV

Let us follow Fig. 13 demonstrating several trajectories comparable in their initial conditions for ${\dot{\gamma }}_0<\varGamma _0$. Recalling Eq. (36a):

$$\begin{aligned} \ddot{\delta } = E(1-\delta )-\omega _0^2\delta (4-3\delta ) \end{aligned}$$

(36a)

for $E= {\dot{\gamma }}_0^2\delta _{c}(2-\delta _{c})+2\omega _0^2\delta _{c}$, cf. Eq. (15), we can see that unlike cases with ${\dot{\gamma }}_0>\varGamma _0$, the first term of the right-hand side loses dominance, and a significant nonlinear character of the equation emerges. Retaining only the second term, we obtain an equation solvable using elliptic functions. This effect is obvious in Fig. 13, where the curves for decreasing IHV become more and more similar to an elliptic sinus. Returning to properties of the rational or irrational spiral form, see also Sect. 4.1, we can conclude that the vertical period length varies in the interval $2T_{\gamma }\in (\pi ,2\pi )$ throughout all $|{\dot{\gamma }}_0|\in (0,\infty )$, and, therefore, obviously no synchronization of the primary type can occur.

4.3.3 Discussion

We shall point out some properties of the length and shift of periods that are specific for the trajectories below the SC. Characteristics of the vertical period length and the lower strip boundary position $\delta _1$ as a function of the IHV and SPT level are demonstrated in Fig. 14; see the solid orange curves below $\varGamma _0=3.63$. Both curves are smooth on the whole interval ${\dot{\gamma }}_0$ including the SPT. It is typical that the period $2{\gamma _T}$ shortens to $\pi$ for ${\dot{\gamma }}_0\rightarrow 0$, as it also corresponds to the character of an elliptical sinus that characterizes the shape of the trajectory; see [1, 21]. The lower boundary of the strip $\delta _1$ evidently tends to zero, i.e., toward the SPC.

With reference to Sect. 6, we can see that the period for ${\dot{\gamma }}_0\rightarrow 0$ does not approach $2\pi$, as could be intuitively expected. However, even cases discussed in this section evince some attributes which are distinctly visible for low IHV values. It is obvious, for instance, that the shape of the relevant curve depends on the SPT level. The phenomenon of period “shift” becomes more prominent as the height of the SPT or the boundary $\delta _1$ level increases. This effect is discussed in more detail for both zero and nonzero ISV for low IHV in Sect. 6.

Fig. 14

A trajectory below the SC. a Spatial period dependent on the IHV $({\dot{\gamma }}_0)$, b lower boundary of the strip $\delta _1$ as a function of ${\dot{\gamma }}_0$

Fig. 15

Example of a trajectory below the SC with no ISV $(\omega _n=0)$: a time history, b top view, c axonometric demonstration

An example of a trajectory below the SC is plotted in Fig. 15. The time history in this domain appears as a simple periodic curve without any higher harmonics and with a weak multiplicative modulation. The SPC always lies inside individual loops. They run round the SPC and never pass it, whatever the SPT and IHV are. The spatial character of the trajectory is obvious in plot (c) of Fig. 15. As in Sect. 4.1, no visible intervention of higher harmonics is observed, although a very light quasi-periodic character can be noticed. The shape in the top view is helical and close to a hypotrochoid form of a prolate type. As no sharp apexes or loop multiple points are detected, we can conclude that no basic cycloid or curtate hypotrochoid is approached.

5 Trajectories influenced by an initial spin of the ball

5.1 Roots of the CE

We again revisit Eq. (14). This time the ISV will be considered nonzero, $\omega _n\ne 0$. The characteristic function Eq. (14a) can be decomposed as in Eq. (31); however, the coefficients L, M will be different:

$$\begin{aligned} f(\delta )=&(\delta -\delta _{c})(K\delta ^2+2L\delta +M)\;=\; 0, \end{aligned}$$

(52a)

$$\begin{aligned} K=&2\omega _0^2,\quad M=\frac{1}{\delta _{c}}(H-\mu \omega _n)^2,\nonumber \\ L=&\omega _0^2\delta _{c}-\frac{1}{2}(E+4\omega _0^2-\mu \omega _n^2(1-\mu )), \end{aligned}$$

(52b)

As in Sect. 4, either the lower or upper boundary of the strip is one of the roots $\delta _i,\;i=1,2$, which are determined by initial condition $\delta _c$, i.e., $\delta _c=\delta _1$ or $\delta _c=\delta _2$ for the strip situated above or below the SC, respectively. The remaining two roots $\delta _1,\delta _3$ or $\delta _2,\delta _3$ can be calculated from the quadratic term in Eq. (52a) as in Eqs. (32a) or (50):

$$\begin{aligned} \delta _{i,3}&=\frac{1}{K}\left( -L\pm \sqrt{L^2-K\cdot M}\right) ,\quad i=1,2, \end{aligned}$$

(32a,50)

$$\begin{aligned} D&=L^2-K\cdot M \nonumber \\&=\left( \omega _0^2\delta _c-\frac{1}{2}\left( E+4\omega _0^2-\mu \omega _n^2(1-\mu )\right) \right) ^2\nonumber \\&\quad -\,\frac{2\omega _0^2}{\delta _c}\left( H-\mu \omega _n\right) ^2. \end{aligned}$$

(53)

The expressions above can be reformulated with reference to Eqs. (9, 12a). Provided E, H are specified with respect to initial conditions, cf. Eq. (15), we obtain

$$\begin{aligned} K&= 2\omega _0^2,\ M=\delta _c\left( {\dot{\gamma }}_0(2-\delta _c)-\mu \omega _n\right) ^2,\nonumber \\ L&=-\frac{1}{2}\left( {\dot{\gamma }}_0^2\delta _c(2-\delta _c)+4\omega _0^2+\mu ^2\omega _n^2\right) , \end{aligned}$$

(54a)

$$\begin{aligned} D&=\frac{1}{4}\left( \mu ^2\omega _n^2+{\dot{\gamma }}_0^2\delta _c(2-\delta _c)+4 \omega _0^2\right) ^2\nonumber \\&\quad -\,2\omega _0^2\delta _c\left( \mu \omega _n-{\dot{\gamma }}_0(2-\delta _c)\right) ^2. \end{aligned}$$

(54b)

Discriminant D is positive with the exception of two cases when $D=0$:

$$\begin{aligned} \omega _n&=0, {\dot{\gamma }}_0\rightarrow \infty&\!\!\implies&&\left\{ \begin{array}{rl} \delta _{1,2}&{}=\, 0,\ \delta _3=2\\ \delta _{1}&{}=\, 0,\ \delta _{2,3}=2 \end{array} \right. \\ \omega _n&=-\frac{{\dot{\gamma }}_0\delta _c}{\mu },\ {\dot{\gamma }}_0^2=2\frac{\omega _0^2}{\delta _c}&\!\!\implies&&\qquad \delta _1=\delta _c,\ \delta _{2,3}=2. \end{aligned}$$

The first option represents two singular cases when $\omega _n=0$ and $\delta _{1,2,3}\in \{0,2\}$, while the latter one describes a particular trajectory which occurs for a given negative spin (assuming positive ${\dot{\gamma }}_0$) and passes through the NPC. This case, however, generally requires IHV ${{\dot{\gamma }}}_0\ne \varGamma _0$; only for $\delta _c=2/3$ does

$$\begin{aligned} {{\dot{\gamma }}}_0=\varGamma _0=\sqrt{3}\omega _0\ \ \text{ and } \ \ \omega _n=-2\frac{\omega _0}{\sqrt{3}\mu }. \end{aligned}$$

(55)

The other particular case occurs when $M=0$, i.e., for

$$\begin{aligned} \omega _n={{\dot{\gamma }}}_0\frac{2-\delta _c}{\mu }. \end{aligned}$$

(56)

Then, $\delta _1=0$ and $\delta _3=2+{{\dot{\gamma }}}_0(2-\delta _c)/\omega _0^2>2$.

In other cases, $K>0,M>0$, and thus, $0<D<L^2$. Hence, both roots are real and positive and fulfill the relation: $0\le \delta _1\le \delta _2<2<\delta _3$. Here, $\delta _c$ is either $\delta _1$ or $\delta _2$, depending on what is considered given. For the rest of this section, we consider the IHV to always be equal to positive velocity $\varGamma _0$, as it corresponds to the IHV of the SC.

5.2 Trajectories above the SC—negative spin

5.2.1 General considerations

We can see that the trajectories influenced by a negative ISV, $\omega _n<0$, lie within a spherical strip above the SC, with the bottom boundary $\delta _c=\delta _1$ and the upper boundary $\delta _2$, which follows from Eq. (53). Note that the third first integral, Eqs. (11c) or (12b), specifies that $\omega _n$ is constant throughout the whole investigated process.

First, we briefly outline the basic character preliminarily inspecting Fig. 19 and taking into consideration also Figs. 16, 17 and 18.

In general, observing the time history of the horizontal and vertical components of the response, it is noticeable that each trajectory consists of two basic components (except for some higher marginal harmonics). They are independent in their frequencies, as the first one is related to the basic spin-free movement and the second follows from the spinning rotation of the ball. The latter one proves to be more or less distinct in its amplitude according to the width of the strip where the particular trajectory is operating. In other words, it is determined by the active area $\delta \in (\delta _c,\delta _2)$. In general, the influence of the spin on the overall shape of the trajectory increases with the value of $|\omega _n|$ and acquires a significant dominance for $\omega _n$ above limit $\omega _{\mathrm{ns}}$, given by Eq. (60), and in particular for $0>\omega _{\mathrm{ns}}\gg \omega _n$ or $\omega _n\rightarrow -\infty$.

Fig. 16

Shape of the trajectory above the SC for various initial spin velocities; colors of curves: $\omega _n=0$—black, $\omega _{\mathrm{ns}}<\omega _n<0$—red, $\omega _n$—bold green (“kings crown” shape—separating case), $\omega _n<\omega _{\mathrm{ns}}$—blue. (Color figure online)

It follows from Eq. (14) that just three types of trajectories can be encountered when $\omega _n\ne 0$. Let us remember that the same equations, presented in Sect. 4.1, conclude that only one type of trajectory can exist within the strip if no spin of the ball is applied. The main reason for this alteration follows from the fact that the right-hand side of Eq. (14b) can have both positive or negative values, when $\omega _n\ne 0$ is considered. The three types of trajectories can be classified with respect to the parameters and initial conditions of the system, namely the ISV. The shapes of the trajectory types differ significantly in the neighborhood of the contact point on the upper boundary of the strip; see Figs. 16 and 17.

Fig. 17

Shapes of trajectories in the neighborhood of the contact point on the upper boundary of the strip $(\delta _2)$; $\omega _n=0$—black, $\omega _{\mathrm{ns}}<\omega _n<0$—red, $\omega _n$—bold green (“kings crown” shape—separating case), $\omega _n<\omega _{\mathrm{ns}}$—blue. The symbol $\varDelta \gamma$ (horizontal axis) means a local coordinate within one period or an increase/decrease of $\gamma$ with respect to $\gamma =\gamma _T$ (position of the tangential point on the $\delta _2$ boundary). (Color figure online)

5.2.2 Wavy trajectories

Let us now discuss some specific details of the individual trajectory types. The general form of the first type is obvious from Fig. 19 (i). The trajectories reflect the ISV in interval $\omega _n\in (0,\omega _{\mathrm{ns}})$, where $\omega _{\mathrm{ns}}$ is the ISV of the separating case, Eq. (60). Roughly observed, these trajectories are not too far from those discussed in Sect. 4, although some influence of the response component caused by the spin is discernible. However, the basic form again resembles an irrational spiral with slightly distorted detailed periods. The difference in the basic shape is rather quantitative, as the frequency of the spin is more or less related with that generated by the basic movement of the ball, and that is why it is hidden in the primary component influencing its amplitude. The trajectory shape is obvious from Fig. 16, where it is plotted as a function of angle $\gamma$ (the two red curves relevant to $\omega _{\mathrm{ns}}<\omega _n<0)$.

Details of the trajectory character near the contact point on the upper boundary are demonstrated in Fig. 17 (the two red curves). For the sake of a better visual comparison of individual trajectory behavior in the neighborhood of the contact point on the upper boundary $\delta _2$, all trajectory graphs have been shifted and concentrated around this point, which serves as the origin of local coordinate $\varDelta \gamma$. The starting and finishing points of one period on the lower boundary $\delta _c$ are denoted $\gamma _{b1},\gamma _{b2}$. The width of the strip increases with descending $\omega _n$ from zero until a maximum width is reached for $\omega _n=\omega _{\mathrm{ns}}$; see Figs. 17 and 18a, b. This corresponds to the total energy conservation principle.

The assessment of the length and duration of one period of the trajectory can be done analogously to Sect. 4.1. The reasoning which brought us to the differential system Eq. (36) is more or less the same, being based on the fact that contact points on the strip boundaries represent points with unavoidable singularity. Hence, with reference to Sect. 4.1, we can deduce the following modified system:

$$\begin{aligned} \ddot{\delta }&= E(1-\delta )-\omega _0^2\delta (4-3\delta ) \nonumber \\&\quad - \mu \omega _n\left( \omega _n(1-\mu )(1-\delta )+H\right) , \end{aligned}$$

(57a)

$$\begin{aligned} {\dot{\gamma }}&= \frac{H-\mu \omega _n(1-\delta )}{\delta (2-\delta )}, \end{aligned}$$

(57b)

where

$$\begin{aligned} E&= \varGamma ^2\delta _{c}(2-\delta _{c})+2\omega _0^2\delta _{c}+\mu \omega _n^2,\\ H&= \varGamma \delta _{c}(2-\delta _{c})+\mu \omega _n(1-\delta _c). \end{aligned}$$

The denominator in Eq. (57b) is identical with that in Eq. (36b) and is positive in the considered interval. The system Eq. (57) is solved for initial conditions: $\delta (0)=\delta _c,{\dot{\delta }}(0)=0$, and, eliminating the time, one obtains $\delta$ as the function of $\gamma$.

Fig. 18

Width of the strip above the SC for descending ISV or $\omega _n<0$ passing throughout all three types of trajectories; a representation as $\alpha _2-\alpha _c$ or b representation as $\delta _2-\delta _c$; c width of the space period along the coordinate $\gamma$ as a function of spin frequency $\omega _n<0$

For the selected $\delta _c=0.3$ and the relevant $\varGamma _0=3.63$, two samples for $\omega _{\mathrm{ns}}<\omega _n<0$ are plotted in Figs. 16 and 17.

Together with the three diagrams in Fig. 18, we can evaluate the character of a single period. Pictures (a) and (b) show the strip width as the difference between relevant angles $\alpha$ and parameters $\delta$, respectively. Picture (c) presents the spatial width of a single period as a function of $\omega _n$. It is obvious that a decrease in ISV $(\omega _n<0)$ leads to an increase in amplitude $\delta _2$ or $\alpha _2$ until a maximum is reached for $\omega _n=\omega _{\mathrm{ns}}$. The length of the period simultaneously decreases. The trajectory for $\omega _n\in (\omega _{\mathrm{ns}},0)$ could only have a periodic character if the ratio $2\pi \rho \sin \alpha _c/2T_{\gamma }$ is a rational number, where $T_{\gamma }$ denotes the half-period in the angular scale.

Fig. 19

Examples of trajectories with negative initial spin velocity ISV $(\omega _n<0)$: 1st row: $\omega _{\mathrm{ns}}<\mathrm{ISV}<0$; 2nd row: $\mathrm{ISV}=\omega _{\mathrm{ns}}$; 3rd row: $\mathrm{ISV}<\omega _{\mathrm{ns}}$. Column a time history, b top view, c axonometric demonstration

5.2.3 The limit case—the pointed trajectory

The second type of trajectory, see Fig. 19 (ii), occurs for the spin frequency $\omega _n=\omega _{\mathrm{ns}}$, see the bold green curve in Figs. 16 and 17. It represents a limit case separating groups below and above $\omega _{\mathrm{ns}}$. The trajectory shape in the axonometric view resembles a “kings crown.” It contains a sharp apex in the midpoint of every period which touches the upper boundary $\delta _2$. The derivative with respect to the circumferential coordinate $\gamma$ is discontinuous in this singular point, jumping from $\infty$ to $-\infty$. All velocity components vanish, and both tangential acceleration components are discontinuous. A spin energy that is proportional to $\omega _n^2$ is retained.

To determine the ISV needed to achieve this special type of trajectory, we take advantage of the fact that values of the first integrals for E and H, Eqs. (9, 12a), are constant throughout the entire time history. Indeed, at the SPT they can be expressed as follows:

$$\begin{aligned} \begin{aligned} E_{\delta _c}&=\;\varGamma ^2\delta _c(2-\delta _c)+\mu \omega _n^2+2\omega _0^2\delta _c,\\ H_{\delta _c}&=\;\varGamma \delta _c(2-\delta _c)+\mu \omega _n(1-\delta _c), \end{aligned} \end{aligned}$$

(58)

while in the sharp apex, at the $\delta _2$ level, it holds that

$$\begin{aligned} E_{\delta _2}= \mu \omega _n^2+2\omega _0^2\delta _2,\qquad H_{\delta _2}= \mu \omega _n(1-\delta _2). \end{aligned}$$

(59)

Evaluating the relevant equivalences $E_{\delta _c}=E_{\delta _2}$ and $H_{\delta _c}=H_{\delta _2}$, after some manipulations one obtains

$$\begin{aligned} \omega _{\mathrm{ns}}=\omega _n=-\frac{2\omega _0^2}{\mu \varGamma }, \quad \delta _2=\delta _c\left( 1+\frac{\varGamma ^2}{2\omega _0^2}(2-\delta _c)\right) . \end{aligned}$$

(60)

Substituting $\varGamma$ from Eq. (23) only provides real results in an unrealistic situation when $\delta _c>1$. However, when $\varGamma =\varGamma _0$, Eq. (24), Eq. (60) transforms to

$$\begin{aligned} \omega _{\mathrm{ns}}=-\frac{2\omega _0^2}{\mu }\sqrt{1-\delta _c}, \qquad \delta _2=\delta _c\, \frac{4-3\delta _c}{2(1-\delta _c)}. \end{aligned}$$

(61)

Relations Eqs. (59–61) are valid only when the apex occurs at $\delta _2$. This is characterized by an infinite curvature of the trajectory, i.e., when $0=({\dot{\gamma }}^2+{\dot{\delta }}^2)^{3/2}$. The necessary condition for ${\dot{\gamma }}=0$, given by Eqs. (9), (15), (24), is $0<\delta _c\le 2/3$, cf. also Eq. (55).

It holds obviously that $\delta _2>\delta _c$, and that the upper boundary of the trajectory can reach into the upper hemisphere of the cavity. This effect occurs when $\delta _c>1-1/\sqrt{3}$. When $\delta _c$ is increased further, the apex may reach the NPC for $\delta _c=2/3$; this case also nullifies the discriminant Eq. (54b). For even larger $\delta _c$, the apex ceases to exist. These remarks are only theoretical, because in the upper hemisphere the contact force becomes negative particularly at the apex point, where all the components of the velocity vector $\mathbf{v}$ identically vanish.

Fig. 20

Top view of one loop of the trajectory: a red: ISV $\omega _n<\omega _{\mathrm{ns}}$, b green: $\omega _n=\omega _{\mathrm{ns}}$ (passing the SPC—“separating case”), c blue: $\omega _n>\omega _{\mathrm{ns}}$. (Color figure online)

Fig. 21

Trajectory shapes below the SC for various initial spin velocities. Colors of curves: $\omega _n=0$—black, $0<\omega _n<\omega _{\mathrm{ns}}$—red, $\omega _n=\omega _{\mathrm{ns}}$—bold green (separating case), $\omega _{\mathrm{ns}}<\omega _{s}$—blue. (Color figure online)

5.2.4 Curly trajectories

The third type of trajectories, see Fig. 19 (iii), emerges provided $\omega _{n}<\omega _{sn}$ and possibly descending as $\omega _n\rightarrow -\infty$. The trajectory has a curly form making a loop every period; see Figs. 16 and 17. As before, one period starts and finishes at the same vertical level, which is given by the root $\delta _c=\delta _1$. The highest point, reached in half a period, lies on the upper circle boundary. It is given by the root $\delta _2$ according to Eqs. (53, 54). Let us turn our attention to the monotonously descending width of the strip, Fig. 18a, b for $\omega _n<\omega _{\mathrm{ns}}$ depicted either with respect to angle $\alpha$ or to parameter $\delta$.

The sharp apex—encountered in the previous type—delimiting the first and third types of trajectories prefigures a formation of two turnabout points $\gamma _{v1},\gamma _{v2}$ with tangents following their relevant meridians and with vanishing velocity ${\dot{\gamma }}$; see Figs. 16, 17 (blue curves).

The vertical position of points $\gamma _{v1},\gamma _{v2}$ above $\delta _c$ can be found with respect to conservation of E and H values along the trajectory and due to the fact that ${\dot{\gamma }}=0$ at turnabout points and ${\dot{\alpha }}=0$ at the SPT. Indeed, using Eqs. (9, 12), we can write the following relations for unknown values ${\dot{\alpha }}_v, \delta _v$ at turnabout points:

$$\begin{aligned} \begin{array}{rcl} {\dot{\alpha }}_v^2+\mu \omega _n^2+2\omega _0^2\delta _v &{}=&{} \varGamma ^2\delta _c(2-\delta _c)+\mu \omega _n^2+2\omega _0^2\delta _c, \\ \mu \omega _n(1-\delta _v) &{}=&{} \varGamma \delta _c(2-\delta _c)+\mu \omega _n(1-\delta _c). \end{array} \end{aligned}$$

(62)

Simple manipulations result in

$$\begin{aligned} \begin{aligned} {\dot{\alpha }}_v=&\left( \varGamma \delta _c(2-\delta _c)(\varGamma +2\frac{\omega _0^2}{\mu \omega _n})\right) ^{\!1/2},\\&\ \text{ where } {\dot{\alpha }}_{v1}={\dot{\alpha }}_v,\;{\dot{\alpha }}_{v2}=-{\dot{\alpha }}_v, \\ \delta _v=&{\delta _c\left( 1-\frac{\varGamma }{\mu \omega _n}(2-\delta _c)\right) ,} \end{aligned} \end{aligned}$$

(63)

which indicates that velocity ${\dot{\alpha }}$ is positive or negative in $\gamma _{v1}$ or $\gamma _{v2}$, respectively, i.e., the trajectory is rising or descending. The level $\delta _v>\delta _c$, but the difference $\delta _v-\delta _c$ descends to zero as $\omega _n\rightarrow -\infty$; see Sect. 5.4. Both expressions remind us that $\omega _n$ should exceed a certain limit in order for the formulae to be meaningful; otherwise, a curly form of the trajectory can exist.

The horizontal position of points $\gamma _{v1},\gamma _{v2}$ can be determined using the same system Eq. (57) as in the case of the first and second trajectory types, where $\gamma _{b1}=0$ is taken as a reference point. The position of the remaining points $\gamma _{v2},\gamma _{b2}$ follows from the symmetry of the second half-period. Evaluation of the period length is already obvious and represents the difference $\gamma _{b2}-\gamma _{b1}$. Note that the point $\gamma _{b1}$ always precedes the point $\gamma _{b2}$ on the advancing coordinate $\gamma$ like in the case of the first and second types of trajectories, and, consequently, the difference $\gamma _{b2}-\gamma _{b1}$ is always positive, but monotonously approaches zero for $\omega _n\rightarrow -\infty$; see Fig. 18c.

Points $\gamma _{v1},\gamma _{v2}$ delimit the boundaries of the upper part of the loop, which is defined in the interval $\gamma \in (\gamma _{v1},\gamma _{v2})$ in the framework of one trajectory period; see Fig. 18. The velocity ${\dot{\gamma }}$ in the upper part of the loop between the points $\gamma _{v1},\gamma _{v2}$ is of the opposite sign than that in the lower part of the loop. The width of interval $\gamma \in (\gamma _{v1},\gamma _{v2})$ starts from zero for $\omega _{n}=\omega _{\mathrm{ns}}$, where the curly form arises, and it increases for descending $\omega _n<\omega _{\mathrm{ns}}$. This width reaches the period length and then becomes larger than the distance between the starting and finishing points $\gamma _{b2}-\gamma _{b1}$ of one period. Although both widths subsequently approach zero, the width of the loop becomes more and more dominant. In addition, the total (curvilinear) length of one loop significantly exceeds the simple distance $\gamma _{b2}-\gamma _{b1}$.

It is useful to define a certain average (or effective) velocity $\varOmega _{\mathrm{av}}$ of the ball advancement along the horizontal SC; it is given by the time and length of one period. This ratio drops as $\omega _n$ increases. Consequently, the average velocity $\varOmega _{\mathrm{av}}$ decreases accordingly. It approaches zero for high values $\omega _n$; see Sect. 5.4. Due to the arrangement of points $\gamma _{b1},\gamma _{b2}$ on the $\gamma$ axis, it holds that $\varOmega _{\mathrm{av}}>0$. This positive sign is maintained throughout the whole interval $\omega _n<0$.

5.3 Trajectories below the SC—positive spin

5.3.1 General considerations

This section will be devoted to trajectories resulting from the same IHV $\varGamma _0$ and a positive ISV. The trajectories lie in a spherical strip below the SC, with the given upper boundary $\delta _c=\delta _2$ and lower boundary $\delta _1$, which is determined by means of Eq. (53), $i=1$. Similarly as before, we go through the main properties of these trajectories inspecting Figs. 20, 21, 22 and 23 together with three examples of the trajectory with positive $\omega _n$, which are presented in Fig. 24.

Trajectories below the SC can also be classified into three types or, in other words, into two groups separated by the special limit case corresponding to a fixed frequency $\omega _n=\omega _{\mathrm{ns}}$ like in the previous section. In general, all trajectories apparently have a spiral form of the prolate hypotrochoid type, as we can see in the top view presented in Fig. 24 (pictures in column (b) for the three types). However, the shape of these spirals differs from those above the SC. The difference between particular types for $\omega _n>0$ consists mainly in their relation to point A (SPC). Let us point out that the time history again has the distinct form of a two-component periodic process resulting from the movement of the ball and its rotation around the moving normal specified by the ISV $(\omega _n)$; see Fig. 24 column (a).

The shape of the trajectory in the neighborhood of the contact point on the lower boundary of the strip and the development of a curly trajectory are obvious in Figs. 21 and 20. The width of the strip $\delta _c-\delta _1\in (0,1)$ changes from zero ($\omega _n=0, \delta _1=\delta _c$) until reaching $\delta _c$ (maximal width) when $\omega _n=\omega _{\mathrm{ns}}$ and $\delta _1=0$. In this latter case, the trajectory passes through the SPC and represents the transition case. A further increase of $\omega _n>\omega _{\mathrm{ns}}$ leads to rising of $\delta _1$ backwards to the SC, which is reached for an infinite ISV; the width of the strip vanishes.

Fig. 22

Shapes of trajectories in the neighborhood of the contact point on the lower boundary of the strip $(\delta _1)$; $\omega _n=0$—black, $0<\omega _n<\omega _{\mathrm{ns}}$—red, $\omega _n=\omega _{\mathrm{ns}}$—bold green (discontinuous—separating case), $\omega _{\mathrm{ns}}<\omega _n$—blue. The symbol $\varDelta \gamma$ (horizontal axis) means a local coordinate within one period or an increase/decrease of $\gamma$ with respect to $\gamma =\gamma _T$ (the position of the tangential point on the $\delta _1$ boundary). (Color figure online)

5.3.2 Trajectories running around the SPC

Let us point out some important features of particular trajectory types. The first type of trajectory, see Figs. 20a, 22 (red curves) and 24 (i), is related to an ISV in the interval $0<\omega _n<\omega _{\mathrm{ns}}$. The trajectory has a spiral shape, where individual loops are prolate and running around the SPC. The influence of the second periodic component is small but still discernible; see Fig. 24(i). This phenomenon becomes more pronounced as the limit case $\omega _{\mathrm{ns}}$ is approached.

Roughly observed (as in Sect. 5.2), the trajectories do not differ significantly from those discussed in Sect. 4.3. The basic form resembles again an irrational spiral with slightly distorted detailed periods. The trajectory shape is obvious from Fig. 21, where three red curves are plotted relevant to $\omega _{\mathrm{ns}}<\omega _n<0$ as functions of angle $\gamma$.

Details of the trajectory character near the contact point on the lower boundary $\delta _1$ are demonstrated in Fig. 22 (three red curves) within a single period. This depiction in Fig. 22 was used in order to facilitate a comparison of the trajectory behavior throughout all $\omega _n$ considered in the neighborhood of the contact point on the lower strip boundary $\delta _1$. The width of the strip increases with ascending $0<\omega _n<\omega _{\mathrm{ns}}$ from zero until a maximum is reached for $\omega _n=\omega _{\mathrm{ns}}$ as it corresponds to the principle of conservation of total energy.

Fig. 23

Width of the strip below the SC for rising ISV or $\omega _n>0$ covering all three types of trajectories; a representation as $\alpha _c-\alpha _1$ or b representation as $\delta _c-\delta _1$; c width of the spatial period along the coordinate $\gamma$ as a function of spin frequency $\omega _n>0$; note a jump in the point $\omega _n=\omega _{\mathrm{ns}}$

The vertical position of the lower boundary $\delta _1$ was determined using Eqs. (53, 54). To assess the angular length and duration of one period of the trajectory and the horizontal position of other important points, the same procedure as in Sect. 5.2 can be applied. The system Eq. (57a) is solved for initial conditions: $\delta (0)=\delta _c,{\dot{\delta }}(0)=0$, and the results of Eq. (a) are put into Eq. (b). This way the time is eliminated, and one obtains $\delta$ as a function of $\gamma$.

For the selected $\delta _c=0.3$ and the relevant $\varGamma _0=3.63$, three samples for $0<\omega _n<\omega _{\mathrm{ns}}$ are plotted in Figs. 21 and 22 (red curves) in $\gamma$ and $\varDelta \gamma$ coordinates. The starting and finishing points of one period on the upper boundary $\delta _c$ are denoted $\gamma _{b1},\gamma _{b2}$ in Figs. 20 and 22. Together with three diagrams in Fig. 23, we can evaluate the character of one period. The graphs in Fig. 23a, b show the strip width expressed in $\alpha$ or $\delta$ variables. Picture (c) presents the spatial width of one period as a function of $\omega _n$. It is obvious that increasing the ISV leads to increasing the amplitude $\delta _1$ (or $\alpha _1$) up to a maximum reached for $\omega _n=\omega _{\mathrm{ns}}$. At the same time, the length of the period reaches its maximum for $\omega _n=\omega _{\mathrm{ns}}$. The periodicity of the trajectory in the interval $0<\omega _{n}<\omega _{\mathrm{ns}}$ is obvious. The half-period is denoted $T_{\gamma }$ (angular scale). Doubtlessly, it can be expected that the spiral is irrational like in Sect. 4.1, except for some special cases.

5.3.3 The limit case—the trajectory passing through the SPC

The second trajectory type represents the trajectory that separates the “lower and upper” groups with respect to frequency $\omega _{\mathrm{ns}}$, see the bold-green curve in Figs. 20b, 21, 22 and also case (ii) in Fig. 24. In this separating case, all loops pass through the SPC (point A). This means that the cubic equation Eq. (14a) possesses one zero root, in particular $\delta _1=0$. The absolute term in Eqs. (14a) or (52) vanishes, a condition which allows us to determine the corresponding ISV $\omega _{\mathrm{ns}}$:

$$\begin{aligned} H-\mu \omega _{\mathrm{ns}}=0,\ \Rightarrow \ \ \varGamma \sin ^2\alpha _c+\mu \omega _{\mathrm{ns}}\cos \alpha _c-\mu \omega _{\mathrm{ns}}=0. \end{aligned}$$

(64)

This means that for the fixed elevation $\delta _c$ of the SPT, and the implicitly defined velocity $\varGamma _0$, the following ISV should be applied:

$$\begin{aligned} \mu \omega _{\mathrm{ns}}\;=\;\varGamma _0(2-\delta _c),\qquad 0<\delta _c<1, \end{aligned}$$

(65)

in order to produce a trajectory the loops of which go through the SPC; see also Eq. (56).

Coordinate $\gamma$ becomes discontinuous when approaching the SPC; see the discontinuity in the green curve and the jump of length $\pi$ in Fig. 22. For the same reason, the derivative of the curve describing the width of the strip in $\omega _{\mathrm{ns}}$ is non-continuous (Fig. 23a). However, the one-sided derivatives with respect to $\omega _n$ at point $\omega _{\mathrm{ns}}$ are finite and equal in absolute value. Therefore, the non-smooth character of the curve is merely a result of maintaining the polar angle $\alpha$ as positive. Similar reasoning regards the discontinuity in the length of the spatial period in Fig. 23c. The $2\pi$ jump originates also from the fact that the relevant radius-vector length is always positive.

Fig. 24

Examples of trajectories with a positive “Initial Spin Velocity ISV” $(\omega _n>0)$: 1st row: $0<\omega _n<\omega _{\mathrm{ns}}$; 2nd row: $\omega _n=\omega _{\mathrm{ns}}$; 3rd row: $\omega _n>\omega _{\mathrm{ns}}$. Column a time history, b top view, c axonometric demonstration

5.3.4 Trajectories passing by the SPC

The third type of trajectory can be observed when $\omega _n>\omega _{\mathrm{ns}}$ and possibly $\omega _n\rightarrow \infty$ as the blue curves in Figs 20c, 21, 22 and also in case (iii) of Fig. 24. These initial conditions result in trajectories that pass by the SPC, which remains outside each loop; see Fig. 20c. The position of the lower boundary $\delta _1$ can be obtained from Eq. (53). It rises monotonously toward $\delta _c$ for an increasing $\omega _n$. Accordingly, as $\omega _n\rightarrow \infty$, we can observe that the width of the spherical strip diminishes; see Fig. 23a, b.

The angular length of one half of the period and the position of the turnabout points $\gamma _{v1},\gamma _{v2}$, see Fig. 22, can be evaluated in a similar way as in the case of the third-type trajectories for the $\omega _n<0$, i.e., to employ system Eq. (57). The analogous deduction from the previous section concerning limitations, singular points and numerical stability remains in force. The results are included in Figs. 21 and 22 (two blue curves).

The turnabout points are characterized by the vertical tangents in Figs. 21 and 22. Two of them, $\gamma _{v1},\gamma _{v2}$, are indicated in Fig. 22. At these points, velocity ${\dot{\gamma }}$ changes direction. Their attributes are similar to those encountered in the previous section, although their geometrical interpretation is slightly different. Also here, we can define a bottom part of the loop within points $\gamma _{v1},\gamma _{v2}$ where velocity ${\dot{\gamma }}$ is of the opposite sign. In contrast to Sect. 5.2 (negative ISV), the loops do not include any twofold point, and the period length $\gamma _{b2}-\gamma _{b1}$ becomes negative for $\omega _{\mathrm{ns}}<\omega _n$; see the two blue curves with an indication of the movement direction in Fig. 22. Both properties of the loops follow from the ordering of starting and finishing points of the period and are typical for trajectories with a positive ISV. It is obvious that as $\omega _n\rightarrow \infty$, the particular loops approach a zero amplitude; see also Sect. 5.4. The same also applies to the length of period $T_{\gamma }$.

With reference to points $\gamma _{v1},\gamma _{v2}$ and the starting and finishing points of the period $\gamma _{b1},\gamma _{b2}$, we can define the average (or effective) angular velocity $\varOmega _{\mathrm{av}}<0$ of the ball movement forward along the SC for $\omega _{\mathrm{ns}}\ll \omega _n$. It is obvious that the clockwise movement of points $\gamma _{b1},\gamma _{b2},\gamma _{v1},\gamma _{v2}$ between consecutive loops gets slower for increasing $\omega _n$, i.e., as $\varOmega _{\mathrm{av}}\rightarrow 0$. The limit case of $\omega _n\rightarrow \infty$ is discussed in detail in Sect. 5.4; see also Fig. 26.

Fig. 25

Shape of the trajectory in the neighborhood of the contact point $\gamma _{T}$; (i) upper boundary $\delta _2$ of the spherical strip for very low ISV: columns a–c—descending $\omega _{n}\ll \omega _{\mathrm{ns}}$; (ii) lower boundary $\delta _1$: columns a–c—rising $\omega _{n}\gg \omega _{\mathrm{ns}}$

Fig. 26

Example of a trajectory with negative ISV $(\omega _n\ll 0)$: a time history, b top view, c axonometric demonstration

5.4 High initial spin values

In this section, we will investigate trajectories starting from the SC, when an extremely high ISV is applied in either the positive or negative sense. We have seen in Sects. 5.2 and 5.3 that the shape of individual trajectories differs significantly if a positive or negative velocity of initial spin is introduced, even though some analogy can be noticed.

Significantly increasing the ISV, the trajectory properties for $\omega _n\ll 0$ and $\omega _n\gg 0$ become more and more related and finally produce a symmetric image with respect to the SC. Both of them are far from the transition case, in which the curly form trajectories start. For $\omega _n\ll 0$, the curly trajectory is located on the upper side of the SC. The upper boundary $\delta _2$ descends from a level above the SC to root $\delta _c=\delta _1$ as $\omega _n\rightarrow -\infty$. Provided $\omega _n\gg 0$, the lower boundary (root $\delta _1$) moves upwards to the SC, represented by $\delta _c=\delta _2$, when $\omega _n\rightarrow \infty$.

Let us outline approximately this process. Either roots $\delta _2,\delta _3$ or $\delta _1,\delta _3$ (depending on which trajectory group is being analyzed) can be determined using the quadratic equation Eq. (52) or its symbolic solution Eq. (53), respectively. For a high $|\omega _n|$, the term L becomes dominant due to the higher power of $|\omega _n|$, so it holds: $L^2\gg K\cdot M$, and, consequently, it can be approximately written for $i=1,2$:

$$\begin{aligned} \delta _{i,3}\approx \frac{1}{K}\left( -L\pm \left( L-\frac{1}{2}\frac{K M}{L}- \frac{1}{8}\frac{(K M)^2}{L^3}-\cdots \right) \right) . \end{aligned}$$

(66)

It is obvious that $\delta _3\rightarrow \infty$ and is therefore of no interest. The other root can be approximated:

$$\begin{aligned} \delta _i= -\frac{1}{2}\frac{M}{L}- \frac{1}{8}\frac{K\cdot M^2}{L^3}-\cdots \end{aligned}$$

(67)

The second and higher terms can be neglected, as they vanish for $|\omega _n|\rightarrow \infty$. Substituting now for K, M from Eq. (54b), and neglecting terms with the zero degree of $\omega _n$, we obtain

$$\begin{aligned} \delta _i\;=\;\displaystyle {\delta _c\Big (1-2\frac{\varGamma (2-\delta _c)}{\mu \omega _n}\Big ), \qquad i=1,2. } \end{aligned}$$

(68)

Formula Eq. (68) shows that $\delta _i$ converges to $\delta _c$ from above or below depending on the sign of the ISV. The width of the spherical strip decreases for increasing $|\omega _n|$ and so does the width of the curly trajectory period, as it has been defined in Sects. 5.2 and 5.3; see the plots in Fig. 25. Graph (i) represents the case $\omega _n\ll 0$ advancing from left to right, whereas (ii) regards $\omega _n\gg 0$ from right to left. The shape of the individual loops for $|\omega _n|\gg 0$ approaches a circle, which is followed in a negative or positive angular sense (with respect to the normal n), when $\omega _n\ll 0$ or $\omega _n\gg 0$, respectively. This approximative circle decreases in diameter and moves slowly along the upper or lower part of the SC.

The increment of the distance which the circle performs during one loop decreases rapidly with growing $|\omega _n|$ because the distance between the starting and finishing points $\gamma _{b1},\gamma _{b2}$ of one period (positive or negative) decreases faster than the circle diameter. In general, the effective distance passed along the SC within one period of a curly trajectory is mostly exploited by the approximating circle. Therefore, the effective velocity $\varOmega _{\mathrm{av}}$ of the ball advancement along the SC decreases. Let us remind the reader, referring to Sects 5.2 and 5.3, that the sense of this slow rotation around the z axis is either positive $(\varOmega _{\mathrm{av}}>0, \omega _n\ll 0)$—above the SC, or negative $(\varOmega _{\mathrm{av}}<0, \omega _n\gg 0)$, but at any time it holds that $\varOmega _av\rightarrow 0$, when $|\omega _n|\rightarrow \infty$.

This deduction concludes with a phenomenon which seems paradoxical at first glance. As we have seen, as $\omega _n\rightarrow \pm \infty$, both the diameter and effective horizontal velocity of the approximating circle degenerate to zero. Therefore, the ball does not evince any movement along the SC, the loops reduce to a single point, and the ball is apparently at rest, only spinning around its normal at the SPT with an infinite spin velocity. Its position, which seems to be out of static equilibrium, is fixed by an infinitely strong gyroscopic effect, whatever the IHV is. This type of trajectory is illustrated as a numerical simulation in Fig. 26: (a) displacement time history, (b) top view, (c) axonometric demonstration. The results correspond to negative $\omega \ll 0$; compare with Fig. 25(i).

6 Vertical plane-related trajectories

6.1 Transformation of the governing system

We examine trajectories which emerge for a very small initial horizontal velocity ($0\le |{\dot{\gamma }}_0|\ll \varGamma$) and a zero or small initial spin. If both IHV and ISV are zero, the trajectory of the ball defines a vertical plane passing through the initial point, ball center and point A of the cavity (SPC). With an injection of a small IHV or ISP, the trajectory mildly declines from this vertical plane. However, the distance from the vertical plane remains small. To describe the character of this spatial curve, it is worthwhile to consider the difference from the vertical plane as a small parameter. This simplification enables us to get into some special properties of this trajectory family.

With reference to Fig. 3, root $\delta _1$ will either coincide with the origin, i.e., $\delta _1=0$, or get a small positive value. The position of root $\delta _2=\delta _c$ depends on the initial position of the ball.

To conveniently describe the movement of the ball in a narrow strip adjacent to plane yz, where the azimuthal angle $\gamma$ exhibits a jump of $2\pi$ when $\alpha$ passes the zero value, it is advisable to rotate the coordinate system around the y axis.

Fig. 27

Arrangement of coordinates for investigation of cases with low IHV

The polar angle of the new coordinate system is denoted by $\xi$, and it represents an almost constant value perturbed by a small parameter, $\xi =\pi /2\pm \varepsilon$ (Fig. 27). The azimuthal angle $\zeta$ then describes movement in the vertical plane parallel to axis y. This way the formula for kinetic energy T [Eq. (7a)] remains identical when we formally substitute $\alpha =\xi$ and $\gamma =\zeta$, and the potential energy V [Eq. (7b)] is modified as follows:

$$\begin{aligned} V_{\mathrm{mod}}\;=\; mg\varrho (1-\sin \xi \cos \zeta ). \end{aligned}$$

(69)

The constant term $mg\varrho \cdot 1$ does not influence the dynamic equilibria, and, therefore, the total energy Eq. (9) can be rewritten in the form

$$\begin{aligned} {\dot{\xi }}^2+{\dot{\zeta }}^2\sin ^2\xi + \mu \omega _n^2-2\omega _0^2 \sin \xi \cos \zeta \;=\;E. \end{aligned}$$

(70)

Revisiting Eq. (1), one can write Lagrangian equations for coordinates $\xi$ and $\zeta$:

$$\begin{aligned}&\xi {:}\quad \ddot{\xi }-\left( {\dot{\zeta }}^2\cos \xi -\mu \omega _n{\dot{\zeta }}\right) \sin \xi -\omega _0^2 \cos \xi \cos \zeta =\, 0, \end{aligned}$$

(71a)

$$\begin{aligned}&\zeta {:}\quad \ddot{\zeta }\sin \xi + 2{\dot{\xi }}{\dot{\zeta }}\cos \xi -\mu \omega _n{\dot{\xi }}+\omega _0^2 \sin \zeta = 0. \end{aligned}$$

(71b)

6.2 Approximated governing system and relevant trajectories

The small parameter can be introduced as follows: $\xi \approx \pi /2-\varepsilon$; see Fig. 27. Hence, Eq. (71) gets a modified form:

$$\begin{aligned}&\xi {:} \quad \ddot{\varepsilon }+{\dot{\zeta }}^2 \varepsilon -\mu \omega _n{\dot{\zeta }} +\omega _0^2 \cos \zeta \cdot \varepsilon = 0, \end{aligned}$$

(72a)

$$\begin{aligned}&\zeta {:} \quad \ddot{\zeta }- 2{\dot{\varepsilon }}\varepsilon {\dot{\zeta }} +\mu \omega _n{\dot{\varepsilon }}+\omega _0^2 \sin \zeta = 0. \end{aligned}$$

(72b)

Provided no IHV or ISV is applied, $\varepsilon$ vanishes. Equation (72a) is fulfilled identically, and Eq. (72b) degenerates into a nonlinear pendulum equation. This equation can be solved in elliptic functions. Depending on the initial velocity ${{\dot{\zeta }}}(0)$, the solution is either periodic—for small ${{\dot{\zeta }}}(0)$—or continuously increasing in time with periodically variable velocity. Detailed discussion can be found, for instance, in [29, 34].

On the other hand, cases for $\zeta _0\ge \pi /2$ are physically meaningless due to a negative contact force. Note that the case considering $\zeta =0$ and $\varepsilon \ne 0$ represents only a perturbation, which can be treated using linearized expressions.

Under the assumption that the amplitudes of $\zeta$ are small, the nonlinear terms in Eq. (72) may be approximated as

$$\begin{aligned} \sin \zeta \approx \frac{y}{\rho },\ \cos \zeta \approx 1\ \text{ and } \ {\dot{\zeta }}^2\varepsilon \approx 0,\ {\dot{\varepsilon }}\varepsilon {\dot{\zeta }}\approx 0, \end{aligned}$$

and the system can be linearized. The following reduced system can be formulated:

$$\begin{aligned}&\xi {:} \quad \ddot{z}-\mu \omega _n\dot{y} +\omega _0^2 \cdot z = 0, \end{aligned}$$

(73a)

$$\begin{aligned}&\zeta {:}\quad \ddot{y} +\mu \omega _n\dot{z}+\omega _0^2 \cdot y = 0. \end{aligned}$$

(73b)

Fig. 28

Trajectory at a low IHV: line (i) linear approach—low level SPT, line (ii) nonlinear approach; column a no initial spin, column b initial spin included

For a nonzero ISV, equations Eq. (73) are coupled by a pair of gyroscopic forces which cause the spatial character of the trajectory even if the IHV vanishes. For a zero ISV, the system is characterized by two independent linear oscillators with identical eigenfrequencies. Consequently, no rosette form trajectory can occur for a non-trivial IHV, and a simple ellipse-like curve is observed. Indeed, considering initial conditions: IHV: $\dot{z}(0)=z_0^*$, and an initial deviation along the meridian in the plane $xy{:}\;\; y(0)=y_0$, a solution to linear system Eq. (73) can be expressed as follows:

$$\begin{aligned} w= & {} \displaystyle \frac{1}{2d} \left( (y_0(d-\varOmega _v)+z_0^*)\exp ({\text{ i }}(d+\varOmega _v)t)\right. \nonumber \\&\left. +(y_0(d+\varOmega _v)-z_0^*)\exp (-{\text{ i }}(d-\varOmega _v)t)\right) ,\nonumber \\ w= & {} y+{\text{ i }}z,\quad { \text{ d }}^2=\displaystyle {\varOmega _v^2+\omega _0^2,\quad \varOmega _v=\;-\frac{\mu \omega _n}{2},} \end{aligned}$$

(74)

which produces a pair of independent components: $w\;=\;y_0\cos \omega _0 t +{\text{ i }}z_0^*/\omega _0 \sin \omega t$, if the ISV $\omega _n=0$. Provided $\omega _n\ne 0$, but is rather small, the top view of the trajectory obviously has the form of a strongly prolate hypotrochoid. Therefore, it is worthwhile to define an affine space which meets a sub-manifold at a point in such a way as to have a second order of contact at this point. From a geometrical point of view, this case is similar to the one we discussed in Sect. 4.2. Indeed, we can assume that the shape of the spiral does not deviate much from the vertical plane during one cycle. The osculating plane rotates slowly around the vertical axis with an angular velocity $\varOmega _v=-\mu \omega _n/2$. Several plots of trajectories following this linear approach are demonstrated in Fig. 28, row (i), without or with an ISV, pictures (a) or (b), respectively.

When the $\zeta$ amplitude is large and the nonlinear character of the system must be respected, the second terms in Eq. (72) remain in force. In a certain sense, they also have the character of gyroscopic forces and lead to a rosette character of the trajectory when a small IHV is applied, see Fig. 28, row (ii), including details near the upper boundary of the strip. Observing picture (b) (a nonlinear approach), we can see that $\varOmega _v>0$ when only IHV is considered and no ISV is applied (counterclockwise). If $\omega _n>0$, then $\varOmega _v<0$ (clockwise).

As a demonstration, we show a trajectory resulting from a numerical simulation; see Fig. 29. The scheme of this figure is the same as that in Fig. 6. Although an ISV was introduced in this initial setting, we can see in the top view that a slightly counterclockwise rotating spiral emerges. Subsequent simulations for the descending SPT level showed a decreasing value of this rotation velocity and vice versa, just as we have seen above for small initial amplitudes (the level of the SPT) [Eq. (74)], and higher amplitudes [Eq. (72)], where the nonlinear character was respected. Hence, the top view of the trajectory is similar to that of the Foucault pendulum even if no ISV is applied. Note that the individual loops go around point A (SPC) regardless of the initial setting if no ISV is applied, i.e., for $\omega _n=0$.

Fig. 29

Example of the trajectory below the SC without ISV $(\omega _n=0)$ for low IHV: a time history, b top view, c axonometric demonstration

7 Conclusions

The dynamic behavior of a non-holonomic system represented by a ball moving inside a spherical cavity has been investigated using an analytical approach based on the Lagrangian governing system. The rolling of the ball is considered to be slipping-less and free of damping at the contact of the ball and cavity. The cavity is assumed to be fixed. Energy is introduced into the system by means of appropriate settings of initial conditions. The reason for this system layout consists in the possibility to differentiate individual groups of trajectory types and to investigate each one separately. This method of investigation indicated a set of trajectory types and limit cases that either separate them or demonstrate a process of limitation of some parameters to certain special values (infinity, zero, etc.). The results were compared to those obtained numerically by the principally different method following from the Appell–Gibbs approach. It can be concluded that the results coincide perfectly. The two sets of results complement one another well with regard to their strengths and weaknesses. The positions of limit cases match with respect to their relevant parameter settings. Using the analytical approach, the transition zones were also qualitatively and quantitatively examined and interpreted, which was impossible using solely numerical simulations.

For the investigation of particular properties of the system, two types of basic relations were used: (1) a Lagrangian governing system with incorporated Pfaff-type non-holonomic constraints and (2) its three first integrals corresponding to the total energy of the ball, the constant angular momentum with respect to fixed vertical axis and the constant angular momentum with respect to the common normal at the contact point of the ball and cavity. The last one showed that spin velocity is constant throughout the whole period the ball is moving along its trajectory.

A cubic algebraic equation makes up the backbone, characterizing the energy flow within the system. This equation always possesses three real roots, where the two lower roots are physically meaningful. They delimit a spherical strip on the spherical surface of the cavity inside of which the relevant trajectory emerges. The properties of these roots enabled us to define individual trajectory types, indicate transition, limit and other cases as well as perform parametric analyses within each group of trajectories.

• The “Separation Circle” (SC), a horizontal circular trajectory at a certain level of the lower hemisphere of the cavity, was revealed to be the main classification element. The selected level on a meridian of the cavity determines the particular critical value of “Initial Horizontal” and/or “Spin Velocity” (IHV or ISV, respectively) necessary to maintain this trajectory. Trajectories within the close neighborhood of the SC were examined, confirming the stability of the SC trajectory with respect to perturbation in initial conditions and to a small cross-impulse imparted to the ball.

• Trajectory types above the SC, assuming a positive IHV, are characterized by an IHV that is higher than the critical velocity and/or a nonnegative ISV. As a rule, spin-free trajectories have a form of irrational spirals similar to a prolate hypotrochoid occurring within the spherical strip. Increasing IHV beyond all limits, one approaches the limit state, which is represented by a planar trajectory that is symmetrically distributed with respect to the cavity center.

• Settings with a positive ISV can be classified into two sub-groups separated by the case that is represented by a “kings-crown” shaped trajectory. In this particular case, every loop contains a singular point (apex) which lies on the upper boundary of the spherical strip. A closed loop expands from this point for higher spin velocities, imparting a curly form to the trajectory. Lower spin velocities produce spirals without loops analogously to cases without spin, despite the fact that the shape of the spiral itself is different.

• Trajectories below the SC have a significantly different character. Depending on the decreasing IHV, the lower limit of spin-free trajectories successively descends down to the “Southern Pole of the Cavity” (SPC), where it degenerates into a point. Regarding cases with negative spin, they can also be classified into two sub-groups. The separating case is represented by a curly trajectory, where every loop passes through the SPC. Higher spin velocities result in trajectories the individual loops of which go round the SPC, while lower spin velocities produce loops that miss this point. Thus, the lower boundary of the strip again rises toward the SC.

• An interesting process is encountered when the ISV velocity tends to $\pm \infty$. The width of the spherical strip where the trajectory is being traced decreases to zero; the trajectory has a shape approaching a small circle slowly moving along the SC. For infinite spin velocity, the ball is seemingly fixed at the initial point of the trajectory due to infinite energy concentrated in the spin of the ball producing an infinite gyroscopic effect.

• The last group includes associated cases characterized by low initial energy of the ball, i.e., the IHV and/or ISP of the ball is small or vanishing when the trajectory starts at a certain height. For zero ISV and low-level initial position, the problem can be linearized. In the top view, a simple ellipse emerges, degenerating into a line segment for zero IHV. If nonlinear terms are respected, then this ellipse-like curve slowly rotates counterclockwise and remotely resembles the Foucault pendulum trajectory. Introducing a slight ISV, the elliptical trajectory shape changes into a significantly prolate hypotrochoid.

Ball-type vibration absorbers are equally effective regardless of the particular direction of excitation, yet they are mostly designed for unidirectional cases, sometimes with the assumption of small angles. It is rare to take into account the auto-parametric effects that stem from the nonlinear nature of the system, and which cause the spatial movement of the sphere in the cavity even under a unidirectional load. Although the movement of the sphere in the presence of external excitation differs from the trajectories described in this work, their character is broadly similar. As shown by numerical simulations of an externally excited case and also by experiments in analogous pendulum absorbers, the spatial motion of the damping mass along paths related to those described above shows an undesirable stability against fluctuations of possible ambient excitation. In the case of pendulum absorbers, such phenomena can be effectively mitigated by applying sufficient damping. However, this can be a problem in the case of ball-type absorbers. The spatial resonance movement of the ball in a limit cycle may negatively affect the structure. Thus, it seems necessary to include a detailed auto-parametric analysis of the complete system during the design stage.

As regards future investigation in this area, non-conservative systems should be paid attention on the basis of analytical processes. Specific difficulties should be expected related to first integrals, which will need to be generalized significantly in order to overcome variable total energy, angular momentum and spin velocity. The mutual penetration of groups and sub-groups when moving along one particular trajectory will have to be taken into account. Chaotic response processes will doubtlessly emerge, as has already been shown when doing numerical simulations. It will be necessary to reconcile with the fact that some steps used in this paper will become inapplicable. Nevertheless, some alternative methods of analysis are emerging and seem promising.