1 Introduction

Some dynamical systems are very sensitive to small changes of initial conditions leading the system to different responses. Consequently, even subtle changes of parameters may cause huge deviations in the response of such a system making the long-term predictions of the system response impossible. Such systems are called chaotic, and their mathematical model, i.e., their governing differential equations of motion, has no analytical solutions. Besides, several solutions can coexist for various initial conditions. The influence of initial conditions on coexisting solutions is most often presented in the form of multi-colored maps of basins of attraction. However, the basins of attraction are difficult to implement in higher-dimensional phase space. Bifurcation diagrams and Lyapunov exponents (which shows divergence in close trajectories defined in phase space) are most often used to distinguish a chaotic response from a periodic response. Another popular numerical tool used to analyze nonlinear systems is the Poincare cross section. In the case of periodic responses, there are single points, the number of which corresponds to the periodicity. On the other hand, in the case of chaotic solutions, on Poincare cross sections, complex geometric structures, called chaotic attractors, are observed [1, 2].

Since Galileo’s time pendulums have fascinated physicists to learn nature and to study natural phenomena. Moreover, this is the simplest nonlinear system manifesting a chaotic behavior. A physical model of such a pendulum may show a transition from regular to chaotic dynamics. Complex dynamics of a spherical pendulum has been attracting attention of scientists for a long time. Miles [3] was one of the first who considered the dynamics of a lightly damped spherical pendulum subjected to an in-plane harmonic excitation. Neglecting terms of order three and higher, the approximate analytical analysis has shown that a planar harmonic motion of the spherical pendulum is unstable. Later, Miles [4] returned to this problem and reported a number of bifurcation diagrams on the stability of planar and non-planar motion. Braynt studied analytically the transition of chaos [5], whereas [6] observed these transitions experimentally. Krasnopolskaya and Shvets [7] considered the chaotic motion of deterministic parametric oscillations of a kinematically driven spherical pendulum. Later, Shvets [8] investigated the appearance, development, and vanishing of deterministic chaos in a “spherical pendulum–electric motor of limited power” dynamical system. He described discovered chaotic attractors. Aston [9] studied the bifurcations that breaks the reflectional symmetry of the problem and lead to non-planar oscillations. Dynamics of the spherical pendulum under harmonic excitation in both horizontal and vertical directions was investigated in [10] using Melnikov approach. Later, the spherical pendulum excited vertically was studied experimentally and theoretically by Naprstek and Fischer [11, 12], Pospíšil et al. [13], and Fischer et al. [14]. They provided a large number of solutions with some approximate considerations. These solutions were classified by the Lyapunov exponent. In this paper, we continue this line of research.

As a spherical pendulum can perfectly model a behavior of a crane with a suspended payload, the dynamics of the spherical pendulum can help in understanding the payload dynamics for hoisting as well as other operations [15]. Ghigliazzaa and Holmes [16] considered the dynamics of a spherical pendulum with a moving support, studying the payload motion of a tower crane under linear acceleration and moving on a circular support. With a similar application in mind, the dynamics of a spherical pendulum under stochastic excitation and the first passage problem was studied in [17]. Perig et al. concentrated on small oscillations of a spherical pendulum with a uniformly rotating suspension center for modeling slewing crane motion [18]. La and Nguyen presented numerical simulations and laboratory tests of spherical pendulum, modeling rope using radial spring-damper, to reduce the load-carrying structure vibration [19].

The subject of the research contained in this paper is a spherical pendulum, whose mathematical model can be used to simulate the movement of payload of cranes and other reloading machines. From the point of view of the dynamics of the handling machines, such as overhead cranes or jib cranes, the spatial movement of the payload takes place when at least two drives of the winch and crane mechanisms are activated. In this paper, the phenomenological model of the spherical pendulum, which is affected by kinematic excitations along the x and y axes, has been formulated. Such model can simulate overhead-traveling cranes, which are characterized by negligible small deflections, and the impact of rail unevenness is neglected. (Kinematic vertical excitations are absent.)

2 Formulation of the mathematical model

Regarding the clearness of the schematics, the dissipating energy elements were not mapped (Fig. 1). Two kinds of elements were put into account, the first damps oscillations of the generalized coordinate ψ whereas the second one maps aerodynamic resistance. It is also worth mentioning that in the context of load positioning, the model based on a spherical pendulum with a rigid rope is an over-idealization of the problem. Therefore, in the case of this type of analysis, the elastic properties of the hoisting rope should always be taken into account. In the analyzed case, the elastic properties of the rope were not taken into account, as the positioning case was not considered.

Fig. 1
figure 1

Phenomenological model of the spherical pendulum being examined

The formal basis for derivation of differential equations are Cartesian coordinates, which uniquely define the position of the moving point A in three-dimensional space. The position of point A in the classic spherical pendulum is determined as the result of the assembly of two rotations, relative to the z and y axes and the translation of the local coordinate system along the z axis, by the value equal to the length of the rope. In the case of kinematically enforced, a column matrix of excitation q should be added to such equations:

$$ \varvec{P}_{\varvec{A}} = \varvec{R}_{\varvec{Z}} \left( \varphi \right) \cdot \varvec{R}_{\varvec{Y}} \left( \psi \right) \cdot \varvec{T}_{\varvec{Z}} + \varvec{q}, $$
$$ \varvec{P}_{\varvec{A}} = \left[ {\begin{array}{*{20}l} {a_{X} \cos \left( {\omega_{X} t} \right) - L\sin \left( \psi \right)\cos \left( \varphi \right)} \hfill \\ {a_{Y} \cos \left( {\omega_{Y} t} \right) - L\sin \left( \psi \right)\sin \left( \varphi \right)} \hfill \\ { - L\cos \left( \psi \right)} \hfill \\ \end{array} } \right], $$
(1)

where

$$ \varvec{R}_{\varvec{Z}} \left( \varphi \right) = \left[ {\begin{array}{*{20}c} {\cos \left( \varphi \right)} & { - \sin \left( \varphi \right)} & 0 \\ {\sin \left( \varphi \right)} & {\cos \left( \varphi \right)} & 0 \\ 0 & 0 & 1 \\ \end{array} } \right],\quad \varvec{R}_{\varvec{Y}} \left( \psi \right) = \left[ {\begin{array}{*{20}c} {\cos \left( \psi \right)} & 0 & {\sin \left( \psi \right)} \\ 0 & 1 & 0 \\ { - \sin \left( \psi \right)} & 0 & {\cos \left( \psi \right)} \\ \end{array} } \right],\quad \varvec{T}_{\varvec{Z}} = \left[ {\begin{array}{*{20}c} 0 \\ 0 \\ { - L} \\ \end{array} } \right], $$
$$ \varvec{q} = \left[ {\begin{array}{*{20}c} {f_{X} } \\ {f_{Y} } \\ 0 \\ \end{array} } \right] = \left[ {\begin{array}{*{20}l} {a_{X} \cos \left( {\omega_{X} t} \right)} \hfill \\ {a_{Y} \cos \left( {\omega_{Y} t} \right)} \hfill \\ 0 \hfill \\ \end{array} } \right]. $$

Whereby the kinematic excitation are harmonic functions, potential energy and Rayleigh’s dissipation function were given by equations:

$$ V = mgL\left( {1 - \cos \left( \psi \right)} \right),\quad R = \frac{1}{2}\left( {b_{1} \dot{\psi }^{2} + b_{2} v_{A}^{2} } \right). $$
(2)

In the examined model, the velocity of the mass suspended on flexible connector is given by the following mathematical expression:

$$ \begin{aligned} v_{A}^{2} & = L^{2} \left( {\sin^{2} \left( \psi \right)\dot{\varphi }^{2} + \dot{\psi }^{2} } \right) + a_{X}^{2} \omega_{X}^{2} \sin^{2} \left( {\omega_{X} t} \right) + a_{Y}^{2} \omega_{Y}^{2} \sin^{2} \left( {\omega_{Y} t} \right) \\ &\quad \quad + \;2L\omega_{Y} a_{Y} \sin \left( {\omega_{Y} t} \right)\left( {\cos \left( \varphi \right)\sin \left( \psi \right)\dot{\varphi } + \sin \left( \varphi \right)\cos \left( \psi \right)\dot{\psi }} \right) \\ &\quad \quad + \;2L\omega_{X} a_{X} \sin \left( {\omega_{X} t} \right)\left( {\cos \left( \varphi \right)\cos \left( \psi \right)\dot{\psi } - \sin \left( \varphi \right)\sin \left( \psi \right)\dot{\varphi }} \right), \\ \end{aligned} $$
(3)

After applying Lagrange equations of the second kind and taking into account adequate functions, we obtain the system of two nonlinear differential equations of the form:

$$ \left\{ \begin{aligned} & \sin^{2} \left( \psi \right)L^{2} m\ddot{\varphi } + 2L^{2} m\dot{\varphi }\dot{\psi }\cos \left( \psi \right)\sin \left( \psi \right) + L^{2} \sin^{2} \left( \psi \right)b_{2} \dot{\varphi } \\ & \; - \omega_{X} L\sin \left( \psi \right)\sin \left( \varphi \right)a_{X} \left( {b_{2} \sin \left( {\omega_{X} t} \right) + m\omega_{X} \cos \left( {\omega_{X} t} \right)} \right) \\ & \; + \omega_{Y} L\sin \left( \psi \right)\cos \left( \varphi \right)a_{Y} \left( {b_{2} \sin \left( {\omega_{Y} t} \right) + m\omega_{Y} \cos \left( {\omega_{Y} t} \right)} \right) = 0, \\ & L^{2} m\ddot{\psi } + \dot{\psi }\left( {b_{1} + b_{2} L^{2} } \right) - L^{2} m\dot{\varphi }^{2} \cos \left( \psi \right)\sin \left( \psi \right) + Lmg\sin \left( \psi \right) \\ & \; + \omega_{X} L\cos \left( \varphi \right)\cos \left( \psi \right)a_{X} \left( {b_{2} \sin \left( {\omega_{X} t} \right) + m\omega_{X} \cos \left( {\omega_{X} t} \right)} \right) \\ & \; + \omega_{Y} L\cos \left( \psi \right)\sin \left( \varphi \right)a_{Y} \left( {b_{2} \sin \left( {\omega_{Y} t} \right) + m\omega_{Y} \cos \left( {\omega_{Y} t} \right)} \right) = 0. \\ \end{aligned} \right. $$
(4)

Taking into consideration the possibility of making quantity and quality investigations at this stage, the mathematical notation of the motion equations was simplified, which mainly brings to dividing by sides the first equation by \( L^{2} m\sin \left( \psi \right) \) and the second by \( L^{2} m \). We obtain the mathematical model given by the equations:

$$ \left\{ \begin{aligned} & \sin \left( \psi \right)\ddot{\varphi } + 2\cos \left( \psi \right)\dot{\varphi }\dot{\psi } + \sin \left( \psi \right)\beta_{2} \dot{\varphi } \\ & \; - P\mu \omega \sin \left( \varphi \right)\left( {\beta_{2} \sin \left( {\mu_{\omega } \omega \tau } \right) + \mu_{\omega } \omega \cos \left( {\mu_{\omega } \omega \tau } \right)} \right) \\ & \; + P\omega \cos \left( \varphi \right)\left( {\beta_{2} \sin \left( {\omega \tau } \right) + \omega \cos \left( {\omega \tau } \right)} \right) = 0, \\ & \ddot{\psi } + \left( {\beta_{1} + \beta_{2} } \right)\dot{\psi } - \cos \left( \psi \right)\sin \left( \psi \right)\dot{\varphi }^{2} + \sin \left( \psi \right) \\ & \; + P\mu \omega \cos \left( \varphi \right)\cos \left( \psi \right)\left( {\beta_{2} \sin \left( {\mu_{\omega } \omega \tau } \right) + \mu_{\omega } \omega \cos \left( {\mu_{\omega } \omega \tau } \right)} \right) \\ & \; + P\omega \cos \left( \psi \right)\sin \left( \varphi \right)\left( {\beta_{2} \sin \left( {\omega \tau } \right) + \omega \cos \left( {\omega \tau } \right)} \right) = 0. \\ \end{aligned} \right. $$
(5)

To so transformed system of differential Eqs. (5), there were applied given transformations:

$$ \begin{aligned} & \omega_{0}^{2} = \frac{g}{L},\quad \mu_{A} = \frac{{a_{X} }}{{a_{Y} }},\quad \mu_{\omega } = \frac{{\omega_{X} }}{{\omega_{Y} }},\quad \mu = \mu_{\omega } \mu_{A} , \\ & \beta_{2} = \frac{{b_{2} }}{{\omega_{0} m}},\quad \beta_{1} = \frac{{b_{1} }}{{\omega_{0} L^{2} m}},\quad P = \frac{{a_{Y} }}{L},\quad \tau = \omega_{0} t,\quad \omega = \frac{{\omega_{Y} }}{{\omega_{0} }}. \\ \end{aligned} $$

The presented mathematical model of the spherical pendulum excited kinematically in relation to the x and y axes, constitutes a formal basis for conducting quantitative and qualitative research.

3 The results of modeling investigation

For quality assessment of the nature of solutions, one can apply the so-called Lyapunov Exponent [2, 20,21,22,23]. This index is a pivotal concept in chaos theory, which provides a tool for distinguishing an unpredictable chaotic behavior from a predictable periodic. It is considered to be a reliable and authoritative criterion of investigating nonlinear equations of motion. On this basis, the rate of separation of the initially infinitely close trajectories on the phase plane is determined ε(0):

$$ \lambda = \mathop {\lim }\limits_{\varepsilon \left( 0 \right) \to 0,n \to \infty } \frac{1}{n \tau }\mathop \sum \limits_{i = 1}^{n} \ln \left( {\frac{{\varepsilon_{i} \left( \tau \right)}}{\varepsilon \left( 0 \right)}} \right). $$
(6)

In the above equation, εi(t) represents the vectors connecting the examined trajectory of motion with the reference one. In practical applications, the computing procedure leads to averaging after multiple iterations within the adequate embedded space. Positive values of the Lyapunov exponent λ >0 indicate a chaotic behavior of the system. In the opposite case λ <0, the trajectories tend to existing stable points or periodic orbits. The results of the numerical simulation are multi-color contour maps of the solution of the largest Lyapunov exponent. The procedure of the simulations does not differ in fact from the calculations carried out for a single control parameter. The generated multi-color maps allow determining the regions of changeability characterizing the external load with irregular chaotic motion. The exemplary results of the largest Lyapunov exponent for the given initial conditions and parameters characterizing the dynamical system being investigated:\( \psi \left( 0 \right) = 0.01, \;\dot{\psi }\left( 0 \right) = 0,\;\varphi \left( 0 \right) = 0, \;\dot{\varphi }\left( 0 \right) = 0, \) µA=1, β1= 0.7, β2=0.7, \( \omega_{0}^{2} \) = 0.981 are shown in Fig. 2. The results of the included numerical simulations were obtained from the generalized coordinate φ and depict the influence of the parameter µω, on localization of the regions where the motion of the pendulum is chaotic. The use of harmonic functions (fX and fY) for numerical simulations results in the trajectories that take the shape of Lissajous curves.

Fig. 2
figure 2

Multi-color maps of the largest Lyapunov’s exponent generated for generalized coordinate φ under the assumption of the frequency ratio: a µω = 2, b µω = 0.5 for nodal initial conditions

The multi-color distributions of the largest Lyapunov exponent were calculated with the resolution of 500 × 500 points over the range of each control parameter P and ω. The map enables us to distinguish two predominant bands in which the motion of the pendulum becomes unpredictable. The first band comprises the range of parameter P∈ (0.4, 0.7), the beginning of the second one is determined by the value P =1.8. It is worth mentioning that the predictable motion of the pendulum occurs for the value ω =0.5 over the entire range of values P. If the frequency ratio is µω=0.5 (Fig. 2b), a similar distribution of the largest Lyapunov’s exponent is observed but the range of the control parameter p is slightly smaller and within the boundaries P∈ (0.5, 0.7). In reference to the observed second band, it is slightly wider because it starts around P ≈ 1.7. If the frequency ratio takes the value µω=0.5, the predicted motion is limited by ω =1. The above-mentioned values establish the starting point and the width of the bands is given for the non-dimensional frequency ω =10. Let us assume that the multi-color maps of the largest Lyapunov exponent (Fig. 2) were plotted on elastic surfaces. Their detailed comparison suggests that the map (Fig. 2b) is generated as a result of the uniform stretching of the map from Fig. 2a along the ω frequency axis so the side of a map determined by the value ω =0 is anchored along the entire range of parameter P. Mathematically, this stretching comes down to rescaling the x coordinate using the scale coefficient which in our case takes the value of 2.

To confirm the compatibility of stretching, multi-colored maps of the extended Lyapunov exponent (Fig. 2), Poincaré cross sections and amplitude–frequency spectra (Fig. 3) were plotted. In the amplitude–frequency spectrum (Fig. 3a), which corresponds to the case µω = 2, the dominant subharmonics and components are the multiple of the excitation frequency. In regard to the spectrum (Fig. 3b), the harmonic component representing the frequency of kinematic excitations plays a dominant role in the remaining harmonics being excited.

Fig. 3
figure 3

Poincare cross sections generated for the generalized coordinate φ, on the assumption of the frequency ratio: a µω = 2, b µω = 0.5 for nodal initial conditions. Dc denotes correlation dimension whose fractional value indicates chaotic solutions

Graphical images of geometric structures of chaotic attractors are presented against the background of a multi-colored density map of Poincare cross-sectional distribution points, corresponding to the intersection of the phase trajectory with the control plane. This was done because a multi-colored distribution map provides information on where the phase stream most often crosses the control plane. In the classic terms, the Poincare cross section most frequently is shown as the set of points laid on the phase plane. In order to improve results interpretation, a real nonnegative function f(x) defines the probability of a given random event:

$$ {\text{PDF}}\left( B \right) = \mathop \int \limits_{B} f\left( x \right){\text{d}}x, $$
(7)

where the integral of x represents integration over the phase space.

As the probability density function defined in Eq. 7 can take values greater than one, it is normalized to unity by integrating the surface area of available system states in the state space. This paper uses the Wolfram Mathematica software function, which, by default, generates a multi-colored map of the probability density function.

The direct comparison of the attractors and the amplitude–frequency spectrum (Fig. 3) does not indicate important discrepancies in reference to: attractors’ geometrical structures, density resolution maps of the points of Poincaré cross section, excited harmonic components of Fourier-type spectrum.

Graphically shown in the Cartesian three-dimensional space of xyz, the load trajectories were drawn for a steady-state response over seven excitation periods, whereas the steady-state motion was assumed to occur after 100 periods of excitation. The red color (Fig. 4) indicates the motion of the anchor point of the spherical pendulum. The differences in graphic images of Cartesian trajectories result from the characteristics of the excitation and in particular from the trajectory of the pendulum suspension point.

Fig. 4
figure 4

Motion trajectories of the load for the frequency ratio: a µω = 2, b µω = 0.5. Red curves indicate the excitation trajectory of the anchor point of the spherical pendulum. (Color figure online)

The bifurcation diagram is considered one of the basic tools to study the dynamics of nonlinear systems. Theoretically, it can be generated in several ways, based, among others, on time sequence or phase stream. In the majority of cases, these algorithms fulfill their task, nevertheless, in relation to our model, the obtained results were burdened with an error because they did not accurately reflect the periodicity of the system response in a wide frequency range. For this reason, in our modeling studies, bifurcation diagrams were generated using an algorithm based on Poincaré cross-sectional points. These points are uniquely determined by means of two coordinates defining displacement φ and velocity \( \dot{\varphi } \) so it is sufficient to know one of them to draw a diagram. The coordinate in our case was the quantity determining the angular displacement \( \varphi \). The results of numerical simulations are shown as two- and three-dimensional parametric plots for six different values of dimensionless frequencies ω. The obtained results of computer simulations were carried out for the frequency ratio µω = 2 (Fig. 5) and µω = 0.5 (Fig. 6), illustrating the influence of ω on the shape of the trajectory of the payload movement.

Fig. 5
figure 5

Trajectories of the payload motion referred to the bifurcation diagram (see the inset) with respect to the coordinate \( \psi \) against ω, for parameter values: µω= 2 and P = 0.75 and nodal initial conditions. The central red lines indicate the excitation Lissajous route of the pendulum support, while Poincaré points are denoted in black. (Color figure online)

Fig. 6
figure 6

Trajectories of the payload motion referred to the bifurcation diagram (see the inset) µω= 0.5 and P = 0.75. The central red lines indicate the excitation Lissajous route of the pendulum support, while Poincaré points are denoted in black. (Color figure online)

In the range of predictable periodic solutions (points I, II, III on the bifurcation diagrams in Figs. 5, 6), the identified motion trajectories of the pendulum do not differ in the geometrical shape. The solutions recorded in the zones with period doubling (points IV and V on the bifurcation diagrams Figs. 5, 6) substantially do not indicate the important differences in the geometry of the trajectory. Nevertheless, in the case of points IV and V, the calculated indicators defined as the average values of squared velocity of generalized coordinate of mathematical model show differences of observed solutions. The further increase in the frequency of kinematic excitation intensifies the phenomenon of period doubling, which results in the rapid increase of the average value of squared velocity of the generalized coordinate ψ (points V on the bifurcation diagrams Figs. 5, 6). Geometrically, the observed differences in the trajectories in the points from I to V come down only to their rotation of the constant angular value equal to − π/2 and mirroring relative to the horizontal axis unlike in chaotic solutions where the differences in the trajectory course are directly noticeable (point VI on the bifurcation diagrams Figs. 5, 6). In such cases, estimated indicators representing the average value of the squared velocity of the generalized coordinate ψ were essentially limited.

Despite the trajectories corresponding to the solutions in the points from I to V do not differ topologically, the generated on their basis Poincaré cross sections indicate important differences. It should be explained here that the intersection points of the load trajectory with the control plane are marked with crosses to represent predictable solutions. In the case of chaotic motion, the transition places through the control plane were mapped with blacked out points to make the diagrams legible. In the case of solutions obtained for µω = 2, on the cross sections there are single constant points. Whereas in the case of solutions obtained µω=0.5, the constant points are double. In the points VI, the geometrical structure of the Poincaré cross sections is shaped in the form of chaotic attractors. Accordingly, the attractor placed on Fig. 6 is a result of the rotation and mirror reflection of the attractor on Fig. 5.

3.1 Impact of kinematic excitation on the dynamics of the spherical pendulum

The dynamics of the spherical pendulum model, considered in the work, is significantly determined by the external excitation forced into a system. For this reason, the influence of external kinematic excitations on the trajectory of the load movement was examined. In particular, the focus was on the parameter characterizing the relationships between the frequencies of horizontal excitations μω. The presented bifurcation diagrams, plotted for the φ coordinate responsible for the rotation of the payload in the xy plane, were generated assuming the initial conditions: φ(0) = 0.00001, \( \dot{\varphi }\left( 0 \right) \) = 0, ψ(0) = 0, \( \dot{\psi }\left( 0 \right) \) = 0.

For µω= n, n integer (Fig. 7), the chaotic motion zones are shifted toward lower frequencies. Until the unpredictable chaotic solutions appear, the time responses of the generalized coordinates reflect 1T-periodic vibrations. It is worth noting that in the considered range of the control parameter variability ω, the geometric structure of bifurcation diagrams does not change significantly because periodic solutions with low periodicity dominate outside the chaotic motion zones ω >1. The case below is examined when the frequency ratio of horizontal kinematic excitations is given by the relation μω = 1/n.

Fig. 7
figure 7

Bifurcation diagrams plotted for different values of the frequency ratio of horizontal kinematic excitations μA = 1: a μω =3, b μω =4, c μω =5, d μω =6

The bifurcation diagrams (Fig. 8) clearly show that in the case of µω equal to the reciprocal of an integer, the system’s solution is represented by vibrations whose periodicity is determined by the denominator of the defining relation µ. This situation occurs until the dimensionless frequency ω does not reach the chaotic motion zone. It is also worth noting that in the range interval ω∈ < 1,2 > , the fixed points of the Poincaré cross sections are closer to each other. In addition, regardless of the value of the parameter µω = 1/n, unpredictable system responses appear roughly at ω ≈ 2.1. If the dimensionless frequency ω goes outside the chaotic motion zone, the periodicity of the solution doubles. As in the case analyzed previously, the situation occurs when the ratio of kinematic excitation frequency is given in the form of a rational number and defined by the relation µω = n1/n2. The periodicity of the solution for this case is also determined by the number appearing in the denominator, the relation defining the μω factor (Fig. 9). In contrast to the cases considered so far, the situation is presented in relation to chaotic motion zones because on the bifurcation diagrams drawn out the zones of unpredictable solutions increased in width. The width of the zones of unpredictable solutions increases if the fractional part of the ratio of kinematic excitations tends to unity.

Fig. 8
figure 8

Bifurcation diagrams plotted for different values of the frequency ratio of horizontal kinematic excitations μA = 1: a μω =1/3, b μω =1/4, c μω =1/5, d μω =1/6

Fig. 9
figure 9

Bifurcation diagrams plotted for different values of the frequency ratio of horizontal kinematic excitations μA = 1: a μω  = 2/3, b μω =3/2, c μω = 2/5, d μω =5/2, e μω =3/4, f μω =4/3, g μω =3/5, h μω =5/3

So far, cases have been considered when the parameter μω was given in the form of an integer or rational number. The last possibility of the frequency of external excitations affecting the spherical pendulum is chance when the frequency ratio \( \mu_{\omega } = \frac{{\omega_{X} }}{{\omega_{Y} }} \) is an irrational number. Sample bifurcation diagrams plotted for irrational values of the parameter μω are presented in the charts (Fig. 10).

Fig. 10
figure 10

Bifurcation diagrams plotted for different values of frequency ratio of horizontal kinematic excitations μA = 1: a μω = π, b μω=1, c μω =\( \sqrt 2 \), d μω =\( 1/\sqrt 2 \)

The bifurcation diagrams, plotted for the sample values of irrational numbers, clearly indicate an increase in the range of unpredictable behavior of the system. In the graphs presented, it is difficult to identify zones with periodic solutions of low periodicity. Quasi-periodic solutions can occur only in the range of very low values of the control parameter ω. The presented results illustrate the impact of the ratio of the frequency of external excitations on the periodicity of solutions of the spherical pendulum.

At the same time, the results of model tests will be limited only to the factor μω based on the number π, because for the case presented in the previous section μω = \( 1/\sqrt 2 \), we are dealing with solutions with very high periodicity. Such high periodicity suggests the presence of quasi-periodic or chaotic solutions in the whole range of control parameter variability ω. It is reasonable to ask whether the solutions obtained in relation to μω given by irrational number are quasi-periodic or periodic with high periodicity. The results of numerical calculations were depicted in the form of Poincaré cross sections, which were plotted on the basis of the Cartesian trajectory of tip mass movement projected on the xy plane (Fig. 11).

Fig. 11
figure 11

Poincaré cross sections of 3d Cartesian trajectories projected to the xy plane μA= 1: a μω =\( 1/\pi \), b μω =\( 1/\sqrt 2 \). The shadow corresponds to the PDF (as defined in Eq. 7)

In order to unequivocally state which character is assumed by the system response, DC correlation dimensions were calculated for the plotted Poincaré cross sections. In the case of periodic solutions, the correlation dimension takes values close to zero. With regard to quasi-periodic solutions, the value of the DC indicator is on the level of one.

3.2 Identification of multiple solutions in the vicinity of chaotic motion zones

We limit the test results included in this section to the case of μω = 2 because the similar results are obtained for the state of μω = 0.5, shifted with respect to the frequency axis. The modeling tests were carried out for the considered ranges of changes in the control parameter ω to identify coexisting solutions. Coexisting multiple periodic solutions were identified by searching a four-dimensional phase space (φ, \( \dot{\varphi } \), ψ, \( \dot{\psi } \)). Each periodic solution was assigned a four-element set of initial conditions which was still the basis for plotting the three-dimensional Cartesian trajectory along which the inertial element of the spherical pendulum travels. If the dimensionless frequency of external forcing ω is in the range from 0 to about 0.725, there is a single 1T-periodic solution. On the other hand, we deal with two coexisting Cartesian trajectories of payload movement in the range ω∈ < 0.725, 0.818 >, and their periodicity is correlated with the bifurcation diagram (Fig. 5). Examples of their graphic images are illustrated in the charts (Fig. 12). It should be signaled here that the symmetry of coexisting solutions, relative to the xz plane, is not caused by the symmetrical location of the initial conditions because this fact leads to the same Cartesian trajectory. If we interact the system with the frequency ω = 0.775, we deal with two symmetrical 1T-periodic solutions, which are mapped with the initial conditions: the first solutions φ = − 3.0032, \( \dot{\varphi } \) = − 0.8383, ψ = − 0.3474, \( \dot{\psi } \) = 0.2775 and the second ones: φ = 2.8072, \( \dot{\varphi } \) = − 0.9641, ψ = − 0.9973, \( \dot{\psi } \) = 0.2899 (Fig. 12a), whereas the initial conditions characterizing the 2T-period solutions graphically illustrated (Fig. 12b) take the values: first φ = − 2.7958, \( \dot{\varphi } \) = − 0.9903, ψ = − 0.2187, \( \dot{\psi } \) = 0.2809 and second: φ = − 2.7514, \( \dot{\varphi } \) = − 1.0662, ψ = − 1.0969, \( \dot{\psi } \) = 0.2457.

Fig. 12
figure 12

Coexisting Cartesian trajectories μA = 1, µω= 2: a ω =0.775, b ω =0.825

Sample detailed results of the simulation tests to illustrate coexisting multiple solutions related to specific dimensionless frequencies ω are shown in the diagram (Fig. 13). Much more multiple solutions are observed in the vicinity of zones where chaotic motion occurs. Sample results of the modeling tests were referred to the bifurcation diagram. If the system is affected by external excitation with a frequency of ω = 0.83, two symmetrical chaotic solutions coexist. At this point, we assume that the numerical values defining the initial conditions of individual system responses will be limited to the cases of periodic solutions only. The relatively small increase in frequency ω to 0.833 means that motion of load in the Cartesian space can take four different forms and two of them take the form of periodic motion with a periodicity of 3T. In the diagram (Fig. 13), these Cartesian trajectories are marked in blue (φ = − 2.8078, \( \dot{\varphi } \) = − 0.8874, ψ = − 0.4644, \( \dot{\psi } \) = 0.3292) and red (φ = − 2.7163, \( \dot{\varphi } \) = − 1.0899, ψ = − 1.1293, \( \dot{\psi } \) = 0.24). The other two data are as chaotic motion. It is worth noting here that both periodic and chaotic symmetrical trajectories are relative to the xz plane and, moreover, chaotic trajectories are “limited” by periodic trajectories. At this frequency of extortion, the bifurcation diagram does not show the presence of chaotic motion. This situation occurs because the bifurcation diagram reproduces the dynamics of the spherical pendulum with zero initial conditions.

Fig. 13
figure 13

Coexisting Cartesian trajectories related to the bifurcation diagram μA = 1, μω = 2

A similar situation occurs when ω = 0.837, but here periodic solutions represented in blue (φ = − 2.8402, \( y \) = − 0.9037, ψ = − 0.5129,\( \dot{\psi } \) = 0.3349) and red (φ = 2.7064, \( \dot{\varphi } \) = − 1.0995, ψ = − 1.1398, \( \dot{\psi } \) = 0.2377) have a periodicity of 6T. Only two symmetrical chaotic solutions occur when the frequency of the external motion source is 0.846. Two coexisting periodic solutions with 4T periodicity occur at a frequency of ω = 0.85. The first solution mapped in blue is determined by the initial conditions (φ = − 2.9636, \( \dot{\varphi } \) = − 0.9605, ψ = − 0.6554, \( \dot{\psi } \) = 0.3354). The second one, represented in red, sets the initial conditions (φ = 2.6722, \( \dot{\varphi } \) = − 1.1561, ψ = − 1.1705, \( \dot{\psi } \) = 0.2223). Multiple solutions that are not symmetrical to each other most often occur in the vicinity of chaotic motion zones.

With regard to the frequency of external excitations, given as irrational numbers, coexisting solutions for selected values of dimensionless frequency are presented. In the case of the frequency of kinematic excitations μω=1 interacting in the xy plane, coexisting solutions are observed in the range ω∈ < 1.95, 2.65 >. At this point, the two cases will be considered only. We claim, however, that in the range of low values of the control parameter ω, the load swing from the horizontal axis is much smaller than the one in the solutions obtained when the parameter μω is defined by a rational number. However, this situation occurs for μω=1 because for the cases illustrated in the diagrams (Fig. 10a, b), the “working” space of the load trajectory is comparable or smaller than the space for the coefficient determined by a rational number. The “working” space is understood as points located on the surface of the sphere with a radius L.

Three solutions coexist at the dimensionless frequency ω = 2: two of them are quasi-periodic with a periodicity of 20T and the third one corresponds to chaotic motion (Fig. 14a). It is worth noting that the chaotic solution is characterized by a significantly smaller pendulum deviation from the z axis. If the dimensionless frequency ω increases to 2.4, we also deal with three coexisting solutions, i.e., two quasi-periodic ones with a 24T periodicity and one chaotic (Fig. 14b).

Fig. 14
figure 14

Coexisting solutions for P = 0.75 μω =1 and μA = 1: a ω = 2, b ω = 2.4

4 Summary and conclusions

The oscillations of a spherical pendulum with a horizontal excitation were analyzed by means of phase portraits, Poincare cross sections with their probability densities, bifurcation diagrams and Lyapunov exponents. In particular, we studied sub- and super-harmonic (Lissajous) with respect to x and y horizontal excitation. The response amplitude was growing with the frequency of excitation leading to series of period doubling bifurcations and finally the transition to the chaotic solution. As expected, the variation in the phase of the Lissajous type of kinematic excitation does not change much the solutions. Furthermore, for the same system parameters the topology of periodic trajectories with different phases in excitation coincides. The chaotic solution is realized for a fairly high frequency. Note that this solution is characterized by the nonzero winding (which refers to the node x = 0 and y = 0) number which reflects the nonlinear composition of rotation and oscillation modes.

In handling machines such as cranes, the executive mechanisms should be controlled in such a way that the positioning procedure of the load goes smoothly. Bearing in mind the obtained results of the modeling tests, it is possible to formulate a conclusion about engineering nature:

  • The executive mechanisms of cranes should be controlled in such a way that the ratio of the frequency of dynamic interactions takes the integer value. Parameters of mechanisms settings should be received so that the working point of the machine is outside the chaotic motion zone.

  • It is absolutely necessary to avoid the case that the coefficient \( \mu_{\omega } \) assumes irrational numbers, because such frequency association favors the induction of unpredictable phenomena.